基于体素互信息的多模式和多时相图像中股骨近端配准

朱文杰 丁 辉 黄 青 彭 江 王广志#*

1(清华大学生物医学工程系，北京 100084)2(解放军总医院全军骨科研究所，北京 100853)



基于体素互信息的多模式和多时相图像中股骨近端配准

朱文杰1丁 辉1黄 青1彭 江2王广志1#*

1(清华大学生物医学工程系，北京 100084)2(解放军总医院全军骨科研究所，北京 100853)

在股骨近端骨质疏松进程以及股骨头坏死状况评估方法中，图像分析是常用的工具，通过不同时相以及不同模式的多组影像可以对病人病情进行更全面的综合评估。然而，在综合评估过程中，由于病人多次在不同系统中成像，体位的差异使不同图像组之间的解剖点位置无法一一对应，因此分析之前需要将多组图像对齐，才能观察同一感兴趣区在不同模式或不同时间骨组织状况的差异。针对这个问题，设计一种多模式、多时相图像配准的解决方案，通过图像的前处理、双阈值分类并结合贝叶斯分类的股骨分割得到股骨体素，然后通过基于归一化互信息的图像配准获得各组图像中股骨之间的三维空间刚性变换矩阵，其中CT与MR图像的配准误差在4 mm以下，CT与CT图像的配准误差在2 mm以下。利用矩阵传递关系，以CT-CT多时相的配准矩阵为基础，可获取任何两组图像间的变换矩阵。在此基础上，再进行任意两组图像的融合、点对点的分析以及骨质状况和血供状况的定量评估。通过该方案，可以对多时相、多模式图像分析中相同感兴趣的区域进行对比。

股骨近端；多时相；多模式；配准；刚性变换矩阵

引言

股骨近端骨折、坏死是临床常见的疾病。其中，股骨近端骨折包括了股骨颈、转子间、粗隆骨折等，坏死一般指股骨头局部坏死，表现为局部骨小梁的骨折、局部血供异常等。股骨近端骨折的诱因分外力压迫以及病理转移，而坏死通常是病理性原因造成。为研究局部骨质变化，跟踪病情的发展，以及评估生长因子注射等新的微介入治疗的效果，常采用多模式图像对局部骨折、小梁塌陷、血供异常等状况进行评估。

在图像分析评估中，CT图像能清晰地显示局部骨折的区域、骨质变化以及骨质缺损情况，而不同模式的磁共振(MR)图像能清晰显示出局部血供以及松质骨内组织代谢情况。因此，通过同一时相、不同模式的图像联合分析，可以对该时间点病人股骨近端的状况进行更全面的评估。而通过不同时相、同一模式的图像，也可以在不同时间点对病人股骨近端的状况进行跟踪评估。然而，由于多次成像的图像组的尺寸、空间分辨率、空间坐标参考系以及病人的体位都可能不完全一样，造成了不同组图像中股骨的空间位置不是一一对应的，无法直接对不同组图像中的同一感兴趣区进行对比，因此需要对不同时相、不同模式的多组三维图像数据集进行配准。

关于股骨近端的配准，国内外学者已经开展了研究。Janet Blumenfeld和Colin Studholme等人为分析股骨近端松质骨的各种参数，基于互信息的配准方法，对三维(3D)的高分辨率MR图像进行配准[1]；Wenjun Li、Irina Kezele、D. Louis Collins等人为分析股骨近端骨质在药物治疗后，外部压力作用后以及随年龄增长而发生的改变，基于区域生长分割以及互信息配准的方法对股骨近端进行了配准[2-3]；Ferdinando等人采用有界的最近点迭代方法，采用4自由度(4-DOF)变换矩阵对CT图像进行配准[4]；张建国和杨庆铭为解决股骨头缺血性坏死时病灶的准确定位问题，采用MC算法构造了特征曲面，再利用特征点与特征面进行配准[5]。

在前人的研究基础上，笔者设计了一套股骨近端多时相、多模式的图像配准方案。首先利用图像分割技术提取股骨近端的主体信息，再基于归一化互信息的方法实现多时相或者多模式图像中股骨近端的配准，接着将得到的变换矩阵通过矩阵变换传递通路(下文会提到)计算得到任意2组图像之间的变换矩阵。基于这些矩阵，可将各组图像的股骨近端对齐，以便定量观察不同图像组之间相同的解剖位置的骨组织状况，实现图像组的融合分析。

1 方法

1.1 变换矩阵传递通路

由于股骨是刚性结构，因此采用6自由度(6-DOF)刚性变换矩阵构建两组图像之间的空间变换关系。特别需要注意的是，一组图像中有左股骨与右股骨，由于成像中体位变化，人体姿态不可能完全一致。因此，需要2个刚性变换矩阵来分别表示两组图像之间左、右股骨的空间变换关系，同时需要将图像数据分为左、右侧分别配准。

关于图像配准时参考基准图像的选择问题，笔者采取下述方案：方案1，当不同时相的CT图之间进行配准时，以最早时间点的CT图作为参考基准；方案2，当同一时相的CT与MR图进行配准时，由于股骨在CT图中的结构最为清晰，因此以CT图像为参考基准；方案3，当不同时相的MR图之间进行配准以及不同时相的MR图与CT图进行配准时，通过方案1、2得到的变换矩阵可计算出方案3实际的变换矩阵，该计算的过程需要变换矩阵传递通路的辅助。两个时相下两种模式图像的变换矩阵传递通路如图1所示(以左股骨为例)。

图1 变换矩阵传递通路Fig.1 The transfer pathways of transformation

为得到任意两组图像之间的配准关系，需要求出图1中所有的变换矩阵。在配准情况1与2下，可以得到变换矩阵TMR1→CT1、TCT2→CT1和TMR2→CT2，对于情况3的其他变换矩阵，可以按照矩阵传递通路由前3个变换矩阵获得。注意到6DOF-刚性变换矩阵是可逆的，可以得到

(1)

TMR2→CT1=TMR2→CT2×TCT2→CT1

(2)

(3)

用这种矩阵传递的思路获取任何两组图像之间的变换关系，主要基于以下考虑：首先，股骨在MR图像中的对比度比在CT图像中的对比度要弱，同时MR图像在Z轴上(层片方向)的空间分辨率一般都比较低，因此MR图像直接与MR图像配准的精度较低。此外，同一时相的股骨形态变化很小，因此多模图像配准时选取同一时相的图像进行配准效果会更好。而对于不同时相的图像序列来说，由于CT图像的骨结构很清晰，用不同时相的CT图像进行配准可以得到更准确的结果。因此，对于不同时相、不同模式的图像组，应以不同时相的CT图之间的变换矩阵作为桥梁推导得出。

1.2 变换矩阵的求解方案

上述的思路需要求出N-1个相互独立的变换矩阵，而求解变换矩阵则采用下述思路：首先，通过图像的前处理以及分割，获得左股骨(或右股骨)的主体结构信息，然后基于图像体素的归一化互信息实现配准，得到相应的变换矩阵，其具体流程如图2所示。

图2 变换矩阵求解流程Fig.2 Flow chart of solving transformation matrix

1.2.1 图像前处理

首先是规范化坐标。从不同成像设备获得的图像坐标参考系以及图像原点可能不同，笔者将坐标统一为LPS坐标系统。

接着是图像的重采样。采用最近邻的方法进行体素的重采样，将图像统一到一致的分辨率。Ires表示重采样的图，Iori表示原图，Sres表示重采样空间分辨率，Sori表示原图空间分辨率，Round是最近邻函数。其中，Pres∈Ires和Pori∈Iori分别表示重采样后的点以及原图的点。那么，重采样后点坐标为

(4)

然后，为了初步去除大部分干扰信息，需要对图像进行剪裁，分别截取包含左、右股骨和部分髋骨的图像数据子集。在图像剪裁后，图像的原点已经发生改变，假设S表示空间分辨率，剪裁的起点坐标为(x0,y0,z0)，那么相应的变换矩阵为

(5)

最后，根据图像的噪声类型进行图像去噪处理。采用高斯平滑去除高斯噪声，用中值滤波去除椒盐噪声。三维高斯平滑的模板函数和中值滤波函数分别为

(6)

(7)

式中，I(x,y,z)表示图像在(x,y,z)处的灰度值，d是滤波半径，M是求中值函数。

1.2.2 图像中股骨的分割

股骨在CT图像中表现为高亮度信号，在MR的T2图像中表现为低亮度信号，采用双阈值的方法能粗略获取股骨的大部分结构信息，基本上排除其他的不相关组织。在双阈值分割后，部分结构噪声仍存留在股骨结构中，采用数学形态学[6]的开运算去除散点，再用闭运算连接股骨边缘缺口，有

(8)

式中，Imath是数学形态学运算后图像，Ithre是双阈值分割后的结果，SE是运算模板。

股骨与髋骨在生理解剖上是不连接的，中间隔了一层髋关节软骨。但CT图像分割时在双阈值宽度较大的情况下，股骨、髋关节和髋骨经常会连到一起。由于股骨是多运动自由度的结构，成像时其相对髋骨的位置往往又是变化不定的。因此，如果保留髋骨部分，会对后面的股骨配准带来很大的误差，髋骨的去除显得尤为重要，而去除髋骨的关键就是去除分割出来的髋关节部分。本研究采用Cheng等提出的贝叶斯分类的方法[7]解决了这个问题，很好地去除了股骨与髋骨的连接部分，得到独立的左右股骨。

接着通过交互可以将股骨提取出来，并对提取的股骨进行填充，以减少配准时边缘的权重。

通过以上处理，可以从一组三维图像数据集中得到左侧股骨和右侧股骨的两组独立的体素数据，用于进一步的配准。

1.2.3 股骨的图像配准

1.2.3.1 平移粗配

在统一坐标系、重采样以及图像剪裁等操作后，需要配准的股骨与固定的股骨位置可能相距还比较远，为了防止配准时陷入局部最优以及提高配准效率，把两组股骨体素数据通过平移实现重心的重合，在此基础上再进行精细配准。假设基准图像重心坐标为(x1,y1,z1)，配准图像的重心坐标为(x2,y2,z2)，空间分辨率为S，那么平移粗配的变换矩阵Tregist-roughly为

(9)

1.2.3.2 降采样配准

图3 配准框图Fig.3 chart of registration

在本研究中，精细配准的相似性测度采用归一化互信息，空间搜索策略采用Powell算法，收敛优化采用降采样策略，最后的配准效果与迭代次数有关，处理的流程如图3所示。

互信息可用于描述两个系统间的相关性，通常用“熵”表示，常作为多模图像配准的相似性测度[8]。H(A)表示A的熵值，H(B)表示B的熵值，H(A,B)表示联合熵，那么归一化互信息(NMI)[9]为

(10)

Powell算法[10]是一种空间迭代搜索算法，核心思想是在每个方向依次找到极值点，然后进行迭代，直至满足迭代停止条件。

为优化配准中的优化速度，采用降采样策略进行处理，先对降采样的低分辨率图像粗配，然后逐级提高分辨率达到精细的配准。假定Ifix,D和Imove,D分别表示基准图像和配准图像经过降采D倍后得到的图像，Tregist,D表示在图像在降采D倍后经过配准得到的变换矩阵。那么，D倍降采迭代得到的配准后变换矩阵为

(11)

在新一轮降采迭代配准中，降采系数降到D/2，那么初始的变换矩阵Tregist,D/2为

(12)

当降采系数降到1时，得到最终的精细配准矩阵Tregist-fine为

(13)

综合配准过程中的所有变换矩阵，两组图像之间的实际变换矩阵Tregist为

(14)

2 结果

2.1 配准结果

所采用的图像数据来自临床，其中CT图像由SIEMENSEmotion16成像系统得到。多序列的磁共振图像由GEDiscoveryMR750 3.0TMRI成像系统得到，采集了T1、T2、动态增强的掩模图像，以及随后按一定时间间隔采集的对比图像等多个序列的三维图像数据集，分别用于骨结构与组织代谢功能的分析。其中，注射造影剂后的动态增强图像，共有1 200张图片，分30个时间点采集，反映股骨近端血供的情况。通过配准可以对CT图像中感兴趣区域对应的血供状况进行跟踪评估。各组图像的成像时间、成像范围和空间分辨率信息如表1所示。

表1 图像组信息

图4列出了没有对齐的原图像(以第一次成像为例)、不同时相CT图像的配准结果、同一时相不同模式图像的配准结果，以及通过变换矩阵传递通路后其他任意2组图像的配准结果。其中，配准的结果采用棋盘格的方法展示两图像的对比。

在图像采集时，CT1图像与其对应的MR图像时间上相差1天，可以认为是同一时相的图像组，其股骨形态变化很小。从图4可以明显看到，不同时相的CT图像之间的配准程度最好，其次是同一时相CT与MR的配准。CT与MR图像的多模配准准确度取决于MR图像的质量，比如MR12与MR11图像中股骨的对比度高，配准精度较高。对于MR13与MR14，股骨在图像中相对周围组织的对比度较低，而且空间分辨率低，Z轴方向的空间体素点间隔达6.0 mm，因此直接进行配准的误差会比较大。在实际中，由于所有磁共振图像是在一次扫描中得到的，在扫描中要求病人保持体位不变，因此可以通过一系列磁共振图像中最容易配准的一个序列的配准，得到所有序列的空间变换矩阵。

图4(c)～(f)是时相1的各组图像经过矩阵传递通路后得到的配准结果，都能达到了配准的要求，由此验证了矩阵通路的思路是可行的。在此基础上可以做更多时相以及多模式的股骨近端病情评估。

为了定量衡量配准结果，按照文献[3]的方法，选取CT1中比较有特征的解剖点(见图5)，从左到右分别为方形结节(quadrate tubercle)、转子窝(trochanteric fossa)、股骨头中央凹(fovea capitis)。然后在CT2、MR11、MR12、MR13、MR14图中选择相应的解剖点，并比较通过配准后解剖点的位置误差，误差分析见表2。可以看出，配准前的平均误差都在45 mm以上，而配准后的平均误差在4 mm以下，其中同一时相的CT图之间的配准误差最小，在2 mm以下。

图4 各图像组配准前后比较。(a)CT1、CT2、MR11、MR12、MR13、MR14没有对齐时的原图；(b)CT1与 CT2、MR11、MR12、MR13、MR14的配准结果；(c)不同时相CT2与MR11、MR12、MR13、MR14在经过矩阵传递通路后的配准结果；(d)同一时相MR12与MR11、MR13、MR14在经过矩阵传递通路后的配准结果；(e)同一时相MR11和MR13、MR14在经过矩阵传递通路后的配准结果；(f)相同时相M1R3与MR14在经过矩阵传递通路后的配准结果Fig.4 The comparison of images before and after registration.(a)original images of CT1,CT2,MR11,MR12,MR13 and MR14 without alignment (b) registration results of CT1 and CT2,MR11,MR12,MR13,MR14. (c) registration results of CT2 and MR11,MR12,MR13,MR14 with different phases through matrix transfer pathways,(d)registration results of MR12 and MR11,MR13,MR14 with the same phase through matrix transfer pathways,(e)registration results of MR11 and MR13,MR14 with the same phase through matrix transfer pathways,(f)registration result of MR13 and MR14 with the same phase through matrix transfer pathways

图5 定量衡量配准结果所选取的解剖点。 (a)方形结节；(b)转子窝；(c)股骨头中央凹Fig.5 Anatomical landmarks for quantitative measurement of registration results. (a)Quadrate tubercle; (b)Trochanteric fossa; (c)fovea capitis

配准前后图像组方形结节误差/mm转子窝误差/mm股骨头中央凹/mm平均误差/mm配准前配准后CT2911900792868MR111876193020771961MR12506559648571MR132257223023402276MR14388415604469CT227240820MR11374026835MR1224192723MR1318455439MR1423202422

2.2 图像融合结果

通过变换矩阵传递通路，即可得到任意两组图像之间股骨的变换关系，借分割得到股骨的体素数据集并以此作为掩模，可以实现任何两组图像组之间股骨的融合。图6展示了CT2数据与对应的MR12图像的融合结果，从横断面和冠状面进行观察，其中(a)～(d)从横断面进行观察，(e)～(j)从冠状面进行观察，各图中左边高亮的部分来自CT2的股骨，其他部分来自MR T2图像。通过融合图像，可以看到股骨与周围组织边缘衔接性很好，从多角度、多层面看是对齐的，在此基础上可以对股骨骨质及其周围的组织结构做综合评估。

图6 图像融合结果。 (a)～(d)横断面下的融合结果； (e)～(j)冠状面下的融合结果Fig.6 The results of image fusion. (a)～(d) Fusion results in axial view; (e)～(j) are fusion results in coronal view

2.3 对应点配准结果

在得到所有图组之间的变换矩阵后，以其中一组图像作为基准，其他所有图像组的点通过变换矩阵与基准图像组的点一一对应。以时相1的CT、MR图像以及时相2的CT图像配准为例，选取了4个特征明显的点的位置，从横断面和冠状面对配准后的点的位置来观察对应点的配准结果，其对应结果如图7所示。图中第1行是从横断面进行观察，第2行是从冠状面进行观察，从左到右每列依次为CT1、CT2、MR11、MR12、MR13和MR14。虽然各组图像在尺寸、空间分辨率以及病人体位都不一致，但经过上述多时相和多模式的配准处理后，各组图像之间的生理解剖点能够对应到一起，这是后续定量评估的基础。

2.4 观察血供状况

在实现配准以后，采用CT1与MR12、MR22(动态增强像)图像对两侧股骨进行血供情况分析。其中，MR12和MR22图像均为30个时相、共2 400张图，该30个时相可以表征股骨近端血供随时间的变化。CT1与MR12的点配准结果见图8(a)，CT1与MR22的点配准结果见图8(b)，相应点的增强曲线见图8(c)、(d)。其中，标号为11、15的点是CT图像中骨质缺损明显的解剖点。

通过多模配准，可以观察股骨近端在CT图像中感兴趣区域的血供状况。其中，点11、15显示出了该部位明显增强所反映的血供异常。通过MR22与MR12的曲线比较，可以发现15号解剖点在时相2时相对时相1的血供量减少。

通过多模配准，可以观察股骨近端在CT图像中感兴趣区域的血供状况。其中，点11、15显示出了该部位明显增强所反映的血供异常。

图8 CT1与MR12、MR22的配准结果和感兴趣区域的血供分析。 (a)CT1与MR12配准后18个点的对应关系； (b)CT1与MR22配准后18个点的对应关系； (c)MR12的18个点的血供曲线； (d)MR22的18个点的血供曲线Fig.8 Registration result of CT1，MR12 and MR22 and analysis of blood supply from interesting places. (a)Correspondence of 18 points after registration of CT1 and MR12; (b) Correspondence of 18 points after registration of CT1 and MR22; (c)Blood supply curve of 18 points from MR12; (d) blood supply curve of 18 points from MR22

3 讨论

本研究的核心思想是：首先通过对三维图像数据的分割，提取股骨体素的信息，进一步通过体素互信息计算将股骨进行配准，得到不同模式或不同时相数据集中股骨间的变换矩阵，接着通过变换矩阵的传递通路得到所有图像之间的空间变换关系，在此基础上再进行后续的图像融合、点对点分析以及定量评估。

本研究的分割采用双阈值的分割方法而没有采用基于梯度或能量信息的分割方法(如Canny[11]、GVF、Snake[12]分割算法)，主要基于以下考虑：首先，本研究涉及多模图像，在一些MR图像中，股骨内部以及边缘信息与图中其他部分组织结构的边缘信息相当丰富，股骨边缘确实检测出来，但其他组织的边缘也同时检测出来，在处理中就成为了干扰信息。其次，本解决方案应该尽量减少需要调整的参量，对于Snake等一些基于能量的算法，在图像分割领域中的分割效果比较好，但需要调整弹性系数、刚性系数、迭代次数等参数，分割的灵活性不能令人满意。此外，本研究最后的目的是配准，虽然利用双阈值分割的结果会引入一些难以剔除的误差和干扰信息，特别是MR图，但由于配准是采用基于归一化互信息对分割得到的三维体素间的关系进行计算，所以在股骨内部以及边缘信息足够的情况下，配准的结果仍能达到要求。

4 结论

在股骨近端多模式、多时相图像组分析时，股骨空间位置不对应，难以进行关联分析。针对该问题，笔者提出了一种解决方案，借此能获得任意两组图像之间的左、右股骨配准的变换矩阵，从而使股骨对齐，便于进行后续的分析工作，如图像融合、不同图像组之间同一解剖点的观察对比，以及通过体素强度值定量分析不同时相图组中感兴趣区的体素与组织结构的变化。

[1]BlumenfeldJ,StudholmeC,Carballido-GamioJ, et al.Three-dimensionalimageregistrationofMRproximalfemurimagesfortheanalysisoftrabecularboneparameters[J].MedicalPhysics, 2008, 35(10): 4630-4639.

[2]LiWenjun,SodeM,SaeedI, et al.AutomatedregistrationofhipandspineforlongitudinalQCTstudies:integrationwith3Ddensitometricandstructuralanalysis[J].Bone, 2006,38:273-279.

[3]LiWenjun,KezeleI,CollinsDL, et al.Voxel-basedmodelingandquantificationoftheproximalfemurusinginter-subjectregistrationofquantitativeCTimages[J].Bone, 2007, 41(5): 888-895.

[4]BaenaFR,HawkeT,JakopecM.Aboundediterativeclosestpointmethodforminimallyinvasiveregistrationofthefemur[J].ProceedingsoftheInstitutionofMechanicalEngineers,PartH:JournalofEngineeringinMedicine, 2013,227(10):1135-1144.

[5] 张建国，杨庆铭. 股骨头坏死病灶清除计算机辅助系统配准模块和导航模块的建立[J]. 中国组织工程研究, 2011, 15(26): 4768-4772.

[6]SerraJ.Introductiontomathematicalmorphology[J].Computervision,Graphics,andImageProcessing, 1986, 35(3): 283-305.

[7]ChengYuanzhi,ZhouShengjun,WangYadong, et al.AutomaticsegmentationtechniqueforacetabulumandfemoralheadinCTimages[J].PatternRecognition, 2013, 46(11): 2969-2984.

[8]MaesF,CollignonA,VandermeulenD, et al.Multimodalityimageregistrationbymaximizationofmutualinformation[J].IEEETransactionsonMedicalImaging, 1997, 16(2): 187-198.

[9]StudholmeC,HillDLG,HawkesDJ.Anoverlapinvariantentropymeasureof3Dmedicalimagealignment[J].PatternRecognition, 1999, 32(1): 71-86.

[10]PowellMJD.Onsearchdirectionsforminimizationalgorithms[J].MathematicalProgramming, 1973, 4(1): 193-201.

[11]CannyJ.Acomputationalapproachtoedgedetection[J].IEEETransactionsonPatternAnalysisandMachineIntelligence, 1986 (6): 679-698.

[12]XuChenyang,PrinceJL.Snakes,shapes,andgradientvectorflow[J].IEEETransactionsonImageProcessing, 1998, 7(3): 359-369.

Registration for Multi-Modal and Multi-Phase Images of Proximal Femur Based on Voxel Mutual Information

Zhu Wenjie1Ding Hui1Huang Qing1Peng Jiang2Wang Guangzhi1#*

1(DepartmentofBiomedicalEngineering,TsinghuaUniversity,Beijing100084,China)2(InstituteofOrthopedicsinPLAGeneralHospital,Beijing100853,China)

Image analysis is a common method for evaluation of proximal femur osteoporosis and femoral head necrosis. This method evaluates patient’s condition based on images of proximal femur with different modes and phases. However, because the images are generated in different systems and the positions of the patient relative to systems are different, the positions of anatomical points in different images are not one-to-one consistent. We need to align points of images with different modes and phases before we can research on interesting area. To solve the problem, we proposed a solution to get spatial rigid transformation through image pre-processing, voxels segmentation of femur based on dual-threshold combined with Bayes decision rule, and femurs registration based on normalized mutual information. The error of CT-MR registration and CT-CT registration were below 4 mm and 2 mm respectively. Using the matrix transfer relationship, rigid transformation between any two images based on multi-phase CT-CT registration matrix was obtained. On this basis, image fusion of any two images, point-to-point analysis and quantitative evaluation of osseous and blood supply condition were processed. The interesting area in images with different modes and phases through the solution was compared.

proximal femur; multi-phase; multi-modal; registration; rigid transformation

10.3969/j.issn.0258-8021. 2015. 05.001

2015-06-15， 录用日期:2015-08-18

清华大学实验室创新基金；北京市科技项目(Z131100006413027)；国家自然科学基金(81127003，51361130032)

R318

A

0258-8021(2015) 05-0513-09

# 中国生物医学工程学会会员(Member, Chinese Society of Biomedical Engineering)

*通信作者(Corresponding author)， E-mail: wgz-dea@mail.tsinghua.edu.cn