非刚性医学图像配准算法的设计与实现

1、非刚性医学图像配准算法的设计与实现    【关键词】  医学图像；非刚性；图像配准；匹配矩阵；薄板样条摘要：非刚性图像匹配问题已成为医学图像分析中一个非常具有挑战性的问题。基于薄板样条插值方法,引入实匹配矩阵,并给出相应配准变换算法，该算法将薄板样条参数表示成仿射分量和非仿射分量，并分别进行求解。与其它非刚性匹配算法相比，该算法不仅保证了对应特征点的双向对应，也实现了自动特征点选择，实验结果令人满意。    关键词：医学图像；非刚性；图像配准；匹配矩阵；薄板样条1引言    在医学诊断

2、和治疗过程中，常需要对比分析多幅图像，以获得更为精确和全面的信息。图像分析大都要求多幅图像的几何位置一致，因此，配准是医学图像分析的一个重大课题。医学图像配准是指对于一幅医学图像寻求一种(或一系列)空间变换，使它与另一幅医学图像上的对应点达到空间上的一致。这种一致是指人体上的同一解剖点在两张匹配图像上有相同的空间位置。配准的结果应使两幅图像上所有的解剖点，或至少是所有具有诊断意义的点及手术感兴趣的点都达到匹配。图像配准不仅可以校正病人多次成像间的位置变化，也可以校正由于成像模式本身导致的畸变。对同一个病人的不同时间的图像进行配准，可以了解发育过程及肿瘤病变的病情；对不同人的图像进行配准，去除种

3、族、年龄等临床及遗传差异，从而形成疾病或人群特异性图谱，可用于正常与否的分析；对不同成像模式进行配准，可以获得互补信息。    医学图像配准可分为刚性配准和非刚性配准两类。刚性配准在许多情况下不能满足临床的需要，因为很多形变的性质是非刚体、非线性的。比如为了精确定位MR图像左心室，常常伴有组织磁化系数差异、非水分子的化学位移以及血流流动等因素导致的几何畸变以及由于磁场不均匀、磁场梯度非线性及涡流等导致的探测畸变，因此在放疗计划制定中，将MR图像配准时，不能单纯地使用刚性配准，必须使用非刚性配准。    非刚性配准算法可分为灰度驱动、模

4、型驱动及混合算法三种13。灰度驱动方法基于数学或统计尺度将一个灰度模式与另一个对准。典型情况下，需要定义源系统与目标系统之间的灰度相似性的数学量度。灰度相似性测度包括象素灰度的均方差、相关或互信息。模型驱动方法首先建立明确的几何模型，以此表示解剖标志。这些解剖标志包括有重要功能的表面、曲线和点。将源系统的解剖标志参数化，与目标系统的对应部分对准，以这种对应关系引导系统其余部分的变换。模型驱动算法包括点约束法、线约束法和面约束法。混合算法是结合使用以上两种算法的方法。薄板样条插值方法是非刚体变换中的一种特殊的变换，它允许局部调整，并符合某种连续性或平滑性要求。第2节讨论刚性能量函数；第3节给出非

5、刚性能量函数；第4节设计并实现一个非刚性配准算法；最后给出实验结果。    2刚性能量函数    本研究之所以采用薄板样条，是因为它的独特性质，就是能够将空间变换分解为一个全局仿射变换和一个局部非仿射变换。Booksteein4首先将薄板样条函数应用于标志点的匹配，结果证明它是一个非常有用的形状分析工具。假设在二维空间，已知两个具有N对对应点的点集，Q=Qi，i=1，2，n和P=Pi，i=1，2，n，将点集Q，P表示为：    Q=1   x1   y1 

6、   1   x2   y2              1   xn   ynP=1   x1   y1    1   x2   y2              1  

7、xn   yn    下面我们建立从点集P到点集Q的薄板样条映射f(Pi)，由于薄板样条是不对称的，因此从P到Q的映射不能简单地反转为从Q到P的映射。通过最小化下面的能量函数，可以得到一个刚性能量函数：    Etps(f)=n   i=1Q-f(P)2+J(f)(1)    其中，J(f)=R22 f   x22+22 f   xy2+2 f   2y2dxdy    (1)

8、式第一项代表经过变换的源标志点与目标标志点之间的距离和；第二项代表了获得的变换的不平滑度，也叫惩罚函数。使该式最小化的变换既满足变换后源标志点与目标标志点间接近(近似)的要求，同时也加入了足够的平滑。系数(>0)表征了近似和平滑之间的相对关系:当较小时，获得的变换表现了很好的近似效果;当较大时，就获得了比较平滑的变换，对较大的局部畸变进行了调整。薄板函数计算如下：    设z(x,y)=-U(r)=-r2logr2，其中，r=x2+y2，U(r)是构建薄板样条的基函数，设rij=|Pi-Pj|为点Pi与点Pj的欧几里德距离。对分散点数据集Pi进行薄板样条弹性

9、插值后可以得到曲面。插值过程形象地模拟为一个薄金属板在点约束下的扭曲变形，要使金属板在点(xi，yi)处高度为zi，并且该板具有最小弯曲能量，即薄板函数f(x,y)使罚函数J(f )最小。定义n×n矩阵：    K=0   U(r12)      U(r1n)    U(r21)   0      U(r2n)          

10、0;      U(rn1)   U(rn2)      0    V=(z(x1,y1),z(x2,y2),z(xn,yn)T    通过解线性方程组（2）可以得到W=(w1,w2,wn)T和T=(a1,ax,ay)T    KW+PT=V    PTW=0(2)    W是n×3的非仿射变换形变参数矩阵，T是3×3的仿射形变

11、参数矩阵，K是薄板样条的核，为n×n矩阵。    然后构造函数：     f(x,y)=a1,axx+ayy+n   i=1wiU(|(xi,yi)-(x,y)|)（3）    此时该函数对于所有i，有f(xi,yi)=zi，并使罚函数J(f)最小。    事实上，直接解方程组（2）是困难的，也不现实，我们将通过迭代求解点集之间的匹配矩阵来求方程（2）的参数W和T。    3非刚性能量函数  &#

12、160; 由刚性能量函数推导表明，只要已知两个点集之间的对应点，就可以得到它们之间的薄板样条映射参数。但是当对应点未知时，该如何处理呢?传统的方法往往都是手动选点,这种方法费时费力,同时在结构不清的情况下,很难选择到足够多的精确对应点。而且其准确性也只是相对的,误差是不可避免的。文献10定义两个点集之间的匹配矩阵M=Mij：    Mij=1，若点Qi对应于点Pi    0，其他    由于两个点集之间是双向一一对应的，即一个点集中的每个点在另一个点集中至多有一个对应点，反之亦然。匹配矩阵Mij具有下面约束：    j，N1   i=1Mij=1，i，N2   j=1Mij=1（4）    N1和N2分别是两个点集的点数，将匹配矩阵考虑到式(1)中，得到基于薄板样条映射的非刚性