CN102831599A - 一种明暗不均匀的医学图像的配准方法 - Google Patents

Abstract

本发明涉及一种明暗不均匀的医学图像的配准方法，该方法包括以下步骤：分别读入待配准的参考图像和浮动图像；用参考图像中每一像素点的灰度值减去浮动图像对应的像素点的灰度值，得到残差图像；将局部方差代入式（Ⅰ）算出与参考图像相对应的残差图像的每一像素点的权重；将残差图像与残差图像的每一像素点的权重点乘后替换式（Ⅱ）中的r，得到如式（Ⅲ）所示的求解配准图像的目标函数；采用梯度下降法迭代求解最小化的（Ⅲ）式所示的目标函数，当相邻两次迭代满足预设的终止条件时，即完成配准。本发明所述的配准方法具有平滑性好、配准精度高和鲁棒性强的优点。

Description

一种明暗不均匀的医学图像的配准方法

技术领域

本发明涉及通用图像数据处理,具体涉及在图像平面内的图形图像转换方法，该方法适用于明暗不均匀的医学图像配准。

背景技术

迄今为止，针对明暗不均匀的医学图像配准问题，主要有以下三类解决方法：

一是利用局部相似性测度。这种方法的理论依据可以直观的理解为：在小邻域范围内，可将图像视为明暗均匀的。Ardecani等人提出在以每个像素点为中心的小窗口内，使用局部相关系数作为相似性测度。为了克服陷入局部极值，使用多分辨率的块匹配。Hermosillo等人通过局部化相关比和互信息量，提出一个局部相似性测度的框架。Yi和Soatto考虑到局部归一化互信息对明暗不均匀不敏感，全局归一化互信息不容易陷入局部极值的特点，从而将二者结合起来并加权。这类基于局部相似性测度的方法虽能对明暗不均匀的图像产生较好的配准结果，但由于在小窗口内进行局部统计，将会使目标函数产生更多的局部极值。另外，此类方法对噪声和出格点比较敏感。

第二种方法是基于贝叶斯理论，使用更复杂的概率模型来解决明暗不均匀图像的配准。El-Baz等人提出使用马尔可夫-吉布斯随机场学习一幅图像的先验表观模型，再对第二幅图像进行变换使其概率最大。Wyatt和Noble提出使用马尔可夫随机场模型迭代地分割和配准图像的类标签。Zheng和Zhang提出使用基于灰度相似性测度的最大后验概率-马尔可夫随机场框架。这些基于马尔可夫随机场的方法严重依赖于局部灰度的关系和初始参数的选择。

第三种解决方法是同时进行灰度校正和图像配准。Frison提出使用非线性灰度变换和卷积滤波对图像进行灰度校正，再使用最小平方差方法进行图像配准，但针对不同的问题，必须手动选择不同的卷积滤波。Ashburner和Frison提出关于联合配准、灰度归一化和分割的概率框架，这种方法虽然可以产生较准确的配准结果，但需要手动选择参数且计算量大。AndriyMyronenko提出基于残差复杂度的相似性测度（Andriy Myronenko*and Xubo Song,Intensity-Based Image Registration by Minimizing Residual Complexity，IEEE Trans.Med.Imag.,vol.29,no.1,pp.1882-1891,Nov.2010），通过解析方法求取灰度校正场并推导出残差复杂度。与前边方法的不同之处在于，残差复杂度的表达式中不含灰度校正场，从而自适应地约束了校正场，保证了需要配准的两幅图在灰度空间上的一致。该方法根据给定的基函数－离散余弦变换，度量残差图像的稀疏性。当使用少量的基函数可以对残差图像进行稀疏编码时，残差复杂度达到最小。但是，Andriy Myronenko利用残差复杂度构造的目标函数

中的r包含了由噪声或其它偶然因素所产生的突变值，因此所输出的目标图像的平滑性不好，鲁棒性较差。

发明内容

本发明所要解决的技术问题的提供一种明暗不均匀的医学图像的配准方法，该方法具有平滑性好、配准精度高和鲁棒性强的优点。

本发明解决上述问题的技术解决方案是：

一种明暗不均匀的医学图像的配准方法，该方法包括以下步骤：

（1）分别读入待配准的参考图像和浮动图像；

（2）计算残差图像和残差图像的每一像素点的权重，其中，

所述残差图像的计算方法为：用参考图像中每一像素点的灰度值减去浮动图像对应的像素点的灰度值，得到残差图像；

所述残差图像的每一像素点的权重的计算方法为：分别以参考图像的每个像素点为中心，逐一求得每一像素点灰度值的局部方差，然后将所得到的局部方差代入下式（Ⅰ）算出与参考图像相对应的残差图像的每一像素点的权重：

$\mathrm{\ω}=e^{-\frac{V_i}{\underset{1}{\overset{N}{\mathrm{\Σ}}}{V}_i}/\mathrm{\ϵ}}---\left(I\right)$

式（Ⅰ）中，ω表示权重，V表示所述的局部方差，i为自然数，它表示参考图像中的任一像素点的编号，N表示参考图像的局部方差对应的局部块的数量，ε为等于0.01的常数，e为自然对数的底；

（3）将步骤（2）所得到的残差图像与残差图像的每一像素点的权重点乘后替换下式（Ⅱ）中的r，得到如下式（Ⅲ）所示的求解配准图像的目标函数：

$E\left(\mathrm{\τ}\right)=-\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\log \frac{1}{{\left(q_n^{T}r\right)}^2/\mathrm{\α}+1}---\left(\mathrm{II}\right)$

$E=-\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\log \frac{1}{{\left(q_n^T\left(\mathrm{r\ω}\right)\right)}^2/\mathrm{\α}+1}---\left(\mathrm{III}\right)$

式（Ⅱ）和（Ⅲ）中，q表示正向离散余弦变换，r表示步骤（2）所述的残差图像，α为等于0.0001的常数，T表示转置变换，n为浮动图像中像素点的数量，其余与式（Ⅰ）相同；

（4）采用梯度下降法迭代求解最小化的（Ⅲ）式所示的目标函数，当相邻两次迭代满足预设的终止条件时，即完成配准。

上述方案中的局部块为大小是3×3、5×5或7×7的图块。本发明人实验显示所述局部块的大小对配准后所得到的图像质量影响不大，因此，为了减小运算量，取3×3局部块为佳。

上述式（Ⅰ）中ε的取值应避免权值过于接近而违背实际情况，选择ε=0.01是本发明人多次实验的优化结果。

上述式（Ⅱ）中a越小，残差复杂度越易陷入局部极值，因此选择α=0.0001可以验证本发明算法的鲁棒性。

本发明人反复试验，上述步骤（4）所述的终止条件为相邻两次迭代的能量函数的改变率小于10^-5。

为了提高计算的准确性并加快程序收敛，上述技术方案中的步骤（1）所述的参考图像和浮动图像在进行后续处理前，最好先进行下述预处理：将参考图像和浮动图像由无符号整数类型转化为双精度浮点类型；将参考图像和浮动图像的每一像素点的灰度值线性变换为0~1。

将本发明所述式（Ⅲ）与式（Ⅱ）[Andriy Myronenko所构建的目标函数]比较可见，本发明所述式（Ⅲ）相当于用rω替换了Andriy Myronenko所构建的目标函数中的r。再由式（Ⅰ）可见，ω的值在0~1范围内变化，而且与局部方差成反比例变化。由此可知，本发明所述的技术方案在变换过程中将图像残差值比较大的地方给予小的权重约束，而将图像残差值比较小的地方给予大的权重约束，即相当于滤除了残差图像中r－rω的部分，从而使输出的目标图像平滑性好，明暗程度相对均匀。此外，所滤除掉的r－rω的部分都是一些由噪声或其它偶然因素所产生的突变，因此显著提高了抗噪性能和鲁棒性。

附图说明

图1为本发明所述配准方法的流程图。

图2为曾患肝肿瘤瘘的病人经治疗后的腹部CT增强扫描，其中，（b）图和（a）图为相隔数十秒钟依次扫描的增强图像。

图3为图2中（a）图的局部方差图；

图4为用图2（a）中每一像素点的灰度值减去图2（b）对应的像素点的灰度值所得到的残差图像。

图5为以图2（a）为参考图像，图2（b）为浮动图像采用本发明所述方法配准后输出的变换后的浮动图像。

图6为以图2（b）为参考图像，图2（a）为浮动图像采用本发明所述方法配准后输出的变换后的浮动图像。

图7为图2（a）与图5的绝对残差图。

图8为图2（b）与图6的绝对残差图。

图9为以图2（a）为参考图像，图2（b）为浮动图像采用Andriy Myronenko所述方法配准后输出的变换后的浮动图像。

图10为图2（a）与图9对应的绝对残差图。

具体实施方式

例1

本实施例以曾患肝肿瘤瘘的病人经治疗后进行的腹部CT增强扫描图像为例描述本发明所述配准方法。参见图1，具体配准过程如下所述：

步骤1，分别读入图2（a）和图2（b），并将图2（a）作为参考图像，图2（b）作为浮动图像。为了提高计算的准确性，并加快程序收敛，接着先将所读入的图2（a）和图2（b）分别由无符号整数类型转化为双精度浮点类型，再将它们的灰度值线性变换为0~1。

步骤2，分别以图2（a）中每个像素点为中心，逐一求得每一像素点8邻域内灰度值的局部方差，即局部块大小为3×3的局部方差，所得到的局部方差图如图3所示。将所得到的局部方差代入式

[即发明内容中所述的式（Ⅰ）]，算出与参考图像相对应的残差图像的每一像素点的权重，便得到下式所示的权重矩阵 $W=\left[\begin{array}{ccccc}{\mathrm{\ω}}_{11}& {\mathrm{\ω}}_{12}& {\mathrm{\ω}}_{13}& ...& {\mathrm{\ω}}_{1n}\\ {\mathrm{\ω}}_{21}& {\mathrm{\ω}}_{22}& {\mathrm{\ω}}_{23}& ...& {\mathrm{\ω}}_{2n}\\ {\mathrm{\ω}}_{31}& {\mathrm{\ω}}_{32}& {\mathrm{\ω}}_{33}& ...& {\mathrm{\ω}}_{3n}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ {\mathrm{\ω}}_{m1}& {\mathrm{\ω}}_{m2}& {\mathrm{\ω}}_{m3}& ...& {\mathrm{\ω}}_{\mathrm{mn}}\end{array}\right].$

步骤3，用参考图像的灰度矩阵减去浮动图像的灰度矩阵得到残差矩阵R，该残差矩阵R所对应的残差图像如图4所示。用r_ij表示残差矩阵R中的每个元素，R可表达成如下所示的m行×列n的矩阵 $R=\left[\begin{array}{ccccc}{r}_{11}& r_{12}& r_{13}& ...& r_{1n}\\ r_{21}& r_{22}& r_{23}& ...& r_{2n}\\ r_{31}& r_{32}& r_{33}& ...& r_{3n}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ r_{m1}& r_{m2}& r_{m3}& ...& r_{\mathrm{mn}}\end{array}\right].$

步骤4，将步骤2和3所得到的残差图像与残差图像的每一像素点的权重点乘后替换式（Ⅱ）[即发明内容中所述的式（Ⅱ）]中的r，得到式（Ⅲ）[即发明内容中所述的式（Ⅲ）]所示的目标函数。上述式（Ⅲ）所示残差图像与残差图像的每一像素点的权重点乘结果为，

$\mathrm{r\ω}=R\×W=\left[\begin{array}{ccccc}{r}_{11}{\mathrm{\ω}}_{11}& r_{12}{\mathrm{\ω}}_{12}& r_{13}{\mathrm{\ω}}_{13}& ...& r_{1n}{\mathrm{\ω}}_{1n}\\ r_{21}{\mathrm{\ω}}_{21}& r_{22}{\mathrm{\ω}}_{22}& r_{23}{\mathrm{\ω}}_{23}& ...& r_{2n}{\mathrm{\ω}}_{2n}\\ r_{31}{\mathrm{\ω}}_{31}& r_{32}{\mathrm{\ω}}_{32}& r_{33}{\mathrm{\ω}}_{33}& ...& r_{3n}{\mathrm{\ω}}_{3n}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ r_{m1}{\mathrm{\ω}}_{m1}& r_{m2}{\mathrm{\ω}}_{m2}& r_{m3}{\mathrm{\ω}}_{m3}& ...& r_{\mathrm{mn}}{\mathrm{\ω}}_{\mathrm{mn}}\end{array}\right].$

步骤5，采用梯度下降法迭代求解最小化的式（Ⅲ）所示的目标函数，当相邻两次迭代的能量函数的改变率小于10^-5时终止迭代，完成配准，并输出变换后的浮动图像，如图5所示。

例2

为了验证本发明所述方法的有效性，本实验例将例1所使用的参考图像与浮动图像对换一下，即以图2（b）为参考图像，图2（a）为浮动图像按例1所述的方法进行配准，输出变换后的浮动图像如图6所示。

在此使用绝对残差图像观察实验结果，具体过程如下所述：用图2（a）中每一像素点的灰度值减去图5对应的像素点的灰度值，其差值取绝对值得到如图7所示的绝对残差图像。用图2（b）中每一像素点的灰度值减去图6对应的像素点的灰度值，其差值取绝对值得到如图8所示的绝对残差图像。

将图7与图8比较可见，将参考图像与浮动图像对换后，二者的效果相当，本发明仍然适用。

例3

为了验证本发明所述方法的效果，本实施例将例1所得到的图5与采用现有技术所得到的变换后浮动图像进行视觉对比。

将图2（a）作为参考图像，图2（b）作为浮动图像，采用Andriy Myronenko在ntensity-BasedImage Registration by Minimizing Residual Complexity一文中所述的方法进行配准，所输出的图像如图9所示。为了便于观察实验结果，在此使用常用绝对残差色图进行视觉对比，具体过程如下所述：

用图2（a）中每一像素点的灰度值减去图9对应的像素点的灰度值，其差值取绝对值得到如图10所示的绝对残差图像。然后，将图10与图7进行视觉比较。比较的结果是，在图10中白色框内患者的肝、脾组织边缘处的像素点的灰度值较高，表现为白色，而图7中白色框内患者的肝、脾组织的边缘的白色像素点数量急剧减少，这说明本发明所述配准方法所输出的变换后的浮动图像的效果明显优于现有技术。

Claims (3)

1.一种明暗不均匀的医学图像的配准方法，该方法包括以下步骤：

（1）分别读入待配准的参考图像和浮动图像；

（2）计算残差图像和残差图像的每一像素点的权重，其中，

所述残差图像的计算方法为：用参考图像中每一像素点的灰度值减去浮动图像对应的像素点的灰度值，得到残差图像；

所述残差图像的每一像素点的权重的计算方法为：分别以参考图像的每个像素点为中心，逐一求得每一像素点灰度值的局部方差，然后将所得到的局部方差代入下式（Ⅰ）算出与参考图像相对应的残差图像的每一像素点的权重：

$\mathrm{\ω}=e^{-\frac{V_i}{\underset{1}{\overset{N}{\mathrm{\Σ}}}{V}_i}/\mathrm{\ϵ}}---\left(I\right)$

式（Ⅰ）中，ω表示权重，V表示所述的局部方差，i为自然数，它表示参考图像中的任一像素点的编号，N表示参考图像的局部方差对应的局部块的数量，ε为等于0.01的常数，e为自然对数的底；

（3）将步骤（2）所得到的残差图像与残差图像的每一像素点的权重点乘后替换下式（Ⅱ）中的r，得到如下式（Ⅲ）所示的求解配准图像的目标函数：

$E\left(\mathrm{\τ}\right)=-\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\log \frac{1}{{\left(q_n^{T}r\right)}^2/\mathrm{\α}+1}---\left(\mathrm{II}\right)$

$E=-\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\log \frac{1}{{\left(q_n^T\left(\mathrm{r\ω}\right)\right)}^2/\mathrm{\α}+1}---\left(\mathrm{III}\right)$

式（Ⅱ）和（Ⅲ）中，q表示正向离散余弦变换，r表示步骤（2）所述的残差图像，α为等于0.0001的常数，T表示转置变换，n为浮动图像中像素点的数量，其余与式（Ⅰ）相同；

（4）采用梯度下降法迭代求解最小化的（Ⅲ）式所示的目标函数，当相邻两次迭代满足预设的终止条件时，即完成配准。

2.根据权利要求1所述的一种明暗不均匀的医学图像的配准方法，其特征在于，步骤（1）还包括将所读入的参考图像和浮动图像进行下述预处理的步骤：将参考图像和浮动图像由无符号整数类型转化为双精度浮点类型；将参考图像和浮动图像的每一像素点的灰度值线性变换为0~1。

3.根据权利要求1或2所述的一种明暗不均匀的医学图像的配准方法，其特征在于，步骤（4）所述的终止条件为相邻两次迭代的能量函数的改变率小于10^-5。

Images