声纳图像配准中的Demons算法
王达，卞红雨
(哈尔滨工程大学 水声技术重点实验室，黑龙江 哈尔滨，150001)
摘要：针对传统基于光流场模型的Demons算法变形方向不易确定，同时考虑声纳图像序列中连续帧图像缺乏梯度信息，提出一种结合梯度互信息的改进Demons算法。该方法在原有图像变形力的基础上，增加两幅图像间梯度互信息作为驱动图像变形的附加力。在图像配准的同时，使两幅图像的梯度互信息达到最大，避免Demons算法仅依靠图像灰度梯度变形，从而得到更为精确的配准变换。通过声纳图像配准实验，研究结果表明：该算法在互信息量上提高5%以上，验证该方法具有配准精度高，鲁棒性强，是一种有效的非刚性配准方法。
关键词：图像配准；互信息；图像梯度；声纳图像；Demons算法
中图分类号：TP317.4 文献标志码：A 文章编号：1672-7207(2013)10-4074-07
A variational Demons algorithm on sonar image registration
WANG Da, BIAN Hongyu
(Science and Technology on Underwater Acoustic Laboratory, Harbin Engineering University, Harbin 150001, China)
Abstract: Considering that a Demons algorithm based on the model of the optical flow field deformation direction can not be determined, and the sonar image sequence lack of gradient information, an improved “Demons” algorithm was proposed. On the basis of the original image deformation force, the additional force of the mutual information between two images gradient was increased as the driving image deformation. The maximum gradient mutual information of the two images was obtained at the same time of the image registration, preventing the Demons algorithm from relying solely on the image gray gradient deformation, thus resulting in a more accurate registration transformation. Experiments of sonar image are compared with classic Demons algorithm, and the results show that the proposed algorithm increases more than 5% on the mutual information. It has high accuracy, robustness, and is an effective method of non-rigid registration.
Key words: image registration; mutual information; image gradient; sonar images; Demons algorithm
声纳图像配准是声纳图像融合的基础，同时也是声纳图像处理领域的一个重要分支。配准的目的是使待配准图像之间的对应点达到空间上的一致。图像的刚性配准只能对整幅图像进行全局变换，当需要对图像进行精确局部变形时，如配准同一目标图像序列中连续帧的声纳图像等，则需要对图像进行非刚性配准。由于声纳图像具有成像物体边缘弱化，图像对比度低、强度非均匀、分辨率不高等特点[1-2]，给配准过程带来很大的困难，较光学图像更易产生误配准的现象。目前，非刚性图像配准方法主要可以分为基于图像共有特征的配准方法和基于像素点的配准方法[3-4]。基于像素点的配准方法由于不受特征提取质量的影响，以其精度高、鲁棒性强、自动化程度高等特点已经成为图像配准过程中的重要方法。其中基于光流场模型的Demons算法[5-6]是一种有效的基于灰度驱动的非刚性配准方法。该模型表达了图像的变化，由于其包含了目标运动的信息，可以补偿时序图像两幅连续帧之间物体和视点的相对运动。经典Demons算法采用参考图像的灰度梯度信息来决定浮动图像中每个像素的移动方向。Wang等[7-8]提出了一种改进的“Active Demons”算法，将浮动图像的梯度信息也加入到光流场方程中，实现了更准确和高效的配准。但在声纳图像处理领域中，灰度梯度不明显是普遍存在的，从而造成图像变形力方向不能确定，容易导致错误的配准结果。Vercauteren等[9-11]提出了一种将图像灰度均方差作为相似性测度的配准模型，通过对目标函数进行优化寻找出一个最优的变形场，但是没有考虑到随着图像配准过程的进行，待配准图像之间互信息对配准过程的影响。考虑到最大互信息配准方法在刚性配准中的成功应用[12-14]，本文作者提出了一种基于互信息的改进Demons算法，使声纳图像在缺乏梯度信息时能够正确的变形。经声纳图像配准实验验证了本文所提方法能够产生精确的配准变换。
1 Demons算法
经典Demons算法是由Thirion提出的，即“Demons-base”算法[5]，在概念上，它和19世纪Maxwell的实验原理很相似。该算法的基本思想是：假设运动目标的灰度不随时间改变，那么非刚性配准可以看作是浮动图像中各个像素向参考图像逐步扩散的过程，浮动图像各个像素的扩散速度由参考图像的灰度梯度信息决定。
假设参考图像r，浮动图像f都是图像函数I(x(t)，y(t)，t)某一时刻的“快照”，图像函数的灰度值保持常数，即
(1)
(2)
(3)
在初始时刻t0，图像函数I等于f，经过一定时间到达t1时刻后，图像函数被完全变形为r。图像配准的过程就是要得出一个能驱动f中的各个像素点向r中对应像素点移动的向量场。为了得出驱动力，将式(1)偏微分，得到
(4)
式(4)可以简化为：
(5)
其中：
为从浮动图像向参考图像变化的速度场；
为r的梯度向量。
进一步可以得到
(6)
但是，当||||→0时，式(6)的定义变得很不稳定，导致较大的扩散速度v，为了避免这一问题，在分母增加一个分量，得到
(7)
式(7)使得图像的形变得到了一定程度的控制。但是该方法仅采用了参考图像的梯度信息来驱动图像的形变，当参考图像梯度信息不明显时，容易产生误配准现象。因此，Wang等[7-8]提出了将浮动图像的梯度信息作为一种正内力，参考图像的梯度信息作为负内力，利用这2种力同时驱动形变，得到的扩散速度为
(8)
由||
||2+α2(r-f)2≥2α||||·||r-f||可知：通过调整α可以控制扩散速度的大小。
由于Demons算法是利用图像局部信息来驱动形变的，为了在全局范围内使形变连续，从而保证图像的拓扑结构，通常的做法是在每一次迭代后，使用高斯滤波来平滑所得到的偏移，其扩散速度定义为
(9)
其中：Gσ为高斯滤波器。
在Demons模型中，浮动图像各个像素都可以自由移动，使得在浮动图像中具有某个特定灰度值的所有像素有可能映射到参考图像上的同一像素点，从而改变了图像的拓扑结构，导致错误的配准结果，为了有效解决这一问题，本文结合了梯度互信息对Demons模型进行改进。
2 基于梯度互信息的Demons方法
在Demons算法中，当||||为0时，像素点的扩散速度不能确定，为了使图像像素能够沿着正确的方向移动，本文提出了一种基于梯度互信息的改进Demons算法，该算法考虑到配准结束时，配准结果之间的梯度互信息最大这一特点，在原有扩散速度的基础上，增加了2幅图像之间梯度互信息的作用，使两幅图像配准结束的同时梯度互信息也达到最大。改进后的Demons算法扩散速度模型为
(10)
其中：Max(IGMI(vn))为2幅图像的梯度互信息；β为正常数，表示此项的权重。
2.1 梯度互信息相似性的计算方法
Pluim方法[15]比传统没有考虑空间信息的互信息配准方法效果要好，但是在计算梯度项时，忽略了噪声和梯度模值较小的像素对梯度项的影响，因此本文在Pluim方法的基础上，针对以上两方面问题对其进行了改进，提出了新的梯度方向角和梯度模值的计算方法。
2.1.1 改进的梯度项
(1) 改进的梯度方向角。由于梯度方法本身容易受到噪声的影响，因此，本文考虑先对梯度幅值做阈值处理，若像素点的梯度幅值大于阈值，则认为该像素点处于灰度变化剧烈的边缘区域；否则，则认为该像素点处于灰度变化不剧烈的非边缘区域或者属于噪声点。其中阈值的定义为
(11)
其中：M和N为参考图像大小；
为参考图像中像素点(x，y)梯度的模值。
为了克服噪声的影响，将小于阈值的梯度方向角设为0，具体表达为
(12)
与Pluim方法相同，对应像素点的梯度权函数定义为
(13)
由式(13)可以看出：梯度权函数越接近最大值1时，梯度方向角越接近于0或π，待处理的两幅图像越接近于完全配准。
(2) 改进的梯度模值。在参考图像和浮动图像中，对应像素点的梯度模值越相近，则像素点对在模值上越相似，因此可以通过对应像素点对梯度模值的比值来衡量像素点对的相似性，定义为
(14)
其中：当max(，
)=0时，意味着对应像素点对的梯度模值相等，并且等于0，表明这对像素点在模值上最相似，但出于比值的考虑，将其设置成1。
2.1.2 结合互信息的联合相似性测度
对于参考图像和浮动图像，分别计算各对应像素点的改进梯度权函数和改进梯度模值，则新的梯度项定义为
(15)
此外，为了消除重叠区域大小对计算互信息量的影响，将改进后的梯度项系数进行归一化，定义为
(16)
其中：count(r,f)为r；f重叠区域像素的总个数。
最后，将改进以后的归一化梯度项与互信息相结合，得出相似性测度。
(17)
将IGMI作为图像配准中新的相似性测度，利用适当的空间变换和寻优策略，对声纳图像进行配准。
2.1.3 基于梯度互信息的Demons算法的分析
基于图像灰度互信息的配准方法，最大的不足是忽略了图像的空间信息。香农熵的定义是建立在一个独立连续信号的基础上，但是这种假设在声纳图像中一般是不成立的。通过熵的定义可知：香农熵的大小并不依赖于灰度本身，而是取决于图像中灰度出现的概率。显然变换(包括平移、旋转等)只是改变了像素点的空间位置，并没有对像素点的灰度值产生影响，也就没有改变灰度出现的概率，进而不会影响香农熵的大小。这样容易造成作为图像配准的相似性测度，在整个搜索空间中，同时在2个变换位置出现峰值，因此将空间信息引入到互信息配准中是有必要的。
在声纳图像中，梯度较大的像素点通常是组成目标边缘的像素点，也是包含信息较丰富的地方。就声纳图像配准过程而言，对相同目标所成图像可能具有不同的灰度，不同的分辨率，甚至图像本身大小也会不同，但是对于同一目标，其边界是确定的，不会发生很明显的变化；当两幅图像配准时，对应位置像素点的梯度矢量将会具有相同或者相反的指向[15]。
基于以上分析，本文提出的梯度互信息采用Parzen窗方法[16]计算。设r为参考图像，f为浮动图像，g(x，y；μ)为形变，其中参数μ=(μ1，μ2，…，μN)，N为像素数目。同时，假设Lr和Lf分别为参考图像和浮动图像的灰度级，ω为离散Parzen窗，通过3次B样条函数定义。
因此，联合概率密度定义为
(18)
其中：lf∈Lf；lr∈Lr；f(g(x,y;μ))和r(x,y)为灰度强度；△bf和△br分别服从Lf和Lr，以约束灰度映射在窗口函数范围内；α为系数。
根据互信息理论[17-18]，本文提出的梯度互信息可以表示为
(19)
其中：p(lf，lr；μ)为联合概率密度；pf(lf；μ)和pr(lr；μ)分别为浮动图像和参考图像的边缘概率密度。
梯度互信息第i个参数μ的偏微分为
(20)
其中：
p(lf，lr；μ)/ μi满足如下约束条件
(21)
其中：N为像素数；ω(ξ)/ ξ为Parzen窗函数的派生形式；f(t)/ t为浮动图像梯度；g(x,y; μ)/ μi实际可以看作为雅可比行列式，在Demons算法中，g(x，y；μ)/ μi=1。
最后，形变可以表示为
(22)
2.2 算法流程
具体的算法步骤如下。
Step 1：将浮动图像每一位置(x，y)像素点i的初始扩散速度vi设置为对图像进行刚性配准得到的初始偏移量。刚性配准过程中可以采用本文提出的梯度互信息作为相似性测度。
Step 2：在第n次迭代下，对于每一个浮动图像像素点，其位置向量为(xin，yin)=(xi，yi)+vin-1，扩散速度vin通过式(22)计算。
Step 3：对vin进行评价，若两幅图像之间的互信息变化量小于某个预先设定的阈值△，则认为迭代已经收敛，转入Step 4；否则，继续进行n+1次迭代。
Step 4：当迭代结束后，用所得到的扩散速度v作用于浮动图像，使用双线性插值得到配准后的图像。
3 实验结果与分析
为了验证本文所提配准方法的有效性，对几组不同声纳序列图像中连续帧图像进行实验，实验平台为处理器Pentium Dual-Core 2.6 GHz，内存为2 GB，程序在Windows操作系统下，使用Matlab 7.10语言编写。
首先对如图1所示的海底沉物钢板图像进行配准实验，图1(a)，图1(b)和图1(c)所示分别为250×390像素的78帧连续声纳图像中的第1帧、第30帧和第50帧图像。其中将第30帧图像作为参考图像，将第1帧和第50帧图像分别作为浮动图像，使用经典Demons算法和本文提出的算法进行配准，得到的配准结果分别为图1(d)，图1(e)和图1(f)，图1(g)。图1(h)和图1(i)所示为使用本文算法得到的扩散速度向量场。
由图1可以看出采用本文算法的配准结果较为准确。在采用经典Demons算法中，将浮动图像的每一个像素点全部作为Demons点，这样从浮动图像到参考图像的形变过程演变为完全自由形变(free-form)。由于自由形变本身的特点，使浮动图像中的每个像素点都可以移动到图像的任何位置，这就造成了在配准过程中原来连续的结构变得不连续或者扭曲。尽管在经典Demons算法中采用高斯滤波以保证拓扑结构的不变性，显然在声纳图像配准中没有发挥良好的效果，而本文提出的算法很好地克服了这一缺点，在声纳图像序列连续帧配准问题中取得了良好的效果。
图1 海底沉物钢板配准
Fig. 1 Steel plate registration
图2(a)，图2(b)和图2(c)所示分别为钢板焊缝49帧连续声纳图像中第10帧、第30帧和第49帧图像，图像均为250×415像素。将第30帧图像作为参考图像，第10帧和第49帧图像作为浮动图像，同样采用经典Demons算法和本文提出的算法进行配准，得到的配准图像分别为图2(d)，图2(e)和图2(f)，图2(g)。
通过对比图2的配准结果，可以看出本文提出的改进Demons算法在钢板焊缝图像中也取得了良好的效果。
图2 钢板焊缝配准
Fig. 2 Weld sonar image registration
根据式(10)并结合实验发现，对Demons算法主要4个参数进行不同的组合，配准结果会有明显的不同。
当α较小时，扩散速度过大，可能影响拓扑结构；反之，扩散速度过小，驱动浮动图像形变的能力变弱，收敛速度变慢。再综合考虑高斯滤波器尺寸h和标准差σ的影响，得出参数取值趋势是：当α，h和σ取值过小时，图像拓扑结构损坏严重；当α，h和σ取值过大时，配准结构没有明显改善，但程序运行时间加大。经过反复实验得出最佳经验值，即海底沉物钢板配准实验中α=0.8，h=70，σ=15；钢板焊缝实验中α=1，h=70，σ=10。
表1和表2所示分别为图1和图2配准实验中经典Demons算法和本文提出的改进算法的性能比较。
从表1和表2可以看出：改进后的Demons算法配准结果之间的互信息都有明显的提高，但是程序运行时间开销较经典Demons算法加大许多，经过分析得知，这主要是因为在进行梯度互信息计算时耗费了大量的运行时间。
表1 海底沉物钢板配准结果比较
Table 1 Steel plate registration results compare
表2 钢板焊缝配准结果比较
Table 2 Weld sonar image registration results compare
4 结论
(1) 基于Demons模型的算法是一种图像灰度的非刚性配准算法，但是当待配准图像灰度梯度信息很小甚至为零时，浮动图像中像素点的扩散速度不能确定，容易造成误配准；而且经典Demons算法在声纳图像序列连续帧配准过程中效果很不理想。
(2) 本文考虑到梯度互信息在刚性图像配准中的成功应用，针对梯度方向角和梯度模值两方面进行了改进，有效克服了噪声和梯度变化剧烈等问题对梯度互信息的影响，最后结合Demons算法，将待配准图像间的梯度互信息作为驱动图像变形的附加力，从而使得两幅图像配准的同时，梯度互信息达到最大，有效克服了在灰度梯度不明显时，图像仅靠灰度梯度信息驱动变形容易造成误配准的问题，使声纳图像序列在连续帧配准过程中能够得到正确的配准变换。
(3) 实验结果对比验证了改进的Demons算法对声纳图像配准非常有效，是一种准确，鲁棒的非刚性配准算法。此外，由于声纳图像自身的特点造成算法时间开销较大，如何在保证配准精度和鲁棒性的情况下，缩短算法时间的开销，是需要进一步研究和解决的问题。
参考文献：
[1] Waite A D. 实用声纳工程[M]. 王德石, 译. 北京：电子工业出版社, 2004: 79-97.
Waite A D. Sonar for practising engineers[M]. WANG Deshi, trans. Beijing: Publishing House of Electronics Industry, 2004: 79-97.
[2] 沈郑燕. 声纳图像去噪与分割技术研究[D]. 哈尔滨: 哈尔滨工程大学水声工程学院, 2010: 18-81.
SHEN Zhengyan. Research on techniques of sonar image de-noising and segmentation[D]. Harbin: Harbin Engineering University. College of Underwater Acoustic Engineering, 2010: 18-81.
[3] 林相波, 邱天爽, 阮素, 等. Demons非刚性配准算法拓扑保持性的研究[J]. 自动化学报, 2010, 36(1): 179-183.
LIN Xiangbo, QIU Tianshuang, RUAN Su, et a l. Research on the topology preservation of the Demons non-rigid registration algorithm[J]. Acta Automatica Sinica, 2010, 36(1): 179-183.
[4] 韩雨, 王卫卫, 冯象初. 基于迭代重加权的非刚性图像配准[J]. 自动化学报, 2011, 37(9): 1059-1066.
HAN Yu, WANG Weiwei, FENG Xiangchu. Iteratively reweighted method based nonrigid image registration[J]. Acta Automatica Sinica, 2011, 37(9): 1059-1066.
[5] Thirion J P. Image matching as a diffusion process: An analogy with Maxwell’s Demons[J]. Medical Image Analysis, 1998, 2(3): 243-260.
[6] 江济良, 屠大维, 周许超, 等. 复杂光流场运动分析与特征提取[J]. 电子测量与仪器学报, 2011, 25(3): 285-290.
JIANG Jiliang, TU Dawei, ZHOU Xuchao, et al. Motion analysis of complex optical flow field and feature extraction[J]. Journal of Electronic Measurement and Instrument, 2011, 25(3): 285-290.
[7] Wang H, Dong L, O’Daniel J, et al. Validation of an accelerated ‘Demons’ algorithm for deformable image registration in radiation therapy[J]. Physics in Medicine and Biology, 2005, 50(12): 2887-2905.
[8] Gu X, Pan H, Liang Y, et al. Implementation and evaluation of various Demons deformable image registration algorithms on a GPU[J]. Physics in Medicine and Biology, 2010, 55(1): 207-219.
[9] Vercauteren T, Pennec X, Perchant A, et al. Diffeomorphic Demons: Efficient non-parametric image registration[J]. NeuroImage, 2009, 45(1): 61-72.
[10] Vercauteren T, Pennec X, Malis E, et a l. Insight into efficient image registration techniques and the Demons algorithm[J]. Information Processing in Medical Imaging, 2007, 16(2): 495-506.
[11] Kroon D J, Slump C H. MRI modality transformation in demon registration[C]//Proceedings of the IEEE International Symposium on Biomedical Imaging. Boston, USA: IEEE, 2009: 963-966.
[12] 朱冰莲, 田学隆, 宋维杰. 基于人工免疫系统的医学图像配准[J]. 仪器仪表学报, 2009, 30(7): 1416-1419.
ZHU Binglian, TIAN Xuelong, SONG Weijie. Medical image registration algorithm based on artificial immune system[J]. Chinese Journal of Scientific Instrument, 2009, 30(7): 1416-1419.
[13] 刘兆英, 周付根, 白相志, 等. 基于感兴趣区域互信息的多模图像配准方法[J]. 航空兵器, 2011(4): 7-12.
LIU Zhaoying, ZHOU Fugen, BAI Xiangzhi, et al. Multi-mode image registration based on mutual information of region of interests[J]. Aero Weaponry, 2011(4): 7-12.
[14] Anthony A. Mutual information for image registration[J]. Computer Technology and Application, 2011, 2(1): 9-14.
[15] Pluim J P W, Maintz J B A, Viergever M A. Image registration by maximization of combined mutual information and gradient information[J]. IEEE Transactions on Medical Imaging, 2000, 19(8): 809-814.
[16] Thevenaz P, Unser M. Optimization of mutual information for mutual information for multiresolution image registration[J]. IEEE Transactions on Image Processing, 2000, 9(12): 2083-2099.
[17] Maes F, Collignon A, Vandermeulen D, et al. Multimodality image registration by maximization of mutual information[J]. IEEE Transactions on Medical Imaging, 1997, 16(2): 187-198.
[18] Tsimpiris A, Vlachos L, Kugiumtzis D. Nearest neighbor estimate of conditional mutual information in feature selection[J]. Expert Systems with Applications, 2012, 39(16): 12697-12708.
(编辑 杨幼平)
收稿日期：2012-09-18；修回日期：2012-12-05
基金项目：中央高校基本科研业务费资助项目(HEUCF120501)
通信作者：王达(1985-)，男，天津人，博士研究生，从事水下图像处理研究；电话：0451-82589511；E-mail：wangda19600913@163.com
