CN102024256B - 一种基于梯度场的变约束图像变形配准方法 - Google Patents

Abstract

本发明公开了一种基于梯度场的变约束图像变形配准方法，在来自参考图象与目标图像两个变形力的驱动下更快速的收敛，通过配准过程中变约束的策略，使算法更少的陷入局部最小，配准结果更准确。

Description

一种基于梯度场的变约束图像变形配准方法

技术领域

本发明涉及一种基于梯度场的变约束图像变形配准方法。

背景技术

现代放射治疗技术，比如调强适形放疗，可以将高度适形的剂量分布投射到静止靶区。然而，放射治疗过程通常采用分割照射技术，即放疗需要分若干次进行。分次放疗过程中，肿瘤及其周围重要器官的位置和形状可能产生放疗分次间变化(如肿块的增大或退缩)和分次内变化(如呼吸等生理运动引起的肿块位置和形状变化)。基于对这些运动变化的考虑，常规的做法是将肿瘤区外放一定范围作为照射靶区，从而造成更多的正常组织受到放射损伤。

放疗前，通常需要获取图像质量高的扇形束CT(Fan-Beam CT，FBCT)用来制定放疗计划。分次放疗实施前，在线采集的锥形束CT(Cone-beam CT，CBCT)影像包含了当前肿瘤和周围器官的位置和形状信息。通过对FBCT和在线CBCT图像进行快速准确的变形配准，可以获取肿瘤及其周围器官的运动和变形信息，并据此修改先前的放疗计划，从而减少考虑到运动变化而增加的不必要的照射范围，减少肿瘤周围正常组织的放射损伤。

变形配准问题可以表达为，已知参考图像R(x)，目标图像M(x′)，求参考图像到目标图像的变形场u(x)，使得在变形场作用下的目标图像M′(x)＝M(x+u(x)与参考图像R(x)中的对应对象相吻合。

当前的FBCT和CBCT图像的变形配准方法主要分为两类，一类方法基于目标的轮廓，这类方法需要在配准前对CBCT中的目标进行人工或自动的分割以获取目标的轮廓。配准的准确性也高度依赖于图像分割的准确性。另一类方法基于图像的灰度值，但是由于电子散射，CBCT图像的CT值并不十分准确，受其影响这些基于图像灰度值的配准算法得出的结果也不十分准确。本发明提出了一种基于梯度场的快速准确的图像变形配准算法，在来自于参考图像和目标图像两个变形力的驱动下更快速的收敛，通过在配准过程中实施变约束的策略，使算法更少的陷入局部最小，配准结果更准确。

发明内容

因此，本发明的目的在于提供一种配准结果更准确的基于梯度场的变约束图像变形配准方法。

本发明采用以下技术方案：

该发明基于梯度场的变约束图像变形配准方法，其包括以下步骤：

1)给定参考图象，记为R(x)，并给定该参考图象的目标图像，记为M(x′)；

2)对所述参考图象和目标图像中带有噪声的图像进行平滑滤波；

3)在参考图象与目标图像之间进行刚体配准；

4)提取参考图象和目标图像的边缘信息并计算各自边缘信息的梯度场或者直接计算参考图象和目标图像的梯度场；

5)基于所获得的所述梯度场，通过在不同的配准阶段对不同的配准对象设置不同的约束值进行参考图象到目标图像变形场的迭代优化；

6)迭代优化完成后，将得到的参考图像到目标图像的变形场输出。

依据本发明技术方案的基于梯度场的变约束图像变形配准方法，在来自参考图象与目标图像两个变形力的驱动下更快速的收敛，通过配准过程中变约束的策略，使算法更少的陷入局部最小，配准结果更准确。

上述配准方法，该配准方法、该配准方法中的参考对象和目标对象都是基于向量描述的，该配准方法对向量所包含的所有分量进行配准。当向量为一维向量时，该配准算法进行一维图像变形配准。当向量为二维向量时，该配准算法进行二维图像变形配准。当向量为更高维向量时，该配准算法进行对应的更高维的图像变形配准。

上述配准方法，所述步骤2)中的平滑滤波采用平滑滤波器。

上述配准方法，所述平滑滤波器优选高斯滤波器。

上述配准方法，所述刚体配准包括包括基于图像特征的刚体配准和基于优化图像相关性或互信息指标的刚体配准。

上述配准方法，所述步骤4)中提取参考图像和目标图像边缘信息的方法采用Canny边缘提取算子。

上述配准方法，在步骤5)中进行迭代优化生成变形场的迭代公式为：

x′＝x+u(x，i)                (1)

${\mathrm{GD}}_n=\frac{\∂\left|\right|\▿R\left(x\right)-\▿M\left(x^{\′}\right)\left|\right|}{\∂{x^{\′}}_n}---\left(2\right)$

${\mathrm{GS}}_n=\frac{\∂\left|\right|\▿R\left(x\right)-\▿M\left(x^{\′}\right)\left|\right|}{\∂x_n}---\left(3\right)$

$u_n(x,i+1)=u_n(x,i)+\mathrm{\α}\×\frac{\mathrm{\λ\Δ}{u}_n(x,i)-\left|\right|\▿R\left(x\right)-\▿M\left(x^{\′}\right)\left|\right|\×({\mathrm{GD}}_n-{\mathrm{GS}}_n)}{\mathrm{\λ}+{({\mathrm{GD}}_n-{\mathrm{GS}}_n)}^2}---\left(4\right)$

式中R(·)表示参考图像中的像素值；M(·)表示运动图像中的像素值；x是参考图像中像素的坐标；x′是迭代优化过程中参考图像上像素x在运动图像上对应像素的坐标；n表示向量的第n维分量；u表示从参考图像到运动图像的变形场；u_n(x，i)表示参考图像上像素x在第i次迭代后第n维的变形量；符号Δ表示拉普拉斯算子；符号

表示梯度；其中项λΔu_n(x，i)是约束平滑变形场的项；约束图像变形后保持像素之间的相对位置关系；α是步长因子；λ是变形约束参数；项

是驱动变形场变形的项；GD_n是驱动变形场的动态变形力；GS_n是驱动变形场的静态变形力。

上述配准方法，迭代优化采用象素值从粗到精的多分辨率方法。

上述配准方法，所述步骤5)中在迭代优化过程中采用先设置数值较小的变形约束参数λ进行迭代优化，然后使用较大的变形约束参数λ进行迭代优化的变约束策略。

上述配准方法，在迭代优化过程中所述优选变约束策略可以重复多次使用。

附图说明

图1为依据本发明技术方案的配准方法的流程图。

具体实施方式

参照说明书附图1，一种基于梯度场的变约束图像变形配准方法，其包括以下步骤：

1)给定参考图象，记为R(x)，并给定该参考图象的目标图像，记为M(x′)；

2)对所述参考图象和目标图像中带有噪声的图像进行平滑滤波，以减少噪声对配准结果造成的干扰；

3)在参考图象与目标图像之间进行刚体配准；由于变形图像配准属于局部优化算法，对于目标图像有较大总体位移和旋转的情况，需要先进行刚体配准；

4)提取参考图象和目标图像的边缘信息并计算各自边缘信息的梯度场或者直接计算参考图象和目标图像的梯度场；

5)基于所获得的所述梯度场，通过在不同的配准阶段对不同的配准对象设置不同的约束值进行参考图象到目标图像变形场的迭代优化；

6)迭代优化完成后，将得到的参考图像到目标图像的变形场输出。

该配准方法和该配准方法中的参考对象和目标对象都是基于向量描述的。由于是基于向量的描述，所以不再限定参考图象与目标图像的维数，参考图象和目标图像可以是一维、二维、三维或者更高维的图像。

梯度场可以直接计算参考图像和目标图像的梯度场。也可以先对参考图像和目标图像提取边缘信息，比如可以采用canny边缘提取算子等提取边缘信息。然后再对边缘信息图像计算梯度场。

关于目标图像与参考图象中带有噪声的图象的平滑滤波采用普通的平滑滤波器就可以满足要求，仅用于消除噪音，对平滑滤波器没有特别的要求。

优选地，所述平滑滤波器优选高斯滤波器，可以反映更真实的信号。

刚体配准在本方案中要求不高，但又必不可少，即使存在明显的刚体配准误差，但只要不是太大就不会对最终的配准结果产生有意义的影响。在后面的FBCT与CBCT间三维变形图像配准的实施例中不大于7个像素的误差通常不会对最终的配准结果产生有意义的影响，因此，许多刚体配准算法都可以满足本方案的要求。优选地，所述刚体配准包括基于图像特征的刚体配准和基于优化图像相关性或互信息指标的刚体配准。

对于提取边缘信息，可以采用微分滤波器，比如一阶微分滤波器之梯度算子，如Prewitt和Sobel算子；二阶微分滤波器，如LoG算子；或ratio边缘检测算子；也可以采用已经非常成熟的Canny算子，即在所述步骤4)中提取参考图像和目标图像边缘信息的方法采用Canny边缘提取算子，该算子最早见于Canny，J.，AComputational Approach To Edge Detection，IEEE Trans.Pattern Analysis and Machine Intelligence，8(6)：679-698，1986，是一种多级边缘检测算法，该算法适用于不同场合。它的参数允许根据不同实现的特定要求进行调整以识别不同的边缘特性。

在步骤4)中计算出梯度场后，进行迭代优化生成变形场：x′＝x+u(x，i)

(1)

${\mathrm{GD}}_n=\frac{\∂\left|\right|\▿R\left(x\right)-\▿M\left(x^{\′}\right)\left|\right|}{\∂{x^{\′}}_n}---\left(2\right)$

${\mathrm{GS}}_n=\frac{\∂\left|\right|\▿R\left(x\right)-\▿M\left(x^{\′}\right)\left|\right|}{\∂x_n}---\left(3\right)$

$u_n(x,i+1)=u_n(x,i)+\mathrm{\α}\×\frac{\mathrm{\λ\Δ}{u}_n(x,i)-\left|\right|\▿R\left(x\right)-\▿M\left(x^{\′}\right)\left|\right|\×({\mathrm{GD}}_n-{\mathrm{GS}}_n)}{\mathrm{\λ}+{({\mathrm{GD}}_n-{\mathrm{GS}}_n)}^2}---\left(4\right)$

其中R(·)表示参考图像中的像素值，M(·)表示运动图像中的像素值，x是参考图像中像素的坐标，x′是迭代优化过程中参考图像上像素x在运动图像上对应像素的坐标，n表示向量的第n维分量，u表示从参考图像到运动图像的变形场，u_n(x，i)表示参考图像上像素x在第i次迭代后第n维的变形量。符号Δ表示拉普拉斯算子，符号

表示梯度。其中项λΔu_n(x，i)是约束平滑变形场的项，用来保证变形场的平滑，连续，约束图像变形后保持像素之间的相对位置关系。α是步长因子；λ是变形约束参数，值越大，此平滑项的作用越大，图象越不易变形。项是驱动变形场变形的项，用来驱动变形场变形。其中GD_n是驱动变形场的动态变形力，GS_n是驱动变形场的静态变形力。

上面给出了公式和部分参数与本方案相关的作用，可以根据需要进行选择。此迭代优化策略通过约束平滑项λΔu_n(x，i)，来保证变形场的平滑、连续，约束图像变形后保持像素之间的相对位置关系；通过在目标图像中的动态变形力GD_n和在参考图像中的静态变形力GS_n的共同作用，驱动变形场变形。

迭代优化采用象素值从粗到精的多分辨率方法。多分辨率策略可以有效地减少整个配准过程的计算量，对于较大的变形还可以避免陷入某些局部最小点。

所述步骤5)中在迭代优化过程中采用先设置数值较小的变形约束参数λ进行迭代优化，然后使用较大的变形约束参数λ进行迭代优化的变约束策略。初始较小的变形约束参数虽然会导致图像严重扭曲变形，但可以使优化更少的陷入局部最小点。在后期重新设置较大的变形约束参数，使图像重新恢复原有的位置约束。

整个优选的变约束策略的过程是：在配准的起始阶段设置较小的λ值，参考图象在较小约束的变形场作用下严重变形，图象看起来变得混乱无序，甚至相邻的像素会失去先前的相对位置关系。在配准的最后阶段设置较大λ值，图像在较强约束变形场的作用下渐渐重新变的清晰起来。

在迭代优化过程中所述优选变约束策略可以重复多次使用。

本发明提出的图像变形配准算法，不限制具体的实施硬件与环境，本算法可以通过串行指令在串行处理器中实施；也可以通过并行编程在并行处理器和计算机网络中实施。

下面提供了此方法的两个具体实施例。分别为此方法应用于FBCT影像与CBCT影像间的三维变形图像配准的一个具体实施例和此方法应用于FBCT影像与磁共振(MR)影像间的三维变形图像配准的一个具体实施例。

此方法应用于FBCT影像与CBCT影像间的三维变形图像配准的一个具体实施例：

此实施例中FBCT图像作为参考图像，CBCT图像作为目标图像。

首先应用高斯滤波函数对包含较多伪影和噪声的CBCT图像进行滤波。

然后采用基于密度的刚体配准，输出三个刚体平移量。目标函数为：

$f\left(u\right)=\underset{\mathrm{pi}=1}{\overset{\mathrm{Number}}{\mathrm{\Σ}}}{(R\left(x\right)-M(x+u))}^2---\left(5\right)$

这里u是沿各个坐标轴的刚体平移向量，pi表示当前像素的编号，Number表示总共像素的个数。采用贪婪算法优化目标函数到最小值，输出平移向量u。

在刚体配准完成后，直接计算FBCT和CBCT的梯度场，用于变形图像配准。变形图像配准根据公式(1)(2)(3)(4)迭代生成变形场。变形场每个像素的初始值为刚体配准中求得的平移向量u。在迭代优化过程中图像平面内x，y轴方向降采样率依次分别设为8、4、2、1。当降采样率为8时，约束参数λ设置为20000；当降采样率为4时，约束参数λ设置为40000；当降采样率为2时，约束参数λ设置为80000；当降采样率为1时，约束参数根据FBCT图像中当前像素CT值设置不同的值，当CT值大于1150时，认为是骨性物质，约束参数λ设置为500000，当CT值小于1150大于400时，认为是软组织，约束参数λ设置为250000，CT值小于400时，约束参数设置为120000。与图像平面垂直的z轴方向没有降采样。整个过程中步长因子α恒为0.4，每个采样率迭代32次后转入下一级采样率的迭代直到采样率为1的所有迭代完成后，输出变形场。

另一个实施例，此方法应用于FBCT影像与MR影像间的三维变形图像配准的一个具体实施例：

此实施例中FBCT图像作为参考图像，MR图像作为目标图像。

由于FBCT图像和MR图像的噪声都较小，没有进行平滑滤波。

首先采用基于像素值互信息的刚体配准，输出三个刚体平移量和绕三个坐标轴的旋转量。参考文献”Multimodality Image Registration byMaximization of Mutual Information”IEEE Transactions on Medical ImagingVol.16.No.2.1997p187-198。

在刚体配准完成后，采用Canny边缘检测算子计算FBCT和MR的边缘图像。

计算FBCT和MR边缘图像的梯度场，用于变形图像配准。变形图像配准根据公式(1)(2)(3)(4)迭代生成变形场。由刚体配准的结果确定变形场每个像素的初始值。在图像平面内x，y轴方向降采样率依次分别为8，4，2，1。在此实施例中降采样率为8、4、2时，反复使用强约束和弱约束进行优化迭代。当降采样率为8，迭代次数为奇数时，约束参数λ设置为200；当降采样率为8，迭代次数为偶数时，约束参数λ设置为250000；当降采样率为4，迭代次数为奇数时，约束参数λ设置为800；当降采样率为4，迭代次数为偶数时，约束参数λ设置为250000；当降采样率为2，迭代次数为奇数时，约束参数λ设置为1000；当降采样率为2，迭代次数为偶数时，约束参数λ设置为250000；当降采样率为1时，约束参数λ设置为10000。与图像平面垂直的z轴方向没有降采样。整个过程中步长因子α恒为0.4，每个采样率迭代32次后转入下一级采样率的迭代直到采样率为1的所有迭代完成后，输出变形场。

Claims (9)

1.一种基于梯度场的变约束图像变形配准方法，其特征在于其包括以下步骤：

1)给定参考图像，记为R(x)，并给定该参考图像的目标图像，记为M(x′)；

2)对所述参考图像和目标图像中带有噪声的图像进行平滑滤波；

3)在参考图像与目标图像之间进行刚体配准；

4)提取参考图像和目标图像的边缘信息并计算各自边缘信息的梯度场或者直接计算参考图像和目标图像的梯度场；

5)基于所获得的所述梯度场，通过在不同的配准阶段对不同的配准对象设置不同的约束值进行参考图像到目标图像变形场的迭代优化；

6)迭代优化完成后，将得到的参考图像到目标图像的变形场输出；

其中，在步骤5)中进行迭代优化生成变形场的迭代公式为：

x′＝x+u(x，i)                      (1)

式中R(·)表示参考图像中的像素值；M(·)表示目标图像中的像素值；x是参考图像中像素的坐标；x′是迭代优化过程中参考图像上像素x在目标图像上对应像素的坐标；n表示向量的第n维分量；u表示从参考图像到目标图像的变形场；u_n(x，i)表示参考图像上像素x在第i次迭代后第n维的变形量；符号Δ表示拉普拉斯算子；符号 

表示梯度；其中项λΔu_n(x，i)是约束平滑变形场的项；约束 图像变形后保持像素之间的相对位置关系；α是步长因子；λ是变形约束参数；项 

是驱动变形场变形的项；GD_n是驱动变形场的动态变形力；GS_n是驱动变形场的静态变形力。

2.根据权利要求1所述的配准方法，其特征在于：该配准方法、该配准方法中的参考对象和目标对象都是基于向量描述的，该方法对向量所包含的所有分量进行配准。

3.根据权利要求1所述的配准方法，其特征在于：所述步骤2)中的平滑滤波采用平滑滤波器。

4.根据权利要求3所述的配准方法，其特征在于：所述平滑滤波器优选高斯滤波器。

5.根据权利要求1所述的配准方法，其特征在于：所述刚体配准包括基于图像特征的刚体配准和基于优化图像相关性或互信息指标的刚体配准。

6.根据权利要求1所述的配准方法，其特征在于：所述步骤4)中提取参考图像和目标图像边缘信息的方法采用Canny边缘提取算子。

7.根据权利要求1所述的配准方法，其特征在于：迭代优化采用像素值从粗到精的多分辨率方法。

8.根据权利要求7所述的配准方法，其特征在于：所述步骤5)中在迭代优化过程中采用先设置数值较小的变形约束参数λ进行迭代优化，然后使用较大的变形约束参数λ进行迭代优化的优选变约束策略。

9.根据权利要求8所述的配准方法，其特征在于：在迭代优化过程中所述优选变约束策略可以重复多次使用。 

Images