CN1975781A - 基于有限采样全局优化的图像弹性配准方法 - Google Patents

Abstract

本发明公开一种基于有限采样全局优化的图像弹性配准方法，该方法的配准步骤为：先读入变形图像，计算定义的目标函数的最高频率，再以目标函数的最高频率的2倍对变形系数进行均匀采样计算，得到变形图像系列，然后读入参考图像，以最小目标函数值所对应的变形函数对变形图像进行变形，完成粗配准，进而以粗配准中的变形系数值为起始点，用梯度下降法搜索目标函数的全局最小值的点，完成图像的配准。本发明依据采样定理，采用有限的采样点就可以得到目标函数的全局最小值的点，即可实现图像的弹性配准，因此具有计算工作量小，配准时间短的显著技术效果。此外，本发明还克服了随机采样需要操作者根据经验进行人工干预，难以保证配准效果一致性的不足。

Description

基于有限采样全局优化的图像弹性配准方法

技术领域

本发明涉及图像平面内的图形图像转换，具体涉及一种全局优化的图像弹性配准方法，适用于图像对比、数据融合、变化分析和目标识别等领域。

背景技术

图像配准是指对于一幅图像寻求一种空间变换，使它与另一幅图像上的对应点达到空间上的匹配。这种匹配是指内容相同的点在两张匹配图像上有相同的空间位置。

常用的图像配准方法有很多种，但它们的配准步骤基本是相同的，主要包括以下步骤：

(1)构造对图像进行变形的变形函数。现有的配准方法可以分为刚性配准方法和弹性配准方法两种。刚性配准方法中使用的变形函数是一种刚性变形函数，它只能对图形进行平移、旋转等相对简单的变换。如果图像畸变严重，配准时就需要引入弹性变形函数，相应的配准方法也就是弹性配准方法。弹性变形函数具有足够的通用性，可以逼近任意的非线性变换。弹性变形函数基本可以分为两类：一类是非参数化模型。在这种模型中，图像被看成是一片有弹性的薄膜，在外力和内力的作用下变形，并最终达到平衡。外力由参考图像和变形图像的差异确定，内力由薄膜的强度和平滑程度确定。这种模型中最有名的是流体模型，此外还有扩散模型、光流模型等。另一类是参数化模型，模型需要使用一些参数和基函数来表示，变形函数的计算过程也就是参数的计算过程。基函数的种类很多，常见的基函数包括多项式、谐波函数、分级基函数、小波以B样条等。

(2)构造需要配准的参考图像和变形图像的相似度准则。相似度准则有时也称为目标函数，它是衡量参考图像和变形图像之间的相似性的。常用的目标函数有互信息量、均方差等。

(3)对目标函数进行优化计算，直到目标函数达到它的极值。该过程就是对变形图像用变形函数不断进行变形，当目标函数达到极值时，即参考图像和变形后的变形图像之间的相似性最大时，配准过程也就完成了。上述配准过程中常用的优化方法多是一些局部优化方法，如牛顿法、半牛顿法、共轭向量法等。这种方法的缺点是不能保证找到目标函数的全局极值，从而不能保证配准结果的正确性。法也可以使用全局优化算法，相应的的全局优化方法有模拟退火、遗传方法、隧道法、随机穷尽法、随机多起点法等方法。其中，模拟退火和遗传方法计算时间太长，并且不能保证最后得到全局解；隧道法的缺点就是它的优化对象有一定的局限性，通用性不好；而各种基于随机穷尽的方法在随机点的数量的计算上也没有明确的计算方法，要得到准全局解，只有加大随机点的数量，这就造成计算时间过长，同时也不能保证最后配准结果的正确性。

发明内容

本发明所要解决的技术问题是减小计算目标函数的全局极值点的计算量，缩短图像的弹性配准的时间，提高图像配准的鲁棒性。

本发明解决上述的技术问题的技术解决方案如下所述。

一种基于有限采样全局优化的图像弹性配准方法，该方法由以下步骤组成：

(1)读入变形图像，用快速傅立叶变形计算它的频谱，由该频谱计算它的最高频率ω_m。

(2)用傅立叶变换方法计算目标函数E(C)相对于变形系数C的频谱，再由该频谱计算得到目标函数的最高频率。该最高频率为ω_max＝2ω_m·max(g′(x，a))，其中ω_m为步骤(1)中变形图像的最高频率，g′(x，a)为变形函数的导数。当样条次数n＝1时，g′(x，C)为1。

(3)按步骤2得到的目标函数的最高频率ω_max的2倍对变形系数C对应的空间进行均匀采样计算，得到一系列不同的变形系数C的值c_i所对应的变形函数g(x，c_i)，再用所得到的这些变形函数分别对变形图像进行变形，得到变形图像系列；

(4)读入参考图像，分别计算参考图像与步骤(3)得到的变形图像系列的均方差，得到目标函数值数据组E(c_i)，比较这些目标函数值的大小，找到其中最小的目标函数值E(c^*)，用该最小目标函数值所对应的变形函数g(x，c^*)对变形图像进行变形，得到粗配准的变形图像。

(5)以步骤(4)所找到的最小的目标函数值E(c^*)所对应的变形系数值c^*为起始点，用梯度下降法搜索目标函数E(C)的最小值

再用该最小值所对应的变形函数 对变形图像进行变形，便得到精配准的变形图像。

根据现有技术可知，图像配准实质上就把待配准的两幅图像，其中一幅作为参考图像，另一幅作为变形图像，以参考图像作为基准对变形图像进行变形，使二者匹配。据此，配准后的变形图像可理解为是由参考图像变形得到的。因此，本发明与现有技术一样，在配准前先要定义联系参考图像和变形图像的变形函数和目标函数。变形函数和目标函数的定义方法有很多种，本发明采样B样条来定义变形函数和目标函数，具体定义如下：

(1)用B样条构造的变形函数为：

$g(x,C)=x+\underset{i\∈I_c}{\mathrm{\Σ}}{C}_i{\mathrm{\β}}_n(x/2^w-i)$

其中C为变形系数，β_n是n次B样条的向量积，即 ${\mathrm{\β}}_n\left(x\right)={\mathrm{\Π}}_{k=1}^m{\mathrm{\β}}_n\left(x_k\right),$ x为m维向量，2^*为B样条的节点间距，w表示层次，即相邻层次间的伸缩尺度为2，I_c为B样条的节点空间集合，i为节点在节点空间I_c中的位移。x为原图像的坐标，g(x)为经变换后相应的坐标值。

(2)设参考图像和变形图像分别为f_r(x)和f_t(x)，以参考图像和用变形函数变形后的变形图像(它可以表示为f_t(g(x，C)))的均方差作为目标函数，则目标函数可以定义为：

$E\left(C\right)=\frac{1}{\left|\right|I\left|\right|}\underset{x\∈I}{\mathrm{\Σ}}{e}_i^2=\frac{1}{\left|\right|I\left|\right|}\underset{x\∈I}{\mathrm{\Σ}}{(f_t\left(g(x,C)\right)-f_r\left(x\right))}^2$

其中I是图像坐标所对应的空间，g(x，C)是变形函数；

根据采样定理，当采样频率大于或等于一个函数最高频率的两倍时，这个函数可以用这些采样点完整的恢复出来。本发明所述方法，先读入变形图像，计算定义的目标函数的最高频率，再以目标函数的最高频率的2倍对变形系数进行均匀采样计算，得到变形图像系列，然后读入参考图像，以最小目标函数值所对应的变形函数对变形图像进行变形，完成粗配准，进而以变形函数系列中每一具体值为起始点，用梯度下降法搜索目标函数的全局最小值的点，完成图像的配准。由此可见，本发明依据采样定理，采用有限的采样点就可以得到目标函数的全局最小值的点，即可实现图像的弹性配准，因此具有计算工作量小，配准时间短的显著技术效果。此外，本发明还克服了随机采样需要操作者根据经验进行人工干预，难以保证配准效果一致性的不足。

附图说明

图1为本发明所述方法的流程图；

图2是从一组心脏灌注的MRI图中挑选出的序号为5和3的相邻两幅图像图，分别冠名为图2a和图2b，其中图2a作为参考图像，图2b作为变形图像；图2c是图2b所示图像采用本发明配准方法变形后的图像；图2d是图2b所示图像与图2c所示图像的差值图像。

图3为噪声对配准准确率的影响曲线图。

图4为当噪声方差为0.1时，本发明配准方法与参考配准方法对同一变形图像的配准结果对比图，其中，图4a为参考图像；图4b为图4a在加入方差0.1的噪声的变形图像；图4c为以图4a为参考图像，图4b为变形图像时，采用本发明方法配准得到的图像；图4d为以图4a为参考图像，图4b为变形图像时，采用参考方法配准得到的图像。

具体实施方式

以下通过四个非限制性实施例对发明进行详细说明。

例1(二维图像的配准)：

为医学研究的需要，估计人体心脏的运动参数，本发明人从天津医科大学总医院采集了一组心脏灌注的MRI图像，共20幅，大小为256×256象素，从中选择的图像是序号为5和3的相邻两幅图像为参考图像和变形图像(参见图2a和图2b)。配准的试验平台是PC机，以WindowXP为操作***，Matlab6.5为编程语言。具体配准过程如下所述：

配准前先定义变形函数和目标函数，具体步骤如下：

(1)设参考图像和变形图像分别定义为f_r(x)和f_t(x)，其中变量x为二维向量，则用B样条构造的联系参考图像和变形图像的变形函数为：

$g(x,C)=x+\underset{i\∈I_c}{\mathrm{\Σ}}{C}_i{\mathrm{\β}}_n(x/2^w-i)$

其中C为变形系数，β_n是n次B样条的向量积，即 ${\mathrm{\β}}_n\left(x\right)={\mathrm{\Π}}_{k=1}^m{\mathrm{\β}}_n\left(x_k\right).$ 在这里B样条取为一次样条。

(2)同样设参考图像和变形图像分别定义为f_r(x)和f_t(x)，以参考图像和用变形函数变形后的变形图像(它可以表示为f_t(g(x，C)))的均方差作为目标函数，则目标函数定义为：

$E\left(C\right)=\frac{1}{\left|\right|I\left|\right|}\underset{x\∈I}{\mathrm{\Σ}}{e}_i^2=\frac{1}{\left|\right|I\left|\right|}\underset{x\∈I}{\mathrm{\Σ}}{(f_t\left(g(x,C)\right)-f_r\left(x\right))}^2$

其中x为二维向量。使目标函数E(C)取得极值的变形函数g就对应最佳变形，即求g＝argmin_g∈GE(g(x，C))，其中G是变形函数g所对应的空间；

作好配准前的准备后，即可进行配准，具体步骤如下：

(1)读入图2b所示变形图像，该图像可表示为f_t(x)，其中x为二维向量，用快速傅立叶变形计算出f_t(x)的频谱为：

$F_t\left(\mathrm{\ω}\right)={\∫}_{-\∞}^{+\∞}{f}_t\left(x\right)\·e^{-\mathrm{j\ωC}}\mathrm{dC}$

由该频谱计算f_t(x)的最高频率ω_m

${\mathrm{\ω}}_m=\underset{F\left(\mathrm{\ω}\right)\≠0}{\max }\left(\mathrm{\ω}\right)$

对于实际的医学图像，ω_m不会很大，经过计算统计ω_m一般为3～10之间。对于图2b，经计算ω_m为5。

(2)用傅立叶变换方法计算目标函数E(C)相对于变形系数C的频谱，由该频谱计算得到目标函数的最高频率为：

ω_max＝2ω_m·max(g′(x，C))

该式中ω_m为步骤(1)中变形图像的最高频率，g′(x，C)为变形函数的导数。本发明中B样条的次数n取为1时，这时目标函数E(C)的最大频率ω_max为2ω_m。

目标函数的最高频率ω_max＝2ω_m·max(g′(x，C))的具体推导过程如下所述：

目标函数E(C)对于变形系数C作傅立叶变换：

$E\left(\mathrm{\ω}\right)={\∫}_{-\∞}^{+\∞}\underset{x\∈I}{\mathrm{\Σ}}{(f_r\left(x\right)-f_t\left(g(x,C)\right))}^2\·e^{-\mathrm{j\ωC}}\mathrm{dC}$

$={\∫}_{-\∞}^{+\∞}(f_r^2\left(x\right)-2f_r\left(x\right)f_t\left(g(x,C)\right)+f_t^2\left(g(x,C)\right))\·e^{-\mathrm{j\ωC}}\mathrm{dC}$

$=\underset{x\∈I}{\mathrm{\Σ}}{\∫}_{-\∞}^{+\∞}(f_r^2\left(x\right)-2f_r\left(x\right)f_t\left(g(x,C)\right)+f_t^2\left(g(x,C)\right))\·e^{-\mathrm{j\ωC}}\mathrm{dC}$

$=\underset{x\∈I}{\mathrm{\Σ}}({\∫}_{-\∞}^{+\∞}{f}_r^2\left(x\right)e^{-\mathrm{j\ωC}}\mathrm{dC}-2{\∫}_{-\∞}^{+\∞}{f}_r\left(x\right)f_t\left(g(x,C)\right)e^{-\mathrm{j\ωC}}\mathrm{dC}+{\∫}_{-\∞}^{+\∞}{f}_t^2\left(g(x,C)\right)e^{-\mathrm{j\ωC}}\mathrm{dC})$

上式共有三项，其中第一项

$\underset{x\∈I}{\mathrm{\Σ}}{\∫}_{-\∞}^{+\∞}{f}_r^2\left(x\right)e^{-\mathrm{j\ωC}}\mathrm{dC}=2\mathrm{\π}\·\mathrm{\δ}\left(\mathrm{\ω}\right)\·\underset{x\∈I}{\mathrm{\Σ}}{f}_r^2\left(x\right)$

第二项为

$\underset{x\∈I}{\mathrm{\Σ}}{\∫}_{-\∞}^{+\∞}{f}_r\left(x\right)f_t\left(g(x,C)\right)e^{-\mathrm{j\ωC}}\mathrm{dC}=\underset{x\∈I}{\mathrm{\Σ}}{f}_r\left(x\right)\·{\∫}_{-\∞}^{+\∞}{f}_t\left(g(x,C)\right)e^{-\mathrm{j\ωC}}\mathrm{dC}$

上式的也就是 ${\∫}_{-\∞}^{+\∞}{f}_t\left(g(x,C)\right)e^{-\mathrm{j\ωC}}\mathrm{dC}$ 的。现在计算 ${\∫}_{-\∞}^{+\∞}{f}_t\left(g(x,C)\right)e^{-\mathrm{j\ωC}}\mathrm{dC}$ 的。

把g(x，C)在任意点C＝a处泰勒展开，得

$g(x,C)=G(x,a)+g^{\′}(x,a)\·(C-a)+\frac{g^{\′\′}(x,a)}{2!}{(C-a)}^2+\·\·\·\·\·\·$

用一阶近似，得

g(x，C)＝g(x，a)+g′(x，a)·(C-a)＝(g(x，C)-g′(x，a)·a)+g′(x，a)·C

在步骤(1)中，f_t(x)频谱为F(ω)，它的最高频率为ω_m，则f_t(g′(x，a)·C+b)的频谱的绝对值为

它的为ω_m·g′(x，a)。因为a为任意点，所以 ${\∫}_{-\∞}^{+\∞}{f}_t\left(g(x,C)\right)e^{-\mathrm{j\ωC}}\mathrm{dC}$ 的为ω_m·max(g′(x，C))。

第三项为 ${\∫}_{-\∞}^{+\∞}{f}_t^2\left(g(x,C)\right)e^{-\mathrm{j\ωC}}\mathrm{dC}.$ 因为不需要求上式的频谱，只是需要求它的最高频率，而f_t(g(x，C))的最高频率为ω_m·max(g′(x，C))，所以f_t ²(g(x，C))的最高频率为2ω_m·max(g′(x，C))。

比较以上三项的最大频率，我们可以看出目标函数E(C)的最高频率是2ω_m·max(g′(x，C))。

由g(x，C)的定义可以知道，g(x，C)中的变形系数为C，所以g′(x，C)＝β_n(x/2^w-i)β_n(y/2^w-j)，β_n(x/2^w-i)的最大值和n的大小有关。当n＝1时，β_n(x/2^w-i)的最大值为1，所以目标函数E(C)的最大频率ω_max为2ω_m。

(3)目标函数E(C)的最高频率是ω_max＝2ω_m，其中ω_m为图像2b的最高频率。那么以2ω_max进行采样就是以4ω_m对变形系数C进行均匀采样，这时采样间隔为 $T=\frac{\mathrm{\π}}{2{\mathrm{\ω}}_m\mathrm{}},$ 对图像2而言，T为0.314.。在C的定义域(定义域的大小和对图像变形大小的估计有关，一般在10以内，本例取10)内每隔0.314采一个点，采样后得到一系列的具体的值c_i(i∈M，M为采样的点集)，对应的也就有一系列的具体的变形函数g(x，c_i)，用这些变形函数对图2b所示图像进行变形，得到变形图像系列。

(4)读入图2a所示参考图像，分别计算它与步骤(3)得到的变形图像系列的均方差，就得到一系列的E(c_i)值，比较这些E(c_i)的大小，从中找出它们的最小值E(c^*)，这时对应的变形函数g(x，c^*)就已经是一个比较好的对图2b所示图像的变形了(粗配准后的)。

(5以步骤(4)所找到的最小的目标函数值E(c^*)对应的变形系数值c^*作为起始点，用梯度下降法搜索目标函数E(C)的全局最小值点 。具体过程描述如下：

以在步骤(4)里得到的变形函数g(x，C^*)中的变形系数C^*为初始值。在迭代过程的第i步，我们根据现有的变形系数c_i计算更新的步长Δc_i＝-μ_i__cE(c_i)，如果在变形系数为c＝c_i+Δc_i的情况下，E的值更小，那么这步更新就是成功的，以新的参数值替换以前的参数值，即c_i+1＝c_i+Δc_i，μ_i+1＝μ_i；如果E的值反而增大了，就取μ_i+1＝μ_i/2，再次计算梯度值Δc_i＝-μ_i+1__cE(c_i)；这个迭代过程一直继续下去，直至Δc小于事先设定的某一阈值。在c_i确定后，变形函数也就确定了。用该变形函数对步骤(4)得到的粗配准的变形图像进行变形，便得到精配准的变形图像。

在上述优化过程中，需要计算E(C)的偏导数向量__cE(c_i)。E(C)的一阶偏导数为

$\frac{\∂E}{\∂c_{j,m}}=\frac{1}{\left|\right|I\left|\right|}\underset{i\∈I}{\mathrm{\Σ}}\frac{\∂e_i}{\∂f_w\left(i\right)}\frac{\∂f_t^c\left(x\right)}{\∂x_m}{|}_{x=g\left(i\right)}\frac{\∂g_m\left(i\right)}{\∂c_{j\·m}}$

其中 $\frac{\∂e_i}{\∂f_w\left(i\right)}=2{(f_w\left(i\right)-f_r\left(i\right))}^2$

$\frac{\∂g_m}{\∂c_{j,m}}={\mathrm{\β}}_{n_m}(x/h-j)$

$\frac{\∂f_t^c}{\∂x_m}\left(x\right)=\underset{k\∈I}{\mathrm{\Σ}}{b}_k{\mathrm{\β}}_n^{\′}(x_m-k_m)\underset{l=1,l\≠m}{\overset{N}{\mathrm{\Π}}}{\mathrm{\β}}_n(x_l-k_l)$

用目标函数E(C)取得全局最小值点

时所对应的变形函数

对变形图像进行变形，得到精配准的变形图像，配准后的图2b所示图像如见图2c所示。图2b所示图像与图2c所示图像的差值图像如图2d所示，它表示所估计的心脏的运动场，可以看出图像中心脏的运动场非常清晰。

例2(三维图像配准)：

本实施例所用的三维图像是从南方医科大学生物医学工程学院医学信息研究所图像数据库中选取的MRI图像，图像大小为128×128×128象素。配准的试验平台是PC机，以WindowXP为操作***，Matlab6.5为编程语言。具体配准过程如下所述：

本实施例中，变形函数和目标函数的定义方法与例1相同，具体为：

(1)设参考图像和变形图像分别定义为f_r(x)和f_t(x)，其中x为三维向量。用B样条来构造联系参考图像和变形图像的变形函数，得到参数化弹性模型：

$g(x,C)=x+\underset{i\∈I_c}{\mathrm{\Σ}}{C}_i{\mathrm{\β}}_n(x/2^w-i)$

其中C为变形系数，β_n是n次B样条的向量积，即 ${\mathrm{\β}}_n\left(x\right)={\mathrm{\Π}}_{k=1}^m{\mathrm{\β}}_n\left(x_k\right).$ 在这里B样条取为一次样条。

(2)同样设参考图像和变形图像分别定义为f_r(x)和f_t(x)，以参考图像和用变形函数变形后的变形图像(它可以表示为f_t(g(x，C)))的均方差作为目标函数，则目标函数定义为：作好配准前的准备后，即可进行配准，具体步骤如下：

$E\left(C\right)=\frac{1}{\left|\right|I\left|\right|}\underset{x\∈I}{\mathrm{\Σ}}{e}_i^2=\frac{1}{\left|\right|I\left|\right|}\underset{x\∈I}{\mathrm{\Σ}}{(f_t\left(g(x,C)\right)-f_r\left(x\right))}^2$

(1)读入三维变形图像，该图像可表示为f_t(x)，其中x为三维向量，用快速傅立叶变形计算它的频谱 $F_t\left(\mathrm{\ω}\right)={\∫}_{-\∞}^{+\∞}{f}_t\left(x\right)\·e^{-\mathrm{j\ωC}}\mathrm{dC},$ 由该频谱计算它的最高频率 ${\mathrm{\ω}}_M=\underset{F\left(\mathrm{\ω}\right)\≠0}{\max }\left(\mathrm{\ω}\right).$

(2)用傅立叶变换方法计算目标函数E(C)相对于变形系数C的频谱，由该频谱计算得到目标函数的最高频率ω_max＝2ω_m·max(g′(x，a))。该式中ω_m为步骤(1)中变形图像的最高频率，g′(x，a)为变形函数的导数。具体推导与例1相同。

由g(x，C)的定义可以知道，g(x，C)中的变形系数为C，所以本实施例中g′(x，C)＝β_n(x/2^w-i)β_n(y/2^w-j)β_n(z/2^w-k)。当n＝1时，β_n(x/2^w-i)的最大值为1，所以目标函数E(C)的最大频率ω_max为2ω_m。

(3)目标函数E(C)的最高频率是ω_max＝2ω_m，其中ω_m为变形图像的最高频率。那么以2ω_max进行采样就是以4ω_m对变形系数C进行均匀采样(相对于二维图像而言，三维图像的不同点是采样是在三维中进行)，采样后得到一系列的具体的值c_i(i∈M，M为采样的点集)，对应的也就有一系列的具体的变形函数g(x，c_i)，用这些变形函数变形图像进行变形，得到变形图像系列。

(4)读入参考图像，分别计算它与步骤(3)得到的变形图像系列的均方差，就得到一系列的E(c_i)值，比较这些E(c_i)值的大小，从中找出它们的最小值E(c^*)，完成粗配准。

(5)以步骤(4)所找到的最小的目标函数值E(c^*)对应的变形系数值c^*作为起始点，用梯度下降法搜索目标函数E(C)的全局最小值点

完成精配准。具体过程同例1。

例3(配准准确率对比实验)：

由于在变形已知的条件下可以更好的检验本发明所述方法的试验效果，所以本发明人采用人工变形的方法来产生变形图像，随机从本实验室图库中选取了16幅图像作为参考图像。所选取的参考图像其中，CT图像8幅，MR图像8幅，并分别对它们进行相同的变形处理。原图像的尺寸是256×256象素。变形的内容包括沿X和Y轴分别平移10各象素、顺时针旋转10度、以及用软件photoshop实现的扭曲弹性变形。与本发明对比的参考方法为Kybic J和Unser M于2003年在《IEEE Trans.Image Processing》第12卷11期1427-1442页发表的文章“Fast Parametric Elastic Image Registration”中基于B样条的弹性配准方法。比较结果如表1所示。由表1可见，本发明算法实现了总体93.75％的准确率，而参考算法配准效果较差，准确率只有31.25％。本发明算法的配准效果明显好于参考算法。

表1 16组试验的结果对比

	CT图像	MR图像	总和
				试验样本数量	8	8	16
配准正确数量(本发明方法)	8	7	15
				准确率(本发明方法)	100％	87.5％	93.75％
配准正确数量(参考方法)	3	2	5
				准确率(参考方法)	37.7％	25％	31.25％

例4(抗噪声试验)：

1、鲁棒性比较试验

本发人从例3的16组图像中选择了5组，分别采用本发明方法和例3所述的参考方法进行配准实验。在实验的5组图像中，加入不同强度的噪声来检验两种方法对噪声的鲁棒性。加入的噪声是高斯噪声，噪声方差分别为0.1和0.5。配准准确率和噪声方差的关系如图3所示。从图3可以看出，本发明方法在不同噪声水平上的配准准确率都要高于参考算法，本发明方法对噪声的鲁棒性较强。

2、准确性比较实验

为了更清楚的显示配准的实际效果，本发明人又选取了一CT图像(如图4a所示)，并加了方差为0.1的噪声(见图4b)，然后分别用本发明方法和例3所述的参考方法进行配准实验。本发明方法的配准结果如图4c所示，参考方法的配准结果如图4d所示。从试验结果可以清楚地看到，本发明算法实现了正确的配准，而参考算法对图像的变形在左下角是错误的。

Claims (1)

1、一种基于有限采样全局优化的图像弹性配准方法，该方法由以下步骤组成：

(1)读入变形图像，用快速傅立叶变换计算它的频谱，再由该频谱计算它的最高频率；

(2)用傅立叶变换方法计算目标函数相对于变形函数中的变形系数的频谱，再由该频谱计算得到目标函数的最高频率；

(3)先按步骤2得到的目标函数的最高频率的2倍对变形系数对应的空间进行均匀采样计算，得到变形函数系列，再用得到的变形函数系列分别对变形图像进行变形，得到变形图像系列；

(4)读入参考图像，分别计算参考图像与步骤(3)得到的变形图像系列的均方差，得到目标函数值数据组，比较这些目标函数值的大小，找到其中最小目标函数值，用该最小目标函数值所对应的变形函数对变形图像进行变形，得到粗配准的变形图像；

(5)以步骤(4)所找到的最小的目标函数值所对应的变形系数值为起始点，用梯度下降法搜索目标函数的全局最小值，再用该最小值所对应的变形函数对变形图像进行变形，便得到精配准的变形图像；

其中，所述变形函数和目标函数的定义如下所述：

(1)用B样条构造变形函数：

$g(x,C)=x+\underset{i{\∈I}_C}{\mathrm{\Σ}}{C}_i{\mathrm{\β}}_n(x/2^w-i)$

其中C为变形系数，β_n是n次B样条的向量积，即 ${\mathrm{\β}}_n\left(x\right)={\mathrm{\Π}}_{k=1}^m{\mathrm{\β}}_n\left(x_k\right);$

(2)设参考图像和变形图像分别定义为f_r(x)和f_t(x)，以参考图像和用变形函数变形后的变形图像(它可以表示为f_t(g(x，C)))的均方差作为目标函数：

$E\left(C\right)=\frac{1}{\left|\right|I\left|\right|}\underset{x\∈I}{\mathrm{\Σ}}{e}_i^2=\frac{1}{\left|\right|I\left|\right|}\underset{x\∈I}{\mathrm{\Σ}}{(f_t\left(g(x,C)\right)-f_r\left(x\right))}^2$

其中I是图像坐标所对应的空间，g(x，C)是变形函数。

Images