CN104637061A - 一种二维三维医学图像配准方法 - Google Patents

Abstract

本发明提供能够快速对术中二维X光图像和术前三维CT图像配准的一种二维三维医学图像配准方法，包括：对三维CT图像进行滤波预处理，根据最大密度投影MIP算法，将所述三维CT图像在互相正交的两平面系内进行投影获得对应的两个数字重建影像DRR集；获取另外两个正交平面上的两个二维X光图像，并对所述二维X光图像进行滤波预处理；遍历所述两个DRR集与相应的所述另外两个正交的二维X光图像进行配准，确定对应的最相似的所述DRR及其平面内位置坐标、平面内旋转角度和平面外旋转角度；将所述坐标和角度转换到三维CT图像坐标系下得到六个配准参数。本发明所述方法能够对X光图像与预先生成的CT图像准确快速的进行实时配准。

Description

一种二维三维医学图像配准方法

技术领域

本发明涉及医学成像技术领域，特别涉及一种二维三维医学图像配准方法。

背景技术

为实现图像引导治疗(Image Guided Therapy，IGRT)，图像引导放射治疗(Image Guided Radiosurgery,IGRS)和图像引导微创治疗(Image GuidedMinimally Invasive Therapy，IGMIT)，现有的手术病灶定位过程主要包括：在术前进行三维(Three Dimension，3D)的计算机断层放射扫描成像(Computed Tomography，CT)；在术中实时扫描病灶二维(Two-Dimension，2D)的X光图像，并从CT图像中根据光线投影算法生成2D的数字重建影像(Digitally Reconstructed Radio-graphs，DRR)，通过测DRR与X光图像相似度对CT图像与X光图像进行配准，从而获得病灶(X光图像)当前在CT图像中的准确位置和方向，进而实时引导治疗，这种方法在配准精度要求较高时则需消耗长时间生成大量DRR来进行配准，因此难以实现快速准确的CT与X光间的配准，使CT图像不能快速准确的引导治疗。

发明内容

本发明提供一种二维三维医学图像配准方法，以实现快速准确的三维CT图像与二维X光图像间的配准，使三维CT图像能够快速准确的引导治疗。

本发明一种二维三维医学图像配准方法，包括：

对病灶的三维CT图像进行滤波预处理，根据最大密度投影(maximumintensity projection，MIP)算法，将所述三维CT图像在互相正交两平面系内进行投影获得对应的两个数字重建影像DRR集；

获取所述病灶在最接近所述互相正交两平面系的另外两个正交平面上的两个二维X光图像，并对所述二维X光图像进行滤波预处理；

遍历所述两个DRR集与相应的所述另外两个正交二维X光图像进行配准，确定对应的最相似的所述DRR图像及其平面内位置坐标、平面内旋转角度和平面外旋转角度；

将所述平面内位置坐标、平面内旋转角度和平面外旋转角度转换到所述三维CT图像坐标系下得到六个配准参数。

本发明的有益效果为：

本发明一种二维三维医学图像配准方法通过对三维CT图像和二维X光图像进行曲率各向异性弥散(Curvature Anisotropic Diffusion)滤波预处理并采用MIP算法生成DRR集，采用近似方法将同样进行曲率各向异性弥散滤波预处理的二维X光图像与投影生成的DRR集分别在平面内和平面外配准，再将平面的四个配准参数变换到三维CT图像坐标系得到六个配准参数，解决了对实时生成二维X光图像与预先生成的三维CT图像进行2D/3D配准时由于重建DRR数量大而不能快速配准的技术问题，从而实现了快速准确的CT图像与X光间的配准，使CT图像能够快速准确的引导治疗。

附图说明

图1是本发明一种二维三维医学图像配准方法实施例一的流程图；

图2是本发明一种二维三维医学图像配准方法实施例一中DRR集与二维X光图像配准流程图；

图3是本发明一种二维三维医学图像配准方法实施例一中两正交方向X光图像的效果图；

图4是本发明一种二维三维医学图像配准方法实施例一中平面到三维配准参数转换的原理图；

图5是本发明一种二维三维医学图像配准方法实施例二中DRR集与二维X光图像配准流程图。

具体实施方式

图1是本发明一种二维三维医学图像配准方法实施例一的流程图，如图1所示，本发明一种二维三维医学图像配准方法，包括：

S101、对病灶的三维CT图像进行滤波预处理，根据最大密度投影算法，将所述三维CT图像在互相正交两平面系内进行投影获得对应的两个数字重建影像DRR集；

所述三维CT图像还可以用CT、MRI、PER或3DRA图像替换；优选的，所述x轴为人体三维CT图像的纵轴，所述三维CT图像为术前也即在手术计划阶段生成的计算机断层放射扫描成像CT图像；

优选的，所述滤波预处理为曲率各向异性弥散(Curvature AnisotropicDiffusion，CAD)滤波预处理；

各向异性弥散算法通过创建尺度空间图像集，并与二维各项同性高斯滤波器进行卷积，从而达到上述效果,其定义方程为：

$\frac{\∂I}{\∂t}=\mathrm{div}(c(x,y,t)\▿I)=\▿c\·\▿I+c(x,y,t)\mathrm{\ΔI}---\left(1\right)$

其中，I为尺度空间图像集，Δ为拉普拉斯算子，为图像梯度，div为散度算子，c(x,y,t)为弥散系数，是一个用来控制平滑程度的单调递减函数，满足随梯度增大而减小。这样保证了在区域内部以较快速度弥散，在边缘处不再弥散，进而起到边缘增强的作用；时间步长最大不超过0.5/2^N，N为图像维度，对三维CT图像取时间步长为0.12，对于二维X光图像，取时间步长为0.06，以保证两种图像不但得到平滑效果，和保证能保留图像的解剖结构细节。

需要说明的是，曲率各向异性弥散滤波是开源医学图像算法库(InsightToolkit，ITK)内置的实现各向异性弥散滤波程序，不但保证图像得到平滑效果，而且保证能边缘增强的作用，从而保留图像的解剖结构细节；

优选的，所述根据最大密度投影MIP算法，将三维CT图像沿互相正交的方向进行投影获得两组数字重建影像DRR集包括：

将三维CT图像直角坐标系xyz的xy平面和xz平面绕x轴转动第一角度分别形成第一参考平面和第二参考平面，在将所述第一参考平面和第二参考平面同时在小于采样角阈值范围内等间隔或不等间隔绕x轴转动形成的投影平面上对所述三维CT图像执行MIP算法获得对应的两组DRR后组成第一DRR集和第二DRR集；

优选的，所述采样角阈值为5°；所述在将第一参考平面和第二参考平面同时在小于采样角阈值范围内等间隔或不等间隔绕x轴转动形成的投影平面上对所述三维CT图像执行MIP算法获得对应的两组DRR包括：

在将所述第一参考平面和第二参考平面同时在0°～3°内以第一采样密度等间隔绕x轴转动形成的投影平面上、在3°～5°内以第二采样密度等间隔绕x轴转动形成的投影平面上对所述三维CT图像执行所述MIP算法获得对应的两组DRR，所述第一采样密度大于所述第二采样密度；所述第一采样密度可以为0.05°，所述第二采样密度可以为0.02°；在本步骤中，通过对参考平面邻域内的投影平面进行深度采样，对偏转较大的投影平面进行稀疏采样，可以进一步提高生成DRR集的效率；

S102、获取所述病灶在最接近所述互相正交两平面系的另外两个正交平面上的两个二维X光图像，并对所述二维X光图像进行滤波预处理；；

优选的，获取所述病灶在最接近所述互相正交两平面系的另外两个正交平面上的两个二维X光图像包括：

将三维直角坐标系x₀y₀z₀的x₀y₀平面和x₀z₀平面绕x₀轴转动所述第一角度分别形成第一图像平面和第二图像平面、所述三维直角坐标系的位置与取向趋近所述三维CT图像直角坐标系，获取分别在所述第一图像平面和第二图像平面内的第一二维X光图像和第二二维X光图像；

优选的，所述x'轴可以为人体的纵轴；所述二维X光图像可以为在手术进行过程中实时生成的X光图像；图3是本发明一种二维三维医学图像配准方法实施例一中两正交方向X光图像的效果图；

优选的，所述对二维X光图像进行滤波预处理的过程与S101中所述对三维CT图像预处理的过程类似，不再赘述；具体来说，这里分别对术前CT图像和术中X光图像进行曲率各向异性弥散(Curvature AnisotropicDiffusion,CAD)滤波预处理，保证能保留图像的解剖结构细节。

S103、遍历所述两个DRR集与相应的所述另外两个正交平面上的两个二维X光图像进行配准，确定对应的最相似的所述DRR的平面内位置坐标、平面内旋转角度和平面外旋转角度；图2是本发明一种二维三维医学图像配准方法实施例一中DRR集与二维X光图像配准流程图，如图2所示，优选的，所述遍历所述两个DRR集与相应的所述另外两个正交平面上的两个二维X光图像进行配准，确定对应的最相似的所述DRR图像及其平面内位置坐标、平面内旋转角度和平面外旋转角度包括：

遍历所述第一DRR集和第二DRR集分别与所述第一二维X光图像和第二二维X光图像进行配准并确定对应的最相似的所述DRR图像及其相对于所述第一参考平面的第一平面内位置坐标、第一平面内旋转角度和所述绕x轴转动的第一x向平面外旋转角度，及其相对于所述第二参考平面的第二平面内位置坐标、第二平面内旋转角度和所述绕x轴转动的第二x向平面外旋转角度，包括：

S1031、对所述第一二维X光图像与所述第一DRR集中第一个DRR根据多分辨率搜索和均方差和的相似性测度进行配准，确定所述第一个DRR相对于所述第一参考平面的第一平面内位置坐标、第一平面内旋转角度；对所述第二二维X光图像与所述第二DRR集中第一个DRR根据多分辨率搜索和均方差和的相似度测度进行配准，确定所述第二参考平面的第二平面内位置坐标、第二平面内旋转角度；

S1032、根据黄金分割法，对所述第一二维X光图像和所述第一DRR集进行配准，确定所述第一DRR集中最相似的DRR相对于第一参考平面的所述绕x轴转动的第一x向平面外旋转角度，和对所述第二二维X光图像和所述第二DRR集进行配准，确定所述第二DRR集中最相似的DRR相对于第二参考平面的所述绕x轴转动的第二x向平面外旋转角度；

S104、将所述平面内位置坐标、平面内旋转角度和平面外旋转角度转换到所述三维CT图像坐标系下得到六个配准参数；图4是本发明一种二维三维医学图像配准方法实施例一中平面到三维配准参数转换的原理图，如图4所示，为描述方便，图中坐标系x'y'z'为所述三维CT图像坐标系转动第一角度后、x轴与共x轴的所述第一参考平面和所述第二参考平面组成的虚拟三维CT图像坐标系；

优选的，所述将平面内位置坐标、平面内旋转角度和平面外旋转角度转换到所述三维CT图像坐标系下得到六个配准参数包括：

根据所述第一角度，所述第一角度可以为45度，将所述第一平面内位置坐标、第一平面内旋转角度和所述绕x轴转动的第一x向平面外旋转角度，或第二平面内位置坐标、第二平面内旋转角度和所述绕x轴转动的第二x向平面外旋转角度平面内位置坐标、平面内旋转角度和所述绕x轴转动的x向平面外旋转角度转换成所述三维CT图像直角坐标系下的六个配准参数，包括：

$x=({\∂}_B{x}_B-{\∂}_A{x}_A)/2---\left(2\right)$

$y=({\∂}_A{y}_A-{\∂}_B{y}_B)/\sqrt{2}---\left(3\right)$

$z=({\∂}_A{y}_A+{\∂}_B{y}_B)/\sqrt{2}---\left(4\right)$

${\mathrm{\θ}}_x=({\mathrm{\θ}}_{x_A}+{\mathrm{\θ}}_{x_B})/2---\left(5\right)$

${\mathrm{\θ}}_y=({\mathrm{\θ}}_B-{\mathrm{\θ}}_A)/\sqrt{2}---\left(6\right)$

${\mathrm{\θ}}_z=({\mathrm{\θ}}_B-{\mathrm{\θ}}_A)/\sqrt{2}---\left(7\right)$

其中，分别表示所述在投影平面上对所述三维CT图像执行所述MIP算法时的所述DRR对所述三维CT图像在所述第一参考平面和第二参考平面内的放大因子；x_A、y_A分别表示所述最相似的DRR相对于所述第一参考平面的第一平面内位置坐标，x_B、y_B分别表示所述最相似的DRR相对于所述第二参考平面的第二平面内位置坐标；θ_A、θ_B分别表示所述对应的最相似的DRR相对于所述第一参考平面和第二参考平面的第一平面内旋转角度和第二平面内旋转角度；表示所述绕x轴转动的第一x向平面外旋转角度或所述绕x轴转动的第二x向平面外旋转角度。

本发明一种二维三维医学图像配准方法通过对三维CT图像和二维X光图像进行曲率各向异性弥散(Curvature Anisotropic Diffusion)滤波预处理并采用MIP算法生成DRR集，采用近似方法将同样进行曲率各向异性弥散滤波预处理的二维X光图像与投影生成的DRR集分别在平面内和平面外配准，再将平面的四个配准参数变换到三维CT图像坐标系得到六个配准参数，解决了对实时生成二维X光图像与预先生成的三维CT图像进行2D/3D配准时由于重建DRR数量大而不能快速配准的技术问题，从而实现了快速准确的CT图像与X光间的配准，使CT图像能够快速准确的引导治疗。

图5是本发明一种二维三维医学图像配准方法实施例二中DRR集与二维X光图像配准流程图，如图5所示，本发明一种二维三维医学图像配准方法实施例二中S201～S202与实施例一中S101～S102对应相同，S204与实施例中S104对应相同，不同之处在于，S203改变为：

S203、遍历所述两个DRR集与相应的所述另外两个正交平面上的两个二维X光图像进行配准，确定对应的最相似的所述DRR图像及其平面内位置坐标、平面内旋转角度和平面外旋转角度；

优选的，所述遍历所述两个DRR集与相应的所述另外两个正交平面上的两个二维X光图像进行配准，确定对应的最相似的所述DRR的平面内位置坐标、平面内旋转角度和平面外旋转角度包括：

遍历所述第一DRR集和第二DRR集分别与所述第一二维X光图像和第二二维X光图像进行配准并确定对应的最相似的所述DRR图像及其相对于所述第一参考平面的第一平面内位置坐标、第一平面内旋转角度和所述绕x轴转动的第一x向平面外旋转角度，以及相对于所述第二参考平面的第二平面内位置坐标、第二平面内旋转角度和所述绕x轴转动的第二x向平面外旋转角度，包括：

S2031、对所述第一二维X光图像与所述第一DRR集中第一个DRR根据多分辨率搜索和均方差和的相似性测度进行配准，确定所述第一个DRR相对于所述第一参考平面的第一平面内位置坐标、第一平面内旋转角度；对所述第二二维X光图像与所述第二DRR集中第一个DRR根据多分辨率搜索和均方差和的相似度测度进行配准，确定所述第二参考平面的第二平面内位置坐标、第二平面内旋转角度；

S2032、根据黄金分割法，对所述第一二维X光图像和所述第一DRR集进行配准，确定所述第一DRR集中最相似的DRR图像及其相对于第一参考平面的所述绕x轴转动的第一x向平面外旋转角度，和对所述第二二维X光图像和所述第二DRR集进行配准，确定所述第二DRR集中最相似的DRR图像及其相对于第二参考平面的所述绕x轴转动的第二x向平面外旋转角度；

S2033、根据陡降最小优化算法，通过利用所述第一x向平面外旋转角度和所述第二x向平面外旋转角度执行所述相似性测度的迭代，获得优化后的所述第一平面内位置坐标、第一平面内旋转角度、第二平面内位置坐标和第二平面内旋转角度。

可以看到，本实施例在一维搜索得到初步的配准结果后对配准参数进一步优化，本领域技术人员容易理解，所述优化过程是在最后一步变换得到配准参数之前进行的，且本步骤可以反复迭代以得到更优的结果，从而使得图像配准的准确度更高。

本领域普通技术人员可以理解：实现上述各方法实施例的全部或部分步骤可以通过程序指令相关的硬件来完成。前述的程序可以存储于一计算机可读取存储介质中。该程序在执行时，执行包括上述各方法实施例的步骤；而前述的存储介质包括：ROM、RAM、磁碟或者光盘等各种可以存储程序代码的介质。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (10)

1.一种二维三维医学图像配准方法，其特征在于，包括：

对病灶的三维CT图像进行滤波预处理，根据最大密度投影MIP算法，将所述三维CT图像在互相正交的两平面系内进行投影获得对应的两个数字重建影像DRR集；

获取所述病灶在最接近所述互相正交两平面系的另外两个正交平面上的两个二维X光图像，并对所述二维X光图像进行滤波预处理；

遍历所述两个DRR集与相应的所述另外两个正交的术中X光图像进行配准，确定对应的最相似的所述DRR及其平面内位置坐标、平面内旋转角度和平面外旋转角度；

将所述平面内位置坐标、平面内旋转角度和平面外旋转角度转换到所述三维CT图像坐标系下得到六个配准参数。

2.根据权利要求1所述的一种二维三维医学图像配准方法，其特征在于，所述滤波预处理为曲率各向异性弥散滤波预处理；时间步长最大不超过0.5/2^N；

其中，N为图像维度。

3.根据权利要求1所述的一种二维三维医学图像配准方法，其特征在于，所述根据最大密度投影MIP算法，将三维CT图像沿互相正交的方向进行投影获得两组数字重建影像DRR集包括：

将三维CT图像所在直角坐标系xyz的xy平面和xz平面绕x轴转动第一角度分别形成第一参考平面和第二参考平面，在将所述第一参考平面和第二参考平面同时在小于采样角阈值范围内等间隔或不等间隔绕x轴转动形成的投影平面上对所述三维CT图像执行MIP算法获得对应的两组DRR后组成第一DRR集和第二DRR集；

相应的，遍历所述两个DRR集与相应的术中两个正交平面上的X光图像进行配准，确定对应的最相似的所述DRR图像及其平面内位置坐标、平面内旋转角度和平面外旋转角度包括：

遍历所述第一DRR集和第二DRR集分别与所述两张正交的X光图像进行配准并确定对应的最相似的所述DRR图像和相对于所述第一张X光图像的第一平面内位置坐标、第一平面内旋转角度和所述绕x轴转动的第一x向平面外旋转角度，和相对于所述第二张X光图像的第二平面内位置坐标、第二平面内旋转角度和所述绕x轴转动的第二x向平面外旋转角度；

相应的，所述将平面内位置坐标、平面内旋转角度和平面外旋转角度转换到所述三维CT图像坐标系下得到六个配准参数包括：

根据所述第一角度，将所述第一平面内位置坐标、第一平面内旋转角度和所述绕x轴转动的第一x向平面外旋转角度，和第二平面内位置坐标、第二平面内旋转角度和所述绕x轴转动的第二x向平面外旋转角度平面内位置坐标、平面内旋转角度和所述绕x轴转动的x向平面外旋转角度转换成所述三维CT图像直角坐标系下的六个配准参数。

4.根据权利要求2所述一种二维三维医学图像配准方法，其特征在于，所述滤波预处理的时间步长为0.12。

5.根据权利要求3所述一种二维三维医学图像配准方法，其特征在于，所述遍历第一DRR集和第二DRR集分别与所述第一X光图像和第二X光图像进行配准确定对应的最相似的所述DRR及其相对于所述第一参考平面的第一平面内位置坐标、第一平面内旋转角度和所述绕x轴转动的第一x向平面外旋转角度，以及，相对于所述第二参考平面的第二平面内位置坐标、第二平面内旋转角度和所述绕x轴转动的第二x向平面外旋转角度包括：

对所述第一X光图像与所述第一DRR集中第一个DRR根据多分辨率搜索和均方差和的相似性测度进行配准，确定所述第一个DRR相对于所述第一参考平面的第一平面内位置坐标、第一平面内旋转角度；对所述第二X光图像与所述第二DRR集中第一个DRR根据多分辨率搜索和均方差和的相似度测度进行配准，确定所述第二参考平面的第二平面内位置坐标、第二平面内旋转角度；

根据黄金分割法，对所述第一X光图像和所述第一DRR集进行配准，确定所述第一DRR集中最相似的DRR图像及其相对于第一参考平面的所述绕x轴转动的第一x向平面外旋转角度，和对所述第二二维X光图像和所述第二DRR集进行配准，确定所述第二DRR集中最相似的DRR图像及其相对于第二参考平面的所述绕x轴转动的第二x向平面外旋转角度。

6.根据权利要求3所述一种二维三维医学图像配准方法，其特征在于，在所述对第一X光图像和所述第一DRR集进行配准，确定所述第一DRR集中最相似的DRR相对于第一参考平面的所述绕x轴转动的第一x向平面外旋转角度，和对所述第二X光图像和所述第二DRR集进行配准，确定所述第二DRR集中最相似的DRR图像及其相对于第二参考平面的所述绕x轴转动的第二x向平面外旋转角度之后，还包括：

根据陡降最小优化算法，通过利用所述第一x向平面外旋转角度和所述第二x向平面外旋转角度执行所述相似性测度的迭代，获得优化后的所述第一平面内位置坐标、第一平面内旋转角度、第二平面内位置坐标和第二平面内旋转角度。

7.根据权利要求3所述的一种二维三维医学图像配准方法，其特征在于，所述第一角度为45度，相应的，所述根据第一角度，将所述第一平面内位置坐标、第一平面内旋转角度和所述绕x轴转动的第一x向平面外旋转角度，和第二平面内位置坐标、第二平面内旋转角度和所述绕x轴转动的第二x向平面外旋转角度平面内位置坐标、平面内旋转角度和所述绕x轴转动的x向平面外旋转角度转换成所述三维CT图像直角坐标系下的六个配准参数包括：

$x=({\∂}_B{x}_B-{\∂}_A{x}_A)/2---\left(2\right)$

$y=({\∂}_A{y}_A-{\∂}_B{y}_B)/\sqrt{2}---\left(3\right)$

$z=({\∂}_A{y}_A-{\∂}_B{y}_B)/\sqrt{2}---\left(4\right)$

${\mathrm{\θ}}_x=({\mathrm{\θ}}_{x_A}+{\mathrm{\θ}}_{x_B})/2---\left(5\right)$

${\mathrm{\θ}}_y=({\mathrm{\θ}}_B+{\mathrm{\θ}}_A)/\sqrt{2}---\left(6\right)$

${\mathrm{\θ}}_Z=({\mathrm{\θ}}_B+{\mathrm{\θ}}_A)/\sqrt{2}---\left(7\right)$

其中，分别表示所述在投影平面上对所述三维CT图像执行所述MIP算法时的所述DRR对所述三维CT图像在所述第一参考平面和第二参考平面内的放大因子；x_A、y_A分别表示所述最相似的DRR相对于所述第一参考平面的第一平面内位置坐标，x_B、y_B分别表示所述最相似的DRR相对于所述第二参考平面的第二平面内位置坐标；θ_A、θ_B分别表示所述对应的最相似的DRR相对于所述第一参考平面和第二参考平面的第一平面内旋转角度和第二平面内旋转角度；表示所述绕x轴转动的第一x向平面外旋转角度或所述绕x轴转动的第二x向平面外旋转角度。

8.根据权利要求3所述一种二维三维医学图像配准方法，其特征在于，所述采样角阈值为5°；所述在将第一参考平面和第二参考平面同时在小于采样角阈值范围内等间隔或不等间隔绕x轴转动形成的投影平面上对所述三维CT图像执行MIP算法获得对应的两组DRR包括：

在将所述第一参考平面和第二参考平面同时在0°～3°内以第一采样密度等间隔绕x轴转动形成的投影平面上、在3°～5°内以第二采样密度等间隔绕x轴转动形成的投影平面上对所述三维CT图像执行所述MIP算法获得对应的两组DRR，所述第一采样密度大于所述第二采样密度。

9.根据权利要求8所述一种二维三维医学图像配准方法，其特征在于，所述第一采样密度为0.05°，所述第二采样密度为0.02°。

10.根据权利要求1～9所述一种二维三维医学图像配准方法，其特征在于，所述x轴为病灶三维CT图像的纵轴，所述x'轴为病灶的纵轴。