工程科学学报，第41卷，第7期：955-960.2019年7月 Chinese Journal of Engineering,Vol.41,No.7:955-960,July 2019 DOI:10.13374/j.issn2095-9389.2019.07.015;http://journals.ustb.edu.cn 一种改进的非刚性图像配准算法 何凯，魏颖，王阳，黄婉蓉 天津大学电气自动化与信息工程学院，天津300072 区通信作者，E-mail:hekai(@tju.edu.cn 摘要非刚性图像配准一直是计算机视觉领域的研究重点.为解决上述问题，提出一种改进的光流场模型算法，以提高光 流估计的准确度.算法首先对原始变分光流模型进行了改进，提出利用新的各向异性正则项来代替原来的同向扩散函数，以 避免图像模糊，保留图像的边缘特征与细节特征：此外，通过引入包含邻域信息的非局部平滑项来去除光流噪点，同时增加了 一个结合图像结构与光流运动信息的权函数，以减少过平滑所造成的细节丢失，提高算法的鲁棒性.最后，利用交替最小化与 金字塔分层迭代策略相结合的方法求解位移场，实现非刚性图像的自动配准.仿真实验结果表明，与传统方法相比，本文算法 对不同类型的非刚性图像均具有较高的鲁棒性，取得了理想的图像配准效果 关键词非刚性图像配准；光流场模型；各向异性正则项；非局部平滑项：交替最小化 分类号TP391.41 An improved non-rigid image registration approach HE Kai,WEI Ying,WANG Yang,HUANG Wan-rong School of Electrical and Information Engineering,Tianjin University,Tianjin 300072,China XCorresponding author,E-mail:hekai@tju.edu.cn ABSTRACT With the rapid development of image registration technology,it is being widely used in the fields of medical image pro- cessing,remote sensing image analysis,computer vision,and others.Image registration involves two or more images that contain the same object that are obtained under different conditions.Geometric mapping between images is realized by spatial geometric transforma- tion,so that the points in one image can be related to their corresponding points in the other.Compared with rigid transformations, non-rigid transformations usually have severe local distortions and obvious nonlinear characteristics.So,it is difficult to describe non- rigid transformations using a unified transformation model.For this reason,non-rigid image registration has always been an issue and a source of difficulty in the field of computer vision.To solve this problem,an improved optical-flow-model algorithm was proposed to more accurately estimate the optical flow field.First,the original variational optical flow model was improved.To prevent blurring and preserve the edge and detail features of images,a new anisotropic regular term was proposed to replace the original homologous diffusion term.Then,to remove optical flow outliers,a non-local smoothness term was introduced that contained neighborhood information.Mo- reover,a weight function that combines image-structure and optical-flow information was added to reduce the loss of detail caused by over-smoothing and to improve robustness.Finally,to solve the displacement field and realize the automatic registration of non-rigid im- ages,an alternating minimization method and pyramid hierarchical iteration strategy were utilized.To verify the effectiveness of the pro- posed algorithm,subjective and objective evaluation values such as the peak signal-to-noise ratio PSNR)and normalized mutual infor- mation (NMI)were adopted to analyze the registration results.Compared with state-of-the-art methods,experimental results reveal the robustness and ideal registration effects of the proposed method on different types of non-rigid images. KEY WORDS non-rigid image registration;optical flow model;anisotropic regularization term;non-local smoothness term;alternate minimization 收稿日期：2018-05-16 基金项目：国家自然科学基金资助项目(61271326)
工程科学学报,第 41 卷,第 7 期:955鄄鄄960,2019 年 7 月 Chinese Journal of Engineering, Vol. 41, No. 7: 955鄄鄄960, July 2019 DOI: 10. 13374 / j. issn2095鄄鄄9389. 2019. 07. 015; http: / / journals. ustb. edu. cn 一种改进的非刚性图像配准算法 何 凯苣 , 魏 颖, 王 阳, 黄婉蓉 天津大学电气自动化与信息工程学院, 天津 300072 苣通信作者, E鄄mail: hekai@ tju. edu. cn 摘 要 非刚性图像配准一直是计算机视觉领域的研究重点. 为解决上述问题,提出一种改进的光流场模型算法,以提高光 流估计的准确度. 算法首先对原始变分光流模型进行了改进,提出利用新的各向异性正则项来代替原来的同向扩散函数,以 避免图像模糊,保留图像的边缘特征与细节特征;此外,通过引入包含邻域信息的非局部平滑项来去除光流噪点,同时增加了 一个结合图像结构与光流运动信息的权函数,以减少过平滑所造成的细节丢失,提高算法的鲁棒性. 最后,利用交替最小化与 金字塔分层迭代策略相结合的方法求解位移场,实现非刚性图像的自动配准. 仿真实验结果表明,与传统方法相比,本文算法 对不同类型的非刚性图像均具有较高的鲁棒性,取得了理想的图像配准效果. 关键词 非刚性图像配准; 光流场模型; 各向异性正则项; 非局部平滑项; 交替最小化 分类号 TP391郾 41 收稿日期: 2018鄄鄄05鄄鄄16 基金项目: 国家自然科学基金资助项目(61271326) An improved non鄄rigid image registration approach HE Kai 苣 , WEI Ying, WANG Yang, HUANG Wan鄄rong School of Electrical and Information Engineering, Tianjin University, Tianjin 300072, China 苣Corresponding author, E鄄mail: hekai@ tju. edu. cn ABSTRACT With the rapid development of image registration technology, it is being widely used in the fields of medical image pro鄄 cessing, remote sensing image analysis, computer vision, and others. Image registration involves two or more images that contain the same object that are obtained under different conditions. Geometric mapping between images is realized by spatial geometric transforma鄄 tion, so that the points in one image can be related to their corresponding points in the other. Compared with rigid transformations, non鄄rigid transformations usually have severe local distortions and obvious nonlinear characteristics. So, it is difficult to describe non鄄 rigid transformations using a unified transformation model. For this reason, non鄄rigid image registration has always been an issue and a source of difficulty in the field of computer vision. To solve this problem, an improved optical鄄flow鄄model algorithm was proposed to more accurately estimate the optical flow field. First, the original variational optical flow model was improved. To prevent blurring and preserve the edge and detail features of images, a new anisotropic regular term was proposed to replace the original homologous diffusion term. Then, to remove optical flow outliers, a non鄄local smoothness term was introduced that contained neighborhood information. Mo鄄 reover, a weight function that combines image鄄structure and optical鄄flow information was added to reduce the loss of detail caused by over鄄smoothing and to improve robustness. Finally, to solve the displacement field and realize the automatic registration of non鄄rigid im鄄 ages, an alternating minimization method and pyramid hierarchical iteration strategy were utilized. To verify the effectiveness of the pro鄄 posed algorithm, subjective and objective evaluation values such as the peak signal鄄to鄄noise ratio (PSNR) and normalized mutual infor鄄 mation (NMI) were adopted to analyze the registration results. Compared with state鄄of鄄the鄄art methods, experimental results reveal the robustness and ideal registration effects of the proposed method on different types of non鄄rigid images. KEY WORDS non鄄rigid image registration; optical flow model; anisotropic regularization term; non鄄local smoothness term; alternate minimization
956. 工程科学学报，第41卷，第7期 图像配准)是针对含有同一对象于不同条件 于以下两个基本假设：1)瞬时灰度值保持不变，即 下获得的2幅或多幅图像，通过空间几何变换来实 灰度不变原理：2)光流的变化是连续平滑的：基于 现图像之间的几何映射，使得不同图像的对应点在 上述假设，在光流约束方程的基础上增加一个速度 空间位置上保持一致.该技术在医学图像处理、遥 平滑约束，将光流场的求解转换成一个变分问题. 感图像分析和计算机视觉等领域都有着广泛的应 能量函数定义如下： 用2-).图像配准大致可以分成刚性图像配准与非 刚性图像配准2大类.其中，非刚性图像配准以其 E(,)=J[(1,u+1,+l,)2+ 局部性、非线性、不满足统一的变换模型等性质，一 a(IVul2+IVvl2)]dxdy (1) 直是计算机视觉领域的研究难点与热点，目前常用 其中，2代表空间图像域：“和：分别为光流水平方 的非刚性图像配准方法主要包括：正则化、变分、光 向与垂直方向的速度分量；I,I,和I,分别为光流在 流优化等. x,y和13个方向上的梯度；V=(8,a,)是空间梯 基于光流场模型的方法是实现非刚性图像配准 度算子；α为权重系数 的一种典型方法.光流场指的是图像灰度模式的表 式(1)可以简写为： 观运动，它包含了各像素的瞬时运动信息.Thii- E(u,v)=Ep(u,v)+aEs(u,v) (2) oms)最早引入光流方法到图像配准领域.Hom和 其中，E。(u,)称为数据项；Es(u,v)称为正则项；a Schunck[6于1981年首先把速度场与图像灰度相结 为正则项的权重系数，其数值越大，对图像的约束能 合，提出了光流约束方程，给出了计算光流的基本方 力越强 法.在此基础上，大量光流方法不断涌现，例如：Na- gel与Enkelmann]提出各项异性正则项，改善了同 2本文算法 向正则项的过平滑问题.Brox与Malik[s]提出包括 研究表明，可以通过最小化能量函数来求解图 灰度一致与梯度一致的数据项，提高了光流模型对 像之间的运动位移场，进而实现图像的自动配准. 光照变化的鲁棒性，并且将描述子信息引入到模型 传统Hom-Schunck光流模型的数据项采用了平方 中，对运动细节起到了保护作用.Amiaz等[]提出 了由粗到精的分层迭代框架，抑制了迭代求解过程 函数，它会放大溢出点的位移估计误差：为了提高数 中的光流误差.Bao等[1o]利用快速保边缘区域匹配 据项对奇异值的鲁棒性，提高光流估计精度，常采用 策略处理大位移问题，同时提高了光流估计速度. 一个非平方形式的惩罚函数山(s2),如式(3)所示： Lu等)采用特征向量守恒作为约束条件，可保护 E=(IL.(X+w)-4,(X12)dK(3) 空间不连续性.Sun等]在优化过程中对中间获得 的光流场进行中值滤波以减少奇异值，提高了光流 本文选用惩罚函数山(s2)=√2+e(ε= 算法的精度.刘洪彬与常发亮]结合模糊c-均值 0.001);其中，L1和12为待配准的2幅图像：X=(x, (FMC)聚类自适应更新光流法的权重系数，使得运 y)T代表图像空间域2中的某点；W=(u,)T代表 动目标更加明显.Revaud等[4]提出一种由稀疏到 图像11和12之间的运动位移场. 稠密的保持边缘的插值方法，能够处理遮挡问题. 2.1正则项的改进 李帅等5]结合机器学习和生物模型提出一种运动 传统的光流场模型采用各向同性扩散的正则 自适应的光流估计算法，能处理光照变化与大位移 项，在图像的边缘区域容易出现过平滑，这会导致边 光流.Dosovitskiy等16]利用神经网络来实现光流 缘特征信息的丢失.为得到稠密位移场，同时尽量 预测. 保留图像纹理和边缘信息，本文采用各项异性正则 针对传统光流模型在非刚性图像配准中出现的 项来代替原来的同性扩散正则项： 边缘模糊、细节丢失和精度不足等现象，本文对传统 Es=(a。+a'g(VI))(17u2+17l2)dx 光流场模型做出了改进，提出利用各向异性正则项 来保护物体的边缘轮廓：同时通过引入非局部平滑 (4) 项，来去除噪点、保护图像细节信息.最后给出了配 g(VI)exp(-BII) (5) 准算法的整体流程，实验结果验证了算法的效果. 其中，α。代表全局平滑因子，4代表局部平滑因子； 若a为0，则该正则项与式(1)中的原正则项相同， 1传统Horm-Schunck光流场模型 g(I)是关于图像梯度VI的单调递减函数：在图像 传统的Hom-Schunck(H-S)光流场模型，是基 边缘区域可减小平滑力度，抑制边缘过平滑现象，避
工程科学学报,第 41 卷,第 7 期 图像配准[1] 是针对含有同一对象于不同条件 下获得的 2 幅或多幅图像,通过空间几何变换来实 现图像之间的几何映射,使得不同图像的对应点在 空间位置上保持一致. 该技术在医学图像处理、遥 感图像分析和计算机视觉等领域都有着广泛的应 用[2鄄鄄4] . 图像配准大致可以分成刚性图像配准与非 刚性图像配准 2 大类. 其中,非刚性图像配准以其 局部性、非线性、不满足统一的变换模型等性质,一 直是计算机视觉领域的研究难点与热点,目前常用 的非刚性图像配准方法主要包括:正则化、变分、光 流优化等. 基于光流场模型的方法是实现非刚性图像配准 的一种典型方法. 光流场指的是图像灰度模式的表 观运动,它包含了各像素的瞬时运动信息. Thiri鄄 on [5]最早引入光流方法到图像配准领域. Horn 和 Schunck [6]于 1981 年首先把速度场与图像灰度相结 合,提出了光流约束方程,给出了计算光流的基本方 法. 在此基础上,大量光流方法不断涌现,例如:Na鄄 gel 与 Enkelmann [7]提出各项异性正则项,改善了同 向正则项的过平滑问题. Brox 与 Malik [8] 提出包括 灰度一致与梯度一致的数据项,提高了光流模型对 光照变化的鲁棒性,并且将描述子信息引入到模型 中,对运动细节起到了保护作用. Amiaz 等[9] 提出 了由粗到精的分层迭代框架,抑制了迭代求解过程 中的光流误差. Bao 等[10]利用快速保边缘区域匹配 策略处理大位移问题,同时提高了光流估计速度. Liu 等[11]采用特征向量守恒作为约束条件,可保护 空间不连续性. Sun 等[12]在优化过程中对中间获得 的光流场进行中值滤波以减少奇异值,提高了光流 算法的精度. 刘洪彬与常发亮[13] 结合模糊 c鄄均值 (FMC)聚类自适应更新光流法的权重系数,使得运 动目标更加明显. Revaud 等[14] 提出一种由稀疏到 稠密的保持边缘的插值方法,能够处理遮挡问题. 李帅等[15]结合机器学习和生物模型提出一种运动 自适应的光流估计算法,能处理光照变化与大位移 光流. Dosovitskiy 等[16] 利用神经网络来实现光流 预测. 针对传统光流模型在非刚性图像配准中出现的 边缘模糊、细节丢失和精度不足等现象,本文对传统 光流场模型做出了改进,提出利用各向异性正则项 来保护物体的边缘轮廓;同时通过引入非局部平滑 项,来去除噪点、保护图像细节信息. 最后给出了配 准算法的整体流程,实验结果验证了算法的效果. 1 传统 Horn鄄鄄Schunck 光流场模型 传统的 Horn鄄鄄 Schunck(H鄄鄄 S)光流场模型,是基 于以下两个基本假设:1) 瞬时灰度值保持不变,即 灰度不变原理;2) 光流的变化是连续平滑的;基于 上述假设,在光流约束方程的基础上增加一个速度 平滑约束,将光流场的求解转换成一个变分问题. 能量函数定义如下: E(u,v) = 乙 赘 [(Ixu + Iy v + It) 2 + 琢( | 驻 u | 2 + | 驻 v| 2 )]dxdy (1) 其中,赘 代表空间图像域;u 和 v 分别为光流水平方 向与垂直方向的速度分量;Ix,Iy 和 It 分别为光流在 x,y 和 t 3 个方向上的梯度; 驻 = ( 鄣x,鄣y)是空间梯 度算子;琢 为权重系数. 式(1)可以简写为: E(u,v) = ED(u,v) + 琢ES (u,v) (2) 其中,ED(u,v)称为数据项;ES ( u,v)称为正则项;琢 为正则项的权重系数,其数值越大,对图像的约束能 力越强. 2 本文算法 研究表明,可以通过最小化能量函数来求解图 像之间的运动位移场,进而实现图像的自动配准. 传统 Horn鄄鄄 Schunck 光流模型的数据项采用了平方 函数,它会放大溢出点的位移估计误差;为了提高数 据项对奇异值的鲁棒性,提高光流估计精度,常采用 一个非平方形式的惩罚函数 鬃(s 2 ),如式(3)所示: ED = 乙 赘 鬃( |I2 (X + W) - I1 (X) | 2 )dX (3) 本文 选 用 惩 罚 函 数 鬃 ( s 2 ) = s 2 + 着 2 ( 着 = 0郾 001);其中,I1 和 I2 为待配准的 2 幅图像;X = (x, y) T 代表图像空间域 赘 中的某点;W = (u,v) T 代表 图像 I1 和 I2 之间的运动位移场. 2郾 1 正则项的改进 传统的光流场模型采用各向同性扩散的正则 项,在图像的边缘区域容易出现过平滑,这会导致边 缘特征信息的丢失. 为得到稠密位移场,同时尽量 保留图像纹理和边缘信息,本文采用各项异性正则 项来代替原来的同性扩散正则项: ES = 乙 赘 (琢g + 琢l·g( 驻 I))·鬃( | 驻 u | 2 + | 驻 v| 2 )dX (4) g( 驻 I) = exp( - 茁 | 驻 I| k ) (5) 其中,琢g 代表全局平滑因子,琢l 代表局部平滑因子; 若 琢l 为 0,则该正则项与式(1)中的原正则项相同, g( 驻 I)是关于图像梯度 驻 I 的单调递减函数;在图像 边缘区域可减小平滑力度,抑制边缘过平滑现象,避 ·956·
何凯等：一种改进的非刚性图像配准算法 ·957· 免出现边缘模糊.其中B和k为正常数，可用于控 行实验，图1给出了序列的真实光流，以及分别采用 制光流的扩散速率，值越大，平滑程度越高 各向同性与各向异性正则项获得的光流场估计结 为了验证各向异性正则项的效果，选取明德光 果.从图中可以看出，各向异性正则项有助于更好 流数据库中富含复杂性边缘的Grove3图像序列进 地保留非刚性物体的边缘信息. (a) (c) 图1采用各向同性与各向异性正则项的光流场对比.(a)真实光流：(b)各向同性正则项：(c)各向异性正则项 Fig.1 Comparison of optical flow fields with isotropic and anisotropic regularization terms:(a)ground truth;(b)isotropic regularization term;(c) anisotropic regularization term 2.2非局部平滑项的引入 其中，I(X,)为彩色向量，W(X)为位移场向量. 传统光流算法在迭代估计位移场时，经常会产 综合考虑图像结构与光流边界信息，本文实验中选 生误差累积.为了减少光流噪点对估计精度的影 取参数：01=7,02=7,03=0.5. 响，尝试通过光流的邻域信息来构建光流的空间相 该非局部项通过最小化位移矢量与其邻域值之 关信息，在传统能量函数基础上引入一个非局部平 间的平方误差，将较大的空间邻域信息整合到光流 滑项E,通过借助空间邻域信息来减少误差累积， 模型中，有助于对光流进行滤波去噪：同时，权函数 减少图像细节丢失.非局部平滑项E定义如下： 的引入又可避免出现过平滑，有助于提高整体模型 Ea=∑∑n0w[(u-u)2+ 的鲁棒性与光流估计精度.图2给出了对明德光流 ij (i')) 数据库中RubberWhale图像序列进行实验的效果 (u-")2]dK= 图，其中图(a)是序列间的真实光流场，图(b)与图 0[(u-u)2+(n-")》P]dK (6) (©)分别为未增加和增加了局部平滑项的光流场估 其中，N∈2代表邻域空间，权函数0反映了点 计结果.从图中可以看出，未增加局部平滑项，其光 X,一与点X处的像素运动矢量保持一致的可能性， 流结果图存在不少噪点，导致光流图像对比度降低， 如下式所示： 些运动细节估计也不够准确：而增加局部项后，光 流估计精度得到了显著提高. jexp 1i-12+j-12 20 2.3本文图像配准算法 I(X)-I(X:)2 IW(X:)-W(X) 综上所述，为了实现更有效的配准，本文的光流 2σ3 23 场模型的能量函数可归结为： 7 E(W)=Ep(W)+Es(W)+E(W)(8) (b) (c) 图2增加非局部平滑项前后的光流场结果对比.(a)真实光流：(b)未加局部平滑项：(c)加局部平滑项 Fig.2 Comparison of optical flow fields before and after adding non-local smoothing term:(a)ground truth;(b)without non-local smoothing term; (c)with non-local smoothing term
何 凯等: 一种改进的非刚性图像配准算法 免出现边缘模糊. 其中 茁 和 k 为正常数,可用于控 制光流的扩散速率,值越大,平滑程度越高. 为了验证各向异性正则项的效果,选取明德光 流数据库中富含复杂性边缘的 Grove3 图像序列进 行实验,图 1 给出了序列的真实光流,以及分别采用 各向同性与各向异性正则项获得的光流场估计结 果. 从图中可以看出,各向异性正则项有助于更好 地保留非刚性物体的边缘信息. 图 1 采用各向同性与各向异性正则项的光流场对比. (a)真实光流;(b)各向同性正则项;(c)各向异性正则项 Fig. 1 Comparison of optical flow fields with isotropic and anisotropic regularization terms: (a) ground truth;(b) isotropic regularization term;(c) anisotropic regularization term 2郾 2 非局部平滑项的引入 传统光流算法在迭代估计位移场时,经常会产 生误差累积. 为了减少光流噪点对估计精度的影 响,尝试通过光流的邻域信息来构建光流的空间相 关信息,在传统能量函数基础上引入一个非局部平 滑项 Enl,通过借助空间邻域信息来减少误差累积, 减少图像细节丢失. 非局部平滑项 Enl定义如下: Enl = 移i,j (i忆,j忆移 )沂N(i,j) wi,j,i忆,j忆[(uXi,j - uXi忆,j忆 ) 2 + (vXi,j - vXi忆,j忆 ) 2 ]dX = 乙 N wi,j,i忆,j忆[(u - uXi忆,j忆 ) 2 + (v - vXi忆,j忆 ) 2 ]dX (6) 图 2 增加非局部平滑项前后的光流场结果对比. (a)真实光流;(b)未加局部平滑项;(c)加局部平滑项 Fig. 2 Comparison of optical flow fields before and after adding non鄄local smoothing term: (a) ground truth;(b) without non鄄local smoothing term; (c) with non鄄local smoothing term 其中,N沂赘 代表邻域空间,权函数 wi,j,i忆,j忆反映了点 Xi忆,j忆与点 Xi,j处的像素运动矢量保持一致的可能性, 如下式所示: wi,j,i忆,j忆邑exp ( - | i - i忆| 2 + | j - j忆| 2 2滓 2 1 - |I(Xi,j) - I(Xi忆,j忆) | 2 2滓 2 2 - | W(Xi,j) - W(Xi忆,j忆) | 2 2滓 2 ) 3 (7) 其中,I(Xi,j ) 为彩色向量,W(Xi,j ) 为位移场向量. 综合考虑图像结构与光流边界信息,本文实验中选 取参数:滓1 = 7,滓2 = 7,滓3 = 0郾 5. 该非局部项通过最小化位移矢量与其邻域值之 间的平方误差,将较大的空间邻域信息整合到光流 模型中,有助于对光流进行滤波去噪;同时,权函数 的引入又可避免出现过平滑,有助于提高整体模型 的鲁棒性与光流估计精度. 图 2 给出了对明德光流 数据库中 RubberWhale 图像序列进行实验的效果 图,其中图(a)是序列间的真实光流场,图( b)与图 (c)分别为未增加和增加了局部平滑项的光流场估 计结果. 从图中可以看出,未增加局部平滑项,其光 流结果图存在不少噪点,导致光流图像对比度降低, 一些运动细节估计也不够准确;而增加局部项后,光 流估计精度得到了显著提高. 2郾 3 本文图像配准算法 综上所述,为了实现更有效的配准,本文的光流 场模型的能量函数可归结为: E(W) = ED(W) + ES (W) + 酌Enl(W) (8) ·957·
·958· 工程科学学报，第41卷，第7期 其中，y是非局部平滑项的权重系数，代表了邻域信 保持W,通过最小化式(11)求解W.通过交替最小 息对光流的正则化强度 化式(10)和式(11)可最终最小化式(8)所示的能量 为了获得图像之间的运动位移场，需要对式 函数，即可获得光流场W (8)的能量函数进行最小化求解.考虑到直接求解 由于非刚性图像配准往往存在大位移形变，这 上式比较复杂，本文通过引入辅助光流场来进行求 会引起像素点漂移，容易造成误匹配，为此，采用由 解，添加辅助项后，能量函数如下所示： 粗到精的金字塔分层细化策略计算位移场.为了进 E(w,向=(，(X+wW）-L,(X)IP)d+ 一步提高配准精度，采用文献[17]中的方法获得初 始位移场，以代替原来的零初始位移场. (a,+ag(0)(1u2+1)dx+ 算法过程为：首先对待配准的两幅图像构建N 阶金字塔：从0到N层分辨率逐层增加，0层为最高 (w-wI2)dx+yw[w-w,I]dx 层：利用文献[17]的方法求得一个位移场作为第0 (9) 层的初始位移场W;利用初始位移场，对能量函数 其中，W为辅助光流场；入少(1W-?)dx为辅 进行最小化求解获得该层的位移场：对每一层获得 的位移场进行带权值的中值滤波，并将其作为下一 助项，表示W与W间的误差平方和 层的初始光流场，直到获得第N层的位移场为止 上式可分解为2个部分： 3实验结果与分析 E.(W)=((X+W)-1(X))dx+ 为了验证本文算法对不同非刚性图像配准的有 0(a.+ag(D)(1u+1P)dK+ 效性，分别选取医学图像、柔性图像与人脸表情图像 (w-wi)ax 进行仿真实验，图3给出了进行配准实验的图像，其 (10) 中，图(a)为参考图像，图(b)为浮动图像 E2(W)=入(1W-W12)dK+ 分别利用H-S算法[6]、Box算法[8]和SIFT Fow算法[四，以及本文算法获取位移场，并纠正浮 yww-W:r1]dx (11) 动图像，结果分别如图4~图6所示.图4~图6分 首先固定W,通过最小化式(10)求解W:然后 别为一组真实脑部的核磁共振图像(MRI)、柔性图 KANGOL a 图3原始参考图像和浮动图像.(：)参考图像：(b)浮动图像 Fig.3 Original reference and floating images:(a)reference images;(b)floating images
工程科学学报,第 41 卷,第 7 期 其中,酌 是非局部平滑项的权重系数,代表了邻域信 息对光流的正则化强度. 为了获得图像之间的运动位移场,需要对式 (8)的能量函数进行最小化求解. 考虑到直接求解 上式比较复杂,本文通过引入辅助光流场来进行求 解,添加辅助项后,能量函数如下所示: E(W,W^ ) = 乙 赘 鬃( |I2 (X + W) - I1 (X) | 2 )dX + 乙 赘 (琢g + 琢l·g( 驻 I))·鬃( | 驻 u | 2 + | 驻 v| 2 )dX + 姿 乙 鬃( | W - W^ | 2 )dX + 酌 乙 N wi,j,i忆,j忆[ | W^ - W^ i忆,j忆 | ]dX (9) 其中,W^ 为辅助光流场;姿 乙 鬃( | W - W^ | 2 ) dX 为辅 助项,表示 W^ 与 W 间的误差平方和. 图 3 原始参考图像和浮动图像. (a)参考图像;(b)浮动图像 Fig. 3 Original reference and floating images: (a) reference images;(b) floating images 上式可分解为 2 个部分: E1 (W) = 乙 赘 鬃( |I2 (X + W) - I1 (X) | 2 )dX + 乙 赘 (琢g + 琢l·g( 驻 I))·鬃( | 驻 u | 2 + | 驻 v| 2 )dX + 姿 乙 鬃( | W - W^ | 2 )dX (10) E2 (W^ ) = 姿 乙 鬃( | W - W^ | 2 )dX + 酌 乙 N wi,j,i忆,j忆[ | W^ - W^ i忆,j忆 | ]dX (11) 首先固定 W^ ,通过最小化式(10)求解 W;然后 保持 W,通过最小化式(11)求解 W^ . 通过交替最小 化式(10)和式(11)可最终最小化式(8)所示的能量 函数,即可获得光流场 W. 由于非刚性图像配准往往存在大位移形变,这 会引起像素点漂移,容易造成误匹配,为此,采用由 粗到精的金字塔分层细化策略计算位移场. 为了进 一步提高配准精度,采用文献[17]中的方法获得初 始位移场,以代替原来的零初始位移场. 算法过程为:首先对待配准的两幅图像构建 N 阶金字塔;从 0 到 N 层分辨率逐层增加,0 层为最高 层;利用文献[17]的方法求得一个位移场作为第 0 层的初始位移场 W 0 ;利用初始位移场,对能量函数 进行最小化求解获得该层的位移场;对每一层获得 的位移场进行带权值的中值滤波,并将其作为下一 层的初始光流场,直到获得第 N 层的位移场为止. 3 实验结果与分析 为了验证本文算法对不同非刚性图像配准的有 效性,分别选取医学图像、柔性图像与人脸表情图像 进行仿真实验,图 3 给出了进行配准实验的图像,其 中,图(a)为参考图像,图(b)为浮动图像. 分别利用 H鄄鄄 S 算 法[6] 、 Brox 算 法[8] 和 SIFT Flow 算法[11] ,以及本文算法获取位移场,并纠正浮 动图像,结果分别如图 4 ~ 图 6 所示. 图 4 ~ 图 6 分 别为一组真实脑部的核磁共振图像(MRI)、柔性图 ·958·
何凯等：一种改进的非刚性图像配准算法 959· 像和不同表情人脸图像的纠正效果，其中，图4~图 (d)分别是利用传统H-S模型、Brox模型和SIFT 6(a)、图4~图6(b)、图4~图6(c)和图4~图6 FIow模型以及本文算法获得的图像纠正结果 回 c 图4核磁共振图像纠正结果对比.(a)H-S算法；(b)Borx算法：(c)SIFT Flow算法：(d)本文算法 Fig.4 Comparison of aligned results on MRI images:(a)H-S method;(b)Brox method;(c)SIFT flow method;(d)proposed method (a (b d 图5柔性图像纠正结果对比.(a)H-S算法：(b)Bow算法：(c)SFT1ow算法：(d)本文算法 Fig.5 Comparison of aligned results on flexible images:(a)H-S method;(b)Brox method;(c)SIFT flow method;(d)proposed method a (b) (c) d 图6人脸图像纠正结果对比.(a)H-S算法：(b)Box算法；(c)STow算法：(d)本文算法 Fig.6 Comparison of aligned results on human face images:(a)H-S method;(b)Brox method;(c)SIFT flow method;(d)proposed method 众所周知，脑部组织一般有较明显的非刚性. 其结果图像在边缘与内部存在严重模糊：Box算法 从图4能够看到，H-S算法对浮动图像起到了一定 的配准结果图像在中间区域（如数字‘38'处）存在 的纠正作用，但整体边界模糊，留存许多噪点：与之 明显的误匹配：SIFT Flow算法虽然取得了较好的配 相比，Brox算法在保持图像结构完整方面有很大提 准效果，但是仍然存在一些断点与块效应，特别是在 高，但对一些存在严重形变的部位与细节部分矫正 衣服边缘区域，存在明显的阶梯效应：而使用本文提 效果也不够理想：SIFT Flow的配准结果存在明显的 出的算法则达到了最优的配准效果. 断点，且存在明显的块效应：而本文算法由于采用各 从图6中可看出：H-S算法的配准结果图像边 项异性正则项，同时引入非局部项，在图像细节保持 界模糊严重：BOx算法虽然整体配准效果还可以， 及配准精度方面明显优于传统算法 但嘴巴部分矫正效果不佳：SIFT Flow算法结果存在 柔性图像通常具有比较严重的局部非刚性形 较多误匹配及块效应：而本文算法在眉毛、嘴巴、眼 变，从图5中可看出，利用传统H-$算法进行配准， 睛等细节部位均实现了比较精细的矫正，取得了最
何 凯等: 一种改进的非刚性图像配准算法 像和不同表情人脸图像的纠正效果,其中,图 4 ~ 图 6(a)、图 4 ~ 图 6( b)、图 4 ~ 图 6( c) 和图 4 ~ 图 6 (d)分别是利用传统 H鄄鄄 S 模型、Brox 模型和 SIFT Flow 模型以及本文算法获得的图像纠正结果. 图 4 核磁共振图像纠正结果对比. (a)H鄄鄄 S 算法;(b)Brox 算法;(c)SIFT Flow 算法;(d)本文算法 Fig. 4 Comparison of aligned results on MRI images: (a) H鄄鄄 S method;(b) Brox method;(c) SIFT flow method;(d) proposed method 图 5 柔性图像纠正结果对比. (a)H鄄鄄 S 算法;(b)Brox 算法;(c)SIFT Flow 算法;(d)本文算法 Fig. 5 Comparison of aligned results on flexible images: (a) H鄄鄄 S method;(b)Brox method;(c)SIFT flow method;(d) proposed method 图 6 人脸图像纠正结果对比. (a)H鄄鄄 S 算法;(b)Brox 算法;(c)SIFT Flow 算法;(d)本文算法 Fig. 6 Comparison of aligned results on human face images: (a) H鄄鄄 S method;(b) Brox method;(c) SIFT flow method;(d) proposed method 众所周知,脑部组织一般有较明显的非刚性. 从图 4 能够看到,H鄄鄄 S 算法对浮动图像起到了一定 的纠正作用,但整体边界模糊,留存许多噪点;与之 相比,Brox 算法在保持图像结构完整方面有很大提 高,但对一些存在严重形变的部位与细节部分矫正 效果也不够理想;SIFT Flow 的配准结果存在明显的 断点,且存在明显的块效应;而本文算法由于采用各 项异性正则项,同时引入非局部项,在图像细节保持 及配准精度方面明显优于传统算法. 柔性图像通常具有比较严重的局部非刚性形 变,从图 5 中可看出,利用传统 H鄄鄄 S 算法进行配准, 其结果图像在边缘与内部存在严重模糊;Brox 算法 的配准结果图像在中间区域(如数字‘38爷处)存在 明显的误匹配;SIFT Flow 算法虽然取得了较好的配 准效果,但是仍然存在一些断点与块效应,特别是在 衣服边缘区域,存在明显的阶梯效应;而使用本文提 出的算法则达到了最优的配准效果. 从图 6 中可看出:H鄄鄄 S 算法的配准结果图像边 界模糊严重;Brox 算法虽然整体配准效果还可以, 但嘴巴部分矫正效果不佳;SIFT Flow 算法结果存在 较多误匹配及块效应;而本文算法在眉毛、嘴巴、眼 睛等细节部位均实现了比较精细的矫正,取得了最 ·959·
.960· 工程科学学报，第41卷，第7期 佳的图像配准效果 (潘博，贲进，闫佳雯，等.全局约束下超声图像微器械轮廓 为进一步验证本文算法的准确性，采用峰值信 提取方法.哈尔滨工业大学学报，2018,50(1)：24) 噪比(PSNR)与归一化互信息(NMI)对上述图像纠 [3]Dellinger F,Delon J,Gousseau Y,et al.SAR-SIFT:a SIFT- like algorithm for SAR images.IEEE Trans Geosci Remote Sens, 正结果与原始参考图像的差值进行客观评估，结果 2015,53(1):453 如表1和表2所示. [4] Sariyanidi E,Gunes H,Cavallaro A.Automatic analysis of facial 表1不同算法的峰值信噪比对比结果 affect:a survey of registration,representation,and recognition Table 1 Comparison of PSNR values of different methods dB IEEE Trans Pattern Anal Mach Intell,2015,37(6):1113 [5]Thirion J P.Image matching as a diffusion process:an analogy 不同算法 图像 with Maxwell's demons.Med Image Anal,1998,2(3):243 H-S Brox SIFT Flow本文算法 [6]Horn B K P,Schunck B G.Determining optical flow.Artif Intell, 核磁共振 13.30 16.99 16.89 18.09 1981,17(1-3):185 柔性 18.69 21.57 25.23 25.24 [7]Nagel HH,Enkelmann W.An investigation of smoothness con- 人脸 20.27 20.35 22.08 22.82 straints for the estimation of displacement vector fields from image sequences.IEEE Trans Pattern Anal Mach Intell,1986,8(5):565 表2不同算法的归一化互信息对比结果 [8]Brox T.Malik J.Large displacement optical flow:descriptor Table 2 Comparison of NMI values of different methods matching in variational motion estimation.IEEE Trans Pattern A- nal Mach Intell,2011,33(3):500 不同算法 图像 [9] Amiaz T,Lubetzky E,Kirvati N.Coarse to over-fine optical flow H-S Brox SIFT Flow本文算法 estimation.Pattern Recognit,2007,40(9):2496 核磁共振 1.0279 1.1031 1.1255 1.1545 [10]Bao LC,Yang QX,Jin H L.Fast edge-preserving patch match 柔性 1.5076 1.5027 1.1790 1.5216 for large displacement optical flow /Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition.Colum- 人脸 1.4180 1.4324 1.4060 1.4554 bus,2014:3534 从表中可以看出：本文算法在峰值信噪比与归 [11]Liu C.Yuen J,Torralba A.SIFT flow:dense correspondence 一化互信息这两项评价指标上，数值结果均高于传 across scenes and its applications.IEEE Trans Pattern Anal Mach Intell,2011,33(5):978 统算法，表明本文算法取得了最优的配准效果 [12]Sun D Q,Roth S,Black M J.A quantitative analysis of current 4结论 practices in optical flow estimation and the principles behind them.Int J Comput Vision,2014,106(2):115 (1)针对传统光流模型在非刚性图像配准中存 [13]Liu H B,Chang F L.Moving object detection by optical flow 在的问题，提出利用各向异性正则项，根据梯度进行 method based on adaptive weight coefficient.Opt Precis Eng, 自适应平滑，保护了图像的边缘特征：此外，通过在 2016,24(2):460 (刘红彬，常发亮.权重系数自适应光流法运动目标检测. 光流场模型中引入一个非局部平滑项，可充分利用 光学精密工程，2016,24(2)：460) 邻域信息来正则化光流矢量场，通过自适应邻域加 [14] Revaud J,Weinzaepfel P,Harchaoui Z,et al.EpicFlow:edge- 权，既达到了去除光流噪点的目的，又较好地保持了 preserving interpolation of correspondences for optical flow/ 图像的细节特征，提高了图像配准精度 Proceedings of the IEEE International Conference on Computer Vi- (2)利用交替求解的方法计算位移场，用最终 sion and Pattern Recognition.Boston,2015:1164 [15]Li S,Fan X G.Xu Y L,et al.Bio-inspired motion-adaptive es- 位移场矫正浮动图像实现了图像配准.与传统算法 timation algorithm of sequence image.Chin J Eng,2017,39 的仿真实验对比结果表明：本文算法具有较高的鲁 (8):1238 棒性，能够获得更理想的非刚性图像配准效果. (李帅，樊晓光，许悦雷，等。序列图像运动自适应V1-T 光流估计算法.工程科学学报，2017,39(8)：1238) 参考文献 [16]Dosovitskiy A,Fischer P,Ilg E,et al.FlowNet:learning optical [1]Chen S W,Zhang S X,Yang X G,et al.Registration of visual- flow with convolutional networks /Proceedings of the IEEE In- infrared images based on ellipse symmetrical orientation moment. ternational Conference on Computer Vision.Santiago,2015:2758 Chin J Eng,2017,39(7):1107 [17]He K,Yan JX,Wei Y,et al.Non-rigid image registration using (陈世伟，张胜修，杨小冈，等.基于椭圆对称方向矩的可见 improved optical flow field.J Tianjin Unin Sci Technol,2018,51 光与红外图像配准算法.工程科学学报，2017.39(7)：1107) (5):491 [2]Pan B,Ben J,Yan J W,et al.Method to extract micro device (何凯，闫佳星，魏颖，等。基于改进光流场模型的非刚性图 profile in ultrasound image under global constraints.Harbin Inst 像配准.天津大学学报：自然科学与工程技术版，2018,51 Technol,2018,50(1):24 (5):491)
工程科学学报,第 41 卷,第 7 期 佳的图像配准效果. 为进一步验证本文算法的准确性,采用峰值信 噪比(PSNR)与归一化互信息(NMI)对上述图像纠 正结果与原始参考图像的差值进行客观评估,结果 如表 1 和表 2 所示. 表 1 不同算法的峰值信噪比对比结果 Table 1 Comparison of PSNR values of different methods dB 图像 不同算法 H鄄鄄 S Brox SIFT Flow 本文算法 核磁共振 13郾 30 16郾 99 16郾 89 18郾 09 柔性 18郾 69 21郾 57 25郾 23 25郾 24 人脸 20郾 27 20郾 35 22郾 08 22郾 82 表 2 不同算法的归一化互信息对比结果 Table 2 Comparison of NMI values of different methods 图像 不同算法 H鄄鄄 S Brox SIFT Flow 本文算法 核磁共振 1郾 0279 1郾 1031 1郾 1255 1郾 1545 柔性 1郾 5076 1郾 5027 1郾 1790 1郾 5216 人脸 1郾 4180 1郾 4324 1郾 4060 1郾 4554 从表中可以看出:本文算法在峰值信噪比与归 一化互信息这两项评价指标上,数值结果均高于传 统算法,表明本文算法取得了最优的配准效果. 4 结论 (1)针对传统光流模型在非刚性图像配准中存 在的问题,提出利用各向异性正则项,根据梯度进行 自适应平滑,保护了图像的边缘特征;此外,通过在 光流场模型中引入一个非局部平滑项,可充分利用 邻域信息来正则化光流矢量场,通过自适应邻域加 权,既达到了去除光流噪点的目的,又较好地保持了 图像的细节特征,提高了图像配准精度. (2)利用交替求解的方法计算位移场,用最终 位移场矫正浮动图像实现了图像配准. 与传统算法 的仿真实验对比结果表明:本文算法具有较高的鲁 棒性,能够获得更理想的非刚性图像配准效果. 参 考 文 献 [1] Chen S W, Zhang S X, Yang X G, et al. Registration of visual鄄 infrared images based on ellipse symmetrical orientation moment. Chin J Eng, 2017, 39(7): 1107 (陈世伟, 张胜修, 杨小冈, 等. 基于椭圆对称方向矩的可见 光与红外图像配准算法. 工程科学学报, 2017, 39(7): 1107) [2] Pan B, Ben J, Yan J W, et al. Method to extract micro device profile in ultrasound image under global constraints. J Harbin Inst Technol, 2018, 50(1): 24 (潘博, 贲进, 闫佳雯, 等. 全局约束下超声图像微器械轮廓 提取方法. 哈尔滨工业大学学报, 2018, 50(1): 24) [3] Dellinger F, Delon J, Gousseau Y, et al. SAR鄄鄄 SIFT: a SIFT鄄鄄 like algorithm for SAR images. IEEE Trans Geosci Remote Sens, 2015, 53(1): 453 [4] Sariyanidi E, Gunes H, Cavallaro A. Automatic analysis of facial affect: a survey of registration, representation, and recognition. IEEE Trans Pattern Anal Mach Intell, 2015, 37(6): 1113 [5] Thirion J P. Image matching as a diffusion process: an analogy with Maxwell蒺s demons. Med Image Anal, 1998, 2(3): 243 [6] Horn B K P, Schunck B G. Determining optical flow. Artif Intell, 1981, 17(1鄄3): 185 [7] Nagel H H, Enkelmann W. An investigation of smoothness con鄄 straints for the estimation of displacement vector fields from image sequences. IEEE Trans Pattern Anal Mach Intell, 1986, 8(5): 565 [8] Brox T, Malik J. Large displacement optical flow: descriptor matching in variational motion estimation. IEEE Trans Pattern A鄄 nal Mach Intell, 2011, 33(3): 500 [9] Amiaz T, Lubetzky E, Kiryati N. Coarse to over鄄fine optical flow estimation. Pattern Recognit, 2007, 40(9): 2496 [10] Bao L C, Yang Q X, Jin H L. Fast edge鄄preserving patch match for large displacement optical flow / / Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. Colum鄄 bus, 2014: 3534 [11] Liu C, Yuen J, Torralba A. SIFT flow: dense correspondence across scenes and its applications. IEEE Trans Pattern Anal Mach Intell, 2011, 33(5): 978 [12] Sun D Q, Roth S, Black M J. A quantitative analysis of current practices in optical flow estimation and the principles behind them. Int J Comput Vision, 2014, 106(2): 115 [13] Liu H B, Chang F L. Moving object detection by optical flow method based on adaptive weight coefficient. Opt Precis Eng, 2016, 24(2): 460 (刘红彬, 常发亮. 权重系数自适应光流法运动目标检测. 光学精密工程, 2016, 24(2): 460) [14] Revaud J, Weinzaepfel P, Harchaoui Z, et al. EpicFlow: edge鄄 preserving interpolation of correspondences for optical flow / / Proceedings of the IEEE International Conference on Computer Vi鄄 sion and Pattern Recognition. Boston, 2015: 1164 [15] Li S, Fan X G, Xu Y L, et al. Bio鄄inspired motion鄄adaptive es鄄 timation algorithm of sequence image. Chin J Eng, 2017, 39 (8): 1238 (李帅, 樊晓光, 许悦雷, 等. 序列图像运动自适应 V1鄄鄄 MT 光流估计算法. 工程科学学报, 2017, 39(8): 1238) [16] Dosovitskiy A, Fischer P, Ilg E, et al. FlowNet: learning optical flow with convolutional networks / / Proceedings of the IEEE In鄄 ternational Conference on Computer Vision. Santiago, 2015: 2758 [17] He K, Yan J X, Wei Y, et al. Non鄄rigid image registration using improved optical flow field. J Tianjin Univ Sci Technol, 2018, 51 (5): 491 (何凯, 闫佳星, 魏颖, 等. 基于改进光流场模型的非刚性图 像配准. 天津大学学报:自然科学与工程技术版, 2018, 51 (5): 491) ·960·