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

1、非刚性医学图像配准算法的设计与实现【关键词】医学图像；非刚性；图像配准；匹配矩阵；薄板样条摘要：非刚性图像匹配问题已成为医学图像分析中一个非常具有挑战性的问题。基于薄板样条插值方法,引入实匹配矩阵,并给出相应配准变换算法，该算法将薄板样条参数表示成仿射分量和非仿射分量，并分别进展求解。与其它非刚性匹配算法相比，该算法不仅保证了对应特征点的双向对应，也实现了自动特征点选择，实验结果令人满意。关键词：医学图像；非刚性；图像配准；匹配矩阵；薄板样条1引言在医学诊断和治疗过程中，常需要比照分析多幅图像，以获得更为准确和全面的信息。图像分析大都要求多幅图像的几何位置一致，因此，配准是医学图像分析的一个重

2、大课题。医学图像配准是指对于一幅医学图像寻求一种(或一系列)空间变换，使它与另一幅医学图像上的对应点到达空间上的一致。这种一致是指人体上的同一解剖点在两张匹配图像上有一样的空间位置。配准的结果应使两幅图像上所有的解剖点，或至少是所有具有诊断意义的点及手术感兴趣的点都到达匹配。图像配准不仅可以校正病人屡次成像间的位置变化，也可以校正由于成像形式本身导致的畸变。对同一个病人的不同时间的图像进展配准，可以理解发育过程及肿瘤病变的病情；对不同人的图像进展配准，去除种族、年龄等临床及遗传差异，从而形成疾病或人群特异性图谱，可用于正常与否的分析；对不同成像形式进展配准，可以获得互补信息。医学图像配准可分为

3、刚性配准和非刚性配准两类。刚性配准在许多情况下不能满足临床的需要，因为很多形变的性质是非刚体、非线性的。比方为了准确定位R图像左心室，常常伴有组织磁化系数差异、非水分子的化学位移以及血流流动等因素导致的几何畸变以及由于磁场不均匀、磁场梯度非线性及涡流等导致的探测畸变，因此在放疗方案制定中，将R图像配准时，不能单纯地使用刚性配准，必须使用非刚性配准。非刚性配准算法可分为灰度驱动、模型驱动及混合算法三种13。灰度驱动方法基于数学或统计尺度将一个灰度形式与另一个对准。典型情况下，需要定义源系统与目的系统之间的灰度相似性的数学量度。灰度相似性测度包括象素灰度的均方差、相关或互信息。模型驱动方法首先建立

4、明确的几何模型，以此表示解剖标志。这些解剖标志包括有重要功能的外表、曲线和点。将源系统的解剖标志参数化，与目的系统的对应部分对准，以这种对应关系引导系统其余部分的变换。模型驱动算法包括点约束法、线约束法和面约束法。混合算法是结合使用以上两种算法的方法。薄板样条插值方法是非刚体变换中的一种特殊的变换，它允许部分调整，并符合某种连续性或平滑性要求。第2节讨论刚性能量函数；第3节给出非刚性能量函数；第4节设计并实现一个非刚性配准算法；最后给出实验结果。2刚性能量函数本研究之所以采用薄板样条，是因为它的独特性质，就是可以将空间变换分解为一个全局仿射变换和一个部分非仿射变换。BksteEin4首先将薄板

5、样条函数应用于标志点的匹配，结果证明它是一个非常有用的形状分析工具。假设在二维空间，两个具有N对对应点的点集，Q=Qi，i=1，2，n和P=Pi，i=1，2，n，将点集Q，P表示为：Q=1x1y11x2y21xnynP=1x1y11x2y21xnyn下面我们建立从点集P到点集Q的薄板样条映射f(Pi)，由于薄板样条是不对称的，因此从P到Q的映射不能简单地反转为从Q到P的映射。通过最小化下面的能量函数，可以得到一个刚性能量函数：Etps(f)=ni=1Q-f(P)2+J(f)(1)其中，J(f)=R22fx22+22fxy2+2f2y2dxdy(1)式第一项代表经过变换的源标志点与目的标志点之间

6、的间隔 和；第二项代表了获得的变换的不平滑度，也叫惩罚函数。使该式最小化的变换既满足变换后源标志点与目的标志点间接近(近似)的要求，同时也参加了足够的平滑。系数(0)表征了近似和平滑之间的相对关系:当较小时，获得的变换表现了很好的近似效果;当较大时，就获得了比拟平滑的变换，对较大的部分畸变进展了调整。薄板函数计算如下：设z(x,y)=-U(r)=-r2lgr2，其中，r=x2+y2，U(r)是构建薄板样条的基函数，设rij=|Pi-Pj|为点Pi与点Pj的欧几里德间隔 。对分散点数据集Pi进展薄板样条弹性插值后可以得到曲面。插值过程形象地模拟为一个薄金属板在点约束下的扭曲变形，要使金属板在点(

7、xi，yi)处高度为zi，并且该板具有最小弯曲能量，即薄板函数f(x,y)使罚函数J(f)最校定义nn矩阵：K=0U(r12)U(r1n)U(r21)0U(r2n)U(rn1)U(rn2)0V=(z(x1,y1),z(x2,y2),z(xn,yn)T通过解线性方程组2可以得到=(1,2,n)T和T=(a1,ax,ay)TK+PT=VPT=0(2)是n3的非仿射变换形变参数矩阵，T是33的仿射形变参数矩阵，K是薄板样条的核，为nn矩阵。然后构造函数：f(x,y)=a1,axx+ayy+ni=1iU(|(xi,yi)-(x,y)|)3此时该函数对于所有i，有f(xi,yi)=zi，并使罚函数J(f

8、)最校事实上，直接解方程组2是困难的，也不现实，我们将通过迭代求解点集之间的匹配矩阵来求方程2的参数和T。3非刚性能量函数由刚性能量函数推导说明，只要两个点集之间的对应点，就可以得到它们之间的薄板样条映射参数。但是当对应点未知时，该如何处理呢?传统的方法往往都是手动选点,这种方法费时费力,同时在构造不清的情况下,很难选择到足够多的准确对应点。而且其准确性也只是相对的,误差是不可防止的。文献10定义两个点集之间的匹配矩阵=ij：ij=1，假设点Qi对应于点Pi0，其他由于两个点集之间是双向一一对应的，即一个点集中的每个点在另一个点集中至多有一个对应点，反之亦然。匹配矩阵ij具有下面约束：j，N1

9、i=1ij=1，i，N2j=1ij=14N1和N2分别是两个点集的点数，将匹配矩阵考虑到式(1)中，得到基于薄板样条映射的非刚性匹配的能量函数为：Etps(，T，)=N1i=1N2j=1ijQ-PT-K2+J(f)(5)式(5)的第一项,可以使点集P中的点尽可能近地映射到Q中的点；第二项是平滑性约束,用于映射的调整，调整参数决定了映射的形变程度,当0时，将得到对应点的准确匹配。4非刚性配准算法的设计与实现在保证式(4)的约束下,放松对匹配矩阵的约束,将二值的匹配矩阵转化为连续实数矩阵,即ij0,1ij0,1,允许部分匹配的存在,称这样的匹配矩阵为模糊匹配矩阵。由第下面的算法可以看到，随着时间的

10、推移，ij的值逐渐变大，越来越接近二值矩阵，当时间足够长时，就会得到最终的二值匹配矩阵。根据匹配矩阵元素的性质，令：ij=1Q-PT-K2+1(6)当Q-PT-K0时，ij1。所以，非刚性能量函数4式可改写为：Etps(，T，)=N1i=1N2j=1(ij2(Q-PT-K2+1)+2ij)+J(f)(7)通过Eij=0可以得到使能量函数(7)式极小的匹配矩阵元素ij。在第2节我们知道，直接解方程组2中参数和T是困难的，对于仿射变换T的计算是独立于2式。因为二维欧氏空间上的仿射变换可写为：S(Pj)=TPj+A，其中，A=(x,y)T为平移量，T=kssin-sins,平移、旋转、缩放及反射和剪

11、切等是二维仿射变换的特例。此模型中的参数k、x和y，即为两图像的配准参数。确定这几个参数的步骤为：首先对需配准的两幅图像估计初始值k0、0、x0和y0，建立两个点集的坐标对应关系。计算两幅图像对应点互信息，可得到x0和y0。对k0、0的选取可先确定一个大致范围，然后设定一定的间隔作步长因子，设k=1k0，=20，以互信息最大为原那么进展迭代搜索，自适应获得最正确值。实际中图像经过预处理后，图像之间的旋转角比拟小，取值范围为(-/4，/4)，就能保证找到正确的值。按照最大相关原那么迭代搜索，以获得最正确值，对k0的处理方法与此类似。获得了最正确值k和后，再对x0和y0按最大互信息原那么沿图像两个

12、正交方向逐像素搜索，以获得最正确值x和y。确定了配准参数K、x和y，就可对图像进展平移、旋转和缩放。将得到的T代入2式，求得参数。的作用是将图像经T变换后坐标值不落在像素点上的点调整到像素点上。本研究的算法主要包括以下几个步骤：给定特征点集Q和P；初始化：ij=1(全1矩阵)，=1(全1矩阵),T=0零矩阵,A=0零向量，=0，N=100，构造初始薄板样条；根据仿射坐标最大相关原那么迭代计算T，A；根据4式计算；根据6式计算ij；根据3式，构造薄板样条；假如满足4式约束或迭代次数大于N，那么转，否那么转；end。需要指出的是,由于互信息量是统计量,因此我们对标准图确定的点数不能太少，保证互信息

13、量的统计有效性。5实验结果图1显示了对一对9696,灰度级为256级的RI图像进展非刚性配准的实例。(a)是一幅标准的RI心脏长轴的一个切片；(b)是一幅有变形的RI图像心脏长轴同一周期另一个切片,此图像是用图像处理软件Phtshp手动变形处理获得；()是标准图像中特征点的选取,图中的小红圆圈为选择的特征点；(d)是最大互信息搜索方法在变形图像上搜索到的对应特征点；(e)是配准结果。由图1(a)、(b)、(e)我们可以看出,经过配准变换后,图1(b)变形为图1(e),与图1(a)已根本完全一致了。通过计算整幅图像之间的互信息量我们得到:图1(a)和图1(b)的互信息量为1.037,图1(a)和

14、图1(e)之间的互信息量为3.44,而图1(a)与自身计算得到的互信息量为3.46。可见,经过非刚性性配准变换,图像到达了较高的匹配精度。abde图1RI非刚性配准的实例6完毕语一直以来，学者们提出了各种非刚性图像匹配方法。文献5将图像分解为许多子图像，估计每个子图像之间仿射变换，这样用多个仿射变换近似反映整个图像之间的非刚性映射。文献6首先从图像中提取出特征点，然后将特征点拟合成曲线或曲面进展匹配。这种方法在拟合的曲线或曲面比拟光滑的情况下效果好，但是当图像所包含的形状很复杂时，曲线或曲面的拟合就变得非常困难，这种方法的好处是曲线的匹配相对于点集的匹配要容易些。最近还有许多非刚性匹配的研究主

15、要集中于非刚性形状的统计特性的学习上7，其方法类似于迭代最近点(IP)算法，是基于部分的启发式搜索。以上这些方法大都存在鲁棒性差，而且通常不能保证图像之间的一一对应和特征点的自动确定。本研究提出了一种通过薄板样条函数来表征特征点集之间的非刚性映射，把该映射分解为仿射变换和非仿射变换，并分别计算求解薄板样条的参数并满足双向对应的约束。参考文献1LikarB,pernuiF.Ahierarhialapprahtelastiregistratinbasednutualinfratin.IagVisinput，2001，19:3344.2KyriauSK，Davatziks.Nnlinearelast

16、iregistratinfbrainiagesithturpathlgyusingabiehanialdelRI.IEEEtransedIag，1999，18:580592.3ButtA，AharyaR，Sibata，etal.Surfaeathingfultidalityiagevluesbyafuzzyelastiregistratintehnique.putedIagGraph，1998，22:1323.4BkstEinFL.Prinipalarps:Thinplatesplinesandthedepsitinfdefratins.IEEETarans.1989，PAI11(6):567

17、585.5PappuS,GldS.AFraerkfrNnRigidathingandrrespndene.AdvaneinNeuralinfratinPressingSystes,ITpress,abridge:1996,8.795801.6SzeliskiR,LavalleeS.athing3Danatisurfaesithnnrigiddefratinsusingtreesplines.Int1JputerVisin,1996,18:171186.7DutaN,JainAK,etal.Learning2Dshapedels.IEEEnf.nputerVisinandPatternRegnitin(VPR).lbia:Frtlins,1999,2：814.8SunDngei,QiuZhen