CN107230223B - 基于不连续运动的肝部三维多模态图像配准方法 - Google Patents

Abstract

本发明公开一种基于不连续运动的肝部三维多模态图像配准方法，其包括以下步骤：对获取的原始肝部序列图像进行预处理和刚性配准；采用结合三阶B样条函数的自由形变模型模拟肝部图像的弹性形变，采用基于全变分的正则项和混合测度来构建代价函数，衡量两幅图像的相似性程度，采用有限内存拟牛顿插值法来优化求解代价函数。本发明提供的基于不连续运动的肝部三维多模态图像配准方法，避免了现有技术中二范数作为正则项导致全局平滑的结果，能够很好的保留肝部器官运动导致的边界不连续性，配准精度高。

Description

基于不连续运动的肝部三维多模态图像配准方法

技术领域

本发明涉及数字图像处理技术领域，更具体地说，本发明涉及一种基于不连续运动的肝部三维多模态图像配准方法。

背景技术

在图像配准当中，通常希望得到的是一个全局平滑的位移场，但是当人体进行呼吸运动时，肺部和肝部这些器官和它们周围组织之间会存在一个相对的滑移运动，此时器官的运动和它们边界的运动是断裂的，因此会在肺和胸膜之间、肝部和腹壁之间会产生一个不连续的位移场，传统的保证全局平滑的代价函数约束项已经不适用于这种特殊的情况。

现有技术中，处理这类不连续位移场的配准问题比较主流的方法是基于分割的运动掩膜法，运动掩膜法仅能保留器官运动的真实性，不能够反应器官外部组织的有效的位移场；运动掩膜法还需要一个初始的分割操作，这对于难以分割的腹部器官来说并不适用。另外，一些基于扩散模型的配准方法被提出，如方向相关的正则项，基于Demons的自适应正则项等方法，这些方法基于非参的形变模型，没有显式的表达式来进行最优化的求解。

发明内容

针对上述技术中存在的不足之处，本发明提供一种基于不连续运动的肝部三维多模态图像配准方法，避免了现有技术中二范数作为正则项导致全局平滑的结果，能够很好的保留肝部器官运动导致的边界不连续性，配准精度高。

为了实现根据本发明的这些目的和其它优点，本发明通过以下技术方案实现：

本发明提供一种基于不连续运动的肝部三维多模态图像配准方法，其包括以下步骤：

对获取的原始肝部序列图像进行预处理和刚性配准；

采用结合三阶B样条函数的自由形变模型模拟刚性配准后肝部图像的弹性形变，采用基于全变分的正则项和混合测度来构建代价函数，衡量两幅图像的相似性程度；采用基于图像的多分辨率金字塔的优化策略，利用有限记忆拟牛顿插值法来优化所述代价函数，求得最优的变换解。

优选的是，模拟弹性形变，包括以下步骤：

定义所述肝部三维CT图像为参考图像，定义所述肝部三维MR图像为浮动图像，分别加载所述参考图像和所述浮动图像；

对所述浮动图像进行弹性形变，并采用结合三阶B样条函数的自由形变模型来模拟所述浮动图像的弹性形变；设所述浮动图像的图像域为Φ，包含n_x×n_y×n_z个间距为δ的控制点，则，

U是所述浮动图像的位移，φ为所述浮动图像的控制点，(x,y,z)为所述浮动图像上的像素点；(p_x,p_y,p_z)是距离像素点最近的控制点，l、m、n分别表示像素点(x，y，z)周围的4*4*4控制点，且

u,v,w是浮动图像上的像素点在x，y，z三个方向上与控制点间的距离，且

B_th表示三阶B样条的基函数。

优选的是，所述浮动图像的控制点网格间距设为7*7*7mm，所述浮动图像上每个像素点的形变由周围64个控制点的形变决定。

优选的是，所述混合测度，指的是基于归一化互信息的相关比，公式如下：

CRMI(M,F；φ)＝(2-NMI(M,F；φ))·(1-CR(M,F；φ))；

其中，M表示所述浮动图像，F表示所述参考图像，NMI表示归一化互信息，CR表示相关比；

在归一化互信息中，p(M；U)表示所述浮动图像的边缘概率密度，p(F)表示所述参考图像的边缘概率密度，p(M,F；U)表示两幅图像的联合概率密度；

在相关比中，X和Y分别表示参考图像和浮动图像的随机变量，Var[Y]为Y的方差，Var[Y-E(Y|X)]为Y独立于X方差。

优选的是，所述基于全变分的正则项，公式如下：

；

其中，

是求导符号，

和

分别表示形变场在三维空间三个方向上的导数，β是一个用于避免全变分函数在零处不可导的常数。

优选的是，所述代价函数为：

C＝D(F(a),M(U(a)))+λR^TV(U)；；

其中，F(x)为所述参考图像，M(x)为所述浮动图像，a为图像上的像素点，R^TV为正则项，λ为正则项的权重因子，D(F(a),M(U(a)))表示混合测度，U(a)为浮动图像的位移。

优选的是，所述λ的值为100。

优选的是，所述预处理和刚性配准，包括以下步骤：

对获取的原始肝部序列图像进行三维重建和尺度归一化处理，输出肝部三维CT图像和肝部三维MR图像；

采用水平集的方法，分别对所述肝部三维CT图像进行去重处理、对所述肝部三维MR图像进行高斯平滑去噪处理；

对去重处理的所述肝部三维CT图像和去噪处理的所述肝部三维MR图像进行刚性配准，获得两幅边界位置大致对齐的肝部图像。

优选的是，利用有限记忆拟牛顿插值法优化所述代价函数，指的是对所述代价函数进行求导，获得最优的变换解；

求导的公式为：

则，全变分的求导展开式为：

；

其中，

φ_i,j,k表示浮动图像的控制点，i，j，k值不同对应不同的控制点。

优选的是，所述基于图像的多分辨率金字塔的层数为三层，包括：

将图像降采样，先降采样为原来图像的八分之一大小，进行配准得到一个最优解，将此次结果作为第二层配准的初始值；

第二层是将图像降采样成原来的四分之一大小，进行配准求得此次最优解，将此作为第一层图像的配准起点，对原始图像进行配准从而得到最终的优化变换解。

本发明至少包括以下有益效果：

1)对原始获取的肝部序列图像进行预处理和刚性配准，用于获得初步粗匹配的肝部图像；

2)将全变分作为代价函数的正则项来约束不连续运动的位移场，此方法能够很好的保留肝部器官运动导致的边界不连续性，避免了传统技术中二范数作为正则项导致全局平滑的结果，不需要初始的分割操作，可显示求解优化函数，配准精度高；

3)采用结合三阶B样条函数的自由形变模型来模拟形变，构造一个有参的形变场以优化求解；

4)采用基于归一化互信息的相关比作为混合测度，此测度对噪声具有较好的鲁棒性，在矫正局部形变的同时也能够考虑图像间函数依赖关系，弥补了使用单一测度的缺陷，提高了配准的精度。

本发明的其它优点、目标和特征将部分通过下面的说明体现，部分还将通过对本发明的研究和实践而为本领域的技术人员所理解。

附图说明

图1(a)-1(b)为本发明所述的单模肺部CT图像配准验证的参考图像和浮动图像；

图2(a)-2(b)为采用二范数作为正则项所得到的位移场及局部边界放大图；

图3(a)-3(b)为采用全变分作为正则项所得到的位移场及局部边界放大图；

图4(a)-4(b)为原始待配准的三维肝部CT和MR图像；

图5为原始肝部三维多模态图像的融合图；

图6(a)-6(b)为分别采用二范数和全变分作为正则项所得的配准后的融合图；

图7(a)-7(b)为采用二范数作为正则项所得的肝部三维多模位移场图及局部边界放大图；

图8(a)-8(b)为采用全变分作为正则项所得的肝部三维多模位移场图及局部边界放大图；

图9为本发明所述的基于不连续运动的肝部三维多模态图像配准方法的示意图。

具体实施方式

下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。

应当理解，本文所使用的诸如“具有”、“包含”以及“包括”术语并不排除一个或多个其它元件或其组合的存在或添加。

如图9所示，本发明提供一种基于不连续运动的肝部三维多模态图像配准方法，其包括以下步骤：

S10，对获取的原始肝部序列图像进行预处理和刚性配准；

S20，采用结合三阶B样条函数的自由形变模型模拟刚性配准后肝部图像的弹性形变，采用基于全变分的正则项和混合测度来构建代价函数，衡量两幅图像的相似性程度；

S30，采用基于图像的多分辨率金字塔的优化策略，利用有限记忆拟牛顿插值法来优化所述代价函数，求得最优的变换解。

作为本发明的一种具体实施方式，模拟弹性形变，包括以下步骤：

S21，定义肝部三维CT图像为参考图像，定义肝部三维MR图像为浮动图像，分别加载参考图像和浮动图像；

S22，对浮动图像进行弹性形变，并采用结合三阶B样条函数的自由形变模型来模拟浮动图像的弹性形变，具体为：设浮动图像的图像域为Φ，包含n_x×n_y×n_z个间距为δ的控制点，则，

U是浮动图像的位移，φ为浮动图像的控制点，(x,y,z)为浮动图像上的像素点；(p_x,p_y,p_z)是距离像素点最近的控制点，l、m、n分别表示像素点(x，y，z)周围的4*4*4控制点，且

u,v,w是浮动图像上的像素点在x，y，z三个方向上与控制点间的距离，且

B_th表示三阶B样条的基函数。

该实施方式中，利用三阶B样条函数与自由形变模型结合来模拟刚性配准后肝部图像的弹性形变，构造一个有参的形变场来优化求解，有利于提高配准精度。每个像素的形变都由周围的若干个控制点的形变来决定。若干个控制点之间的网格间距太小会导致过拟合、太大会导致配准精度不够，作为优选，浮动图像的控制点网格间距设为7*7*7mm，浮动图像上每个像素点的形变由周围64个控制点的形变决定，具有较好的拟合度和配准精度。

作为本发明的一种具体实施方式，混合测度，指的是基于归一化互信息的相关比，公式如下：

CRMI(M,F；φ)＝(2-NMI(M,F；φ))·(1-CR(M,F；φ))；

其中，M表示浮动图像，F表示参考图像，NMI表示归一化互信息，CR表示相关比；在归一化互信息中，p(M；U)表示浮动图像的边缘概率密度，p(F)表示参考图像的边缘概率密度，p(M,F；U)表示两幅图像的联合概率密度；在相关比中，X和Y分别表示参考图像和浮动图像的随机变量，Var[Y]为Y的方差，Var[Y-E(Y|X)]为Y独立于X方差。

该实施方式中，采用归一化互信息与相关比相结合的混合测度，此测度对噪声具有较好的鲁棒性，在矫正局部形变的同时也能够考虑图像间函数依赖关系，弥补了使用单一测度的缺陷，提高了配准的精度。

作为本发明的一种具体实施方式，基于全变分的正则项，公式如下：

其中，

是求导符号，

和

分别表示形变场在三维空间三个方向上的导数，β是一个用于避免全变分函数在零处不可导的常数。

该实施方式中，将全变分作为代价函数的正则项来求解不连续位移场的配准，能够很好的保留肝部器官运动导致的边界不连续性，避免了传统技术中二范数作为正则项导致全局平滑的结果，不需要初始的分割操作，可显示求解优化函数，配准精度高。

作为本发明的一种具体实施方式，基于全变分的正则项和混合测度来构建代价函数，代价函数公式为：

C＝D(F(a),M(U(a)))+λR^TV(U)；

其中，F(x)为参考图像，M(x)为浮动图像，a为图像上的像素点，R^TV为正则项，λ为正则项的权重因子,D(F(a),M(U(a)))表示混合测度，U(a)为浮动图像的位移。

该实施方式中，代价函数的构建是基于全变分的正则项和混合测度，具有较好的配准精度。当λ的值设置为100，可以获得较好的变换解。

作为本发明的一种具体实施方式，预处理和刚性配准，包括以下步骤：

S11，对获取的原始肝部序列图像进行三维重建和尺度归一化处理，输出肝部三维CT图像和肝部三维MR图像；

S12，采用水平集的方法，分别对肝部三维CT图像进行去重处理、对肝部三维MR图像进行高斯平滑去噪处理；

S13，对去重处理的肝部三维CT图像和去噪处理的肝部三维MR图像进行刚性配准，获得两幅边界位置大致对齐的肝部图像。

该实施方式中，对获取的原始肝部序列图像进行预处理和刚性配准，用于获取粗匹配的肝部图像。

作为本发明的一种具体实施方式，利用有限记忆拟牛顿插值法优化所述代价函数，指的是对代价函数进行求导，获得最优的变换解。

该实施方式中，基于图像的多分辨率金字塔的层数为三层，包括：将图像降采样，先降采样为原来图像的八分之一大小，进行配准得到一个最优解，将此次结果作为第二层配准的初始值；第二层是将图像降采样成原来的四分之一大小，进行配准求得此次最优解，将此作为第一层图像的配准起点，对原始图像进行配准从而得到最终的优化变换解。作为进一步优选，每层的最大迭代次数设为200，如果在200次迭代之前实现函数收敛，即求得本次迭代极值，如果到达最大的迭代次数200函数仍然没有收敛，那将最后一次的结果作为所求的解。对代价函数求导的公式为：

则，全变分的求导展开式为：

其中，

φ_i,j,k表示浮动图像的控制点，i，j，k值不同对应不同的控制点。

为了验证不连续运动的肝部三维多模态图像配准的精确度，进行单模肺部CT图像配准验证，并结合附图进行说明。

单模肺部CT图像配准时，采用呼气时的肺部CT图像作为参考图像，吸气时的肺部CT图像作为浮动图像，如图1(a)-1(b)所示。图2(a)-2(b)、图3(a)-3(b)分别为采用二范数和全变分作为正则项所得到的位移场及局部边界放大图。由图2(a)-2(b)可知，当采用二范数作为正则项时，从位移场的放大图2(b)中可以看到，在边界处可以得到一个连续的位移场并且肋骨和肺部的运动位移大小差别不大，这显然和实际运动时的位移场不相符；当采用本发明提供的全变分作为正则项时，从位移场放大图3(b)中可以看到，肺部是沿着胸腔壁进行运动，肺部的运动位移比较大而边界处位移比较小，边界处存在着一个明显的不连续的位移场，因此，上述对比结果表明，全变分作为正则项时，可以很好的保留器官滑移运动所带来的边界不连续性。

下面进一步给出其他三种应用在单模肺部CT图像行的配准方法和本发明提供的配准方法在单模肺部CT上所得的误差值比较。

表1三种单模肺部CT图像配准方法和本发明提供的配准方法的误差值比较

表格1中：1、Vandemeulebroucke指的是运动掩膜法；2、Pace指的是方向相关的正则项方法；3、Papiz(·)指的是基于Demons的自适应正则项法。这三种方法在单模肺部CT图像的配准上都取得了比较好的结果。表格1中，采用五组实验数据，配准精度指标为参考图像和浮动图像对应300个特征点之间的位移，这些特征点由专家手动提取。配准之后的值越小，说明配准的精度越高。将本发明提供的配准方法和其他三种方法来进行比较，在五组实验中，本发明提供多模的配准方法对应的结果中有三组是最好的，并且平均配准误差也最小。

接下来进行多模CT/MR肝部三维图像配准的实验，用本发明提供的基于全变分作为代价函数的正则项与传统采用二范数作为代价函数的正则项所得的位移场图进行对比。

图4(a)-4(b)为原始待配准的三维肝部CT和MR图像；图5为原始肝部三维多模态图像的融合图；图6(a)-6(b)为分别采用二范数和全变分作为正则项所得的配准后的融合图。图7(a)-7(b)为采用二范数作为正则项所得的肝部三维多模位移场图及局部边界放大图；图8(a)-8(b)为采用全变分作为正则项所得的肝部三维多模位移场图及局部边界放大图。从图6(a)-6(b)可以看出，两种方法所得到的配准结果都比较好。图8(a)-8(b)与图7(a)-7(b)对比可知，当用二范数作为约束项的时候，在肝部和腹壁的边界处存在着一个连续的位移场，位移场的运动穿过滑移边界，这个位移场不能真实的反应肝部的运动。当用全变分作为正则项的时，边界处的位移场可以得到一定程度的矫正，肝部沿着腹壁做滑移运动，肝部和腹壁的运动模态不同，所以应用此方法能够反应边界处运动的真实性。

尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用。它完全可以被适用于各种适合本发明的领域。对于熟悉本领域的人员而言可容易地实现另外的修改。因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。

Claims (7)

1.一种基于不连续运动的肝部三维多模态图像配准方法，其特征在于，其包括以下步骤：

对获取的原始肝部序列图像进行预处理和刚性配准；

采用结合三阶B样条函数的自由形变模型模拟刚性配准后肝部图像的弹性形变，采用基于全变分的正则项和混合测度来构建代价函数，衡量两幅图像的相似性程度；

采用基于图像的多分辨率金字塔的优化策略，利用有限记忆拟牛顿插值法来优化所述代价函数，求得最优的变换解；

所述混合测度，指的是基于归一化互信息的相关比，公式如下：

CRMI(M,F；φ)＝(2-NMI(M,F；φ))·(1-CR(M,F；φ))；

其中，M表示浮动图像，F表示参考图像，φ表示浮动图像的控制点，NMI表示归一化互信息，CR表示相关比；

在归一化互信息中，p(M；U)表示所述浮动图像的边缘概率密度，p(F)表示所述参考图像的边缘概率密度，p(M,F；U)表示两幅图像的联合概率密度；

在相关比中，X和Y分别表示参考图像和浮动图像的随机变量，Var[Y]为Y的方差，Var[Y-E(Y|X)]为Y独立于X方差；

所述基于全变分的正则项，公式如下：

其中，

是求导符号，U是所述浮动图像的位移，

和

分别表示形变场在三维空间三个方向上的导数，β是一个用于避免全变分函数在零处不可导的常数；

基于全变分的正则项和混合测度来构建的代价函数为：

C＝D(F(a),M(U(a)))+λR^TV(U)；

其中，F(a)为所述参考图像，M(a)为所述浮动图像，a为图像上的像素点，R^TV为正则项，λ为正则项的权重因子，D(F(a),M(U(a)))表示混合测度，U(a)为浮动图像的位移。

2.如权利要求1所述的基于不连续运动的肝部三维多模态图像配准方法，其特征在于，所述模拟弹性形变，包括以下步骤：

定义所述肝部三维CT图像为参考图像，定义所述肝部三维MR图像为浮动图像，分别加载所述参考图像和所述浮动图像；

对所述浮动图像进行弹性形变，并采用结合三阶B样条函数的自由形变模型来模拟所述浮动图像的弹性形变：设所述浮动图像的图像域为Φ，包含n_x×n_y×n_z个间距为δ的控制点，则，

U是所述浮动图像的位移，φ为所述浮动图像的控制点，(x,y,z)为所述浮动图像上的像素点；(p_x,p_y,p_z)是距离像素点最近的控制点，且

B_th表示三阶B样条的基函数。

3.如权利要求2所述的基于不连续运动的肝部三维多模态图像配准方法，其特征在于，所述浮动图像的控制点网格间距设为7*7*7mm，所述浮动图像上每个像素点的形变由周围64个控制点的形变决定。

4.如权利要求1所述的基于不连续运动的肝部三维多模态图像配准方法，其特征在于，所述λ的值为100。

5.如权利要求1-4中任一项所述的基于不连续运动的肝部三维多模态图像配准方法，其特征在于，所述预处理和刚性配准，包括以下步骤：

对获取的原始肝部序列图像进行三维重建和尺度归一化处理，输出肝部三维CT图像和肝部三维MR图像；

采用水平集的方法，分别对所述肝部三维CT图像进行去重处理、对所述肝部三维MR图像进行高斯平滑去噪处理；

对去重处理的所述肝部三维CT图像和去噪处理的所述肝部三维MR图像进行刚性配准，获得两幅边界位置大致对齐的肝部图像。

6.如权利要求2-4中任一项所述的基于不连续运动的肝部三维多模态图像配准方法，其特征在于，利用有限记忆拟牛顿插值法优化所述代价函数，指的是对所述代价函数进行求导，获得最优的变换解；

求导的公式为：

则，全变分的求导展开式为：

其中，

7.如权利要求1所述的基于不连续运动的肝部三维多模态图像配准方法，其特征在于，所述基于图像的多分辨率金字塔的层数为三层，包括：

将图像降采样，先降采样为原来图像的八分之一大小，进行配准得到一个最优解，将此次结果作为第二层配准的初始值；

第二层是将图像降采样成原来的四分之一大小，进行配准求得此次最优解，将此作为第一层图像的配准起点，对原始图像进行配准从而得到最终的优化变换解。

Images