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

基于粘滞流体B样条模型的快速非刚体配准方法 荣成城1) 周健2) 曹国刚3) 罗立民4) 1) 3) 4)（东南大学影像科学与技术实验室, 中国南京, 210096.） 2)（法国国家健康研究院，信号与图像处理实验室，法国雷恩，35000.） 摘要 基于粘滞流体模型的非刚体配准算法是一种适合个体差异较大配准场合的方法。算法的核心是求解流体偏微分方程组。通常的求解方法是迭代求解，比如逐次超松弛迭代，但是相应的计算复杂度是非常高的。为了降低计算负担，本文提出一种基于粘性流体模型结合B样条的快速配准方法。首先，对速度场进行B样条建模，将数量众多的未知量转化为数量较少的B样条插值系数；其次，利用B样条的特殊性质以及快速傅里叶变换（FFT），推导出B样条插值系数的逆形式。实验结果说明新方法不仅适用于大尺度形变的场合，同时具有较低的计算负担。 介绍 近年来，非刚体配准方法在医学图像分析领域得到了广泛的应用，尤其是在脑部图像的应用中[1-3]。在典型的脑部图像配准应用中，由于人脑间存在着较大的的差异，通常要求非刚体配准方法能够处理大尺度形变的复杂情况。基于粘性流体模型的配准方法最早由Christensen等人提出，被证明是一种能够适应这种需求的方法[4]。在该算法中，图像像素被当作流体粒子，其运动规律满足特定运动偏微分方程。算法的核心步骤是偏微分方程的求解步骤，Christensen给出的数值解法是传统的有限差分结合逐次超松弛迭代，但是这种解法的计算代价非常大。针对这一缺陷，很多作者给出了相应的改进方法。Freeborough和Fox设计了一种多格结合逐次超松弛的方法，在台式机上得到较好的结果[5]。Wollny和Kruggel使用最小残差算法求解方程，减少了迭代的次数，但是却需要更多的计算时间[6]。Bro-Nielsen等人设计了一种偏微分方程的滤波求解方法[7]，使得算法可以在单核计算上得到应用。Crum等人基于多分辨率策略和多格方法，设计了一种多分辨率求解方法，也得到了较好的结果[8]。 本文中，我们设计了一种基于B样条的快速求解方法。基于偏微分方程的连续表达形式，模型的求解可以进行精确的求解，而不需要依赖于有限差分步骤。与此同时，通过使用B样条表达形式，偏微分方程中的未知量可以转化为B样条插值系数，从而大大降低未知量的数量，进而大大降低计算的复杂度。我们设计了一种基于快速傅里叶变换的方法，来快速算得所需的插值系数。 下面章节组织结构如下：首先给出原始算法的描述，然后在第二部分中给出新的快速算法。第三部分给出人造数据以及真实数据实验的结果及比较。第四部分给出了总结与展望。 2.方法 2.1粘性流体非刚体配准方法 不失一般性，本文考虑两维图像的配准。设分别为待配准的源图像和目标图像函数，其中表示像素空间坐标。如前文所述，配准的目的是计算如下的位移场。 我们可以定义如下的相似能量泛函： 然后通过最小化该能量（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)可以表示为 由