CN100454337C - 一种结合空间信息的图像配准方法 - Google Patents

Abstract

本发明涉及医学图像配准技术领域，提出了一种基于互信息并结合像素邻域空间信息的配准方法，对受噪声污染严重的图像有很好的效果。该方法假定图像像素灰度信息和空间信息构成的二维向量呈正态分布，将像素的邻域平均值作为空间信息。在计算互信息表达式中的概率密度时，根据多维变量正态分布的性质，构造出两幅待配准图像的对应向量矩阵。然后只需要计算该矩阵的均值和协方差即可直接计算出两幅图像的边缘熵和联合熵，从而得到它们的互信息。实验结果证明，该方法对于噪声较大的图像有很好的效果，计算速度快，节省内存空间，可以广泛地用于医学图像配准领域中，有重要的临床应用价值和研究意义。

Description

一种结合空间信息的图像配准方法

技术领域

本发明涉及医学图像配准技术研究领域，特别是一种结合空间信息的图像配准方法。

背景技术

医学图像配准是医学图像分析的重要步骤。为了充分利用各种成像方式提供的信息，医生希望能将不同图像的信息进行融合，这就需要将不同图像的相同部位进行对齐，从而更为全面地掌握病变部位的情况。

自从1995年Collignon等人^[1]和Viola及Wells等人^[2]成功地将互信息用在医学图像配准中，它就作为相似性度量被广泛的应用到多模态图像配准当中。由于互信息不要求待配准图像间存在线性关系，并且不受两幅图像特征的相关性限制，因此它在多模态图像配准中作为相似性度量有其独到的优势。它衡量的是两个变量间的相似性程度，即一个变量包含另一个变量的信息量的多少。换言之，两个变量越相似，则其中一个变量包含到的另外一个变量的信息量越大。

然而，互信息作为相似性度量也很容易陷入由插值误差所造成的局部极值中。它的一个最大不足就是缺乏空间信息。在互信息的表达式中，只考虑了图像灰度的信息而忽略了像素之间在空间中的相互关系。它只是对整幅图像的灰度分布做了统计分析，并没有考虑它们之间的位置关系。于是，很多人在如何在互信息中加入空间信息方面做了相关的工作。Pluim等人^[3]将梯度信息作为空间信息加入到互信息中；Rueckert等人^[4]提出了一种扩展互信息的框架，将二阶熵形式的互信息成功的应用到了非刚体配准中；Daniel B.Russakoff等人^[5]提出了基于区域的互信息，该方法认为每个像素和它邻域像素的灰度值构成一个高维向量，该向量服从高斯分布；与之相似，Dejan Tomazevie等人^[6]提出了基于多维特征的互信息。

在[4]中，作者将每个像素的相邻像素灰度值作为联合概率的第二维变量。由于该相邻像素只是代表了按照某一特定顺序位置的信息，因此加入的空间信息是远远不够的。为此，我们将每个像素的邻域平均值作为第二维变量，能够使得空间信息更为丰富。同时，我们的方法与[5][6]相比，更节省内存空间，具有很好的抗噪声能力。

参考文献：

1.A.Collignon，F.Maes，D.Delaere，D.Vandermeulen，P.Suetens，andG.Marchal，″Automated multi-modality image registration based oninformation theory″，in Information Processing in Medical Imaging，pp.263-274，Norwell，MA：Kluwer，1995.

2.P.Viola and W.M.Wells，″Alignment by maximization of mutualinformation″，in Proc.5th Int.Conf.Computer Vision，pp.16-23，Boston，MA，1995.

3.J.P.W.Pluim，J.B.A.Maintz，and M.A.Viergever，“Image Registration byMaximization of Combined Mutual Information and Gradient Information，”IEEE Transactions on Medical Imaging，vol.19，no.8，pp.809-814，2000.

4.D.Rueckert，M.J.Clarkson，D.L.G.Hill，and D.J.Hawkes，“Non-rigidregistration using higher-order mutual information，”In Proc.SPIE MedicalImaging 2000：Image Processing，pp.438-447，San Diego，CA，2000.

5.D.B.Russakoff，C.Tomasi，T.Rohlfing，and C.R.Maurer，“ImageSimilarity Using Mutual Information of Regions，”In Proc.of the 8th EuropeanConference on Computer Vision(ECCV)，LNCS 3023，pp.596-607，2004.

6.D.Tomazevie，B.Likar，F.Pernus，“Multi-Feature Mutual Information，”Proc.SPIE，vol.5370，pp.143-154，2004.

发明内容

本发明利用像素的邻域平均值作为空间信息，将该信息加入作为相似性度量函数的互信息中。该方法假定图像像素灰度信息和空间信息构成的二维向量呈正态分布，在计算互信息表达式中的概率密度时，根据多维变量正态分布的性质，构造出两幅待配准图像的对应向量矩阵，然后只需要计算该矩阵的均值和协方差即可直接计算出两幅图像的边缘熵和联合熵，从而得到它们的互信息。该方法可以很好地用于多模态医学图像配准中。对于受噪声污染严重的多模态图像配准有很好的效果。

在计算两幅待配准图像的互信息时，概率密度函数的估计是首要的核心步骤。传统的直方图统计方法只考虑了图像像素的灰度统计分布，并未考虑像素之间的空间相互关系。此外，基于高阶熵的互信息需要用到多维直方图，会占用很大的内存空间。本发明中，认为图像像素的灰度信息和空间信息所构成的二维向量服从正态分布。根据多维变量的正态分布性质，可以通过计算对应向量矩阵的均值和协方差来得到概率密度。

在实数域R上，对于均值为u、方差为σ²的正态分布变量，其熵可以用下式直接计算得到：

$-{\∫}_R\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}{e}^{-(1/2){((x-\mathrm{\μ})/\mathrm{\σ})}^2}(-\log \left(\sqrt{2\mathrm{\π}}\mathrm{\σ}\right)-\frac{1}{2}{\left(\frac{x-\mathrm{\μ}}{\mathrm{\σ}}\right)}^2)\mathrm{dx}---\left(1\right)$

类似地，呈正态分布的n维变量的高斯密度为

$G\left(x\right)=\frac{1}{\sqrt{{\left(2n\right)}^n\det \mathrm{\Σ}}}{e}^{-(1/2)(x-\mathrm{\μ})\·{\mathrm{\Σ}}^{-1}(x-\mathrm{\μ})}---\left(2\right)$

其熵为

$\frac{1}{2}(n+\log \left({\left(2\mathrm{\π}\right)}^n\det \mathrm{\Σ}\right))$

(3)

附图说明

图1是对应像素对向量矩阵的构造示意图。

图2是实验中用到的三幅图像。

图3是不同相似性度量函数的配准曲线图。

图4是不同相似性度量下每一步的搜索时间图。

图5是本发明的基于互信息并结合像素邻域空间信息的配准方法流程图。

具体实施方式

图2中，(a)为功能磁共振图像及其伪影图像；(b)为猩猩及其噪声图像；(c)为核磁共振图像及其噪声图像。在每组图像中，左侧的是目标图像，右侧的是浮动图像。

由于配准问题实际上是一个优化问题，因此配准目标函数曲线的质量决定了配准结果的好坏。曲线的光滑度和最大峰值处的尖锐程度是影响搜索结果的重要因素。

图3所示的各曲线图中，横坐标表示图像沿水平(或者竖直)方向平移的像素数目或者图像绕中心旋转的角度；纵坐标表示原图像和变换后的图像之间的相似性度量函数值。其中平移范围从左(或者上)20个像素到右(或者下)20个像素，旋转角度从负20度到正20度。

通过比较可以看出，SMI和AHMI的方法做出的曲线很粗糙，出现了很多局部极值，将SMI用于平滑预处理后的图像仍然也很粗糙。相反，用我们的方法可以得到很平滑的曲线。特别是我们的方法还可以校正对于图3(c)第一行配准曲线出现的最大值位置偏离真值的问题。

图3中，(a)，(b)，(c)分别是对应于图2中三组图像的相似性度量函数曲线。在每一组中，第一行为浮动图像沿水平位置平移后与目标图像的相似性度量曲线；第二行为浮动图像沿竖直位置平移后与目标图像的相似性度量曲线；第三行为浮动图像绕中心旋转后与目标图像的相似性度量曲线。

图4给出了不同相似性度量下每一步搜索所需要的时间。虽然基于区域的互信息(RMI)考虑到了邻域区域内的每一个像素，但是在保证精度的同时我们的计算方法速度更快，并且节省了内存空间。

如图5所示，我们可以按照以下流程来进行计算：

步骤一(S1)：构造两幅待配准图像的对应向量矩阵V；

图中阴影部分分别为像素a₀和b₀在3*3窗口内的相邻像素。I(a_i)和I(b_i)分别是a_i和b_i的灰度值。对于给定的图像A和B，构造所有像素对的对应向量 $\stackrel{\→}{v_{\mathrm{AB}}}(i,j)={(\stackrel{\→}{v_a},\stackrel{\→}{v_b})}_{\mathrm{ij}}$ (见图1)，其中 $\stackrel{\→}{v_a}=(v_{a1},v_{a2}),$ v_a1＝I(a₀)， $v_{a2}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}I\left(a_i\right);$ $\stackrel{\→}{v_b}=(v_{b1},v_{b2}),$ v_b1＝I(b₀)， $v_{b2}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}I\left(b_i\right).$ 由此我们得到一个向量矩阵V。假定两幅待配准图像有相同的大小，为l×m，对于给定的邻域窗口半径r，总的像素对数目为(l-2r)×(m-2r)。

步骤二(S2)：计算矩阵V，V_A和V_B的协方差；

首先计算向量中每个分量的均值，比如第一维分量 $\stackrel{\‾}{v_{A1}}=\frac{1}{(l-2r)\×(m-2r)}\underset{i,j=(r+1,r+1)}{\overset{(l-r,m-r)}{\mathrm{\Σ}}}{v}_{a1}(i,j),$ 然后从V中减去相应的均值得到V′。V的协方差为 $\mathrm{Cov}\left(V\right)=\frac{1}{(l-2r)\×(m-2r)}{V}^{\′}{V}^{\′T}.$ 同样很容易可以分别得到图像A的向量矩阵V_A的协方差阵Cov(A)和图像B的协方差阵Cov(B)。步骤三(S3)：计算图像的熵和目标函数的值；

由公式(3)，图像A的边缘熵为 $H\left(A\right)=1+\frac{1}{2}\log \left({\left(2\mathrm{\π}\right)}^2\right|\mathrm{\ΣA}\left|\right).$ 这里|∑A|为矩阵Cov(A)的行列式的值。同样可以计算出H(B)和H(A，B)以及目标函数值 $\mathrm{MI}(A,B)=H\left(A\right)+H\left(B\right)-H(A,B)=\frac{1}{2}\log \frac{\left|\mathrm{\ΣA}\right|\·\left|\mathrm{\ΣB}\right|}{\left|\mathrm{\ΣV}\right|}.$

下面举例说明本发明中方法的有效性。

对于一些质量很差的图像，传统的互信息不能精确地使这些图像配准(见图2)。我们比较了几种相似性度量函数配准曲线。为了与高阶熵互信息[4]加以比较，我们将某特定排列顺序的相邻像素用邻域平均值来替换作为联合熵中的第二维变量，这里简称为AHMI(average high-ordermutual information)的方法。传统互信息用SMI(standard mutualinformation)来表示。为了排除本发明的方法中涉及到的求邻域均值平滑因素，我们还与SMI方法用于平滑后的图像的情况作了比较(见图3)。

本发明中所述的图像配准方法很适合受噪声污染严重的图像。该方法鲁棒性强、计算速度快，具有重要的临床应用价值和研究意义。实验结果证明，该方法对于噪声较大的图像有很好的效果，计算速度快，节省内存空间，可以广泛地用于医学图像配准领域中，有重要的临床应用价值和研究意义。

Claims (2)

1.一种基于互信息并结合像素邻域空间信息的配准方法，其特征是，将像素的邻域平均值作为空间信息，将该信息加入作为相似性度量函数的互信息中，该方法假定图像像素灰度信息和空间信息构成的二维向量呈正态分布，在计算互信息表达式中的概率密度时，根据多维变量正态分布的性质，构造出两幅待配准图像的对应向量矩阵，然后只需要计算该矩阵的均值和协方差即可直接计算出两幅图像的边缘熵和联合熵，从而得到它们的互信息。

2.根据权利要求1的基于互信息并结合像素邻域空间信息的配准方法，其具体步骤如下：

步骤S1，构造两幅待配准图像的对应向量矩阵V；

步骤S2，计算矩阵V的协方差阵，和图像A的向量矩阵VA的协方差阵及图像B的协方差阵；

步骤S3，计算图像的熵和目标函数的值。

Images