CN102446358A - 基于边缘特征和cs信息的多模态医学图像配准方法 - Google Patents

Abstract

本发明公开一种基于边缘特征和CS信息的多模态医学图像配准方法。本发明基于广义形态滤波器原理，对形态学梯度算子进行改进，采用新的梯度滤波算子来进行图像边缘检测，从而提高图像边缘检测的运算速度和抗噪性能，快速实现图像边缘的粗检测，并确定图像边缘点的坐标和方向，为后继图像配准奠定基础；根据CS不等式，构造新的CS几何信息测度，本测度方法在数学计算过程中不涉及除法和对数运算，可有效降低图像配准的计算复杂度，提高计算效率；另外，本发明采用Powell方法进行配准参数的寻优，不必对目标函数求导，具有二次收敛性，收敛速度快，且便于编程实现。相比现有技术，本发明在计算标准差、信息量以及平均计算时间等方面性能更优，鲁棒性更好。

Description

基于边缘特征和CS信息的多模态医学图像配准方法

技术领域

本发明涉及一种医学图像配准方法，尤其涉及一种基于边缘特征和CS(Cauchy-Schwarz)信息的多模态医学图像配准方法。

背景技术

医学图像配准就是指将来自不同形式的成像设备(如MRI，CT，PET，SPECT等)的医学图像，利用计算机技术实现对于一幅医学图像寻求一种或者一系列的空间变换，使它与另一幅医学图像上的对应点达到空间上的一致。

随着计算机技术和电子技术的发展，MRI、PET、CT等各模态影像设备，逐步成为医院的常规诊断设备，医学图像已成为临床诊断与治疗不可缺少的最重要的手段。而如何最大限度地利用影像数据提供的有用信息，提高医生对疾病诊断的准确性，对临床诊断有重大的意义。医学图像配准的主要目的就是对在临床上广泛使用的诸如CT、MRI、PET等多源医学图像信息进行有机综合，最大限度地挖掘有用信息，以提高医学成像质量、诊断的准确性和治疗精度。另外，随着医学影像新的设备不断涌现，图像种类日益丰富，基于图像配准技术可大大提高医疗单位医学图像资源的共享和有效利用率。此外，医学图像配准作为信息技术的一个重要领域，是医学图像可视化的重要支撑技术之一，还可应用于更广泛的研究和工程领域，例如：机器视觉、模式识别、生物特征识别、信息理论和信息安全等众多领域，具有重要的理论意义和应用价值。

根据医学图像所提供的信息内涵，可将医学图像分为两大类：解剖结构图像(CT、MRI、B超等)和功能图像(SPECT、PET等)。这两类图像各有其优缺点：功能图像分辨率较差，但它提供了脏器功能代谢信息是解剖图像所不能替代的；解剖图像以较高的分辨率提供了脏器的解剖形态信息(功能图像无法提供脏器或病灶的解剖细节)，但无法反映脏器的功能情况。目前这两类成像设备的研究都已经取得了很大进展，图像的空间分辨率和图像质量有很大提高，但由于成像原理不同所造成的图像信息的局限性，使得单独使用某一类图像的效果并不理想，而多种图像的利用，又必须借助医生的空间想象和推测去综合判定他们所要的信息，其准确性受到主观影响，更主要的是一些信息将可能被忽略。因此，解决这个问题的最有效的方法就是以医学图像配准技术为基础，利用信息融合技术，将结构和功能图像结合起来，利用各自的信息优势在一副图像上同时表达来自人体的多方面信息，使人体内部的结构、功能等多方面的状况通过影像反映出来，从而更加直观地提供人体解剖、生理及病理等信息。

按照特征空间，图像配准一般分为两种，一种是基于图像特征的配准方法；另一种是基于灰度的配准方法。基于图像特征的图像配准方法，就是通过提取图像共同的特征，如点、线段和面等，并建立基于图像特征的目标函数，最后通过最优化搜索技术得到配准参数。该方法依赖于特征点的准确提取，易受噪声的影响，使得优化过程极易陷入局部极值。基于灰度的配准方法主要是对图像的灰度进行操作，缺点是计算量过大。

为了寻找图像间更加稳定的相似性测度来改善各种图像配准方法的稳定性，以克服图像配准存在的鲁棒性问题，将几何特征和像素相似性的方法相结合，实现两者优势互补，是一种有效和可行的思路。有的研究者利用互信息匹配形状特征点进行配准，有的结合区域边缘信息进行配准，有的利用互信息与曲线、曲面等方法相结合实现图像匹配，另外还有将图像的梯度信息作为空间信息加入到互信息中的方法，在一定的范围内也可以实现图像配准。但是，医学图像因受成像设备、噪声干扰等因素影响，与普通图像比较，本质上具有模糊性和不均匀的特点，且易受到人体解剖组织结构和形状不确定性知识的影响，普遍视觉效果较差。从配准算法自身的角度看，有其配准环节算法设计的问题。当前多模医学图像配准方法所存在的主要问题是：图像数据分割一直没有针对图像配准机制行之有效的方法，图像局部特征不变性缺乏，图像分割效果偏离生理视觉注意特点，且所构造的医学图像配准测度往往计算量过大、速度慢、鲁棒性差，难以满足图像配准的要求。

发明内容

本发明所要解决的技术问题在于克服现有技术所存在的不足，提供一种基于边缘特征和CS信息的多模态医学图像配准方法，该方法在计算标准差、信息量以及平均计算时间等方面性能更优，鲁棒性更好。

本发明的基于边缘特征和CS信息的多模态医学图像配准方法，包括以下步骤：

步骤A、对源图像和目标图像分别进行边缘特征提取，得到目标图像和源图像的边缘特征点集合X、Y；

步骤B、以源图像和目标图像的边缘特征点集合X、Y之间的距离测度D_CS(X，Y)作为目标函数，把求解距离测度D_CS(X，Y)最大值作为实现配准目标，对配准参数进行寻优，得到最终的配准参数，从而完成图像配准；其中，所述距离测度D_CS(X，Y)如下所示：

$D_{\mathrm{CS}}(X,Y)=\underset{i,j}{\mathrm{\Σ}}{p}_{\mathrm{ij}}^2\underset{i,j}{\mathrm{\Σ}}{\left(p_i{p}_j\right)}^2-{\left[\underset{i,j}{\mathrm{\Σ}}{p}_{\mathrm{ij}}^2\left(p_i{p}_j\right)\right]}^2,$

式中，p_i和p_j分别为X、Y的边缘概率分布，p_ij是X、Y的联合概率分布，i＝1，2L N₁，j＝1，2L N₂，N₁和N₂分别为待配准图像的边缘特征点集合X、Y中的边缘特征点个数。

进一步地，所述边缘特征提取具体按照以下方法：

首先选取两个十字形和交叉形的结构元素，并令其关于原点对称且互补，采用两者形态梯度的加权组合形式构造如下形式的梯度滤波算子：

FGrad(f)＝c₁IGrad1(f)+c₂IGrad2(f)

式中c₁和c₂为权系数，IGrad1和IGrad2为与两个结构元素相对应的改进形态学梯度；所述改进形态学梯度按照下式得到：

$\mathrm{IGrad}\left(f\right)=\left(\mathrm{fog}\right)\⊕g-(f\·g)\mathrm{\Θg}$

式中f(x，y)为原图像的图像灰度分布函数，g(x，y)为结构元素，

表示膨胀运算，Θ表示腐蚀运算；

然后利用上述梯度滤波算子提取图像的边缘特征，得到图像的边缘特征点集合。

优选地，因本发明只对二维平面医学图像作配准，所以只需求X轴方向的偏移量Δx、Y轴方向的偏移量Δy和绕图像中心的旋转角度Δθ三个参数。所述对图像配准参数进行寻优采用Powell方法，具体如下：

步骤B1、给定允许误差ε＞0，选取初始点X₀(0，0，0)和一组线性无关的搜索方向向量e⁽ⁱ⁾(i＝1，2，3)，其中e⁽¹⁾＝(1，0，0)，e⁽²⁾＝(0，1，0)，e⁽³⁾＝(0，0，1)，即令e⁽ⁱ⁾分别为沿坐标轴的方向；设置迭代计数器K＝1。

步骤B2、从X₀出发沿方向e⁽ⁱ⁾(i＝1，2，3)依次进行一维搜索，得到

X_(i)＝X_(i-1)+λ_i×e⁽ⁱ⁾    i＝1，2，3

f(X_(i))＝min f(X_(i-1)+λ_i×e⁽ⁱ⁾)    i＝1，2，3

在完成了这3次一维搜索后，得到X₍₃₎；

步骤B3、检验是否已满足终止准则：若|f(X_(i))-f(X₍₀₎)|≤ε或者迭代次数K大于等于M，则停止迭代，退出；其中，ε为预先设定的阈值e^-4，M是预设的最大迭代次数30；步骤B4、按照下式计算最速下降方向上的函数的变化量：

$\mathrm{\&Delta;f}=\underset{1<i<N}{\max }|f\left(X_{\left(i\right)}\right)-f\left(X_{\left(0\right)}\right)|,$

并记最速下降方向为e_max；

步骤B5、引进e＝X₍₃₎-X₍₀₎，Pcs＝2X₍₃₎-X₍₀₎，计算f(e)＝f(X₍₃₎-X₍₀₎)，f(Pcs)＝f(2X₍₃₎-X₍₀₎)，若f(e)≥f(x₍₀₎)或f(Pcs)＜f(x₍₀₎)，且2(f(x₍₀₎)-2f(x_(n))+f(e))[(f(x₍₀₎)-f(x_(n)))-Δf]²≥(f(x₍₀₎)-f(e))²Δf则将X₍₃₎赋值给X₍₀₎作为新的起点，沿上面的一组旧方向向量e⁽ⁱ⁾(i＝1，2，3)重复步骤B2至B5；若上述条件均不满足，则沿方向e＝X₍₃₎-X₍₀₎，以X₍₃₎作为起始点进行搜索得到目标函数在此方向上的极小值点P_min；将原来的下降方向最大的e_max去掉，而保留其它的2个方向，再加上方向e，仍得到3个方向：e⁽¹⁾，e⁽²⁾，e⁽³⁾，把此时的P_min赋给X₍₀₎作为新起点，重复步骤B2至B5。

优选地，所述距离测度D_CS(X，Y)中的联合概率分布p_ij根据下式计算：

$p_{\mathrm{ij}}(a,b)=\frac{h(a,b)}{\underset{a\&Element;X,b\&Element;Y}{\mathrm{\&Sigma;}}h(a,b)},$

式中，h(a，b)为待配准图像的二维联合直方图，具体表示X、Y重叠部分中灰度值为(a，b)对的像素对总数，其中a，b表示两图像中对应点灰度等级。

进一步地，所述边缘概率分布p_i和p_j分别根据以下公式计算得到：

$p_i=\underset{b}{\mathrm{\&Sigma;}}{p}_{\mathrm{XY}}(a,b),$ $p_j=\underset{a}{\mathrm{\&Sigma;}}{p}_{\mathrm{XY}}(a,b).$

本发明基于广义形态滤波器原理，对形态学梯度算子进行改进，采用新的梯度滤波算子来进行图像边缘检测，从而提高图像边缘检测的运算速度和抗噪性能，快速地实现图像边缘的粗检测，并确定图像边缘点的坐标和方向，为后继图像配准奠定基础；根据CS不等式，构造了新的CS几何信息测度，同已有的互信息测度等方法不同，本测度方法在数学计算过程中不涉及除法和对数运算，可有效降低图像配准的计算复杂度，提高计算效率；另外，本发明采用Powell方法进行配准参数的寻优，Powell法是解无约束优化问题行之有效的直接算法，它几何意义明确，不必对目标函数求导，具有二次收敛性，收敛速度快，可加快搜索最大值的速度，且便于编程实现。相比现有技术，本发明在计算标准差、信息量以及平均计算时间等方面性能更优，鲁棒性更好。

附图说明

图1为本发明图像配准方法的流程图；

图2为本发明梯度计算的原理示意图，其中，图2(a)为像素点(x，y)的对应邻域图，图2(b)为应用点(x，y)对应临域来计算f(x+1，y)-f(x，y)所用模板，图2(c)为相应计算f(x，y+1)-f(x，y)所用模板；

图3为效果验证实验中所使用的目标图像；

图4为效果验证实验中所使用的源图像。

具体实施方式

下面结合附图对本发明的技术方案进行详细说明：

本发明的基于边缘特征和CS信息的多模态医学图像配准方法，包括以下步骤：步骤A、对源图像和目标图像分别进行边缘特征提取，得到目标图像和源图像的边缘特征点集合X、Y。

本发明在进行图像边缘特征提取时，采用了一种新的梯度滤波算子，具体如下：

(1)设计一种新的边缘检测算子。

该形态学梯度算子具有如下形式：

$\mathrm{IGrad}\left(f\right)=\left(\mathrm{fog}\right)\&CirclePlus;g-(f\&CenterDot;g)\mathrm{\&Theta;g}$

式中f(x，y)为原始图像的图像灰度分布函数，g(x，y)为结构元素，表示膨胀运算，Θ表示腐蚀运算；

(2)构造新的梯度滤波算子

选取两个的十字形和交叉形的结构元素，并令其关于原点对称且互补，采用两者形态梯度的加权组合形式构造新的梯度滤波算子，该算子具有如下形式：

FGrad(f)＝c₁IGrad1(f)+c₂IGrad2(f)

式中c₁和c₂为权系数，IGrad1和IGrad2为与两个结构元素相对应的改进形态学梯度。

考虑到图像的灰度值函数f(x，y)是一个离散的二元函数，因此通常使用梯度来表示图像的一阶导数。所谓梯度就是有两个一阶导数组成的向量，像素点梯度G定义为：

$G=\left[\begin{array}{c}{G}_x\\ G_y\end{array}\right]=\left[\begin{array}{c}\frac{\&PartialD;f}{\&PartialD;x}\\ \frac{\&PartialD;f}{\&PartialD;y}\end{array}\right]$

梯度向量的模值G大小为：

$G=\mathrm{mag}\left(G\right)=\sqrt{G_x^2+G_y^2}=\sqrt{{\left(\frac{\&PartialD;f}{\&PartialD;x}\right)}^2+{\left(\frac{\&PartialD;f}{\&PartialD;y}\right)}^2}$

为了减少计算量，梯度向量模值G可用绝对值近似表示为：

$G\&ap;\left|G_x\right|+\left|G_y\right|=\left|\frac{\&PartialD;f}{\&PartialD;x}\right|+\left|\frac{\&PartialD;f}{\&PartialD;y}\right|$

梯度向量的方向角大小α(x，y)为：

$\mathrm{\&alpha;}(x,y)=\arctan \left(\frac{\mathrm{Gy}}{\mathrm{Gx}}\right)$

在医学图像的边缘检测技术中，梯度通常指的是梯度向量的模值。在具体实现过程中，有前面公式可知，根据导数和梯度的定义，可得出像素点(x，y)的梯度值G(x，y)：

$G(x,y)=\sqrt{G_x^2+G_y^2}=\sqrt{{\left(\frac{\&PartialD;f}{\&PartialD;x}\right)}^2+{\left(\frac{\&PartialD;f}{\&PartialD;y}\right)}^2}=\sqrt{{(f(x+1,y)-f(x,y))}^2+{(f(x,y+1)-f(x,y))}^2}$

式中的f(x+1，y)-f(x，y)可用点(x，y)应邻域图2(a)与图2(b)所示的模板进行空间卷积得到，模板中的左上角对应图像中的像素(x，y)。同理图2(c)所示的模板可以用来计算f(x，y+1)-f(x，y)。至此，一个像素点的梯度即可求出。而要计算整幅图像所有像素点的梯度，只需在图像上方移动模板，依次进行操作即可。

(3)对目标图像和源图像提取边缘特征点集

使用以上两步骤设计的梯度滤波算子，分别对目标图像和源图像提取边缘特征，得到包含边缘点E的向量组X＝{X_i，i＝1，2，...N₁}，Y＝{Y_i，i＝1，2，...N₂}，X_i和Y_i表示形状特征点在二维平面中的坐标，X和Y分别表示目标图像和源图像的边缘形状特征点集合，

其中，每个特征点的θ方向信息满足

目标图像特征点集合X中的特征点P_mi(E_i，θ_i)，点P_m(i+1)(E_i+1，θ_i+1)，点P_m(i+1)(E_i+2，θ_i+2)描述了在空间结构中，邻居点P_mi凡与邻近点P_m(i+1)在逆时针方向上邻接，而与邻居点P_m(i+2)在逆时针方向上间隔一个邻居点。同理，在源图像特征点集合Y中P_nj也与其邻居点具有同样的约束关系。因此对边缘特征点顺序进行分析等价于对特征邻居点在空间结构上的约束关系进行分析。

步骤B、以源图像和目标图像的边缘特征点集合X、Y之间的距离测度D_CS(X，Y)作为目标函数，以求解距离测度D_CS(X，Y)最大值作为图像配准目标，对配准参数进行寻优，得到最终的配准参数，从而完成图像配准。

本发明构造了新的CS几何距离测度。为便于公众理解，下面对其测度数学基础Cauchy-Schwarz不等式性质进行简要介绍。

Cauchy-Schwarz不等式定理：已知f，g是区间[a，b]上的连续函数，f，g∈C[a，b]则有， ${\left|{\&Integral;}_b^{a}f\left(x\right)g\left(x\right)\mathrm{dx}\right|}^2\&le;{\&Integral;}_b^a{\left|f\left(x\right)\right|}^2\mathrm{dx}{\&Integral;}_b^a{\left|g\left(x\right)\right|}^2\mathrm{dx}$

当且仅当f(x)与g(x)成比例是，等式成立。

从实变函数论的角度而言，我们仅需要求f，g是平方可积分函数(L²[a，b])，则该不等式仍然成立。N维实数空间Rⁿ的推广空间关系可表示为：

(即：

)

本发明所构造的X和Y的距离测度D_CS(X，Y)形式如下：

$D_{\mathrm{CS}}(X,Y)=\underset{i,j}{\mathrm{\&Sigma;}}{p}_{\mathrm{ij}}^2\underset{i,j}{\mathrm{\&Sigma;}}{\left(p_i{p}_j\right)}^2-{\left[\underset{i,j}{\mathrm{\&Sigma;}}{p}_{\mathrm{ij}}^2\left(p_i{p}_j\right)\right]}^2$

式中，p_i和p_j分别为X、Y的边缘概率分布，p_ij是X、Y的联合概率分布，i＝1，2L N₁，j＝1，2L N₂，N₁和N₂分别为X、Y中的边缘特征点个数。D_CS(X，Y)值的大小表示了目标图像和源图像的配准关联性程度，当两幅图像间的配准关联性增大时，p_ij和p_ip_j之间的差异相应增大，因而D_CS(X，Y)也增大，当D_CS(X，Y)取最大值时，两个随机变量之间的相似程度达到了最大，对两幅图像而言，则已经相互配准重合。与基于信息论的配准测度相比，CS几何信息测度在数学计算上因为不涉及除法和对数运算，其算法更加简单，运算效率更高。

上述两幅待配准图像X、Y的联合概率p_ij可以借助两幅图像的联合直方图h(a，b)进行计算。h(a，b)表示两幅图象重叠部分图像的灰度值为(a，b)对的象素对总数。其中a，b表示两图中灰度等级。h(a，b)可以用下面的矩阵表示：

$h=\left[\begin{array}{cccc}h\left(\mathrm{0,0}\right)& h\left(\mathrm{0,1}\right)& K& h(0,N-1)\\ h\left(\mathrm{1,0}\right)& h\left(\mathrm{1,1}\right)& K& h(0,N-1)\\ K& K& K& K\\ h(M-\mathrm{1,0})& h(M-\mathrm{1,1})& K& h(M-1,N-1)\end{array}\right]$

其中，（a，b)是表示两图中同一个位置的象素点，两幅图像中不重合的部分不参加计算，或是用图像背景灰度填充浮动图像因为几何变换而出现的空白区域。

利用联合直方图，两幅图像的联合概率分布为：

$p_{\mathrm{ij}}(a,b)=\frac{h(a,b)}{\underset{a\&Element;X,b\&Element;Y}{\mathrm{\&Sigma;}}h(a,b)}$

边缘概率分布p_i，p_j也可以通过联合概率密度求得：

$p_i=\underset{b}{\mathrm{\&Sigma;}}{p}_{\mathrm{XY}}(a,b),$ $p_j=\underset{a}{\mathrm{\&Sigma;}}{p}_{\mathrm{XY}}(a,b).$

从所构造的距离测度表达式可以看出，当两个随机变量相互独立时，即p_ij＝p_ip_j，这时D_CS(X，Y)值为0；当两幅图像间的配准关联性增大时，p_ij和p_ip_j之间的差异相应增大，因而D_CS(X，Y)也增大，当D_CS(X，Y)取最大值时，两个随机变量之间的相似程度达到了最大，对两幅图像而言，则已经实现配准。

为简化计算，计算过程中可按照两幅图像边缘特征点在空间结构上的约束关系，分别从目标图像和源图像边缘点E的向量组X，Y中提取点P_mi和P_nj，令配准参数ΔT＝(Δx，Δy，Δθ)，以此为起点计算两幅图像的D_CS(X，Y)值。

将距离测度D_CS(X，Y)作为目标函数，用Powell方法对配准参数ΔT＝(Δx，Δy，Δθ)进行寻优，Δx和Δy分别表示在二维坐标内，源图像沿X轴和Y轴的平移量，Δθ为源图像围绕图像中心的的旋转量，利用共轭方向并以此作为搜索方向，在一定迭代次数后输出最终的配准参数，从而完成图像配准。

Powell算法的基本思想是：把整个计算过程分成若干轮迭代，每轮迭代有n+1次(n为搜索函数的参数个数)一维搜索组成。在每轮迭代中，先从初始点X₀出发以此沿着已知的n个方向进行n次一维搜索，得到一个最好点X_best；接着从点X_best出发沿X₀与X_best连线方向进行1次一维搜索，得到本轮最好点；然后改变初始点，构造新的搜索方向，升始新一轮的迭代。

Powell算法的配准参数寻优具体步骤如下：

因为本发明只对二维平面医学图像作配准，所以只需求X轴方向的偏移量Δx、Y轴方向的偏移量Δy和绕图像中心的旋转角度Δθ三个参数。

步骤B1、给定允许误差ε＞0，选取初始点X₀(0，0，0)和一组线性无关的搜索方向向量e⁽ⁱ⁾(i＝1，2，3)，其中e⁽¹⁾＝(1，0，0)，e⁽²⁾＝(0，1，0)，e⁽³⁾＝(0，0，1)，即令e⁽ⁱ⁾分别为沿坐标轴的方向；设置迭代计数器K＝1。

步骤B2、从X₀出发沿方向e⁽ⁱ⁾(i＝1，2，3)依次进行一维搜索，得到

X_(i)＝X_(i-1)+λ_i×e⁽ⁱ⁾    i＝1，2，3

f(X_(i))＝min f(X_(i-1)+λ_i×e⁽ⁱ⁾)    i＝1，2，3

在完成了这3次一维搜索后，得到X₍₃₎；

步骤B3、检验是否已满足终止准则：若|f(X_(i))-f(X₍₀₎)|≤ε或者迭代次数K大于等于M，则停止迭代，退出；其中，ε为预先设定的阈值e^-4，M是预设的最大迭代次数30；

步骤B4、按照下式计算最速下降方向上的函数的变化量：

$\mathrm{\&Delta;f}=\underset{1<i<N}{\max }|f\left(X_{\left(i\right)}\right)-f\left(X_{\left(0\right)}\right)|,$

并记最速下降方向为e_max；

步骤B5、引进e＝X₍₃₎-X₍₀₎，Pcs＝2X₍₃₎-X₍₀₎，计算f(e)＝f(X₍₃₎-X₍₀₎)，f(Pcs)＝f(2X₍₃₎-X₍₀₎)，若f(e)≥f(x₍₀₎)或f(Pcs)＜f(x₍₀₎)，且2(f(x₍₀₎)-2f(x_(n))+f(e))[(f(x₍₀₎)-f(x_(n)))-Δf]²≥(f(x₍₀₎)-f(e))²Δ则将X₍₃₎赋值给X₍₀₎作为新的起点，沿上面的一组旧方向向量e⁽ⁱ⁾(i＝1，2，3)重复步骤B2至B5；若上述条件均不满足，则沿方向e＝X₍₃₎-X₍₀₎，以X₍₃₎作为起始点进行搜索得到目标函数在此方向上的极小值点P_min；将原来的下降方向最大的e_max去掉，而保留其它的2个方向，再加上方向e，仍得到3个方向：e⁽¹⁾，e⁽²⁾，e⁽³⁾，把此时的P_min赋给X₍₀₎作为新起点，重第步骤B2至B5。

为了验证本发明方法的有效性，选取如图3、图4的CT脑部图像，以图2作为目标图像，图3作为源图像，对这两组图像分别采用现有的同类MI算法，Rangarajan算法(特征点互信息算法)以及本发明方法进行配准。实验结果如表1和表2所示，其中表1是配准变换参数测量值间的比较表，表2是各算法均独立运行100次的标准差和信息测度的统计表。

表1 配准变换参数表

算法	Δx	Δy	Δθ
				MI算法	11.27	29.28	8.14
Rangarajan算法	11.18	28.87	8.22
				本发明方法	11.44	29.15	8.47

表2 3种配准算法对比表

算法	MI	Rangarajan	本发明方法
				标准差	7.567	6.124	5.225
信息量	3.722	3.451	3.824
				平均时间(s)	826.4	571.1	523.2

由表1和表2可以看出，三种算法的变换参数基本接近，但本发明方法在计算标准差、互信息量和平均计算时间等方面相对于同类MI算法和Rangarajan算法性能更优，鲁棒性更好。

Claims (5)

1.一种基于边缘特征和CS信息的多模态医学图像配准方法，其特征在于，包括以下步骤：

步骤A、对源图像和目标图像分别进行边缘特征提取，得到目标图像和源图像的边缘特征点集合X、Y；

步骤B、以源图像和目标图像的边缘特征点集合X、Y之间的距离测度D_CS(X，Y)作为目标函数，把求解距离测度D_CS(X，Y)最大值作为实现配准目标，对配准参数进行寻优，得到最终的配准参数，从而完成图像配准；其中，所述距离测度D_CS(X，Y)如下式所示：

$D_{\mathrm{CS}}(X,Y)=\underset{i,j}{\mathrm{\&Sigma;}}{p}_{\mathrm{ij}}^2\underset{i,j}{\mathrm{\&Sigma;}}{\left(p_i{p}_j\right)}^2-{\left[\underset{i,j}{\mathrm{\&Sigma;}}{p}_{\mathrm{ij}}^2\left(p_i{p}_j\right)\right]}^2,$

式中，p_i和p_j分别为X、Y的边缘概率分布，p_ij是X、Y的联合概率分布，i＝1，2L N₁，j＝1，2L N₂，N₁和N₂分别为待配准图像的边缘特征点集合X、Y中的边缘特征点个数。

2.如权利要求1所述基于边缘特征和CS信息的多模态医学图像配准方法，其特征在于，所述边缘特征提取具体按照以下方法：

首先选取两个十字形和交叉形的结构元素，并令其关于原点对称且互补，采用两者形态梯度的加权组合形式构造如下形式的梯度滤波算子：

FGrad(f)＝c₁IGrad1(f)+c₂IGrad2(f)

式中c₁和c₂为权系数，IGrad1和IGrad2为与两个结构元素相对应的改进形态学梯度；所述改进形态学梯度按照下式得到：

$\mathrm{IGrad}\left(f\right)=\left(\mathrm{fog}\right)\&CirclePlus;g-(f\&CenterDot;g)\mathrm{\&Theta;g}$

式中f(x，y)为原图像的图像灰度分布函数，g(x，y)为结构元素，

表示膨胀运算，Θ表示腐蚀运算；

然后利用上述梯度滤波算子提取图像的边缘特征，得到图像的边缘特征点集合。

3.如权利要求1所述基于边缘特征和CS信息的多模态医学图像配准方法，其特征在于，所述对图像配准参数进行寻优采用Powell方法，具体如下：

步骤B1、给定允许误差ε＞0，选取初始点X₀(0，0，0)和一组线性无关的搜索方向向量e⁽ⁱ⁾(i＝1，2，3)，其中e⁽¹⁾＝(1，0，0)，e⁽²⁾＝(0，1，0)，e⁽³⁾＝(0，0，1)，即令e⁽ⁱ⁾分别为沿坐标轴的方向；设置迭代计数器K＝1；

步骤B2、从X₀出发沿方向e⁽ⁱ⁾(i＝1，2，3)依次进行一维搜索，得到

X⁽ⁱ⁾＝X_(i-1)+λ_i×e⁽ⁱ⁾    i＝1，2，3

f(X_(i))＝minf(X_(i-1)+λ_i×e⁽ⁱ⁾)    i＝1，2，3

在完成了这3次一维搜索后，得到X₍₃₎；

步骤B3、检验是否已满足终止准则：若|f(X_(i))-f(X₍₀₎)|≤ε或者迭代次数大于等于M，则停止迭代，退出；其中，ε为预先设定的阈值e^-4，M是预设的最大迭代次数30；步骤B4、按照下式计算最速下降方向上的函数的变化量：

$\mathrm{\&Delta;f}=\underset{1<i<N}{\max }|f\left(X_{\left(i\right)}\right)-f\left(X_{\left(0\right)}\right)|,$

并记最速下降方向为e_max；

步骤B5、引进e＝X₍₃₎-X₍₀₎，Pcs＝2X₍₃₎-X₍₀₎，计算f(e)＝f(X₍₃₎-X₍₀₎)，f(Pcs)＝f(2X₍₃₎-X₍₀₎)，若f(e)≥f(x₍₀₎)或f(Pcs)＜f(x₍₀₎)，且2(f(x₍₀₎)-2f(x_(n))+f(e))[(f(x₍₀₎)-f(x_(n)))-Δf]²≥(f(x₍₀₎)-f(e))²Δf则将X₍₃₎赋值给X₍₀₎作为新的起点，沿上面的一组旧方向向量e⁽ⁱ⁾(i＝1，2，3)重复步骤B2至B5；若上述条件均不满足，则沿方向e＝X₍₃₎-X₍₀₎，以X₍₃₎作为起始点进行搜索得到目标函数在此方向上的极小值点P_min；将原来的下降方向最大的e_max去掉，而保留其它的2个方向，再加上方向e，仍得到3个方向：e⁽¹⁾，e⁽²⁾，e⁽³⁾，把此时的P_min赋给X₍₀₎作为新起点，重复步骤B2至B5。

4.如权利要求1所述基于边缘特征和CS信息的多模态医学图像配准方法，其特征在于，所述距离测度D_CS(X，Y)中的联合概率分布p_ij根据下式计算

$p_{\mathrm{ij}}(a,b)=\frac{h(a,b)}{\underset{a\&Element;X,b\&Element;Y}{\mathrm{\&Sigma;}}h(a,b)},$

式中，h(a，b)为待配准图像的二维联合直方图，具体表示X、Y重叠部分中灰度值为(a，b)对的像素对总数，其中a，b表示两图像中对应点灰度等级。

5.如权利要求4所述基于边缘特征和CS信息的多模态医学图像配准方法，其特征在于，所述距离测度D_CS(X，Y)中的边缘概率分布p_i和p_j分别根据以下公式计算得到： $p_i=\underset{b}{\mathrm{\&Sigma;}}{p}_{\mathrm{XY}}(a,b),$ $p_j=\underset{a}{\mathrm{\&Sigma;}}{p}_{\mathrm{XY}}(a,b).$

Images