光电工程  2019, Vol. 46 Issue (2): 180301      DOI: 10.12086/oee.2019.180301     
结合灰度信息的压敏漆图像配准
梁诚1 , 蒲方圆1 , 梁磊2 , 高志升1     
1. 西华大学计算机与软件工程学院,四川 成都 610039;
2. 中国空气动力研究与发展中心,四川 绵阳 621000
摘要:压敏漆技术是一种经济性高、速度快的风洞测压前沿技术。在风洞试验中,由于强风影响,模型会发生畸变,造成有风图像和无风图像难以配准,从而严重影响测压精度。针对这一问题,本文创新性的将二维非刚性ICP算法用于此问题,采用点云方式使得图像细节区域有效配准,同时也有利于后续三维重建工作。然而由于二维非刚性ICP算法仅考虑二维坐标位置关系,忽略压敏漆图像像素灰度具有的相关性,使得配准精度不高。直接利用三维非刚性ICP算法又会发生误配准,所以为了进一步提高配准精度,本文提出了一种基于像素关联搜索策略的非刚性ICP算法,算法设计了综合考虑2D坐标与像素灰度值的双目标搜索策略,实现了精确的局部匹配点搜索与双目标优化。在多组压敏漆图像上将本文算法与五种配准算法进行了对比实验分析。实验结果表明,本文所提出的算法具有最好的配准精度。相比次优算法,RMSE提升超过15%,NMI提升在5%左右。
关键词压敏漆图像    图像畸变    非刚性ICP    归一化互信息    图像配准    
Pressure sensitive paint image registration combined with gray level information
Liang Cheng1, Pu Fangyuan1, Liang Lei2, Gao Zhisheng1     
1. School of Computer and Software Engineering, Xihua University, Chengdu, Sichuan 610039, China;
2. China Aerodynamics Research and Development Center, Mianyang, Sichuan 621000, China
Abstract: Pressure-sensitive paint technology is a wind tunnel pressure measurement frontier technology with high economical efficiency and high speed. In the wind tunnel test, due to the strong wind, the model will be distorted, making the wind image and the windless image difficult to register, which will seriously affect the pressure measurement accuracy. In response to this problem, this paper innovatively applies the two-dimensional non-rigid iterative closest point (ICP) algorithm to solve this problem. The point cloud method is used to make the image detail area to be effectively registered, and it is also conducive to the subsequent three-dimensional reconstruction work. However, due to the two-dimensional non-rigid ICP algorithm, only the two-dimensional coordinate positional relationship is considered. The correlation of the pixel grayscales of the pressure-sensitive paint image is neglected, so that the registration accuracy is not too high. However, if the three-dimensional non-rigid ICP algorithm is directly used, misregistration will occur. Therefore, in order to further improve the registration accuracy, this paper proposes a non-rigid ICP algorithm based on pixel-based search strategy. The algorithm designs a dual-target search strategy that takes 2D coordinates and pixel gray values into consideration and achieves accurate local matching, realizing point search and double goal optimization. The algorithm of this paper is compared with five registration algorithms on multiple sets of pressure sensitive paint images. The experimental results show that the proposed algorithm has the best registration accuracy. Compared to the suboptimal algorithm, the RMSE is improved by more than 15% and the NMI is increased by about 5%.
Keywords: pressure-sensitive paint image    image distortion    non-rigid ICP    normalized mutual information    image registration    

1 引言

近年来,随着我国航空产业的快速发展,针对设计制造新型航空器所需的时效与成本提出了更高的要求,作为关键环节的风洞测压试验更是如此。目前流行的方法是在试验模型表面涂上压强敏感涂料(简称压敏漆,pressure sensitive paint,PSP)进行测压。其原理是压敏漆会根据试验模型表面压强不同而进行动态氧化猝灭反应,从而发出不同强度的光线。通过拍摄试验模型在有风与无风条件下的两幅图像,利用Stern-Volmer[1]公式计算对应像素点压强数据,得到整个模型表面的二维压强分布,最后将二维数据映射到三维试验模型上即完成测压工作。根据Stern-Volmer公式计算模型表面压强,关键是必须保证两幅图像像素点对应,一旦出现偏差,非对应点进行相除时会严重影响算法的精度。但实际试验中,由于存在空气动力载荷,会使得试验模型产生形变和位移,特别是强风条件下,试验模型易发生弹性形变。如大展弦比飞机机翼或直升飞机螺旋桨机翼,都呈现出几何非线性弹性形变。这就需要将有风情况下的畸变图像和无风情况下的参考图像进行精确配准。

虽然图像配准技术已经提出并发展了几十年,但实现特定场景或领域的图像配准仍具有挑战性并值得进一步研究。就压敏漆技术而言,提高测量精度仍是其主要研究方向,如利用新型涂料[2],或者环境控制等[3]。2004年后出现了少量的压敏漆图像处理与配准的研究工作[4-8],但无论是国内文献还是国外文献相关研究都很少。

早期研究一般先进行滤波降噪,提高两幅图像相关性,再提取特征点进行刚性配准。如文献[5]通过小波变换降低噪声。文献[7-8]是对畸变不明显图像提取少量特征点,再进行刚性或简单多项式变换,实现图像配准。文献[6]提出了一种较新颖的思路,利用人工标志点和自选边缘点构造Delaunay三角网格剖分,计算浮动图像与参考图像对应三角形仿射变换参数,实现整幅图像的配准。但该方法需要标志点数量足够多,而且每个三角形内部点被视为同参数变换。概括起来目前压敏漆图像配准往往是在畸变不是特别明显的情况下,提取少量特征点,再进行刚性配准,或是整体的非刚性配准,且利用的非刚性配准算法较为简单,如仿射变换、低阶多项式变换等。

非刚性图像配准方法主要有以下几种类型,基于特征的算法,如SIFT与形状信息结合的算法[9];以及将两幅图像同时处理,利用中空间(mid-space)特征配准[10];点云配准算法可以归于基于点特征的配准算法,其中经典的迭代最近点算法(iterative closest point,ICP)[11]是刚性算法,有结合主成分分析的改进算法[12]等。本文非刚性ICP[13]算法的正则项特征点采用人工标志点。另外还有利用特征点作为模型主体结构,非特征点不参与计算,而是同参数变换的方法,如文献[6],但这需要足够多的特征点。对于非刚性ICP算法有将统计形状模型纳入其框架去除异常点的改进算法[14],以及分区域统计的算法[15]等,但都没有将灰度因素考虑在其搜索策略中。另外点云的非刚性配准算法还有薄板样条鲁棒点匹配(TPS-RPM)及其改进算法[16],一致性点漂移算法CPD及其加速算法[17]等。由于压敏漆图像分辨率很高,单幅图像模型所占像素点在20万左右,甚至更大,利用其他非刚性点云配准算法的工作将在后续研究中进行。非刚性图像配准还有基于灰度信息的算法,比如互信息与B样条相结合的算法[18],互信息与空间信息结合的算法[19],互信息与薄板样条结合的算法[20]

基于特征和灰度的配准算法各有优缺点,如基于特征的算法对于图像内部结构、空间对应关系等更加明确,但无论是提取特征点还是人工标志点都易受提取或标定精度影响。而基于灰度的配准算法可全自动配准,人为因素少,但其空间变换不具备微分同胚性,所以本文结合两者特性进行配准。

由于试验模型会发生非线性弹性形变,使所得图像产生畸变,但这种畸变又不能仅仅通过简单的非刚性变换(如仿射变换、投影变换)来对图像配准。原因是相机位于试验模型上方,由于形变和位移,镜头所得图像会发生近大远小的投影变换,但并不能保证模型整体都以一定角度或范围发生投影变换,往往模型局部会发生更为明显的形变,如图 1所示。

图 1 畸变图像(a)与原图(b)对比图 Fig. 1 Comparison between distortion image and original image

针对上述问题,本文提出了一种结合灰度信息的点云配准算法。本文是由粗配准到精配准的过程,首先利用活动轮廓模型提取二维图像模型点云[21],再利用迭代最近点[11]或互信息配准[22]进行刚性配准,实现粗配准,最后将像素灰度因素加入非刚性ICP[13]的搜索算法,迭代搜索新定义的对应点,最终实现精配准。

2 结合灰度信息的配准算法 2.1 基于归一化互信息刚性配准

使用互信息(mutual information, MI)进行图像配准是由Collignon[23]和Viola[24]等人率先使用的。为了克服相应缺点,Studholme[25]等提出了归一化互信息(normalized mutual information, NMI)作为相似性测度指标,表达式如下:

$ NMI(R, F) = \frac{{H(R) + H(F)}}{{H(R, F)}}, $ (1)

其中:H(R)和H(F)分别为参考图像与浮动图像的信息熵,H(R, F)为两幅图像的联合信息熵。根据后文配准算法流程,满足一定条件时,利用文献[22]方法对图像进行刚性互信息配准,可得到旋转变换参数RMI和平移变换参数TMI。在完成刚性粗配准后再进行非刚性高精度配准。

2.2 非刚性ICP(non-rigid ICP)

图像配准可分为迭代算法和非迭代算法,Besl等人提出的迭代最近点(iterative closest point,ICP)[11]算法就是一种经典且应用广泛的算法,其可有效实现二维或三维的点云刚性配准。基本思想是将待配准点云与目标点云之间的欧氏距离均方误差之和作为目标函数,不断搜索对应目标点进行迭代,最终求得点云间的刚性变换参数,即旋转变换参数RICP和平移变换参数TICP。由于采用点云配准,常用“目标”一词,所以本文不区分“参考图像”与“目标图像”,“浮动图像”与“待配准图像”。

而Amberg[13]等人在传统刚性ICP算法基础上,利用Delaunay三角网格剖分,进一步考虑待配准点云中的点间关系和标志点因素,实现了点云的非刚性配准。其基本思想是有两个点云集合VT(Target),首先使用三角剖分将其构成网格,然后利用刚性配准算法得到的变换参数组成非刚性ICP变换的初始变换矩阵X1,对于三维图像,X1为3nx4n的矩阵;对于二维图像,X1为2nx3n的矩阵。非刚性ICP配准的目的,就是通过k次迭代,最终找到使得待配准点云最近似目标点云的变换矩阵Xk。其目标函数可表示为

$ E(\mathit{\boldsymbol{X}}) = {E_{\rm{d}}}(\mathit{\boldsymbol{X}}) + \mathit{\boldsymbol{\alpha }}{E_{\rm{s}}}(\mathit{\boldsymbol{X}}) + \beta {E_{\rm{l}}}(\mathit{\boldsymbol{X}}), $ (2)

其中:Ed(X)是对应点之间欧氏距离均方误差之和,Es(X)和El(X)分别为平滑惩罚项和标志点惩罚项。

平滑惩罚项是约束点云中的点间关系的正则项,其惩罚系数α是为实现模型从刚性到非刚性平滑过渡而人为设定的从大到小的数组。当某αi较大时,待配准点的点间约束较大,模型整体近似刚性变换;当αi较小时,点间约束较小,允许相邻点向不同方向移动,模型呈现非刚性变换。就经验而言,数组α的设置应该与所提取点云的数量级相当,本文实验中所提取点云像素点在20万左右,α为一个最大值50000,最小值100,间距均匀的100维数组。

标志点惩罚项是为进一步提高非刚性配准精度,选取两者点云中的特征点对模型进行的约束。目前模型的特征点需要人工选取,一般选取模型上的角点,以及2~3个边缘点。标志点惩罚项系数β为一个较大的常数,其数量级也应与点云像素数量相当或者更大,本文中系数β设置为30万。式(2)变为矩阵形式如下:

$ \overline E (\mathit{\boldsymbol{X}}) = \left\| {\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{WD}}}\\ {\mathit{\boldsymbol{\alpha M}} \otimes \mathit{\boldsymbol{G}}}\\ {\beta {\mathit{\boldsymbol{D}}_{\rm{L}}}} \end{array}} \right]\mathit{\boldsymbol{X}} - \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{WU}}}\\ 0\\ {{\mathit{\boldsymbol{U}}_{\rm{L}}}} \end{array}} \right]} \right\|_{\rm{F}}^2, $ (3)

其中:W是点权重矩阵,D是待配准点矩阵,U是对应目标点矩阵,构成Ed(X)。M为点间关联矩阵,G是对角权重矩阵,用于权衡变换矩阵中各个分量的比重,一般设为1,构成Es(X)。DL为待配准点云标志点矩阵,UL为目标点云中对应标志点矩阵,构成El(X)。通常设定一个很小的τ,当上式小于等于τ时,求得Xk

2.3 基于像素灰度因素的改进搜索算法

无论是刚性还是非刚性迭代最近点算法定义最近点的依据都是点集之间的欧几里得距离,即仅根据像素坐标距离判断点云是否达到配准,并未考虑点与点之间其他因素,在一般的情况下这是有效的。但通过相关实验发现,结合像素灰度能使有风图与无风图具有更高的相关程度;同时,如果直接将灰度因素作为z轴数据,进行三维点云配准,会发生误配准。受此启发,本文对传统算法最近点的定义加以修改,不再以欧氏距离为唯一依据,而是将欧氏距离最近点与像素间灰度差异最小点综合考虑, 本文改进部分主要在精配准过程,利用了非刚性ICP中的剖分网格。式(3)目标函数,每次迭代会发生改变的项仅有对应目标点矩阵Uα,而α是控制模型从刚性到非刚性平滑过渡的超参数,需要人为设定。而传统算法U是采用在KD-Tree[26]上进行KNN搜索来实现的。往往是1NN,即每次仅搜索到与待配准点vi空间距离最近的一个对应点ui,将找到的所有ui组成对应目标点矩阵U,每次迭代重新搜索组成新的U,最终使目标函数满足条件。但这样并不能利用像素灰度信息,为进一步提高配准精度,本文将传统的U变为UU是结合像素灰度因素后的对应目标点矩阵。

$ \mathit{\boldsymbol{\overline U}} : = {[\overline {{u_1}} , ..., \overline {{u_n}} ]^{\rm{T}}}, $ (4)

其中:$\overline {{u_n}} $是某个待配准点vi考虑灰度因素后的对应点。

后续搜索算法需要满足两个假设条件,第一待配准点已经在最终对应点附近,第二待配准图像与目标图像基本对齐,整体不发生大角度的旋转。由于当前步骤是在粗配准之后的精配准算法中进行的,所以有理由认为,上述假设条件成立。

点云数据利用三角形剖分形成网格,不同剖分结构大致相同,一般一个节点与6到7个邻域节点相连,当然边缘处连接点较少,但本算法同样适用。为更形象地阐述本文算法,截取图像中某点加以呈现,如图 2局部网格节点图。

图 2 局部网格节点图 Fig. 2 Local grid node diagram

红色网格为待配准点网格,黑色网格为目标点网格,两组网格形状相似,其符合非刚性配准的初始条件并满足上文两点假设。以待配准点vi为例进行说明,利用KD-Tree进行KNN搜索,k一般取3~9之间,本文取3,如图 2,点vi在三角形abc中,其距离点abc的欧氏距离分别为$D({v_i}a)、D({v_i}b)、D({v_i}c)$。其中

$ D({v_i}a) = \sqrt[2]{{{{({x_{{v_i}}} - {x_a})}^2} + {{({y_{{v_i}}} - {y_a})}^2}}}, $

$D({v_i}b)$$D({v_i}c)$类似,不失一般性的假设$D({v_i}a) < D({v_i}b) < D({v_i}c) $。如果仅以二维欧氏距离大小为标准,点vi的对应点应是点a。新搜索算法的目的是使待配准点搜索的对应点,不仅满足欧氏距离较近的原则,同时满足像素灰度接近,如图 3 x-pixel intensity三点截面图,当然在y方向或其他方向的截面图同理。

图 3 x-pixel intensity三点截面图 Fig. 3 x-pixel intensity three-point longitudinal section

这里以有风图为待配准图像进行说明,红色部分是待配准点x方向与像素灰度方向的三点截面曲线,黑色部分是目标点三点截面曲线,有风图像像素灰度往往会低于无风状态下灰度。实际为离散型,为了便于观察,显示为曲线,曲线形状任意,以钟型为例。这里以点v2搜索其对应点,仅考虑二维坐标的欧氏距离而言,点v2应向u3方向移动,因为点v2u3坐标间距离D(v2u3)小于点v2u2间距离D(v2u2)。由于是非刚性变换,待配准各点间移动方向可以不同,只要最终满足目标函数条件,所以如果点v2u3方向微量的移动并不影响最终结果。但根据直观经验和相关实验,点v2的对应点应是u2,所以对应点的选择有必要考虑像素灰度差异(intensity difference)D1(v2u2)、D1(v2u3)。上面是只考虑三个节点的简单模型,存在多个节点时,如图 4所示。

图 4 x-pixel intensity多点截面图 Fig. 4 x-pixel intensity multi-point longitudinal section

u′、u″,v′、v″泛指其他节点,仍以点v2搜索其对应点为例,当存在一个点u″时,点v2u″无论是坐标间距离$D{\rm{(}}{v_{\rm{2}}}u′′{\rm{)}}$,还是像素灰度差异${D_{\rm{I}}}{\rm{(}}{v_{\rm{2}}}u′′{\rm{)}}$都要小于点v2与点u2间的距离与差异。但正如前文所说,点v2的对应点应是u2,于是说明了在考虑像素灰度差异时,不能仅以某一个待配准点与目标点云中某一点的灰度差异为标准,而应该以某一点为中心,在一定范围内进行像素灰度差异比较,才能使得像素灰度整体趋势相似。这个范围本文利用三角形剖分后,以某点为中心邻接周围各点形成最小范围,如图 2加粗部分,所谓最小是因为加粗部分是以最小半径环绕住中心点。本文搜索对应点算法流程如下:

1) 以点vi为例,利用KD-Tree进行KNN搜索,得到k个二维空间欧氏距离最近点,即向量$\mathit{\boldsymbol{U}}({v_i}) = [{u_1}, ..., {u_j}, ..., {u_k}]$,各点与vi的欧氏距离组成向量$\mathit{\boldsymbol{D}}({v_i}) = [D({v_i}{u_1}), ..., D({v_i}{u_j}), ..., D({v_i}{u_k})]$,分别搜索得到D(vi)中最大值Dmax(vi),和最小值Dmin(vi),再对D(vi)进行归一化(normalization)处理,得到向量:

$ {\mathit{\boldsymbol{D}}_{\rm{N}}}({v_i}) = [{D_{\rm{N}}}({v_i}{u_1}), ..., {D_{\rm{N}}}({v_i}{u_j}), ..., {D_{\rm{N}}}({v_i}{u_k})], $

其中:

$ {D_{\rm{N}}}({v_i}{u_j}) = \frac{{{D_{{\rm{max}}}}({v_i}) - D({v_i}{u_j})}}{{{D_{{\rm{max}}}}({v_i}) - {D_{{\rm{min}}}}({v_i})}}。$

为简化计算,k往往取3~9,如图 2,搜索到abc三个最近点。

2) 在待配准网格中确定一定范围,如图 2加粗部分节点为最小范围,以vi为中心,以纵坐标等于${y_{{v_i}}}$,横坐标小于${x_{{v_i}}}$的射线为基准,按顺时针旋转,将中心点vi和依次出现的m个邻接点(adjacency)组成向量$\mathit{\boldsymbol{ad}}({v_i}) = [{v_i}, {v_{i1}}, {v_{i2}}, ..., {v_{im}}]$图 2为6个邻接点。三角剖分一个点往往邻接6到7个点,m等于6或7,如果待配准点vi与搜索到的各点邻接节点个数不同,以数量少的为准,算法都以m表示。

3) 在目标网格中确定范围,以向量U(vi)中各点为中心,如图 2中点uj为例,以纵坐标等于${y_{{u_j}}}$,横坐标小于${x_{{u_j}}}$的射线为基准,按顺时针旋转,将中心点uj和依次出现的m个邻接点组成向量$\mathit{\boldsymbol{ad}}({u_j}) = [{u_j}, {u_{j1}}, {u_{j2}}, ..., {u_{jm}}]$

4) 由于满足上文两点假设,所以有理由相信向量ad(vi)、ad(uj)中元素能够对应运算,计算ad(vi)、ad(uj)中对应元素像素灰度差值并取绝对值,得到以点vi和点uj为中心一定范围内节点像素灰度差异,组成向量:

$ \begin{array}{l} {\mathit{\boldsymbol{D}}_{\rm{I}}}{\rm{(}}{v_i}{u_j}{\rm{)}} = [{D_{\rm{I}}}{\rm{(}}{v_i}{u_j}{\rm{)}}, {D_{\rm{I}}}{\rm{(}}{v_{i1}}{u_{j1}}{\rm{)}}, ..., {D_{\rm{I}}}{\rm{(}}{v_{im}}{u_{jm}}{\rm{)}}], \\ {D_{\rm{I}}}{\rm{(}}{v_i}{u_j}{\rm{)}} = \left| {I({v_i}) - I({u_j})} \right|, \\ {D_{\rm{I}}}{\rm{(}}{v_{i1}}{u_{j1}}{\rm{)}} = \left| {I({v_{i1}}) - I({u_{j1}})} \right|, \end{array} $

${D_{\rm{I}}}{\rm{(}}{v_{im}}{u_{jm}}{\rm{)}}$同理。其中I(·)为某点的像素灰度值。然后计算这一范围内平均(Average)像素灰度差异:

$ {D_{{\rm{AvgI}}}}{\rm{(}}{v_i}{u_j}{\rm{)}} = \frac{1}{{m + 1}}\sum {{\mathit{\boldsymbol{D}}_{\rm{I}}}{\rm{(}}{v_i}{u_j}{\rm{)}}}。$

同理点vi与向量U(vi)中其他节点计算出一定范围内的平均像素灰度差异,可组成向量:

$ {\mathit{\boldsymbol{D}}_{{\rm{AvgI}}}}{\rm{(}}{v_i}{\rm{)}} = [{D_{{\rm{AvgI}}}}{\rm{(}}{v_i}{u_1}{\rm{)}}, ..., {D_{{\rm{AvgI}}}}{\rm{(}}{v_i}{u_j}{\rm{)}}, ..., {D_{{\rm{AvgI}}}}{\rm{(}}{v_i}{u_k}{\rm{)}}], $

然后对DAvgI(vi)中元素进行搜索,找到最大值DAvgImax(vi),最小值DAvgImin(vi),再对DAvgI(vi)进行归一化处理,得到:

$ {\mathit{\boldsymbol{D}}_{{\rm{NAvgI}}}}{\rm{(}}{v_i}{\rm{)}} = [{D_{{\rm{NAvgI}}}}({v_i}{u_1}), ..., {D_{{\rm{NAvgI}}}}({v_i}{u_j}), ..., {D_{{\rm{NAvgI}}}}({v_i}{u_k})], $

其中:

$ {D_{{\rm{NAvgI}}}}({v_i}{u_j}) = \frac{{{D_{{\rm{AvgImax}}}}({v_i}) - {D_{{\rm{AvgI}}}}{\rm{(}}{v_i}{u_j}{\rm{)}}}}{{{D_{{\rm{AvgImax}}}}({v_i}) - {D_{{\rm{AvgImin}}}}({v_i})}}。$

5) 利用步骤1)得到点vik个二维欧氏距离最近点归一化后向量DN(vi),利用步骤4)得到以点viU(vi)中各点为中心,一定范围内平均像素灰度差异归一化向量DNAvgI(vi)。综合考虑距离向量DN(vi)和灰度差异向量DNAvgI(vi),两组都是k维向量;并且位置对应,即DN(vi)中DN(viuj)表示点vi与点uj的归一化距离,DNAvgI(vi)中DNAvgI(viuj)表示以这两点为中心一定范围内像素灰度差异。利用当前时刻归一化互信息(normal mutual information, NMI, 用INM表示)作为权重(每次迭代INM可能会改变),进行加权计算,得到最终向量F(vi):

$ \mathit{\boldsymbol{F}}({v_i}) = {\mathit{\boldsymbol{D}}_{\rm{N}}}({v_i}) + {I_{{\rm{NM}}}} \times {\mathit{\boldsymbol{D}}_{{\rm{NAvgI}}}}{\rm{(}}{v_i}{\rm{)}}{\rm{。}} $ (5)

再搜索向量F(vi)中最大值所对应的点,不失一般性的,我们定义为点uj,即点vi考虑灰度因素后的对应点。各个待配准点都按上述算法得到其对应点,组成新的对应目标矩阵U,即式(4)。

构成式(5)基于以下原因,第一是模型需要综合考虑二维距离因素和像素灰度因素,但又不能将灰度因素放到与距离因素同等重要的地位,所以需要在灰度因素前加上权重。第二是其权重决定模型考虑多少比例的灰度因素,由于每次迭代,点云位置可能变化,考虑灰度因素的比例也不能一成不变。比如当前时刻两幅图像的互信息较小,说明此时两幅图像在像素灰度上的统计依赖性不强,信息包含程度较低,所以不能过多考虑灰度因素,而应着重考虑距离因素,反之亦然。同时需要说明的是,之所以距离向量DN(vi)中元素仅是以点vi与各搜索点的归一化距离,而非类似灰度差异向量一样考虑一定范围内节点的平均距离再归一化,原因在于,如果考虑以某点为中心的其他邻接点对应距离,再平均化处理,可能会将最重要的中心点与中心点的距离信息湮没,不利于后续计算。

配准流程:

1) 初始化X1:输入待配准图P和目标图T,计算其归一化互信息NMI,利用刚性ICP得变换矩阵XICP和粗配准图像PICP,计算PICPT归一化互信息INM-ICP

if(INM-ICP < INM)

{利用刚性互信息配准得变换矩阵XMI}X1=XMI;

else  X1=XICP;

2) 设定平滑参数α${\alpha ^i} \in \{ {\alpha ^1}, ..., {\alpha ^n}\} $, ${\alpha ^i} > {\alpha ^{i + 1}}$,标志点参数β,目标函数最小误差τ

for$({\alpha ^i} > {\alpha ^n})$

while$(||{\mathit{\boldsymbol{X}}^{j + 1}} - {\mathit{\boldsymbol{X}}^j}|| > \tau \;\;||{\mathit{\boldsymbol{X}}^j} = = {\mathit{\boldsymbol{X}}^1})$

利用本文搜索算法得对应目标点矩阵U

利用正规方程求解方程(3),解得Xj+1

end  while

αi=αi+1

end  for。

上述流程是本文提出的结合灰度信息的压敏漆图像配准算法。由于非刚性ICP算法,初始位置十分重要,先利用归一化互信息作为判断标准,使得粗配准后图像在坐标距离与像素灰度上接近目标图像,虽可能没有刚性ICP算法在坐标距离上达到更近,但为后面非刚性ICP算法提供了能兼顾两种因素的有利初始位置,同时非刚性ICP算法也是以坐标距离最近为主要目标进行迭代的,所以模型最终在距离因素上也将达到最近。改进的搜索算法考虑了像素灰度因素,但搜索所得的对应点仍以二维坐标参与后续计算,避免直接利用三维坐标时的误配准,兼顾了坐标距离与像素灰度,实现了双目标共同优化的目的。

3 实验结果与对比分析

本实验选择在Windows 7(64 bit)的MATLAB2017a上运行,CPU为core i5-3470(主频3.2 GHz),RAM为4 GB。利用活动轮廓模型提取图像模型点云,如图 5所示。

图 5 点云提取图像。(a)活动轮廓模型提取轮廓;(b)所得模型二维网格图像 Fig. 5 Point cloud extraction image. (a) Extracted contour by active contour model; (b) 2D grid image of the model

利用活动轮廓模型提取轮廓后再对图像进行整体提取,使用Delaunay三角剖分得到对应网格图像,待配准网格用红底绿边表示,目标网格用蓝底蓝边表示,图 5(b)右图存在孔洞。

本文采用配准后图像与目标图像对应标志点像素距离的均方根误差(RMSE)和两者图像归一化互信息(NMI)作为评价指标。为了实验的科学性以及体现明显畸变下的性能,本文采用两组数据,一组是由相机拍摄的0°无风与10°有风图像,并且采用拍摄中有漆料脱落的图像,加大了配准的难度(0off/10)。第二组采用将20°吹风后图像手动进行一次仿真形变,再与其20°未吹风图像进行配准(20off/20D)。

本文算法与目前配准领域较新或较经典的一些算法进行了对比分析。分别是基于特征的MSID(mid-space-independent deformable)[10]算法;基于灰度与B样条结合的MI-Bspline[18]算法;粗(coarse)和精(fine)配准均未利用NMI测度的NRICP算法[13];以及粗配准部分未利用NMI测度,但精配准部分采用像素关联搜索策略的NRICP(CnNMI)算法;粗配准部分采用了NMI测度,但精配准部分未采用本文相关策略的NRICP(FnNMI)算法。另外本文算法配准后图像,采用最邻近插值法将点云配准后出现的间隙填补。算法的配准结果如图 6图 7所示。

图 6 0off/10各算法配准后重叠图像。(a)原始重叠图像;(b) MSID配准后重叠图像;(c) MI-Bspline配准后重叠图像;(d) NRICP配准后重叠图像;(e) NRICP(CnNMI)配准后重叠图像;(f) NRICP(FnNMI)配准后重叠图像;(g)本文算法 Fig. 6 Overlapped images after registration of different algorithms in the case of 0off/10. (a) Raw overlapped images; (b) Overlapped images after MSID registration; (c) Overlapped images after MI-Bspline registration; (d) Overlapped images after NRICP registration; (e) Overlapped images after NRICP (CnNMI) registration; (f) Overlapped images after NRICP (FnNMI) registration; (g) Overlapped images based on our algorithm

图 7 20off/20D各算法配准后重叠图像。(a)原始重叠图像;(b) MSID配准后重叠图像;(c) MI-Bspline配准后重叠图像;(d) NRICP配准后重叠图像;(e) NRICP(CnNMI)配准后重叠图像;(f) NRICP(FnNMI)配准后重叠图像;(g)本文算法 Fig. 7 Overlapped images after registration of different algorithms in the case of 20off/20D. (a) Raw overlapped images; (b) Overlapped images after MSID registration; (c) Overlapped images after MI-Bspline registration; (d) Overlapped images after NRICP registration; (e) Overlapped images after NRICP (CnNMI) registration; (f) Overlapped images after NRICP (FnNMI) registration; (g) Overlapped images based on our algorithm

客观评价指标如表 1表 2所示,从视觉效果上,尤其是通过测压孔边缘配准上看,本文算法明显优于所对比的五种算法,从客观评价指标上看,本文算法也获得最好的效果。

表 1 0off/10各算法配准后指标对比表 Table 1 Comparison of indicators after registration of different algorithms in the case of 0off/10
Raw (a) MSID (b) MI-Bspline (c) NRICP (d) NRICP(CnNMI) (e) NRICP(FnNMI) (f) Ours (g)
RMSE 3.8194 3.0133 1.4077 2.7612 1.5211 1.215 1.0529
NMI 0.4721 0.5564 0.5710 0.4892 0.5587 0.5658 0.5970

表 2 20off/20D各算法配准后指标对比 Table 2 Comparison of indicators after registration of different algorithms in the case of 20off/20D
Raw (a) MSID (b) MI-Bspline (c) NRICP (d) NRICP(CnNMI) (e) NRICP(FnNMI) (f) Ours (g)
RMSE 2.8611 1.8630 1.3591 1.3482 1.2808 1.2460 1.0065
NMI 0.5508 0.5891 0.6156 0.5897 0.5926 0.5954 0.6212

表 1中,由于是0°与10°图像的配准,其像素灰度相关性要弱于20°图像,所以其NMI普遍不高。MSID算法是利用中空间特征的配准,但算法过多关注梯度方向信息,效果不是很理想。基于互信息与B样条的配准算法,其目标函数以图像互信息最大进行迭代,利用B样条的约束,考虑了部分坐标间关系,是一种以像素灰度因素为主,坐标关系为辅的非刚性配准算法,对于压敏漆图像特点而言并不太适用。粗和精配准均未采用NMI测度的NRICP算法实验效果要差于单一利用NMI测度的NRICP(CnNMI)或NRICP(FnNMI)算法。NRICP(CnNMI)算法因为在精配准过程中考虑灰度信息,所得效果有所提升,但由于在粗配准过程中未考虑,其效果要略差于NRICP(FnNMI)算法,同时也说明非刚性ICP算法初始位置十分重要,算法流程中加入NMI判断很有必要。NRICP(FnNMI)只在粗配准中考虑灰度因素,效果仍有待提升。而本文算法兼顾坐标关系与灰度因素,每次迭代都会根据当前状况去调整灰度与坐标的关系,找到合理的对应点,双目标共同优化实现配准。本文算法相对于次优的NRICP(FnNMI)算法,RMSE提升超过15%,NMI提升在5%左右。

4 结论

基于压敏漆图像配准,本文创新性地运用了非刚性点云配准技术。同时改进了非刚性ICP的搜索算法,提出了结合灰度信息的压敏漆图像配准算法。新算法既考虑了空间坐标因素又考虑了像素灰度因素,实验证明,所提算法具有不错表现。

本文算法在考虑像素灰度因素时仅利用图像整体互信息,并未对不同区域分别处理。同时模型的第二惩罚项需要人工选取标志点。为进一步提高配准精度与自动化程度,后续研究主要就是对不同区域测度加以区分处理,并找到科学合理的局部特征描述子。

参考文献
[1]
Friedl F, Krah N, Jähne B. Optical sensing of oxygen using a modified stern-volmer equation for high laser irradiance[J]. Sensors and Actuators B: Chemical, 2015, 206: 336-342. [Crossref]
[2]
Peng D, Chen J W, Jiao L R, et al. A fast-responding semi-transparent pressure-sensitive paint based on through-hole anodized Aluminum oxide membrane[J]. Sensors and Actuators A: Physical, 2018, 274: 10-18. [Crossref]
[3]
Kontis K. A review of some current research on pressure sensitive paint and thermographic phosphor techniques[J]. The Aeronautical Journal, 2007, 111(1122): 495-508. [Crossref]
[4]
Park S H, Sung H J. Correlation-based image registration for applications using pressure-sensitive paint[J]. AIAA Journal, 2005, 43(3): 472-478. [Crossref]
[5]
Fujimatsu N, Tamura Y, Fujii K. Improvement of noise filtering and image registration methods for the pressure sensitive paint experiments[J]. Journal of Visualization, 2005, 8(3): 225-233. [Crossref]
[6]
Chen C, Chen Y, Huang Q. Study of image geometric distortion correction for PSP measurements[J]. Computer Engineering and Design, 2010, 31(24): 5298-5301.
陈畅, 陈勇, 黄琦. 压敏漆技术中的图像几何畸变校正研究[J]. 计算机工程与设计, 2010, 31(24): 5298-5301 [Crossref]
[7]
Gao L M, Wei N, Gao J, et al. Image processing of PSP technique in the internal flow[J]. Journal of Experiments in Fluid Mechanics, 2013, 27(1): 93-97.
高丽敏, 韦楠, 高杰, 等. 基于内流场PSP测量技术的图像后处理[J]. 实验流体力学, 2013, 27(1): 93-97 [Crossref]
[8]
Xiong H L, Xiang X J, Lang W D. Image data processing techniques for pressure sensitive paint[J]. Journal of Experiments in Fluid Mechanics, 2011, 25(2): 77-82.
熊红亮, 向星居, 郎卫东. 压敏漆图像数据处理技术[J]. 实验流体力学, 2011, 25(2): 77-82 [Crossref]
[9]
Lu X S, Zhang S, Yang W, et al. SIFT and shape information incorporated into fluid model for non-rigid registration of ultrasound images[J]. Computer Methods and Programs in Biomedicine, 2010, 100(2): 123-131. [Crossref]
[10]
Aganj I, Iglesias E I, Reuter M, et al. Mid-space-independent deformable image registration[J]. Neuroimage, 2017, 152: 158-170. [Crossref]
[11]
Besl P J, McKay N D. A method for registration of 3-D shapes[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1992, 14(2): 239-256. [Crossref]
[12]
Chi Y, Yu X Q, Luo Z Y. 3D point cloud matching based on principal component analysis and iterative closest point algorithm[C]//Proceedings of 2016 International Conference on Audio, Language and Image Processing, Shanghai, China, 2016: 404-408.
[13]
Amberg B, Romdhani S, Vetter T. Optimal step nonrigid ICP algorithms for surface registration[C]//Proceedings of 2007 IEEE Conference on Computer Vision and Pattern Recognition, Minneapolis, MN, USA, 2007: 1-8.
[14]
Hontani H, Matsuno T, Sawada Y. Robust nonrigid ICP using Outlier-Sparsity regularization[C]//Proceedings of 2012 IEEE Conference on Computer Vision and Pattern Recognition, Providence, RI, USA, 2012: 174-181.
[15]
Cheng S Y, Marras I, Zafeiriou S, et al. Statistical non-rigid ICP algorithm and its application to 3D face alignment[J]. Image and Vision Computing, 2017, 58: 3-12. [Crossref]
[16]
Chen J, Ma J Y, Yang C C, et al. Non-rigid point set registration via coherent spatial mapping[J]. Signal Processing, 2015, 106: 62-72. [Crossref]
[17]
Lu M, Zhao J, Guo Y L, et al. Accelerated coherent point drift for automatic Three-Dimensional point cloud registration[J]. IEEE Geoscience and Remote Sensing Letters, 2016, 13(2): 162-166. [Crossref]
[18]
Rueckert D, Sonoda L I, Hayes C, et al. Nonrigid registration using free-form deformations: application to breast MR images[J]. IEEE Transactions on Medical Imaging, 1999, 18(8): 712-721. [Crossref]
[19]
Woo J, Stone M, Prince J L. Multimodal registration via mutual information incorporating geometric and spatial context[J]. IEEE Transactions on Image Processing, 2015, 24(2): 757-769. [Crossref]
[20]
Qiu Z Y, Tang H H, Tian D S. Non-rigid medical image registration based on the thin-plate spline algorithm[C]//Proceedings of 2009 WRI World Congress on Computer Science and Information Engineering, Los Angeles, CA, USA, 2009: 522-527.
[21]
Chan T F, Sandberg B Y, Vese L A. Active contours without edges for vector-valued images[J]. Journal of Visual Communication and Image Representation, 2000, 11(2): 130-141. [Crossref]
[22]
Xu L, Wan J W L, Bian T T. A continuous method for reducing interpolation artifacts in mutual Information-Based rigid image registration[J]. IEEE Transactions on Image Processing, 2013, 22(8): 2995-3007. [Crossref]
[23]
Collignon A. Automated multimodality image registration using information theory[C]//Proceedings of International Conference Information Processing in Medical Imaging, Ile De Berder, France, 1995: 263-274.
[24]
Viola P, Wells III W M. Alignment by maximization of mutual information[J]. International Journal of Computer Vision, 1997, 24(2): 137-154. [Crossref]
[25]
Studholme C, Hill D L G, Hawkes D J. An overlap invariant entropy measure of 3D medical image alignment[J]. Pattern Recognition, 1999, 32(1): 71-86. [Crossref]
[26]
Friedman J H, Bentley J L, Finkel R A. An algorithm for finding best matches in logarithmic expected time[J]. ACM Transactions on Mathematical Software, 1977, 3(3): 209-226. [Crossref]