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

1非刚性医学图像配准算法的设计与实现【关键词】 医学图像；非刚性；图像配准；匹配矩阵；薄板样条摘要：非刚性图像匹配问题已成为医学图像分析中一个非常具有挑战性的问题基于薄板样条插值方法,引入实匹配矩阵,并给出相应配准变换算法，该算法将薄板样条参数表示成仿射分量和非仿射分量，并分别进行求解与其它非刚性匹配算法相比，该算法不仅保证了对应特征点的双向对应，也实现了自动特征点选择，实验结果令人满意关键词：医学图像；非刚性；图像配准；匹配矩阵；薄板样条1 引言在医学诊断和治疗过程中，常需要对比分析多幅图像，以获得更为精确和全面的信息图像分析大都要求多幅图像的几何位置一致，因此，配准是医学图像分析的一个重大课题医学图像配准是指对于一幅医学图像寻求一种(或一系列)空间变换，使它与另一幅医学图像上的对应点达到空间上的一致这种一致是指人体上的同一解剖点在两张匹配图像上有相同的空间位置配准的结果应使两幅图像上所有的解剖点，或至少是所有具有诊断意义的点及手术感兴趣的点都达到匹配图像配准不仅可以校正病人多次成像间的位置变化，也可以校正由于成像模式本身导致的畸变对同一个病人的不同时间的图像进行配准，可以了解发育过程及肿瘤病变的病情；对不同人的图像进行配准，去除种族、年龄等临床及遗传差异，从而形成疾病或人群特异性图谱，可用于正常与否的分析；对不同成像模式进行配准，可以获得互补信息。

2医学图像配准可分为刚性配准和非刚性配准两类刚性配准在许多情况下不能满足临床的需要，因为很多形变的性质是非刚体、非线性的比如为了精确定位MR 图像左心室，常常伴有组织磁化系数差异、非水分子的化学位移以及血流流动等因素导致的几何畸变以及由于磁场不均匀、磁场梯度非线性及涡流等导致的探测畸变，因此在放疗计划制定中，将 MR 图像配准时，不能单纯地使用刚性配准，必须使用非刚性配准非刚性配准算法可分为灰度驱动、模型驱动及混合算法三种［1~3］灰度驱动方法基于数学或统计尺度将一个灰度模式与另一个对准典型情况下，需要定义源系统与目标系统之间的灰度相似性的数学量度灰度相似性测度包括象素灰度的均方差、相关或互信息模型驱动方法首先建立明确的几何模型，以此表示解剖标志这些解剖标志包括有重要功能的表面、曲线和点将源系统的解剖标志参数化，与目标系统的对应部分对准，以这种对应关系引导系统其余部分的变换模型驱动算法包括点约束法、线约束法和面约束法混合算法是结合使用以上两种算法的方法薄板样条插值方法是非刚体变换中的一种特殊的变换，它允许局部调整，并符合某种连续性或平滑性要求第 2 节讨论刚性能量函数；第 3 节给出非刚性能量函数；第 4 节设计并实现一个非刚性配准算法；最后给出实验结果。

2 刚性能量函数本研究之所以采用薄板样条，是因为它的独特性质，就是能够将空间变换分解为一个全局仿射变换和一个局部非仿射变换Booksteein［4］首先将薄板样条函数应用于标志点的匹配，结果证明它是一个非常有用的形状分析工具假设在二维空3间，已知两个具有 N 对对应点的点集，Q={Qi，i=1，2，…，n}和P={Pi，i=1，2，…，n}，将点集 Q，P 表示为：Q=1 x1 y11 x2 y2… … …1 xn ynP=1 x1 y11 x2 y2… … …1 xn yn下面我们建立从点集 P 到点集 Q 的薄板样条映射 f(Pi)，由于薄板样条是不对称的，因此从 P 到 Q 的映射不能简单地反转为从 Q 到 P 的映射通过最小化下面的能量函数，可以得到一个刚性能量函数：Etps(f)=∑n i=1‖Q-f(P)‖2+λJ(f)(1)其中，J(f)=R22 f x22+22 f xy2+2 f 2y2dxdy(1)式第一项代表经过变换的源标志点与目标标志点之间的距离和；第二项代表了获得的变换的不平滑度，也叫惩罚函数使该式最小化的变换既满足变换后源标志点与目标标志点间接近(近似)的要求，同时也加入了足够的平滑。

系数λ(λ0)表征了近似和平滑之间的相对关系:当 λ 较小时，获得的变换表现了很好的近似效果;当 λ 较大时，就获得了比较平滑的变换，对较大的局部畸变进行了调整薄板函数计算如下：设 z(x,y)=-U(r)=-r2logr2，其中，r=x2+y2，U(r)是构建薄板样条的基函数，设4rij=|Pi-Pj|为点 Pi 与点 Pj 的欧几里德距离对分散点数据集 Pi 进行薄板样条弹性插值后可以得到曲面插值过程形象地模拟为一个薄金属板在点约束下的扭曲变形，要使金属板在点(xi，yi)处高度为 zi，并且该板具有最小弯曲能量，即薄板函数 f(x,y)使罚函数 J(f )最小定义 n×n 矩阵：K=0 U(r12) … U(r1n)U(r21) 0 … U(r2n)… … … …U(rn1) U(rn2) … 0V=(z(x1,y1),z(x2,y2),…,z(xn,yn))T通过解线性方程组（2）可以得到 W=(w1,w2,…,wn)T 和 T=(a1,ax,ay)TKW+PT=VPTW=0(2)W 是 n×3 的非仿射变换形变参数矩阵，T 是 3×3 的仿射形变参数矩阵，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 和 T3 非刚性能量函数5由刚性能量函数推导表明，只要已知两个点集之间的对应点，就可以得到它们之间的薄板样条映射参数但是当对应点未知时，该如何处理呢?传统的方法往往都是手动选点,这种方法费时费力,同时在结构不清的情况下,很难选择到足够多的精确对应点而且其准确性也只是相对的,误差是不可避免的文献［10］定义两个点集之间的匹配矩阵 M={Mij}：Mij=1，若点 Qi 对应于点 Pi0，其他由于两个点集之间是双向一一对应的，即一个点集中的每个点在另一个点集中至多有一个对应点，反之亦然匹配矩阵{Mij}具有下面约束：j，∑N1 i=1Mij=1，i，∑N2 j=1Mij=1（4）N1 和 N2 分别是两个点集的点数，将匹配矩阵考虑到式(1)中，得到基于薄板样条映射的非刚性匹配的能量函数为：Etps(M，T，W)=∑N1 i=1 ∑N2 j=1 Mij‖Q-PT-KW‖2+λJ(f)(5)式(5)的第一项,可以使点集 P 中的点尽可能近地映射到 Q 中的点；第二项是平滑性约束,用于映射的调整，调整参数 λ 决定了映射的形变程度,当 λ→0 时，将得到对应点的精确匹配。

4 非刚性配准算法的设计与实现在保证式(4)的约束下,放松对匹配矩阵的约束,将二值的匹配矩阵转化为连续实数矩阵,即 Mij∈{0,1}→Mij∈［0,1］,允许部分匹配的存在,称这样的匹配矩阵为模糊匹配矩阵由第下面的算法可以看到，随着时间的推移，Mij 的值逐渐变大，越来6越接近二值矩阵，当时间足够长时，就会得到最终的二值匹配矩阵根据匹配矩阵元素的性质，令：Mij=1 ‖Q-PT-KW‖2+1(6)当‖Q-PT-KW‖→0 时，Mij→1所以，非刚性能量函数（4）式可改写为：Etps(M，T，W)=∑N1 i=1 ∑N2 j=1 (Mij2(‖Q-PT-KW‖2+1)+2Mij)+λJ(f)(7)通过 E Mij=0 能够得到使能量函数(7)式极小的匹配矩阵元素 Mij在第 2 节我们知道，直接解方程组（2）中参数 W 和 T 是困难的，对于仿射变换 T的计算是独立于（2）式因为二维欧氏空间上的仿射变换可写为：S(Pj)=TPj+A，其中，A=(Δx,Δy)T 为平移量，T=kcosθ sinθ-sinθ cosθ,平移、旋转、缩放及反射和剪切等是二维仿射变换的特例此模型中的参数 k、θ、Δx 和 Δy，即为两图像的配准参数。

确定这几个参数的步骤为：首先对需配准的两幅图像估计初始值 k0、θ0、Δx0 和 Δy0，建立两个点集的坐标对应关系计算两幅图像对应点互信息，可得到 Δx0 和 Δy0对 k0、θ0 的选取可先确定一个大致范围，然后设定一定的间隔 μ 作步长因子，设 k=μ1 k0，θ=μ2θ0，以互信息最大为原则进行迭代搜索，自适应取得最佳值实际中图像经过预处理后，图像之间的旋转角 θ 比较小，取值范围为(-π/4，π/4)，就能保证找到正确的 θ 值按照最大相关原则迭代搜索，以获得最佳值 θ，对 k0 的处理方法与此类似获得了最佳值 k 和 θ 后，再对 Δx0 和 Δy0 按最大互信息原则沿图像两个正交方7向逐像素搜索，以取得最佳值 Δx 和 Δy确定了配准参数 K、θ、Δx 和 Δy，就可对图像进行平移、旋转和缩放将得到的T 代入（2）式，求得参数 WW 的作用是将图像经 T 变换后坐标值不落在像素点上的点调整到像素点上本研究的算法主要包括以下几个步骤：① 给定特征点集 Q 和 P；② 初始化：Mij=1 (全 1 矩阵)，W=1 (全 1 矩阵), T=0（零矩阵）, A=0（零向量），λ=λ0，N=100，构造初始薄板样条；③ 根据仿射坐标最大相关原则迭代计算 T，A；④ 根据（4）式计算 W；⑤ 根据（6）式计算 Mij；⑥ 根据（3）式，构造薄板样条；⑦ 如果 M 满足（4）式约束或迭代次数大于N ，则转⑧， 否则转③；⑧ end。

需要指出的是,由于互信息量是统计量,因此我们对标准图确定的点数不能太少，保证互信息量的统计有效性5 实验结果8图 1 显示了对一对 96×96,灰度级为 256 级的 MRI 图像进行非刚性配准的实例a)是一幅标准的 MRI 心脏长轴的一个切片；(b)是一幅有变形的 MRI 图像心脏长轴同一周期另一个切片,此图像是用图像处理软件 Photoshop 手动变形处理获得；(c)是标准图像中特征点的选取,图中的小红圆圈为选择的特征点；(d)是最大互信息搜索方法在变形图像上搜索到的对应特征点；(e)是配准结果由图 1(a)、(b)、(e)我们可以看出,经过配准变换后,图 1(b)变形为图 1(e),与图 1(a)已基本完全一致了通过计算整幅图像之间的互信息量我们得到:图 1(a)和图 1(b)的互信息量为 1.037,图 1(a)和图 1(e)之间的互信息量为 3.44,而图 1(a)与自身计算得到的互信息量为3.46可见,经过非刚性性配准变换,图像达到了较高的匹配精度abcde图 1MRI 非刚性配准的实例6 结束语一直以来，学者们提出了各种非刚性图像匹配方法文献［5］将图像分解为许多子图像，估计每个子图像之间仿射变换，这样用多个仿射变换近似反映整个图像之间的非刚性映射。

文献［6］首先从图像中提取出特征点，然后将特征点拟合成曲线或曲面进行匹配这种方法在拟合的曲线或曲面比较光滑的情况下效果好，但是当图像所包含的形状很复杂时，曲线或曲面的拟合就变得非常困难，这种方法的好处是曲线的匹配相对于点集的匹配要容易些最近还有许多非刚性匹配的研究主要集中于非刚性形状的统计特性的学习上［7］，其方法类似于迭代最近点(ICP)算法，是基于局部的启发式搜索以上这些方法大9都存在鲁棒性差，而且通常不能保证图像之间的一一对应和特征点的自动确定本研究提出了一种通过薄板样条函数来表征特征点集之间的非刚性映射，把该映射分解为仿射变换和非仿射变换，并分别计算求解薄板样条的参数并满足双向对应的约束参考文献1Likar B, pernui F.A hierarchical approach to elastic registration based on mutual in formation. Imag Vision Comput，2001，19:33~44.2Kyriacou S K， Davatzikos C. No。