CN106846386B - 基于roi和共形几何代数特征不变量的3d颅部图像配准方法 - Google Patents

Abstract

本发明公开了一种基于ROI和共形几何代数特征不变量的3D颅部图像配准方法，以颅部感兴趣区的轮廓数据作为配准点云集，对参考图像和浮动图像的x,y,z三个方向进行几何校正，并截取两个模态的感兴趣区域；利用有效轮廓点云集计算其“最小投影单位球”，使得点云集到单位球球心投影距离的和最小，以此确定“第一几何不变量”以及参考图像相对于浮动图像的平移量；再依据有效轮廓点云集到“最小单位环”的投影距离的和为最值，确定“第二几何不变量”以及浮动图像相对于参考图像的旋转算子；最后两个几何不变量的统一实现浮动图像到参考图像的平移和旋转并达到配准效果。本发明能够精确的定位组织器官的三维位置，执行效率高，稳定性强。

Description

基于ROI和共形几何代数特征不变量的3D颅部图像配准方法

技术领域

本发明涉及一种颅部图像配准方法。

背景技术

现有的医学成像设备，根据医学图像的不同特点可分为解剖与功能这两种图像。解剖医学图(如：CT、MRI和X射线等)用于提供人体中内脏器官的解剖形态信息，功能医学图(如SPECT和PET等)可以反映内脏器官功能的代谢情况。人体内某些方面的状况均可通过使用相应的成像技术获取，不同的成像特性提供了多类诊断参考。而当医疗人员对病情进行判断时，往往需要同时分析病人多种形式的断层医学图，以便获得疑似病患部位的完整解剖细节和功能状况。这种情况下，对多个单一独立的图像进行判断往往不够准确，并伴随着不可预知的风险。这种背景下，医学图像配准和医学图像融合技术成为了解决上述问题的最佳途径。

3D医学图像配准是医学图像配准研究领域内的一个重要分支，与2D图像相比，其空间复杂度与自由度都有很大程度的提升。已有的配准方法如互信息、最近点迭代等，其庞大的互信息数据量使得整个配准过程需要耗费相当长的时间，计算复杂。因此在保证配准精度的同时，如何提高算法的配准速度是现在的主要研究方向。

发明内容

本发明的目的在于提供一种能够精确的定位组织器官的三维位置，且执行效率高，稳定性强的基于ROI和共形几何代数特征不变量的3D颅部图像配准方法。

本发明的技术解决方案是：

一种基于ROI和共形几何代数特征不变量的3D颅部图像配准方法，其特征是：包括下列步骤：

步骤1：预处理参考图像和浮动图像的轮廓点云，使两个模态的切片层的分辨率和尺寸达到统一；

步骤2：获取参考图像与浮动图像的3D模态数据，从中提取各自的外轮廓离散点云集和其中v是外轮廓的离散点云；r表示该点云在参考图像中，f表示该点云在浮动图像中；i则表示该点云在在点云集中的序号其中n_r，n_f分别为参考模态和浮动模态的轮廓点云总数；

步骤3：通过利用两个模态轮廓数据中耳部上端的“耳尖”特殊点，截取两个模态轮廓的有效区域，并使有效区域内两者切片层(3D医学图像是由数十张切面图层堆叠而成)的厚度达到比例一致，实现两个模态x,y,z轴3个方向的比例一致，这里的比例一致是指统一x，y方向的分辨率后，要将z方向的数据坐标要按照x，y方向变换要求，同比例的缩放；

步骤4：获取参考图像和浮动图像校正后的外轮廓离散点云集和其中n_r'，n_f'分别为参考模态和浮动模态校正后的轮廓点云总数；

步骤5：计算两个模态的点云集各自的“最小投影单位球”的第一几何不变量，得到单位球的球心o_r和o_f，并有浮动模态相对于参考模态的平移量Δ＝o_r-o_f，得到平移量Δ后将两个模态的点云集平移到以各自的球心为原点；

步骤6：计算两个模态的点云集到各自的“最小投影单位环”的第二几何不变量，并分别用符号a_r-min、a_f-min表示其点云集到单位环距离和的最小值特征轴，符号a_r-max、a_f-max表示其点云集到单位环距离和的最大值特征轴；

步骤7：计算上一步骤得到的两个模态对应4个特征轴a_r-min、a_f-min、a_r-max、a_f-max，并构造旋转算子T＝T₂T₁，其中T₁是利用“和最小”确定的参考模态和浮动模态的两根特征轴a_r-min、a_f-min计算得出，其表达式为：

T₁＝a_r-min·η₁ (1)

上式中

T₂是利用“和最大”确定的参考模态和浮动模态的两根特征轴a_r-max、a_f-max计算得出，其表达式为：

T₂＝a_r-max·η₂ (2)

上式中

步骤8：最后，利用旋转算子T实现几何变换，完成浮动模态相对于参考模态的配准。

步骤5和步骤6中第一几何不变量、第二几何不变量的具体获取方法是：

(1)在圈定感兴趣区域并三维重构后，颅部的轮廓体仍可以视为刚体G，G是由n'个质点v_i组成，即在几何代数空间下，n'则表示有效配准区域的轮廓点云总数，任意的轮廓点云可以表示为v_i＝x_ie₁+y_ie₂+z_ie₃；i为点云在点云集中的序号，e₁、e₂、e₃分别表示几何代数空间下的三个基矢量；

(2)设颅部空间中有一个单位球Q，列出点云集到球心距离的平方和表达式；

(3)点云集到“最小单位球”球心的距离和最小，利用多元函数求最值的方式，计算得单位球球心为o_r、o_f，确定“第一几何不变量”，并得到偏移量；

(4)设颅部空间中有一个单位环C，列出点云集到单位环上的距离的平方和表达式；

(5)使用拉格朗日求最优解方式，确定“第二几何不变量”。

“耳尖”特征点定位的方式是从x方向的左右两侧扫描轮廓点，并分别读取到左右两边耳部的z方向坐标最大点以确定两个“耳尖”，以这两个特征点连成直线，在颅部确定一条基准线。两个模态的三维图像再以这条基准线为标准，分别取其颅部上端和下端的最大共有区域，以确定颅部的最大感兴趣配准区域。

本发明在几何代数空间和共形几何代数空间下，以颅部感兴趣区的轮廓数据作为配准点云集，对参考模态图像和浮动模态图像的x,y,z三个方向进行几何校正，并截取两个模态的感兴趣区域(ROI)，即有效配准区域。基于有效区域分析参考图像和浮动图像之间的几何位置关系约束，利用有效轮廓点云集计算其“最小投影单位球”，方法是在颅部的三维空间中构造一个单位球，使得点云集到该单位球球心的距离的和最小，求得单位球的球心便确定了“第一几何不变量”以及浮动图像相对于参考图像的平移量；再利用有效轮廓点云集计算其“最小投影单位环”，该单位环是指在共形几何代数域下，存在一个过“最小单位球”球心的平面，该平面和单位球相交形成的一个单位环，该单位环所在平面的法矢量由离散点云集到单位环上投影距离确定，而平面的法矢量即是两个模态的特征轴。其中点云到环的投影距离是指点云到环的最小距离，其计算方法是，过点云和单位球球心的直线在单位环所在平面上的投影与单位环相交于一点，而点云到交点的距离即是该点云到环的最小距离。利用浮动模态的点云集到单位环投影距离和的最大值、最小值，分别求得特征轴af-max，af-min，利用参考模态的点云集到单位环投影距离和的最大值、最小值，分别求得特征轴ar-max、ar-min。计算得到了两种模态的四个特征轴便确定了“第二几何不变量”以及浮动图像相对于参考图像的旋转算子

本发明在共形几何代数空间下分析参考图像和浮动图像的几何位置关系，并通过全体点云参与相对不变量构造的计算方式，降低了噪声点云的干扰，ROI的圈定不仅减少了计算量，而且屏蔽了大量的干扰点云，进一步降低了配准误差，具有鲁棒性。另外，新发明中几何特征轴的构造仅利用轮廓点云数据进行配准，降低了计算量和运算复杂度，减少了配准的时间。共形几何代数概念的引入，将原欧式空间中的点点运算转变为共形点线距离运算，降低了运算复杂度，并且共形空间变换也一定程度上的缩小了运算规模，提高了配准的效率。总而言之，本发明在一定程度上提高了配准精度，减少了配准的时间，是一种稳定的、高效的、有价值的3D颅部图像配准方法。

附图说明

下面结合附图和实施例对本发明作进一步说明。

图1是图像配准流程图。

图2、图3分别是感兴趣区域圈定前、后示意图。

图4是基于“最小投影单位球”的“第一几何不变量”提取以及基于“最小投影单位环”的“第二几何不变量”提取步骤。图4是图1中步骤5和6提取两个几何不变量的具体实施方法。

图5、图6分别为颅部外轮廓点云数据集三维建模的配准前、后效果图。

具体实施方式

从物理对象上分析，由于受到颅骨的保护，同一病人头部的头颅部外轮廓具有很高的相似度，使得整体的图像数据点云可以视为由离散点云构成的刚体，配准过程也可以视为刚体的运动，包括平移与旋转等。

参照说明书附图对本发明所述的算法做以下详细地说明。

参照图1，本发明所述的配准方法具体步骤如下：

步骤1：预处理参考图像和浮动图像的轮廓点云，使两个模态的切片层的分辨率和尺寸达到统一；

步骤2：获取参考图像与浮动图像的3D模态数据，从中提取各自的外轮廓离散点云集和其中v是外轮廓的离散点云；r表示该点云在参考图像中，f表示该点云在浮动图像中；i则表示该点云在在点云集中的序号其中n_r，n_f分别为参考模态和浮动模态的轮廓点云总数；

步骤3：通过利用两个模态轮廓数据中耳部上端的“耳尖”特殊点，截取两个模态轮廓的有效区域，并使有效区域内两者切片层(3D医学图像是由数十张切面图层堆叠而成)的厚度达到比例一致，实现两个模态x,y,z轴3个方向的比例一致，这里的比例一致是指统一x，y方向的分辨率后，要将z方向的数据坐标要按照x，y方向变换要求，同比例的缩放。

步骤4：获取参考图像和浮动图像校正后的外轮廓离散点云集和其中n_r'，n_f'分别为参考模态和浮动模态校正后的轮廓点云总数；

步骤5：计算两个模态的点云集各自的“最小投影单位球”不变量，得到单位球的球心o_r和o_f，并有浮动模态相对于参考模态的平移量Δ＝o_r-o_f，得到平移量Δ后将两个模态的点云集平移到以各自的球心为原点；

步骤6：计算两个模态的点云集到各自的“最小投影单位环”的投影不变量，并分别用符号a_r-min、a_f-min表示其点云集到单位环距离和的最小值特征轴，符号a_r-max、a_f-max表示其点云集到单位环距离和的最大值特征轴；

步骤7：计算上一步骤得到的两个模态对应4个特征轴a_r-min、a_f-min、a_r-max、a_f-max，并构造旋转算子T＝T₂T₁，其中T₁是利用“和最小”确定的参考模态和浮动模态的两根特征轴a_r-min、a_f-min计算得出，其表达式为：

T₁＝a_r-min·η₁ (1)

上式中

T₂是利用“和最大”确定的参考模态和浮动模态的两根特征轴a_r-max、a_f-max计算得出，其表达式为：

T₂＝a_r-max·η₂ (2)

上式中

步骤8最后，利用旋转算子T实现几何变换，完成浮动模态相对于参考模态的配准。

参照图2、图3，所述的是两种模态统一分辨率和尺寸后的X-Z方向建模截面图，分别为配准前后两张图。配准前参考图像和浮动图像在z方向高度不一致，这是因为成像的颅部区域和厚度不一致。因此在利用颅骨外轮廓点云集进行配准之前，必须要圈定参考模态以及浮动模态的ROI，否则无效区的点云将对配准结果产生极大的干扰，影响配准精度和算法的稳定性。不难发现，耳朵上端的“耳尖”这个特殊点突出且不易受到干扰，因此本发明利用此特殊点校准图像，分别截取参考图像和浮动图像的有效配准区域。可以从圈定感兴趣区后的参考模态及浮动模态的建模剖面图看到，两种模态图像的外轮廓非常的接近，这就保证了配准的有效性和准确性。

参照图4，所述的“最小单位球”第一几何不变量和“最小单位环”第二几何不变量的提取步骤

(步骤31)在圈定感兴趣区域并三维重构后，颅部的轮廓体仍可以视为刚体G，G是由n'个质点v_i组成，即在几何代数空间下，n'则表示有效配准区域的轮廓点云总数，任意的轮廓点云可以表示为v_i＝x_ie₁+y_ie₂+z_ie₃。

(步骤32)对于给定的点云集设有一个单位球Q，其球心坐标为o＝x_oe₁+y_oe₂+z_oe₃，那么点云v_i到单位球Q球面的最小投影距离可以等价为点云到球心o的距离再减半径，Q为单位球，固有半径r＝1。为方便计算，取点云v_i到单位球Q球心的距离的平方，即：

这样的计算方法与向量空间上计算两点距离的方法基本相似，其中(x_o,y_o,z_o)为未知量，表示球心的坐标。因此对于n'个点云组成的数据点集所有点云v_i到球心o的距离的平方和为：

上式的最小值所对应的单位球的球心o的解为：o_min＝argD(o)

(步骤33)实质上，目标函数(4)是关于单位球的球心o的3个未知数x_o,y_o,z_o的函数D(o)最小值的解，属于最值问题。多元函数求最值，这里采用让其一阶偏导数等于0的方法，对式(4)中的未知数x_o,y_o,z_o分别求偏导数，得：

对方程组(5)进行求解，可得x_o,y_o,z_o为：

经验证，上式求得的坐标即为点云集投影距离的和最小的单位球球心坐标，为了区别参考模态数据及浮动模态数据所对应球心，这里分别用符号o_r,o_f表示。由于同一患者不同模态的颅部轮廓具有极大的相似性，因此这里利用轮廓点云求得的“最小投影单位球”的球心可以作为两个模态的“第一几何不变量”，并有浮动模态相对于参考模态的平移量Δ＝o_r-o_f，得到平移量Δ后将两个模态的点云集平移到以各自的球心为原点。

(步骤34)“最小环”构造几何不变量的方法是寻找一个单位环C，使得点云到环的距离的和最小。在共形几何代数空间中，点云的表达式为

其中

另外，球s为球心在原点的单位球，其表达式为

平面A为过原点的单位平面，其表达式为

A＝xe₁+ye₂+ze₃+d_Ae_∞ (9)

其中d_A为平面到原点的距离，所以d_A＝0，且A＝xe₁+ye₂+ze₃。

那么点X_i为到球s和平面A相交的环C的距离的和为

因为平面A为单位平面，||A||²＝1可以作为其约束条件，当上式D(c)为最小时所对应的解c_min的表达式为

(步骤35)此时可以通过在单一等约束下构建拉格朗日函数来对c_min进行求解。该等约束条件为：

以此构造拉格朗日函数为：

环的表达方式为球s和平面A相交，即：

C＝s∧A (14)

所以D(c)可表示为：

综上所述，L(c,λ)的表达式为：

分别对x、y、z和λ这四个未知量求偏导数，并转化为矩阵方程组为：

可以设矩阵M为：

其中

(步骤36)通过式(18)可以看出，矢量(x,y,z)^T符合M特征矢量的定义，则λ这个未知量就是相对应的特征值。并且M所对应的任一特征矢量都满足约束条件方程的要求，这就保证了矩阵M的特征矢量和相应的λ(两者都是非零量)就是(17)方程组的解。经验证，第一特征矢量是投影和最小的解，第三特征矢量是投影和最大的解。至此便求得浮动模态相对于参考模态的第二几何不变量，并计算得到旋转算子T。利用平移量将浮动模态和参考模态平移至原点处对齐，再利用旋转算子实现浮动模态到参考模态的几何变换，即可完成配准。

参照图5、图6，在医学图像使用的实际情况中，3D医学图像实际是将人体组织按三维方向生成切片图，对这些切片图进行三维重建即可生成与人体组织基本一致的三维图。假设浮动模态的医学切片图共有n张，按三维方向可将其重建为G＝{G(x,y,z)}，三维空间中的任意像素点可用p(x,y,z)表示，其中z表示切片图的层数。颅内医学图像有显著的外轮廓特征，且受到颅骨的保护使得经过仪器所采集到的图像不易因为仪器本身或外界其他因素而产生变形，这就构成一个刚体。通过扫描每层切片图的轮廓点，对颅部轮廓进行三维重建，利用三维重建图可以更加方便的分析其几何特征，并可以直观的观察配准前后浮动图像和参考图像的相对位置变化。本发明使用美国Vanderbilt大学的“回顾性图像配准评估”项目(The Retrospective Image Registration Evaluation Project，RREP)的实验数据，以患者代号patient_002为例，CT图作为浮动模态，PD为参考模态，本发明配准的评估结果显示，最小误差为1.54mm，最大误差为2.67mm，平均误差为2.13，配准精度达到亚像素级别，可以进行临床应用。

Claims (1)

1.一种基于ROI和共形几何代数特征不变量的3D颅部图像配准方法，其特征是：包括下列步骤：

步骤1：预处理参考图像和浮动图像的轮廓点云，使两个模态的切片层的分辨率和尺寸达到统一；

步骤2：获取参考图像与浮动图像的3D模态数据，从中提取各自的外轮廓离散点云集和其中v是外轮廓的离散点云；r表示该点云在参考图像中，f表示该点云在浮动图像中；i则表示该点云在点云集中的序号, 其中n_r，n_f分别为参考模态和浮动模态的轮廓点云总数；

步骤3：通过利用两个模态轮廓数据中耳部上端的“耳尖”特殊点，截取两个模态轮廓的有效区域，并使有效区域内两者切片层的厚度达到比例一致，实现两个模态x,y,z轴3个方向的比例一致，这里的比例一致是指统一x，y方向的分辨率后，要将z方向的数据坐标要按照x，y方向变换要求，同比例的缩放；

步骤4：获取参考图像和浮动图像校正后的外轮廓离散点云集和其中n_r'，n_f'分别为参考模态和浮动模态校正后的轮廓点云总数；

步骤5：计算两个模态的点云集各自的“最小投影单位球”的第一几何不变量，得到单位球的球心o_r和o_f，并有浮动模态相对于参考模态的平移量Δ＝o_r-o_f，得到平移量Δ后将两个模态的点云集平移到以各自的球心为原点；

步骤6：计算两个模态的点云集到各自的“最小投影单位环”的第二几何不变量，并分别用符号a_r-min、a_f-min表示其点云集到单位环距离和的最小值特征轴，符号a_r-max、a_f-max表示其点云集到单位环距离和的最大值特征轴；

步骤7：计算上一步骤得到的两个模态对应4个特征轴a_r-min、a_f-min、a_r-max、a_f-max，并构造旋转算子T＝T₂T₁，其中T₁是利用“和最小”确定的参考模态和浮动模态的两根特征轴a_r-min、a_f-min计算得出，其表达式为：

T₁＝a_r-min·η₁ (1)

上式中

T₂是利用“和最大”确定的参考模态和浮动模态的两根特征轴a_r-max、a_f-max计算得出，其表达式为：

T₂＝a_r-max·η₂ (2)

上式中

步骤8：最后，利用旋转算子T实现几何变换，完成浮动模态相对于参考模态的配准；

步骤5和步骤6中第一几何不变量、第二几何不变量的具体获取方法是：

(1)在圈定感兴趣区域并三维重构后，颅部的轮廓体仍可以视为刚体G，G是由n'个质点v_i组成，即在几何代数空间下，n'则表示有效配准区域的轮廓点云总数，任意的轮廓点云可以表示为v_i＝x_ie₁+y_ie₂+z_ie₃；i为点云在点云集中的序号，e₁、e₂、e₃分别表示几何代数空间下的三个基矢量；

(2)设颅部空间中有一个单位球Q，列出点云集到球心距离的平方和表达式；

(3)点云集到“最小单位球”球心的距离和最小，利用多元函数求最值的方式，计算得单位球球心为o_r、o_f，确定“第一几何不变量”，并得到偏移量；

(4)设颅部空间中有一个单位环C，列出点云集到单位环上的距离的平方和表达式；

(5)使用拉格朗日求最优解方式，确定“第二几何不变量”。