基于粘滞流体B 样条模型的快速非刚体配准方法.docVIP

基于粘滞流体B样条模型的快速非刚体配准方法 荣成城1) 周健2) 曹国刚3) 罗立民4) 1) 3) 4)（东南大学影像科学与技术实验室, 中国南京, 210096.） 2)（法国国家健康研究院，信号与图像处理实验室，法国雷恩，35000.）分别为待配准的源图像和目标图像函数，其中表示像素空间坐标。如前文所述，配准的目的是计算如下的位移场。 我们可以定义如下的相似能量泛函： 然后通过最小化该能量（1）来获得相应所需的位移场。在粘性流体配准方法中，u是通过求解偏微分方程来得到的。让我们来考虑单个流体粒子的运动过程。设位于处的粒子时刻所受外力为，它是关于当前位置和位移的函数。设该粒子受力而产生的瞬时速度为，则根据粘滞流体运动学原理[9]，则有如下偏微分方程： 此处为拉普拉斯算子，为梯度算子，为散度运算，为粘滞系数，本文考虑一般情况，取，。该偏微分方程反映了u随v和k变化而变化的关系。此时，k可以取泛函关于u的变分形式： 因此，最优的位移场u可以通过求解偏微分方程（2）来获得。 通过引入时域的离散化操作，第一个偏微分方程可以改写为： 其中表示时间间隔。 在二维空间中展开和，可以得到： 通过离散化v1，v2，k1，k2，Christensen【4】等使用逐次超松弛来求解上述方程组。 2.2引入B样条 求解偏微分方程组（5）的准确性取决于有限差分过程。但是如果假设v1，v2可以通过B样条函数表达成线性分段多项式，则可以设计更快的求解过程。假设t时刻下，二维的速度场可以表达成如下的B样条函数形式： 其中p和q是各方向上所使用的B样条基函数个数，分别为第个二维B样条基函数，为对应的控制顶点系数。 现考虑离散图像情况，其像素位置由一系列坐标：（N为像素总数）确定。相应的离散速度场记为 ，。定义矩阵, ，以及向量，则式(8)可表达为： 类似地，利用B样条表达形式，速度场的导数可以表示如下形式：（以v1(x)为例） ,, 接着定义如下矩阵： 则方程组(5)可转化为以B样条系数为待求未知量的形式： 其中为时刻的外力场向量。 2.3插值系数的快速计算 令 ，则(10)可以表示为 由于(11)通常是超定的，所以我们求其最小二乘意义下的解，比如求解如下的线性方程组： 当上述矩阵的逆不存在时（比如的尺寸很大），通常可以使用迭代的方法进行求解。但是在使用B样条插值函数的情况下，我们可以直接基于傅里叶变换来求出上述解形式。实际上，基于Unser[10]提出的B样条插值理论，我们可以很容易地发现均可用卷积结合降采样运算来快速获得。这里所使用的B样条控制网格点进行均匀地摆放，并且保证是像素间距的整数倍。更多的细节可以参考文献[10]。如果考虑信号的循环扩展，可以发现，这些矩阵都是循环卷积矩阵。另外，如果定义： 不难证明为循环Toeplitz矩阵。 根据Toeplitz矩阵的性质，可以表示为： 其中表示共轭转置矩阵，为对角阵，其对角元素由向量中的元素构成。为矩阵的特征值构成的向量。注意到与图像无关，所以可以提前计算出来。 将（14）代入（13），可以得到： 这里表示向量间对应元素相乘。显然，方程组解的直接计算公式为： (16)中的除法是向量间对应元素相除。实际计算中，两维离散傅立叶变换及其反变换均可以利用FFT来实现，因此求解的运算量将大大降低。 新算法的流程同原始算法类似，唯一的区别在于首先计算出B样条插值系数，然后再根据（7）式算出相应的速度场。 3.实验 3.1人造数据 我们选择部分‘C’字符作为模板图像，而将整个字符作为目标图像（见图一）。两幅图像的大小均为128*128.图二中给出了实验结果。本文提出的算法迭代280次，共耗费14.9s。而同样的实验，文献[8]中的算法需要423s。这个结果验证了本文提出的算法适合大尺度形变的场合，并且具有相当快的计算速度。 (a) 源图像 (b) 目标图像 图1 C字符形变实验 3.2真实数据 为了验证本文算法在真实数据上的效果，我们选取了三组真实数据进行配准实验。所有图像的尺寸均为256*256。 本文中使用源与目标图像的SSD(Sum of Squared Difference)，以及相关系数CC(Correlation Coefficient)来作为相似性的评价指标，SSD值越小表明相似度越高，CC值越大说明相似度越大，计算公式如下： 其中，表示各自的均值，N为像素总数。表1给出了相关的实验结果统计。 表1 算法耗