何巍，魏国栋，师为礼，苗语，何飞，杨华民

（长春理工大学　计算机科学与技术学院，长春　130022）

基于点云的膝关节胫骨三维CT与MRI图像配准

何巍，魏国栋，师为礼，苗语，何飞，杨华民

（长春理工大学计算机科学与技术学院，长春130022）

同时利用CT和MRI图像对膝关节及胫骨疾病进行诊断和治疗，是目前临床上的常规方法。精确配准膝关节胫骨部分三维CT和MRI图像可以更有效地利用这两种图像中的组织信息。目前，基于图像特征或者图像灰度值的方法可以实现图像间的配准，但是都存在一些配准精度低或者时间复杂度高的缺陷。针对膝关节胫骨组织的刚性特征，采用一种基于骨组织表面重建的点云配准方法。首先对数据集进行骨表面重建，通过重建后的结果提取点云数据，对点云数据进行下采样并利用重心进行初始变换，最后利用ICP算法对点云数据进行配准。对该方法进行配准结果分析及配准精度实验，表明这种方法简单可靠并且配准精度较高。

配准；三维重建；ICP；点云

随着科学技术的不断创新和发展，利用不同成像原理的医学成像系统被应用于医学辅助诊断和治疗当中。由于不同的成像方式都有其各自不可替代的特点，以至于在医疗诊断或治疗中通常使用多种成像设备的影像共同协作来提供综合信息。目前，广泛使用的成像模式有计算机断层扫描成像（CT）、核磁共振成像（MRI）、超声成像（US）、正电子发射断层成像（PET）、单光子发射断层成像（SPECT）和数字剪影血管造影成像（DSA）等［1］。它们可以从各自的角度观察人体内的结构或功能代谢情况。其中CT对密度大、非流动性结构表现力强，如：肿瘤、骨等；MRI则对软组织比较敏感，如：膝关节软骨。为此在膝关节部位诊断和治疗当中，为了更准确的判断病情或进行治疗，需要进行CT图像和MRI图像的融合，使医生可以同时获得多方面信息［2］。在对两种模式的图像进行融合时，关键的一步就是要将这两种模式的图像进行精确配准。

三维图像配准就是求使两个三维图像中具有相同特征点的位置在某一坐标空间上达到一致时的空间变换。目前主要有基于灰度［3］和基于特征［4］的两种方法。基于特征配准是通过提取点、线、轮廓、边缘等图像特征来计算变换参数。这种算法完全依赖于图像特征的提取精度，所以配准结果往往不理想。基于灰度的配准则是直接利用图像的灰度信息进行图像间的匹配，而不需要提取任何特征信息。但这种算法需要对三维图像中所有的点进行计算，所以较基于特征匹配，计算代价庞大。

在胫骨CT与MRI图像进行配准中，骨组织在这两种图像上组织区域明显且具有无形变的特性。本文就是利用骨组织的这种刚性特征，通过对两个图像中骨组织表面进行刚性空间变换，来获得两种图像的空间位置对准。首先利用面绘制方法对胫骨CT和MRI图像进行三维可视化重建，通过重建后的三维模型提取待配准点云。之后，通过改进的迭代最近点（Iterative Closest Points，ICP）算法对得到的两个点云进行配准。最后对配准过程进行分析并对配准结果进行精度评价。

1　点云重建

1.1面绘制

目前医学图像的存储格式通常为DICOM （Digital Imaging and Communication in Medicine）格式。CT与MRI数据集为一系有序的二维断层图像构成的三维空间体素（voxel：立体像素）阵列。为了得到两种图像中胫骨组织表面的点云进行空间位置配准，首先对CT和MRI数据采用面绘制进行胫骨三维表面重建，提取重建表面模型的点作为待配准的点云。

面绘制方法基本思想是通过求解三维数据中的等值面，经过表面的特征分割处理，再通各个特征面片的拼接组合，形成三维物体表面。面绘制方法有两类，一种是基于断层轮廓线方法，该方法首先构建二维图像轮廓，然后连接相邻两二维图像的轮廓线形成物体表面，这种方法对与具有复杂轮廓的图形处理不理想，并且相邻二维图像之间的灰度信息没有得到利用，重建表面往往容易产生跳跃性形变；另一种是基于体素的表面重建方法，基本思想是在体素中提取等值面，按照一定规则将这些等值面连接形成物体面绘制模型，比如Marching Cube（移动立方体）方法［5］就是典型的代表。面绘制方法通常包括三个过程。分别是断层图像预处理、图像边缘分割提取和多边形拟合，形成物体表面［5］。

1.2Marching Cubes（MC）

Marching Cubes算法由Lorensen等人于1987年提出［6］，该算法是用于提取具有相同灰度特征的等值面（即物体表面）的一种算法。算法的核心思想是利用三维空间图像中每两层相邻的8个顶点的灰度值与所给的阈值，来计算这8个顶点所构成的立方体中的等值面位置，然后通过连接这些等值面形成的曲面构成物体表面。

MC方法主要分为三步：

（1）根据阈值提取含有等值面体素，确定分割模式。

图1为以相邻两层断层图像k与k+1分别取4个点构成的体素模型，整个三维数图像可看作由多个体素构成。

图1　体素模型

图2　Marching cubes的15种拓扑结构

对与给定的阈值c，若顶点值大于c，标记当前顶点为1，则当前顶点在形成的等值面的外侧，若顶点灰度值小于c，将该顶点标记为0，则该点的位置为等值面的内侧。若立方体内一条棱的两顶点分布于等值面的两边，则这条棱与等值面有交点。通过8个顶点不同状态可产生256种不同的体素状态，根据旋转对称性和正反对称性可精简得到最终15种拓扑结构，如图2所示［5］。

（2）求等值面与体素边缘交点。

为方便计算机查找，建立256个全状态的查找表，查找索引表共8位，每一位记录一个顶点在等值面内侧还是外侧的标志，用0和1表示，如图3所示。这样通过八个顶点标志得到索引值，一个8位的索引值对应已有查找表中的一种状态，这样就得到了交点位置和等值面的构成方式。交点位置可以通过对两个点的插值获得，最后结合交点位置和交点连接方式形成体素内等值面。

图3　状态索引表

（3）计算等值面交点法向量。

等值面交点法向量在三维渲染时用于计算光照，取得更好的视觉效果。因为等值面上的点的在面内的梯度分量为零，所以交点处的梯度矢量即为其法向量，可采用中心差分和线性插值的方法来计算各交点梯度。

1.3生成点云

利用MC算法分别对CT与MRI三维图像进行胫骨和胫骨关节的等值面提取，利用每个体素边缘与等值面的交点的集合作为待配准的CT与MRI重建点云。

2　ICP算法

2.1基本原理

ICP（迭代最近点）算法是目前三维点云配准中使用最多的一种配准方法。最初是1992年由Besl 和Mckay［7］提出的一种基于特征点的配准方法。该方法通过迭代优化变换参数矩阵，每次寻找浮动点云中每个点在目标点云中的最近点，构成对应点对，计算使对应点对配准的变换参数，若根据此参数变换没有达到配准精度则继续寻找对近点对进行变换参数求解，直到满足配准精度，得到最终的变换参数。该方法提出较早，但凭借其逻辑简单、速度快且精度高的特点仍广泛应用在各领域的三维点云配准当中［8，9］。基本原理如图4所示，在两个点云中寻找对应点对，计算变换参数（图4（a）），通过变换得到新点云（图4（b）），如不满足收敛精度则继续迭代直到精确配准（图4（c））。

图4　ICP点云配准

2.2具体步骤

假设目标点云P，浮动点云为Q，Pi、Qi分别为点云P和Q在点i上的坐标，迭代k次Q在P中的对应点设为，利用对应点对计算变换矩阵，并通过变换矩阵对浮动点云进行更新，判断点云间平均距离是否小于给定值。

对Qk中的每个点在Pk中搜索相应最近点，判断最近点的依据是使式（1）中描述的对应点均方误差最小；

采用一种闭合求解方法［9］计算旋转矩阵R；

旋转矩阵确定后，最优平移变换T为经过旋转变换后的重心差异；

通过计算得到的变换得到变换后的新点云Qk′为：

如果dk-dk+1不小于给定的阈值τ（τ＞0），返回第1步继续迭代，直到dk-dk+1＜τ或迭代达到所给最大迭代次数则计算结束［11］。

2.3ICP算法改进

传统的ICP算法可以得到精确的配准效果，在初值较好的情况下收敛性也较好。但这种方法也存在不少问题。首先，算法要求一个点集必须包含在另一个点集当中，这一条件在很多情况下无法满足；其次，要求有相近的初始位置，否则会陷入局部最优得到错误结果；而后，算法在寻找最近点对时计算量非常大；最后，采用距离最近的点作为对应点对有可能产生错误的对应点［11］。对于传统ICP算法存在的问题，通常采用分阶段［12-14］的方法对其中某阶段进行改进，例如可分为控制点的选取、特征度量、空间搜索、点对加权和变换求解［15］。

配准过程中我们针对不同的问题采用不同阶段的改进方法，在保证胫骨关节部位配准精度的同时减少配准时间。由于待配准的点云位置相差较大，首先利用两点云重心对点云进行初配准，使两点云具有相近的初始位置。为了减少计算量、加快运算速度，对点云进行一致采样方法［16］，每隔几点选取一点作为采样点云。在特征度量上采用传统ICP算法的点到点（point-to-point）距离。空间搜索上采用多维二元搜索树（K-D Tree）［17］加快查找速度。

3　实验结果及分析

实验数据采用512*512*627的CT数据集与512*512*62 MRI的数据集。图5中左侧为CT数据胫骨面绘制三维重建渲染结果，图5右侧为MRI数据胫骨关节处面绘制三维重建渲染结果。

图5　CT（左）和MRI（右）的三维点云重建渲染结果

三维重建后得到CT点云中含点210192个，MRI重建后点的个数为106450个。重建后的点云模型轮廓清晰特征明显，为之后的配准提供良好的数据基础。

对于位置相差很大的点云利用两个点云的重心进行初配准是必要的，距离较远会使ICP算法陷入局部最优配准当中，导致配准后的结果错误，如图6 （a）为目标点云与浮动点云的初始位置在z方向相差300mm，此时，若不对点云进行初配准则配准的结果如图6（b）所示，浮动图像在其初始位置陷入局部最优，图6（c）为经过重心初始定位后再进行配准的结果，浮动图像变换到目标图像的对应位置。经过初配准有效减少了局部最优的情况，但在某些极端情况下局部最优还是会出现，这时需要人工介入，确定大致的待配准区域，然后再利用ICP算法进行配准。

为了减少计算量，加快配准速度，对点云进行重采样。选择合适的采样点个数对于配准的精度和时间影响较大。采样点个数少将缺少原始点云的特征，使得经过采样后的点云与目标点云配准时产生大量错误对应点对，从而造成配准得到假的最优值，导致配准失败。采样点多会提高配准的准确性，但准确性的提高代价是配准时间大幅增加。如图7所示，在相同最大迭代次数下（1000次），50和100采样点时不能正确配准（红色），如图8所示，200个及以上采样点时，随着采样点的增加配准时间同步增长，而配准精度不再提高。

图6　初始相差较大点云配准（深色为浮动图像）

图7　采样点个数不同时配准结果

图8　采样点个数不同时配准结果

对配准结果分别采用对应点距离、浮动点到目标点中最近三点构成的平面距离和CT与MRI切片图像中人为标记对应点距离对配准精度进行评价。经过多次试验，对应点距离平均误差为1.6mm，点到平面平均误差为0.9mm，人为标记点平均误差为1.8mm。

4　结论

为了使用CT图像与MR图像融合来提高膝关节部分的诊断和治疗的精确性。提出了一种基于点云的三维CT图像与MR图像膝关节胫骨配准方法。这种方法采用MC算法进行点云重建，并对ICP算法进行改进来进行点云配准。利用多种方法对配准精度进行评价表明该方法配准精度高，计算量较小，为之后进行的融合提供有力的支持，保证膝关节部分诊断和治疗的精确性。但该方法也存在一个缺点，在浮动图像表面特征不明显时，若浮动图像相对与目标图像为倒置的，可能会导致配准后的浮动图像也为倒置的，达不到配准的要求，这也是本方法今后的改进方向。

［1］Petra A，den Elsen，et al.Medical image matching-a review with classification［J］.IEEE Engineering in Medical and Biology，1993（3）：26-39.

［2］ 杨镇源.在体三维光声图像的肿瘤微血管信息提取方法研究［D］.长春：长春理工大学，2014.

［3］ Maes F，Collignon A，Vandermeulen D，et al.Multimodality image registration by maximization of mutual information［J］.IEEE Transactions on Medical Imaging，1997（16）：187-198.

［4］Staring M，Vander-Heide U A，Klein S.Registration of Cervical MRI using multi-feature mutual information［J］.IEEE Transactions on Medical Imaging，2009，28（9）：1414-1421.

［5］ 叶青.三维重建技术在医学图像中的研究与应用［D］.西安：西安电子科技大学，2009.

［6］Lorensen WE，Cline HE.Marching cubes：a high resolution3Dsurfaceconstructionalgorithm［J］. Computer Graphics，1987，21（4）：163-169.

［7］Besl PJ，Mckay ND.Amethod for registration of 3-d shapes［J］.IEEE Transactionson Pattern Analysisand Machine Intelligence，1992，14（2）：239-256.

［8］ Furong Peng.Street view cross-sourced point cloud matching and registration［J］.Image Processing（ICIP），IEEE International Conference on，2014，2026-2030.

［9］NieBner M，Dai A，Fisher M.Combining inertial navigation and ICP for real-time 3D surface reconstruction［J］.Eurographics 2014-Short Papers，2014：13-16.

［10］Berthold K P Horn.Closed-form solution of absolute orientation using unit quaternions［J］.Journal of the Optical Society of America A.1987，4（4）：629-642.

［11］ 刘承香，阮双琛，刘繁明，等.基于迭代最近点算法的地形匹配算法可靠性分析［J］.深圳大学学报理工版，2015（01）：22-26.

［12］Rusinkiewicz S，Levoy M.Efficient Variants of the ICP Algorithm［C］.The Third International Conference on 3D Digital Image and Modeling.Quebee city，Canada，2001：145-152.

［13］ 伍毅.三维扫描信息获取的深度图像配准算法设计及开发［D］.杭州：浙江大学，2005.

［14］Nishino K，Ikeuchi K.Robust simultaneous registrationofmultiplerangeimagescomprisinga large number of points［J］.Electroics and Communications in Japan（PartⅡ：Electronics），2004，87（8）：61-74.

［15］ 李世飞，王平，沈振康.迭代最近点算法研究进展［J］.信号处理，2009，25（10）：1582-1588.

［16］Greg T，Marc L.Zippered polygon meshes form range image［C］.Proceedings of the13th InternationalConferenceonPatternRecognition1996：879-883.

［17］Zhang Z.Iterative point matching for registration of freeform curves and surfaces［J］.International Journal of Computer Vision.1994，13（2）：119-15.

Point Cloud Based Registration of 3D CT and MRI Tibia Image

HE Wei，WEI Guodong，SHI Weili，MIAO Yu，HE Fei，YANG Huamin
（School of Computer Science and Technology，Changchun University of Science and Technology，Changchun 130022）

In order to improve the accuracy of diagnosis or operation on knee and tibia，the registration between CT and MRI 3D images is required in clinical application.Accurate registration can help doctors understand anatomical structures correctly.In order to improve the algorithm for registering CT and MRI images of tibia，this paper proposes a registration method based on surface point cloud of bone.First，surface reconstruction of bone is performed by using the Marching Cubes method，and point clouds are extracted as the reconstructed result.Second，a down sample method is utilized for the point clouds to reduce computational consumption，and the transformation of displacement is estimated by using center of gravity.Finally，ICP method is used for point clouds registration.The experimental results show that the proposed method has high registration accuracy and low computational consumption.

registration；3D reconstruction；ICP；point cloud

TP317.4

A

1672-9870（2015）05-0131-05

2015-08-08

何巍（1978-），女，博士，讲师，E-mail：hw@cust.edu.cn

杨华民（1963-），男，教授，博士生导师，E-mail：yhm@cust.edu.cn