一种基于高阶互累积量的序列医学图像配准方法

摘要

本发明属于医学图像处理领域，具体为一种基于高阶互累积量的序列医学图像配准方法。本发明提出了一种适用于序列图像配准的相似性测度--归一化高阶互累积量系数NHCC。分析表明NHCC测度不仅可以同时敏感多幅图像之间的相关程度，而且能免疫加性高斯噪声对配准的影响。在仿真结果验证了该测度的有效性和优越性的同时，也给出了一个采用序列图像配准后进行数字减影较好的实验结果。

说明书

技术领域

本发明属于医学图像处理领域，具体为一种基于高阶互累积量的序列医学图像配准方法。

背景技术

医学成像技术的蓬勃发展为临床提供了大量实用的解剖及功能方面的影像数据，序列医学图像是指在一定时间间隔内的多幅图像的集合。在临床中，序列医学影像可以反映出病人器官病灶位置的动态发展情况以及在不同时间点器官的状态，而通过比较序列图像之间的状态或位置变化可以用来检测病灶可能发生的潜在危险等。因而研究序列图像具有重要的理论和应用价值。但在时间序列图像的采样过程中，由于呼吸、心跳、吞咽、肌肉收缩、镜头震颤等影响，病人和镜头间不可避免地存在着相对运动，从而造成时间序列图像在空间位置上的差异。为了保证数据的可靠性和分析的精确度，必须对序列图像进行配准。

国内外许多学者致力于序列图像配准的研究，但现有解决序列图像配准问题的方法还主要是基于两幅图像的配准方法，即：选取序列图像中的一幅图像作为参考图像，利用通用的图像配准方法将其它图像通过两两配准到该参考图像上。由于配准中只使用了序列图像中一幅作为参考图像，随机选择的参考图像会使配准的结果出现随机性，这样的配准不具备科学性；此外，这种两两配准方法使得序列图像配准的计算复杂度大幅度增加。因此，人们一直致力于寻找一个能够同时敏感多幅图像之间位置关系的配准方法。

互信息(mutual information，MI)已经被广泛用作两幅医学图像配准的相似性测度，基于该测度的互信息配准方法，目前被公认为是配准精度和鲁棒性较好的配准方法之一。在互信息理论的基础上，Boes和Meyer提出利用高维互信息(High-Dimentional Mutual Information，HDMI，又称高阶互信息或广义互信息)同时配准三幅图像；Studholme等在对序列图像分割的研究过程中，也使用高维互信息来配准多幅图像；Wang和Shen改进了高维互信息的计算方法，并将其用于多幅超声图像的配准中。实践证明，高维互信息可以同时敏感多个变量之间的相关性，是一种性能较好的序列医学图像配准测度，受到了广泛关注。但高维互信息同样存在配准鲁棒性问题：即基于图像灰度统计量，当图像分辨率较低或者存在噪声干扰时，其鲁棒性较差，计算量较大。

在医学序列图像的形成过程中，众所周知，它们会受到加性高斯噪声和/或其它噪声的污染而导致成像的医学序列图像不同程度地降质。高维互信息在配准降质序列图像时，会因为这样的噪声污染而导致配准的鲁棒性下降。因此，如何在图像信噪比(Signal-to-Noise Ratio，SNR)较低的情况下保证相似性测度的鲁棒性，是序列配准研究中的重要课题。

发明内容

本发明的目的是提供一种基于高阶互累积量的序列医学图像配准方法，该方法可以同时敏感多幅图像之间的相关程度，而且能降低加性高斯噪声对配准的影响。

为达到上述发明目的，本发明提出了一种新的对高斯噪声鲁棒的相似性测度，称为归一化高阶互累积量系数(Normalized Higher-order Cross-cumulantCoefficient，NHCC)。

高阶累积量是高阶统计理论中的重要概念，一般指高于二阶的累积量。设{x_k}为实随机变量序列，其二阶、三阶和四阶累积量分别定义为：

cum(x₁，x₂)＝E[x₁，x₂]

cum(x₁，x₂，x₃)＝E[x₁，x₂，x₃]

cum(x₁，x₂，x₃，x₄)＝E[x₁，x₂，x₃，x₄]-E[x₁，x₂]E[x₃，x₄]-   (1)

E[x₁，x₃]E[x₂，x₄]-E[x₁，x₄]E[x₂，x₃]

其中，cum(·)、E[·]分别表示统计量的累积量和矩。由式(1)可知，{x_k}的三阶以下累积量和矩是等价的，但四阶以上两者不同。高阶累积量具有一些重要的性质，与本发明相关的性质包括：

性质1：半不变性：如果随机变量{x_k}与{y_k}彼此统计独立，则有

cum(x₁+y₁，x₂+y₂，...，x_k+y_k)

                                                (2)

＝cum(x₁，x₂，...，x_k)+cum(y₁，y₂，...，y_k)

性质2：高斯过程的N阶累积量(N＞2)恒为零。

高阶累积量的物理意义在于提供了随机过程{x_k}偏离高斯性的测度。

性质3：变元对称性：

$\mathrm{cum}(x_1,x_2,...,x_k)=\mathrm{cum}(x_{i_1},x_{i_2},...,x_{i_k})---\left(3\right)$

其中，(i₁，i₂，...，i_k)是(1，2，...，k)的任意一个排列。

对于含有加性高斯噪声的图像，利用该性质2，从理论上讲我们可以完全抑制加性高斯噪声对图像处理的影响。本发明定义了一种新的适用于K(K≥3)幅序列图像配准的归一化高阶互累积量系数(Normalized Higher-orderCross-cumulant Coefficient，记为NHCC^K)，设待配准的K幅序列医学图像分别记为g₁、g₂、…、g_K，则NHCC^K可表示为：

${\mathrm{NHCC}}^K=\frac{c_K(g_1(m,n),g_2(m,n),...,g_K(m,n))}{c_K(g_1(m,n),g_1(m,n),...,g_1(m,n))},K\ge 3---\left(4\right)$

其中，(m，n)为图像像素的坐标值，c_K(g₁(m，n)，g₂(m，n)，...，g_K(m，n))为K幅待配序列图像的K阶互累积量，c_K(g₁(m，n)，g₁(m，n)，...，g₁(m，n))为第一幅图像的K阶自累积量。此时，视第一幅图像为参考图像，当然我们也可以使用序列图像中的任一幅图像作为参考图像，例如：当选择第n幅图像时(n≤K)，式(4)的分母为c_K(g_n(m，n)，g_n(m，n)，...，g_n(m，n))。

下面分析式(4)对加性高斯噪声的免疫性：

设有K幅图像是医学成像设备对某个医学图像目标s在不同时间并具有不同偏移位置时获得的图像，即满足：

g₁(m，n)＝s(m，n)+w₁(m，n)；

g₂(m，n)＝s(m-Δx₂，n-Δy₂)+w₂(m，n)；            (5)

M

g_K(m，n)＝s(m-Δx_K，n-Δy_K)+w_K(m，n)；

式(5)中，(Δx_i，Δy_i)表示第i幅图像相对于第一幅图像的偏移量(这里第一幅图像被认为是参考图像，以第n幅图像为参考图像的情况可由类似的方法获得)，(m，n)为图像像素的坐标值，w_i(m，n)为零均值的加性高斯噪声，i＝1，2，...，K。根据累积量的定义和性质可知：

c_K(g₁(m，n)，g₁(m，n)，...，g₁(m，n))＝c_K(s(m，n)，s(m，n)，...，s(m，n))

c_K(g₁(m，n)，g₂(m，n)，...，g_K(m，n))                                    (6)

＝c_K(s(m，n)，s(m-Δx₂，n-Δy₂)，...，s(m-Δx_K，n-Δy_K))

则其对应的归一化高阶互累积量系数NHCC^K由式(4)知为：

${\mathrm{NHCC}}^K=\frac{c_K(g_1(m,n),g_2(m,n),...,g_K(m,n))}{c_K(g_1(m,n),g_1(m,n),...,g_1(m,n))}$  (7)

$=\frac{c_K(s(m,n),s(m-\Delta x_2,n-\Delta y_2),...,s(m-\Delta x_K,n-\Delta y_K))}{c_K(s(m,n),s(m,n),...,s(m,n))}$

式(7)表明，NHCC^K的计算结果只与原始图像之间的相对偏移量有关，而与加性高斯噪声无关。所以从理论上，NHCC^K可以完全消除高斯噪声的影响。当K幅图像均完全对准时，即(Δx_i，Δy_i)＝(0，0)，i＝2，...，K时，由式(8)知其NHCC^K取得最大值1。

两幅图像配准是上述序列图像配准的一个特例。当将NHCC^K测度用于两幅图像配准时，即K＝2时，NHCC^K退化为基于二阶累积量的测度，但是由性质2可知，二阶累积量对高斯噪声没有免疫性，为了使其在K＝2时也对加性高斯噪声具有免疫性，本发明单独定义用于两幅图像配准的归一化高阶互累积量系数(记为NHCC²)如下：

${\mathrm{NHCC}}^2=\frac{c_3(g_1(m,n),g_1(m,n),g_2(m,n))}{c_3(g_1(m,n),g_1(m,n),g_1(m,n))}---\left(8\right)$

其中，分母上的c₃(g₁(m，n)，g₁(m，n)，g₁(m，n))为g₁的三阶自累积量，分子上的c₃(g₁(m，n)，g₁(m，n)，g₂(m，n))为待配图像g₁、g₂的三阶互累积量。

令(Δx，Δy)表示g₂相对于g₁的偏移量，则：

${\mathrm{NHCC}}^2=\frac{c_3(s(m,n),s(m,n),s(m-\mathrm{\Delta x},n-\mathrm{\Delta y}))}{c_3(s(m,n),s(m,n),s(m,n))}---\left(9\right)$

式(9)表明，它的计算结果只与两幅图像的相对偏移量有关，NHCC²完全免疫高斯噪声的影响。而当两幅图像完全对准时，即(Δx，Δy)＝(0，0)，由式(10)知其NHCC²取得最大值1。类似的，该定义方法可以推广到大于三阶的累积量。

根据高阶累积量的性质3，即其变元对称性，NHCC²也可以定义为：

${\mathrm{NHCC}}^2=\frac{c_3(g_1(m,n),g_2(m,n),g_1(m,n))}{c_3(g_1(m,n),g_1(m,n),g_1(m,n))}$

或${\mathrm{NHCC}}^2=\frac{c_3(g_2(m,n),g_1(m,n),g_1(m,n))}{c_3(g_1(m,n),g_1(m,n),g_1(m,n))}$等等        (10)

式(10)的计算结果与式(8)相同。

本发明提出了一种全新的适用于序列图像配准的相似性测度--归一化高阶互累积量系数，不仅可以同时敏感多幅图像之间的相关程度，而且能免疫加性高斯噪声对配准的影响，对序列医学图像的配准具有重要的理论和应用价值。

附图说明

图1是相似性测度的性能比较。

图2是脑部序列图像。

图3是脑部序列图像配准的对比结果。

图4是DSA序列图像。

图5是未配准情况下的减影结果。

图6是采用MI配准后的减影结果。

图7是采用NHCC²配准后的减影结果。

图8是序列配准后的减影结果。

具体实施方案

以下结合具体的实施例，对本发明做进一步的阐述。实施例仅用于对本发明做说明而不是对本发明的限制。

实施例1：

两幅医学图像配准的性能对比

本实施例中，采用归一化互信息NMI和本发明提出的NHCC²分别来配准两幅脑部MR图像，以验证本发明在抗噪方面的有效性和优越性。

实验数据取自Mcgill大学的BrainWeb脑部图像数据库，待配准的图像为T2加权的脑部磁共振图像(T1、T2、PD)，是MR图像的三个重要参数，对应三种不同模式的MR图像。图像大小为256×256。

1、以图1(a)为参考图像，将图1(a)平移(10，10)个像素后的图像作为偏移图像，如图1(b)所示。分别采用NMI和NHCC²配准参考图像和偏移图像。

2、对图1(a)和图1(b)分别添加高斯噪声，作为待配准的降质图像，噪声的归一化均值及方差为(0，0.01)，如图1(c)和图1(d)所示。分别采用NMI和NHCC²配准这两幅图像。

对比结果如图1(e)和图1(f)所示，其中实线为无噪声情况下的配准结果，菱形线为有噪声情况下的配准结果。由图1(e)可见，NMI在无噪声情况下峰值尖锐，曲线平滑，但在存在噪声的情况下NMI受影响较严重，NMI数值围绕1.06上下波动，已经没有明显的峰值。而NHCC²在两种情况下的曲线几乎完全重合，有效抑制了高斯噪声的影响，与本发明的理论分析结果一致。

实施例2：

脑部序列图像配准

本实施例用以验证本发明提出的序列图像配准测度NHCC^K测度对噪声的免疫性。

待配原始图像来自于NLM-NIH虚拟人项目，为真实的PD加权的MR图像，因为脑部组织由颅骨固定，所以脑部序列图像之间的位移变换可以视作刚性变换。实验图像如图2所示，为三幅具有垂直偏移量的降质医学图像，故本实施例中K＝3，图像大小为256×256。分别采用NHCC³和HDMI配准该序列图像。实验结果如图3所示，图3(a)和图3(b)为NHCC³测度在不同的噪声水平下的评价曲面，图3(c)和图3(d)为HDMI测度在不同的噪声水平下的评价曲面，X坐标和Y坐标分别表示第二幅图/第三幅图相对于第一幅图的垂直偏移量。

由图3可见，在三幅图像完全对准的时候，即(X＝0，Y＝0)时，NHCC³取到最大值1，任何两幅图像间的偏差都会导致NHCC³迅速下降。同时，对比图3(a)和图3(b)可知，在不同的高斯噪声水平下，NHCC³测度的评价曲面变化很小，对高斯噪声的免疫性较强。而由图3(c)和图3(d)可见，HDMI在有噪声的情况下，鲁棒性较差，评价曲面的峰值不明显，最大值仅达到0.3左右。

实施例3：

数字减影序列图像配准的实验

数字减影血管造影技术中，未注射血管造影剂之前的X光片称为蒙片(maskimage)，注入血管造影剂之后拍摄的X光片称为活片(live image)。数字减影血管造影术就是希望通过蒙片和活片的减法操作，能够得到一幅血管特征清晰的数字减影图像。由于呼吸、心跳、吞咽、镜头震颤等影响，病人和镜头间会有相对运动。因此在进行减法操作之前，需要对蒙片和活片进行配准。

本实施例所采用的真实DSA序列图像如图4所示。图T1为蒙片图像，图T2和图T3为活片图像。由于在不同的时刻拍摄，序列图像在成像过程中会受到不同程度的噪声污染。同时，由于血流脉冲依次到达各处血管，各图呈现的血管信息也有所不同：图T2中区域A的血管成分较清晰，图T3中区域B的血管信息比较明显。那么，如果我们经过配准后结合多幅活片图像的血管信息，可以期待得到一幅更清晰更全面的减影结果。DSA图像整体变形多为非刚性的，但如果待配图像的区域比较小或者图像目标为骨骼部血管时，那么该图像活片相对于蒙片的变形可近似为刚性变换，故可以采用本发明方法进行配准。

本实施例首先给出未配准情况下的实验结果：图5(a)和(b)为未配准情况下对活片和蒙片直接进行数字减影的实验结果，其中图5(a)为图T2与图T1直接相减的图像，图5(b)为图T3与图T1直接相减的图像，图5(c)为将图5(a)和图5(b)通过简单相加融合后的实验结果。可见，相比于仅利用两幅图像进行减影的实验结果，经过简单相加融合后的减影结果能够提供更全面的血管信息，区域A和区域B(区域范围见图4中所标注)的血管成分都比较突出，但由于未进行配准，图5中三幅减影图像中均有较多的运动伪影。

接着，本实施例采用MI测度对实验图像进行两两配准后再进行数字减影，实验结果如图6所示：图6(a)和(b)为采用MI测度分别配准活片和蒙片后再进行减影的实验结果，其中图6(a)为对图T2与图T1配准减影的实验结果，图6(b)为对图T3与图T1配准减影的实验结果，图6(c)为将图6(a)和图6(b)通过简单相加融合后的实验结果。通过比较图5和图6可以看出，经过配准后减影的实验结果图6优于图5未配准减影的结果。不过图6的减影图像中仍然存在较多伪影，例如图6(c)中箭头所指示，MI配准的效果不够理想。

采用本发明提出的NHCC²测度对实验图像配准后再进行数字减影，实验结果如图7所示：图7(a)和(b)为采用本文方法配准活片和蒙片后再进行减影的实验结果，其中图7(a)为对图T2与图T1配准减影的实验结果，图7(b)为对图T3与图T1配准减影的实验结果，图7(c)为将图7(a)和图7(b)通过简单相加融合后的实验结果。由图可见，采用本发明方法配准后的减影图像比较平滑，特别是在图6(c)中箭头指示区域的伪影，在图7(c)中得到了有效抑制。

上述两组实验结果均基于两幅图像的配准方法。最后，本实施例给出了对图4的三幅图像进行序列配准后再进行数字减影的实验结果，实验结果如图8所示：图8(a)为采用HDMI测度配准该序列图像后再进行减影和简单相加融合的实验结果，图8(b)为采用本发明NHCC³测度配准该序列图像后再进行减影和简单融合的实验结果。由图8可见，相对于采用HDMI的减影图像，采用本发明方法的减影图像给出了一个更优的视觉效果，有效抑制了HDMI方法中存在的运动伪影。