基于SIFT特征的肺部非刚性配准应用研究

史明霞，张 旭，张 涛

(1.中国科学院 自动化研究所，北京 100190；2.北京中盾安全技术开发公司，北京 100048；3.中国人民解放军总医院胸外科，北京 100853)

基于SIFT特征的肺部非刚性配准应用研究

史明霞1，张 旭2，张 涛3

(1.中国科学院 自动化研究所，北京 100190；2.北京中盾安全技术开发公司，北京 100048；3.中国人民解放军总医院胸外科，北京 100853)

随着肺癌发病率的增高和影像技术的进步，越来越多的微小肺癌(尤其是肺部磨玻璃结节)得以检出，对于此类病变首先考虑手术治疗。但肺是含气组织，在术中会发生萎陷。巨大的体积变化导致肺内结节也随之发生位移，术前CT确定的位置在术中无法定位。针对这一问题，在术中肺萎陷动物模型基础上，提出了一种融合3D-SIFT特征的B样条配准方法。该方法首先提取局部特征，筛选并匹配较好的特征点；然后依据浮动图像与参考图像匹配的特征点，得到浮动图像的初始位置；最后利用控制网格对浮动图像进行配准。经初步研究，对于不同萎陷状态的肺部图像，应用此方法得到了很好的配准结果。因此，融合3D-SIFT特征的B样条配准方法适合肺部图像配准，将来有望应用于微小肺癌的术中定位。

肺部图像配准；图像配准方法；3D-SIFT；非刚性配准；B样条配准

0 引 言

我国肺癌发病率近10年呈逐年上升趋势，平均每五个癌症患者中就有一个罹患肺癌，已成为癌症“第一杀手”。近年来，高分辨率CT的应用使越来越多的早期肺癌得以检出，其中多数为直径小于1 cm的磨玻璃密度的微小肺癌。早期肺癌的体积小，加之术中肺萎陷产生的巨大形变，导致术中无法定位，成为困扰手术的主要难题。为此，研究应用配准技术，对膨胀与系列萎陷状态的肺进行配准，以期实现肺癌术中定位。

图像配准技术在肺部图像分析中已有广泛应用，如：合并多幅扫描图像的信息，即图像融合[1]；或用于肺组织的运动预测[2]，匹配病人的呼气和吸气相图像。但是，目前对肺在术中发生的巨大形变配准研究相对较少。

一般来讲，医学图像的配准方法可分为基于灰度和基于特征两大类[3-4]：前者信息量大，计算复杂度大，但受灰度的影响小[5-6]；而后者提取图像中的特征点，只需要特征点匹配，对计算性能要求低，但其配准性能是由特征点匹配准确度以及特征点的分布决定的[7-8]。特征点提取越稳定，特征点匹配越精确，特征点分布越均匀，最终的整体配准结果越准确。对于复杂的肺部图像配准，结合基于灰度和基于特征的方法，受到越来越多的关注[9-10]。首先，使用提取特征点的方法找到特征明显的点，将特征点的对应位移信息作为灰度信息的约束；然后使用基于灰度的配准算法进行配准，显著提高了配准精度。

通过SIFT(Scale Invariant Feature Transform)算法选取的特征点具有旋转、变换、尺度不变性，对噪声、亮度的变化等都具有较强的鲁棒性，在二维图像处理中，已经得到了很好的应用[11-12]。虽然存在一些提取三维图像特征的算法，但并没有真正地得到旋转不变性的特征[13]。

在二维图像中，为了得到方向不变性，在进行方向统计之前需要根据主方向对局部图像进行旋转。然而，在三维空间中，从一个向量方向旋转到另一个方向的解不是唯一的。因此，之前的许多SIFT方法不具有旋转不变性。Allaire等为了规避该问题，利用类似于二维的方法分别计算出旋转方向角，但该方法需要的计算量大，并且角度直方图的组距没能均匀地量化。对此，Rister等提出了另一种扩展方法，使用结构张量，不仅能够避免方向角的计算，而且能够均匀地统计方向[14]。

针对早期肺癌的术中定位难题，对膨胀和萎陷状态的肺展开图像配准研究，在前述研究的基础上，提出了一种融合3D-SIFT特征的B样条配准方法。首先利用结构张量的特性对肺部CT图像的特征点进行提取、筛选和匹配[15-16]，然后进行几何形变，计算出待配准图像的初始图像，最后利用灰度图B样条配准方法[17]精细地配准肺部CT图像。将该方法应用于自建的肺渐次萎陷动物模型，对于不同萎陷状态的肺部图像能够得到很好的配准结果。

1 3D-SIFT特征点定位及特征提取

通过SIFT提取局部特征，不仅具有尺度不变性，而且具有旋转角度不变性，对于图像的亮度、角度具有很好的鲁棒性。通过高斯滤波和下采样构建尺度空间，计算连续两层之间差值(高斯差分尺度空间)的极值来寻找候选关键点，再除去那些不稳定的具有边缘响应的点即得到稳定的特征点。

1.1局部方向

为了能够构造方向不变性的特性，通常需要每个特征点具有可复现的方向，一般为局部特征的主方向(直方图中的峰值)，称为该特征点的方向。根据此方向旋转浮动图像的局部数据，从而使其具有旋转不变性。

在SIFT中，特征点的方向角度θ是由局部梯度直方图确定的，并据此确定旋转矩阵Rθ。然而，在高维空间中主方向的计算是复杂而不统一的。通常来说，旋转矩阵是正交的，即|R|=1。对于N维图像来说，旋转矩阵可以根据N维欧拉角的三角函数确定。但N维球面坐标仅知道N-1个角度。因此，不可能由梯度直方图确定一个向量作为特征点的方向。所以之前的很多3D-SIFT特征是不具有旋转不变性的。

一种简单的替代指定方向的方法是使用梯度之间的相关性，该方法具有各向同性，能够很好地应用于任意维度，称之为结构张量。

(1)

结构张量是实对称矩阵。因此，该矩阵具有正交的特征向量分解，即K=QΛQT。如果特征值是不同且有序的，那么特征向量便是唯一的(在不考虑正负的情况下)。显然，当对K进行旋转时，其特征向量也同样进行了旋转。事实上，矩阵Q不能提供鲁棒性的方向，主要由于在每个特征向量上方向是可正负的。为了解决方向上的歧义，假设Q相对于梯度是不变的，采用以图像的局部梯度为参考方向。

(2)

(3)

其中，qi为Q的第i列向量。并且约束每个特征向量均是正方向，那么矩阵R的列向量可表达为：ri=siqi。矩阵R记录了图像中每个特征点的方向。由于旋转矩阵R的行列式必须为1，因此，如果一旦发现|R|=-1，那么将最后一个列向量乘上-1，以保持旋转矩阵的性质。

1.2梯度直方图

梯度直方图用于表示在一个窗口中梯度方向的分布情况，Lowe等认为是一种鲁棒的图像表达方法[18]。在二维SIFT中，梯度方向被划分为8个组距，每一个组距的范围是π/4。在高维空间中，较多地将该方法拓展到N维球面坐标上，并且在每个球角度上分割相同大小的组距。然而该方法对维度大于2的图像是有偏差的。因为对于一个球面来说，纬度较高区域与纬度较低区域，在相同区间中所包围的面积是不同的，高纬度的面积要小于低纬度区域的面积。因此按照该方法进行分组，统计的方向是有偏差的。

梯度直方图的组距可以视为N球面上的区域。通过原点且与梯度向量一致的射线通过的区域即为该向量所在的组。之前提出的方法所构成的平面是不均匀的。为了解决此问题，将N维球面划分为面积、形状相同的多边形。且构成的凸多面体的个数满足维度N的限制。在三维空间中，最大的正多面体是正二十面体。

采用了十二面体，通过线性插值的方法以减少量化梯度对直方图的影响，而不是直接统计各个平面所对应的向量的多少。如图1所示，首先找到一个梯度向量V感兴趣的平面p1p2p3。之后利用线形插值方法确定该向量对该三角形的三个顶点的贡献度λ1,λ2,λ3，将此值分别加到各对应的顶点组上。如图1所示，这等价于12面体的每个顶点构成小三角形，且以各个顶点为重心的平面上的插值(12面体与20面体是对偶图)。

图1 统计直方图示例

类似于SIFT，特征描述为在不同的子区域(4×4×4)中构成的直方图向量，且每个直方图具有12个顶点，最终的描述符的维度为4×4×4×12=768。梯度向量对vi的贡献度为：

f(vi)=wwinwsλi‖‖

(4)

其中，wwin为窗口的权重；ws为子区域插值的权重；λi为相对顶点vi梯度的重心坐标。

最后，对这些特征描述符进行归一化，并清空掉小于常量δ的vi以减少噪音干扰，之后再次进行归一化。

1.3特征点筛选

得到的点中存在较多噪声和不稳定的点，因此需要选择一种方法过滤特征不明显的点。首先，考虑结构张量的特性。对于特征明显的角点区域，其对应的结构张量矩阵的迹必须大于0。因此删除那些小于常量β的特征点。其次，根据梯度与特征向量之间的夹角，如果两个向量接近于垂直，那么该梯度方向是不稳定的。因此抛弃那些最小角度的余弦小于γ的关键点。

(5)

最后，判定特征向量的稳定性。拒绝接受那些特征值之间变化较大的特征点。这三个条件可以过滤掉很大一部分不稳定的特征点。

2 配准变换

3D-SIFT提取的CT图像特征点以及特征描述符，可以利用最近邻方法确立两幅CT图像特征点的对应关系。依据特征描述符之间的L2距离作为相似性度量。为了保证一致性，采用同样的方法逆向匹配目标图像与浮动图像之间的特征。也就是说F1与F2匹配，有且仅有F2也与F1匹配。

(6)

使用B样条配准算法[19-20]对浮动图像进行非刚性变换。B样条是具有紧支持集的基函数，其基本思想是确定控制点即可确定一个变形场，位于变形场中的图像随控制点的形变而变形。B样条具有紧支持集，每个控制点只影响其局部领域，且可以得到C2光滑的连续变换。通过匹配的特征点，可以确定模板图像位移场。

(7)

(8)

由于B样条的计算量与控制点的个数相关，控制点越多，自由度越高，计算量就越大。因此一般采用多分辨率的B样条来配准图像。

由于可能存在成像因素导致的全局运动，因此利用特征点的对应关系拟合两幅图像之间的仿射变换，去除不同时间扫描的全局形变[21]。由于特征点匹配可能存在误判的情况，因此使用RANSAC方法来拟合对应特征点之间的全局变换[22-23]。

肺部的萎陷过程并不是简单的仿射变换，肺部的不同部位会发生不同的变换。如果简单依据仿射变换后的结果对肺部进行B样条配准容易陷入局部解中。依据特征点对应关系，将对应的特征点作为B样条向量场的控制点，由控制点之间的相对位移便可计算出浮动图像的位移场，进而得到浮动图像的初始位置，而不仅是使用仿射变换作为B样条配准的初始位置。初始的浮动图像的位置可表示为：

x'=u(x)+x

(9)

最后，使用控制网格再次对浮动图像进行精细配准。配准的目的是寻找变换函数φ，使浮动图像T通过变换尽量与参考图像R相似。一般将形变函数φ表示为一个位移函数u(x)。

φ(x)=x-u(x)

(10)

那么，该优化问题可以使用能量函数来表达。使用常用的SSD(Sum of Squared Differences)作为代价函数。

E=∑(IT(X')-IR(X))2=

∑(IT(u(x)+x)-IR(X))2

(11)

J=E+λS(u)

(12)

采用通用的变分法来解决该优化问题。该方法使用拉格朗日方程的必要条件确定该式的最小值。即寻找一条通向J(u)最小的路径。该路径上的每一个点都是沿着梯度的负方向移动。其对应的微分方程为：

(13)

而能量函数的梯度可以表示为：

(14)

其中，f=E(u)。

该微分方程最简单、直接的迭代求解办法是显式地将时间离散化。通过一个时间间隔τ在梯度方向变化得到位移的变换，按照公式更新位移。

uk+1=uk-τ(f(uk)-αλuk)

(15)

为了得到更加有效稳定的结果，增加了额外的线性搜索过程，用于确定每一次最好的时间间隔。

(16)

3 配准过程

首先，提取匹配浮动图像与参考图像之间的特征。为了能够集中采集肺内部的特征点，在提取特征点之前，使用区域增长的方法分割出CT中的肺部区域，在分割的区域内提取图像特征。

由图2可知，在肺内部，特征一般存在于血管、气管的末梢及分支处[24-25]，但这些区域有一部分并未被分割结果覆盖，因此采取图像的膨胀腐蚀算法使其尽量覆盖更多的特征明显的区域。

图2 人体肺部图像分割结果

为了提高计算效率，在被分割结果Mask覆盖的区域处，提取3D-SIFT特征点及对应的特征描述符。图3中的白色亮点即为检测到的3D-SIFT特征点。由图3可知，得到的特征点分布较为均匀，且大部分分布在特征明显的区域。

图3 人体肺部特征点提取结果

最后，计算出浮动图像变换到参考图像的全局仿射矩阵。再利用特征点的对应关系采用B样条插值的方法计算出浮动图像的初始位置，根据此初始位置利用控制网格再次使用B样条插值，并根据浮动图像的控制点的位移，迭代变换浮动图像，最终得到配准后的浮动图像。

4 实验及结果分析

为进行配准研究，建立了肺渐次萎陷动物模型，通过控制肺内的含气量，模拟肺的逐步萎陷过程。采用256层螺旋CT扫描机进行全肺扫描。经多次CT扫描，获得了肺在不同萎陷状态下的图像，记录了肺从膨胀到萎陷的全过程，同时记录了肺内大的血管、支气管等特征结构在不同萎缩状态下的位置及变化。CT扫描图像清晰，输出Dicom文件。如图4所示，从左向右分别表示肺从完全膨胀到逐渐萎陷的状态(标示1-4)。通过多次实验，逐渐形成了肺渐次萎陷动物模型数据库。

图4 肺萎陷CT扫描图像

然后，分割出目标图像和浮动图像的肺部区域。由于肺部在萎陷状态下，有一些特征会逐渐消失，因此，在提取特征点之后，增加了一些人工标注的特征点，然后利用B样条配准方法，得到匹配之后的图像。肺部在各阶段萎陷状态下的配准效果如图5所示。

图5 配准效果图

5 结束语

肺癌的高发病率使其成为全社会防治重点，而术中定位却始终困扰着肺癌的手术治疗，因此，应用图像配准技术辅助实现无创定位已成为亟待解决的问题。针对肺术中萎陷的大形变特点，采用融合3D-SIFT特征的B样条配准方法。使用改进的3D-SIFT算法提取肺部CT图像的局部特征，在此基础上利用特征点对应关系计算浮动图像的初始位置，使用B样条插值算法对萎陷状态下的肺部图像进行配准。通过肺渐次萎陷动物模型的实验证明，该方法适合肺术中萎陷过程的图像配准，将来有望应用于微小肺癌的术中定位。但是，对于萎陷过程中肺叶发生的分裂现象等，尚需进一步研究。

[1] 余 霞,葛 红,李 彬,等.基于并行计算和多层次B样条的肺部CT-PET图像配准[J].计算机应用,2009,29(7):1940-1942.

[2] Ehrhardt J,Werner R,Schmidt-Richberg A,et al.Statistical modeling of 4D respiratory lung motion using diffeomorphic image registration[J].IEEE Transactions on Medical Imaging,2011,30(2):251-265.

[3] 杨 健,李若楠,黄晨阳,等.基于局部显著边缘特征的快速图像配准算法[J].计算机应用,2014,34(1):149-153.

[4] 夏 威,高 欣,王 雷,等.基于点集与互信息的肺部CT图像三维弹性配准算法[J].江苏大学学报:自然科学版,2014,35(5):558-563.

[5] 王建文,李 青.基于点和边缘相结合特征提取的图像配准算法[J].计算机工程与设计,2009,30(4):928-930.

[6] Huang X,Paragios N,Metaxas D N.Shape registration in implicit spaces using information theory and free form deformations[J].IEEE Transactions on Pattern Analysis & Machine Intelligence,2006,28(8):1303-1318.

[7] Glaunes J,Qiu A,Miller M I,et al.Large deformation diffeomorphic metric curve mapping[J].International Journal of Computer Vision,2008,80(3):317-336.

[8] Castillo R,Castillo E,Guerra R,et al.A framework for evaluation of deformable image registration spatial accuracy using large landmark point sets[J].Physics in Medicine & Biology,2009,54(7):1849-1870.

[9] Yin Y,Hoffman E A,Ding K,et al.A cubic B-spline-based hybrid registration of lung CT images for a dynamic airway geometric model with large deformation[J].Physics in Medicine & Biology,2011,56(1):203-218.

[10] Cao K,Ding K,Reinhardt J M,et al.Improving intensity-based Lung CT registration accuracy utilizing vascular information[J].International Journal of Biomedical Imaging,2012,2012:285136.

[11] Candemir S,Jaeger S,Palaniappan K,et al.Lung segmentation in chest radiographs using anatomical atlases with nonrigid registration[J].IEEE Transactions on Medical Imaging,2014,33(2):577-590.

[12] Gong M,Zhao S,Jiao L,et al.A novel coarse-to-fine scheme for automatic image registration based on SIFT and mutual information[J].IEEE Transactions on Geoscience & Remote Sensing,2014,52(7):4328-4338.

[13] Goncalves H,Corte-Real L,Goncalves J.Automatic image registration through image segmentation and SIFT[J].IEEE Transactions on Geoscience & Remote Sensing,2011,49(7):2589-2600.

[14] Rister B,Reiter D,Zhang H,et al.Scale-and orientation-invariant keypoints in higher-dimensional data[C]//IEEE international conference on image processing.[s.l.]:IEEE,2015:3490-3494.

[15] Holden M.A review of geometric transformations for nonrigid body registration[J].IEEE Transactions on Medical Imaging,2008,27(1):111-128.

[16] Mattes D,Haynor D R,Vesselle H,et al.PET-CT image registration in the chest using free-form deformations[J].IEEE Transactions on Medical Imaging,2003,22(1):120-128.

[17] Yin Y，Hoffman E A，Lin C L.Mass preserving non-rigid registration of CT lung images using cubic B-spline[J].Medical Physics,2009,36(9):4213-4222.

[18] Lowe D G.Distinctive image features from scale-invariant keypoints[J].International Journal of Computer Vision,2004,60(2):91-110.

[19] 冯兆美,党 军,崔崤峣,等.基于B样条自由形变三维医学图像非刚性配准研究[J].影像科学与光化学,2014,32(2):200-208.

[20] Xie Z,Farin G E.Image registration using hierarchical B-splines[J].IEEE Transactions on Visualization & Computer Graphics,2004,10(1):85-94.

[21] 罗文超,刘国栋,杨海燕.SIFT和改进的RANSAC算法在图像配准中的应用[J].计算机工程与应用,2013,49(15):147-149.

[22] 常 青,张 斌,邵金玲.基于SIFT和RANSAC的特征图像匹配方法[J].华东理工大学学报:自然科学版,2012,38(6):747-751.

[23] 马丽丽,曹春梅,陈金广,等.基于RANSAC的特征点匹配算法[J].计算机工程与设计,2016,37(7):1794-1797.

[24] Bauer C,Krueger M A,Lamm W J,et al.Airway tree segmentation in serial block-face cryomicrotome images of rat lungs[J].IEEE Transactions on Bio-medical Engineering,2014,61(1):119-130.

[25] Dhara A K,Mukhopadhyay S,Dutta A,et al.A combination of shape and texture features for classification of pulmonary nodules in lung CT images[J].Journal of Digital Imaging,2016,29(4):466-475.

ResearchonApplicationofPulmonaryNon-rigidRegistrationMethodwith3D-SIFTFeatures

SHI Ming-xia1，ZHANG Xu2，ZHANG Tao3

(1.Institute of Automation，Chinese Academy of Sciences，Beijing 100190，China;2.Beijing Zhongdun Security Technology Development Co,Beijing 100048,China;3.Department of Thoracic Surgery，Chinese PLA General Hospital,Beijing 100853,China)

With the increased incidence of lung cancer and the progress of medical imaging technology，more and more small pulmonary nodules (especially for ground glass nodules) have been detected.Surgical treatment is the first consideration for such nodules.But lung is a gas containing organ,which will collapse during operation.With the huge volume change,pulmonary nodule will change its position,which makes location difficult in the operation.To solve it,an animal model database simulating lung collapse in operation is established.A B-spline curves registration method with fusion of 3D-SIFT features is proposed.Local feature is extracted firstly,and then the feature points are selected and matched.According to the feature points the initial position of the floating image is obtained.Used the control grid,the floating image is registered finally.For a series of lung CT images in different collapse states,satisfactory registration results are obtained by it which is suitable for pulmonary images registration and can be used for intraoperative localization of small pulmonary nodules in the future.

pulmonary images registration;registration methods;3D-SIFT;non-rigid registration;B-spline based registration

2016-09-25

2016-12-28 < class="emphasis_bold">网络出版时间

时间：2017-07-19

北京市自然科学基金项目(7142152)

史明霞(1980-)，女，硕士，高级工程师，研究方向为计算机视觉；张 涛，博士，通信作者，研究方向为早期肺癌的计算机辅助外科。

http://kns.cnki.net/kcms/detail/61.1450.TP.20170719.1108.014.html

TP391.9；R318.04

A

1673-629X(2017)11-0181-06

10.3969/j.issn.1673-629X.2017.11.039