非刚性医学图像配准算法的设计与实现
【关键词】 医学图像；非刚性；图像配准；匹配矩阵；薄板样条
摘要：非刚性图像匹配问题已成为医学图像分析中一个非常具有挑战性的问题。基于薄板样条插值方法,引入实匹配矩阵,并给出相应配准变换算法，该算法将薄板样条参数表示成仿射分量和非仿射分量，并分别进行求解。与其它非刚性匹配算法相比，该算法不仅保证了对应特征点的双向对应，也实现了自动特征点选择，实验结果令人满意。
关键词：医学图像；非刚性；图像配准；匹配矩阵；薄板样条
1引言
在医学诊断和过程中，常需要对比分析多幅图像，以获得更为精确和全面的信息。图像分析大都要求多幅图像的几何位置一致，因此，配准是医学图像分析的一个重大课题。医学图像配准是指对于一幅医学图像寻求一种(或一系列)空间变换，使它与另一幅医学图像上的对应点达到空间上的一致。这种一致是指人体上的同一解剖点在两张匹配图像上有相同的空间位置。配准的结果应使两幅图像上所有的解剖点，或至少是所有具有诊断意义的点及手术感兴趣的点都达到匹配。图像配准不仅可以校正病人多次成像间的位置变化，也可以校正由于成像模式本身导致的畸变。对同一个病人的不同时间的图像进行配准，可以了解发育过程及肿瘤病变的病情；对不同人的图像进行配准，去除种族、年龄等临床及遗传差异，从而形成疾病或人群特异性图谱，可用于正常与否的分析；对不同成像模式进行配准，可以获得互补信息。
医学图像配准可分为刚性配准和非刚性配准两类。刚性配准在许多情况下不能满足临床的需要，因为很多形变的性质是非刚体、非线性的。比如为了精确定位MR图像左心室，常常伴有组织磁化系数差异、非水分子的化学位移以及血流流动等因素导致的几何畸变以及由于磁场不均匀、磁场梯度非线性及涡流等导致的探测畸变，因此在放疗计划制定中，将MR图像配准时，不能单纯地使用刚性配准，必须使用非刚性配准。
非刚性配准算法可分为灰度驱动、模型驱动及混合算法三种［1~3］。灰度驱动方法基于数学或统计尺度将一个灰度模式与另一个对准。典型情况下，需要定义源系统与目标系统之间的灰度相似性的数学量度。灰度相似性测度包括象素灰度的均方差、相关或互信息。模型驱动方法首先建立明确的几何模型，以此表示解剖标志。这些解剖标志包括有重要功能的表面、曲线和点。将源系统的解剖标志参数化，与目标系统的对应部分对准，以这种对应关系引导系统其余部分的变换。模型驱动算法包括点约束法、线约束法和面约束法。混合算法是结合使用以上两种算法的方法。薄板样条插值方法是非刚体变换中的一种特殊的变换，它允许局部调整，并符合某种连续性或平滑性要求。第2节讨论刚性能量函数；第3节给出非刚性能量函数；第4节设计并实现一个非刚性配准算法；最后给出实验结果。
2刚性能量函数
本研究之所以采用薄板样条，是因为它的独特性质，就是能够将空间变换分解为一个全局仿射变换和一个局部非仿射变换。Booksteein［4］首先将薄板样条函数应用于标志点的匹配，结果证明它是一个非常有用的形状分析工具。假设在二维空间，已知两个具有N对对应点的点集，Q={Qi，i=1，2，…，n}和P={Pi，i=1，2，…，n}，将点集Q，P表示为：
Q=1 x1 y1
1 x2 y2
… … …
1 xn ynP=1 x1 y1
1 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)是构建薄板样条的基函数，设rij=|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) … 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的仿射形变参数矩阵，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非刚性能量函数
由刚性能量函数推导表明，只要已知两个点集之间的对应点，就可以得到它们之间的薄板样条映射参数。但是当对应点未知时，该如何处理呢?传统的方法往往都是手动选点,这种方法费时费力,同时在结构不清的情况下,很难选择到足够多的精确对应点。而且其准确性也只是相对的,误差是不可避免的。［10］定义两个点集之间的匹配矩阵M={Mij}：
Mij=1，若点Qi对应于点Pi
0，其他
由于两个点集之间是双向一一对应的，即一个点集中的每个点在另一个点集中至多有一个对应点，反之亦然。匹配矩阵{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的值逐渐变大，越来越接近二值矩阵，当时间足够长时，就会得到最终的二值匹配矩阵。
根据匹配矩阵元素的性质，令：
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按最大互信息原则沿图像两个正交方向逐像素搜索，以取得最佳值Δx和Δy。
确定了配准参数K、θ、Δx和Δy，就可对图像进行平移、旋转和缩放。将得到的T代入（2）式，求得参数W。W的作用是将图像经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实验结果
图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)算法，是基于局部的启发式搜索。以上这些方法大都存在鲁棒性差，而且通常不能保证图像之间的一一对应和特征点的自动确定。本研究提出了一种通过薄板样条函数来表征特征点集之间的非刚性映射，把该映射分解为仿射变换和非仿射变换，并分别计算求解薄板样条的参数并满足双向对应的约束。
文献
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. Nonlinear elastic registration of brain images with tumor pathology using a biomechanical model［MRI］.IEEE trans Med Imag，1999， 18:580~592.
3Butt A， Acharya R， Sibata C，et al.Surfacematching of multimodality image volumes by afuzzy elastic registration technique.Comput Med Imag Graph， 1998，22:13~23.
4Bookstein F L. Principal warps: Thin plate splines and the decomposition of deformations .IEEE Tarans. 1989， PAMI 11(6):567~585.
5Pappu S, Gold S. A Framework for Non Rigid Matching and Correspondence. Advance in Neural information Processing Systems, MIT press, Cambridge:1996,8.795~801.
6Szeliski R, Lavallee S. Matching 3D an atomic surfaces with non rigid deformations using octree splines. Int 1 J Computer Vision,1996,18:171~186.
7Duta N, Jain A K,et al. Learning 2D shape models. IEEE Conf. On Computer Vision and Pattern Recognition(CVPR). Colombia:Fort Colins,1999,2：8~14.
8Sun Dongmei,Qiu Zhengding. A New Non Rigid Image Matching Algorithm Using Thin Plates Spline. ACTA Electronic Sinica(in China). 2002,30(8):1104~1108.
9Ren Haihua， et al. Medical Image Elastic Registration Based on Thin Plate Splines. Beijing Biomedical Engineering(in Chnia). 2003, 22(4):241~245.
10孙冬梅,裘正定.利用薄板样条函数实现非刚性图像匹配算法.学报，2002,30(8): 1104~1107.
