基于丰度反演及光谱变换的书画霉斑虚拟修复

侯妙乐,王庆民,3,4,谭 丽,武望婷,吕书强

[1. 北京建筑大学测绘与城市空间信息学院,北京 100044;2. 建筑遗产精细重构与健康监测北京市重点实验室(北京建筑大学),北京 100044;3. 青海省基础测绘院,青海西宁 810008;4. 青海省地理空间信息技术与应用重点实验室(青海省地理信息和自然资源综合调查中心),青海西宁 810008;5. 北京市规划和自然资源委员会西城分局,北京 100054;6. 首都博物馆,北京 100045]

0 引 言

书画是民族传统文化重要的艺术表现形式之一,但是存世久远与人为因素的影响以及书画本身较为脆弱的特质使其易遭受各种病害。霉斑是其中常见的一种——主要由于环境湿度大,纸张纤维、浆糊等受潮而产生霉菌,形成霉斑污渍。霉斑修复主要用清水和高锰酸钾、草酸水等化学试剂对霉斑处进行清洗,但清洗程度却主要依靠经验,清洗不足无法达到理想效果,清洗过量则可能损伤书画的纸张纤维,易造成纸质变脆或者粉化以及画体颜料颜色变浅等后果[1]。因此,利用数字图像或光谱成像进行虚拟修复,为书画霉斑清洗提供参考依据,在不损伤书画的同时,尽可能恢复画作原貌,已成为书画保护修复中的研究热点。

2004年Pei等结合颜色对比度增强方法和纹理合成的马尔科夫随机场模型对中国古代绘画的污渍和裂缝病害进行虚拟修复,集成了加权掩膜、环形扫描、邻域搜索方法,保证了修复后形状的完整性,避免了边缘错位[2]。Hedjam等分析了书写材料的物理特性及其在近红外波段的特征,利用全变分模型对退化文档多光谱影像进行了虚拟修复,实现了退化文档的文字增强[3]。Valdiviezo-N等设计了一种用于修复受损的古代文档的自动技术,利用像元线性解混方法得到颜料端元分布的丰度图像,用来提取文档中的孔洞和裂缝,然后利用改进的形态学区域填充算法,实现文档中缺失和被遮挡部分区域的虚拟修复[4]。Hou等提出一种最大噪声正逆变换去除污渍的书画虚拟修复方法,通过变换去除污渍特征波段,执行逆变换得到修复后的高光谱影像,实现了书画污渍的虚拟修复[5]。Hanif等针对古代手稿油墨渗透的退化问题,提出了一种非局部的稀疏影像修复方法,首先利用左右页影像对的文字重叠来识别病害像元,然后利用基于字典学习的稀疏影像修复方法重建病害像元,完成对油墨渗透部分的虚拟修复[6]。Ardizzone等提出了一种修复历史照片的解决方案,对历史照片的病害进行分类,提出了修复的框架和流程,并实现了一个原型系统,能够实现对打印照片和数字照片的修复[7]。Zhou等利用污渍的光谱特性提出了一种色彩约束的泊松编辑古书画污渍虚拟修复方法,通过光谱分析提取最佳特征波段组合,平衡了色彩约束块和特征波段梯度的调节能力,解决了梯度阈融合色彩失真的问题,能够有效地去除污渍,恢复污渍覆盖区域的颜料,但是无法保持书画的内容和结构[8]。周平平等针对油渍区域的修复,提出了一种基于高光谱影像分类线性回归的古画油渍虚拟修复方法,建立各类线性回归方程,利用受污渍影响较小的波段去校正受污渍影响较大的真彩色波段,实现污渍下颜料在真彩色波段的色调显示,但无法连贯纹理结构[9]。本研究的方法实现高光谱影像波长在433.57～949.12 nm之间的851个波段的修复,并且纹理结构和绘画内容修复自然连贯。

综上可知,利用彩色数字影像对书画表面病害的虚拟修复研究比较成熟,而针对高光谱影像的病害虚拟修复研究涉及较少。本研究采用连续最大角凸锥(sequential maximum angle convex cone,SMACC)方法对书画高光谱影像提取端元,进而对污渍端元丰度反演结果进行阈值分割,得到污渍区域,形成病害掩膜图。然后利用主成分分析(principal component analysis,PCA)将书画高光谱影像降维,取信息含量最大的前三主成分合成数字影像,利用针对数字影像的Criminisi算法进行虚拟修复,然后再执行PCA逆变换,实现污渍的虚拟修复,为书画高光谱影像的虚拟修复提供了新的方法。

1 研究数据

1.1 研究对象

中国近现代画家倪田(1855—1919)所作的《捕鱼图》,以青绿及墨黑色调为主,意境清新而富野趣。然而,画面霉斑严重,影响其外观和欣赏性。选择受霉斑影响较为严重的四景高光谱影像进行虚拟修复,其正射影像和研究区域如图1所示。

图1 研究区域Fig.1 Study areas

1.2 数据采集

数据采集使用Themis Vision System公司的VNIR/400H高光谱成像仪。这是一种推扫式扫描高光谱成像系统,其空间分辨率为1 392×1 000,光谱范围为400～1 000 nm,包含了可见光与近红外波段,一共有1 040个波段,采样间隔为0.6 nm,光谱分辨率为2.8 nm,光源为两个250 W卤素灯。

获取古书画的高光谱影像时,古书画放置在面向相机的物体支架上,古书画到相机的距离约为0.8 m,保证高光谱影像的空间分辨率约为100 dpi,每幅采集的影像为1 392像素×1 000像素。为了避免其他光源干扰测量,现场为暗室环境,使用成像仪自带的两盏卤素灯作为照明,卤素灯固定在相机两侧。相机、照明与拍摄画作形成距离与角度固定的作业模式,按照矩阵摄影方式,逐行逐列获取古书画的高光谱影像,相邻影像重叠约30%。在获取的多景高光谱影像中选取了四景霉斑污染较严重的作为研究对象。

2 霉斑提取与虚拟修复

书画霉斑区域自动提取与虚拟修复的方法主要包括高光谱数据预处理、霉斑区域提取以及虚拟修复三步,其流程如图2所示。

图2 霉斑提取与虚拟修复的流程图Fig.2 Flow diagram of the extraction and virtual restoration of stains

2.1 预处理

1) 辐射校正。高光谱成像仪采集的原始数据为辐射亮度,一般需要将其转换为反射率再进行后续处理,其计算公式为:

(1)

式中,R为校正后的反射率影像;Rraw为书画原始高光谱数据;Rwhite为参考白板数据;Rdark为在光源关闭且镜头遮盖情况下获得的暗电流数据。标准白板的反射率为99%。

2) 影像去噪。使用成像光谱仪采集的成像光谱数据通常存在噪声,需要进行去噪处理。经反射率辐射校正后,数据仍然存在一定的噪声。一般成像光谱仪在波长两端的波段其数据噪声比较大,通过观察不同波段数据,人为去除噪声较大的前100个波段和最后100个波段,选取高光谱影像的第100～950波段,即波长在433.57～949.12 nm之间的851个波段。对851个波段进行PCA变换,选择特征值包含原始高光谱图像中99%的信息量的前m个主成分,对前m个主成分进行PCA逆变换,得到去噪后的高光谱影像。

2.2 霉斑区域提取

如图2所示,霉斑区域的提取主要包括特征波段选择、SMACC丰度反演、阈值分割、数学形态学变换等步骤。

1) 特征波段选择。书画中一般包含多种颜料,可能存在部分颜料与病害的光谱曲线差异较小,将会影响病害提取的精度。因此,提出利用包络线去除方法来增强光谱特征。首先在高光谱影像的不同物质上选取感兴趣区,计算该类物质的平均光谱曲线。然后对其进行包络线去除,在处理后的光谱曲线中选择霉斑与其他物质差异最大的n个特征波段。特征波段的选择标准是提取的目标与其他目标反射光谱差异越大越好。以第1景影像为例,在处理后的光谱曲线中选择霉斑、人物皮肤、人物衣服、鱼竿和白色背景,通过比较发现,包络线去除后霉斑区域光谱曲线与其他物质在450～600 nm波长区间内存在明显差异,因此选择450～600 nm之间的244个波段作为特征波段。

2) SMACC丰度反演。SMACC是由Gruninger等基于凸锥模型提出的一种端元提取方法,可以借助约束条件来识别高光谱图像的端元波谱[10]。首先采用极点来确定凸锥,并将其定义为第一个端元。再用一个具有约束条件的斜投影,从现有的锥体中生成下一个端元。然后继续增加锥体生成新的端元。重复这个过程直至生成的凸锥中包括了已有的终端单元(满足一定的容差),或者直至满足了指定的端元波谱类别个数[11]。通过对特征波段选择的n个霉斑高光谱影像特征波段进行SMACC处理,自动提取霉斑端元波谱及其相应的丰度图像。

3) 阈值分割与数学形态学变换。霉斑丰度图像中每个像素点的值表示霉斑的可能性,选择恰当的阈值可以将其分为霉斑区域和正常区域两类,形成霉斑掩膜图。由于霉斑在宣纸上会有轻微晕染且没有十分明确的边界,为了减少因为霉斑提取而产生的影响,对提取结果进行形态学滤波运算。考虑到霉斑区域的连通性和孤立点的去除,对图像先做一次开运算再做一次膨胀运算。这样既能去除一些不需要的零碎像素点,又能使提取到的霉斑区域得到适当扩张。

2.3 高光谱影像虚拟修复

首先利用PCA变换将原始高光谱图像降维,选取PCA变换之后的第一主成分、第二主成分和第三主成分合成特征图像,因为PCA变换之后的前三主成分,其特征值包含了原始高光谱图像中98.73%的信息量。一般选择原始数据95%以上的信息量就能满足要求。再将其反射率经线性拉伸后转换到0～255。然后运用数字图像修复算法对霉斑区域进行修复,最后将修复后的图像进行线性拉伸还原到原始的反射率值域,再进行PCA逆变换,得到修复后的高光谱图像。

数字图像修复选择经典的Criminisi算法,该算法是A. Criminisi基于纹理合成的思想提出的一种非常经典的图像修复算法,在数字图像修复领域应用广泛,其核心思想是计算待修复区域中样本块的优先权[12]。当对破损区域进行修复时,首先对图像缺损区域边缘处的样本块计算它们的优先权,优先权高的样本块优先得到纹理填充。图像中的线性结构一般优先权会比较高,能够优先得到填补,这样可以保证图像结构的正确传播。

3 结果与分析

3.1 数据预处理

以第1景影像为例,选取高光谱影像的第100～950波段进行主成分变换,如表1所示,其前六个主成分的特征值包含了原始高光谱图像中98.88%的信息量。

表1 各主成分所对应的特征值与信息量Table 1 Eigenvalue and information corresponding to each principal component

选取信息量较大的前六个波段进行PCA逆变换。图3为三类不同物质在PCA正逆变换前后的400～1 000 nm光谱曲线图,从中可以看出经过去噪处理后,光谱曲线明显得到了平滑。

图3 PCA正逆变换前后光谱曲线Fig.3 Spectral curves before and after PCA forward and inverse transformation

3.2 霉斑区域提取

3.2.1特征波段选择 对霉斑的光谱曲线和画作中人物皮肤、衣服、鱼竿和白色背景的光谱曲线之间经过比较,发现霉斑的光谱曲线与人物皮肤的光谱曲线相近,而与其他画面颜色差异较大(图4a和图4b)。通过对光谱曲线进行包络线去除,可以发现光谱吸收带特征更加明显,去包络线后霉斑区域光谱曲线在450～600 nm波长区间内与其他曲线存在明显差异(图4c),因此对于第1景影像选取450～600 nm之间共244个波段为特征波段,进行后续处理。经实验,第2景影像、第3景影像和第4景影像也分别选取了450～600 nm之间244个波段作为特征波段。

图4 第1景影像及其光谱曲线Fig.4 The first image and its spectral curves

3.2.2丰度图像反演 丰度图像为特定端元的分布图,其值域为[0,1],在图像上某点的丰度值越大,表明该特定端元在该点的含量越高。分别对四景高光谱影像进行SMACC处理,其中初始端元数量是根据经验设置,通过目视观察国画表面的颜色类别数,考虑到国画复合色一般为2～3种颜料混合而成,端元数一般设置为观察到的颜色类别数的1.5倍到2倍。观察发现画作中有白色背景、蓝色衣服、褐色霉斑、红色皮肤和黑色鱼竿等五种主要颜色类别,因此该处端元数设置为10。SMACC算法运算过程中可以自动合并相似端元,以第1景影像为例,经过SMACC端元提取及相似端元合并后得到了六个端元波谱曲线与相应的丰度图像(图5)。

图5 第1景影像的六个端元波谱曲线及相应的丰度图像Fig.5 Six endmember spectral curves of the first image and their corresponding abundance images

经过SMACC处理之后,将提取的端元波谱曲线与霉斑的波谱曲线进行对比分析,采用光谱角距离(spectral angle distance,SAD)和光谱信息散度(spectral information divergence,SID)[13]来衡量提取端元光谱与霉斑光谱的相似性程度。SAD和SID值越小,说明与霉斑光谱越相似,对应的SAD和SID计算公式为:

(2)

(3)

以第1景影像为例,图4b中的霉斑光谱曲线与图5a中端元光谱曲线经过SAD、SID计算的结果如表2所示。

表2 霉斑光谱与各端元光谱相似度评价Table 2 Similarity evaluation between the spectra of stains and endmembers

以第1景影像为例,通过对比光谱角距离(SAD)和光谱信息散度(SID)发现端元3的光谱曲线与霉斑的光谱曲线最为接近,SAD和SID值也是最小,并结合目视解译,其波形和反射率也十分接近。因此,通过光谱曲线分析和目视解译,选取第三个端元的丰度图像进行后续处理。类似地,分别对第2景影像、第3景影像和第4景影像选择霉斑对应的丰度图像。

3.2.3密度分割与形态学滤波 密度分割的目的是最大限度地将霉斑区域和背景区域分开。以第1景影像为例,图6为不同的分割阈值对霉斑提取结果的影响。最终确定第1景影像的第三端元丰度影像在灰度分割阈值为0.19时,可以最大限度将霉斑与其他物质区分开来。再利用3×3的结构元对图像先做一次开运算再做一次膨胀运算(图6f)。运算核的大小会对膨胀腐蚀的效果产生影响,如5×5的运算核比3×3的运算核产生的膨胀腐蚀效果更加明显。由于提取的霉斑区域比较零散、碎小,为防止出现遗漏,所以采用3×3小运算核,这样既能去除一些不需要的零碎像素点,又能使提取到的霉斑区域得到扩张。

图6 霉斑提取结果Fig.6 Results of stain extraction

3.3 霉斑区域虚拟修复

3.3.1PCA降维与Criminisi虚拟修复 以第1景影像为例,对高光谱影像进行PCA变换,取前三主成分合成假彩色影像(图7a)。然后利用Criminisi算法对霉斑区域进行虚拟修复,得到的结果图像如图7b所示。

图7 PCA前三主成分假彩色合成图像Fig.7 False color images synthesized using PCA first three principal components

3.3.2PCA逆变换 对前三主成分合成的特征影像进行Criminisi修复后,进行反向线性拉伸,使其影像像元灰度值由0～255变换到原来的值域范围(如第一主成分为-1414.547974～854.142273),然后利用PCA逆变换将影像变换到原高光谱波段数,完成高光谱影像的虚拟修复。图8为第1景高光谱影像虚拟修复后部分波段的效果。

图8 虚拟修复后部分波段Fig.8 Some bands after virtual restoration

3.4 结果分析

对第2～4景高光谱影像进行同样的处理。图9为《捕鱼图》第1～4景高光谱影像修复结果。

图9 霉斑区域虚拟修复结果Fig.9 Virtual restoration results of stains of the four images

图9a、图9d、图9g和图9j为四个样本区域修复前的高光谱图像,图9b、图9e、图9h和图9k为运用本研究方法提取霉斑后的掩膜图像(白色区域为待修复区域),图9c、图9f、图9i和图9l为修复后的图像。观察图9可知经过虚拟修复后,样本数据中霉斑区域得到了较好的修复,修复后的霉斑区域与画体其他区域本身融入性较好,修复后的边界也比较自然、平滑。

因为无法得到霉斑区域的真实灰度值,对霉斑虚拟修复的客观量化评价比较困难。理论上,霉斑区与其邻近的同一颜色非霉斑区应具有相似的灰度值,修复效果越好,两个区域对应的灰度值越接近。因此,借鉴遥感影像阴影去除后精度评价方法[14],根据人为判断霉斑覆盖区域的颜色,在其邻域非霉斑区随机选择对应颜色覆盖的20个像元作为参考真值(理论上做参考真值的像元数目越多越好,但是其邻域非霉斑区物质种类复杂,根据人为判断霉斑覆盖区域的颜色,在其邻域非霉斑区随机选择对应颜色覆盖的20个像元,足以代替霉斑覆盖区域作为参考真值),将其平均后与霉斑虚拟修复前后霉斑区域同种颜色的像元进行比较,计算修复前后霉斑区域与参考真值的均方根误差(root mean square error,RMSE),来客观评价影像修复的质量,其单个波段的RMSE计算公式为:

以上述方法,分别统计了虚拟修复前后高光谱数据中的R(640.3051 nm)、G(549.7856 nm)、B(460.2024 nm)三个波段霉斑区域各种颜色与周围非霉斑区域对应颜色像元之间的RMSE,对RGB三个波段的RMSE进行平均后,其结果如表3所示。从表3中可以看出,虚拟修复后的RMSE值相对于修复前的普遍变小,说明虚拟修复后霉斑区域的颜色值更接近无霉斑区域对应的颜色值,虚拟修复效果明显。

表3 虚拟修复前后RGB波段霉斑区与对应邻域像元的平均RMSETable 3 Average RMSEs of RGB band stained areas and corresponding neighboring pixels before and after virtual restoration

(续表3)

为了验证方法的有效性,选取倪田《捕鱼图》的另外两景(第5景和第6景)高光谱数据,按照提出的方法进行了同样的处理,结果如图10所示。

图10 其他两景影像霉斑区域虚拟修复结果Fig.10 Virtual restoration results of stains of other two images

以第1景影像为例,各波段修复效果如图11所示,从图中可以看出,经过虚拟修复后,两景影像中霉斑区域均得到了较好的修复,修复后的病害区域效果自然,无明显修复边界。

图11 虚拟修复后部分波段Fig.11 Some bands after virtual restoration

4 结 论

本研究以古书画中的霉斑为研究对象,设计了高光谱影像虚拟修复方法:首先通过光谱变换降维、SMACC丰度反演提取霉斑特征丰度图,利用灰度分割实现提取污渍掩膜;然后利用PCA变换,采用Criminisi修复PCA前三主成分,再执行PCA逆变换,实现了书画表面污渍高光谱影像的虚拟修复。和现有的虚拟修复方法相比,本方法可以实现高光谱影像的修复,可以为书画污渍的清洗提供参考。同时,仍需对本方法做进一步研究,以提高污渍提取的准确性,并且对高光谱影像修复前后的数据质量评价展开深入的研究。