CN104599268A - 一种联合点配准的局部区域精确变形配准算法 - Google Patents

Abstract

本发明公开了一种联合点配准的局部区域精确变形配准算法，本发明方法通过将基于特征的点配准方法和基于强度的图像配准方法结合在一起，尤其是通过B样条模型将点配准的结果转化成变形场，并以该变形场作为后续基于强度图像配准的初始变形场，从而将两种不同方式的配准方法结合在一起，巧妙地利用了点配准配准精度较高的特点，提高了感兴趣区域的配准精度，从而提高受照累积剂量的精度，实现了对图像中某些局部区域进行配准的精度干涉及调控。本发明作为一种联合点配准的局部区域精确变形配准算法可广泛应用于图像变形配准领域。

Description

一种联合点配准的局部区域精确变形配准算法

技术领域

本发明涉及图像变形配准，尤其是一种联合点配准的局部区域精确变形配准算法。

背景技术

常规的放射治疗是以第一次获取的病人计算机断层图像(CT，Computed Tomography)作为计划图像，并只以该计划CT图像生成治疗计划，对病人实施治疗。但实际上，在每一次分次治疗中，病人的内部解剖结构都会发生改变，如技师摆位的误差、病人内部脏器的蠕动、肿瘤随着放疗的进行体积会缩小等，如果仍基于首次计划实施放疗，必将带来严重的剂量误差，并有可能使正常组织卷入高剂量区域。自适应放射治疗(ART，Adaptive Radiotherapy)技术可以解决上述问题。ART技术在治疗分次中获取当次治疗的病人CT图像，并基于该图像重新设计计划，对病人实施治疗。而图像变形配准技术(DIR，Image Deformable Registration)是ART的一项关键技术，它主要用于将计划CT图像上的轮廓线推衍至当前CT图像，避免物理师重新勾画轮廓线的麻烦；另一方面，变形配准可以将各分次的剂量分布变形到同一解剖空间后进行叠加，从而评价病人所受的总剂量照射。

图像变形配准本质上是要寻找一种点对点的空间变换，将一幅图像(浮动图像)的解剖结构变形到另一幅图像(参考图像)上。目前变形配准算法可分为基于物理模型变形配准算法：如基于线性弹性模型的，基于流体力学模型的，基于光流场模型的；也有基于数学函数的变形配准算法：如基于小波变换的，基于薄板样条函数或B样条函数的，基于调和函数的等等。但这些配准算法的精度只和算法模型有关，配准的精度是全局的，它与待配准的图像中某个具体的解剖结构没有关系，即无法通过调节算法参数提高图像中某个感兴趣区域的配准精度。

但在ART中，我们往往要求某些感兴趣的区域(Region ofInterest，ROI)，例如肿瘤靶区或危及器官(Organ At Risk，OAR)，有较高的配准精度，因为这些区域的配准精度直接关系到剂量计算以及计划重设计的精度，关系到靶区能否接受准确的剂量照射以及危及器官是否受到过量的照射，这与放疗的疗效及愈后有着密切的联系。

发明内容

为了解决上述技术问题，本发明的目的是：提供一种能提高图像中某些感兴趣区域配准精度的CT图像间的变形配准的方法。

本发明所采用的技术方案是：一种联合点配准的局部区域精确变形配准算法，包括有以下步骤：

A、在浮动图像I_m和参考图像I_s中分别记录同一个感兴趣区域的轮廓线的三维浮动点集和三维参考点集

B、对三维浮动点集和三维参考点集利用基于特征点配准算法进行非刚性点配准；

C、根据步骤B中计算得到的非刚性点配准后的浮动点集及配准前的浮动点集通过B样条模型估算出变形场Δ；

D、根据变形场Δ对浮动图像I_m进行形变，得到初步变形后的浮动图像I′_m，然后基于图像强度的配准算法对初步变形后的浮动图像I′_m和参考图像I_s进行变形配准得到配准变形场u；

E、计算最终变形场u+Δ，并得到最终变形的浮动图像I_m ^o(Δ+u)。

进一步，所述步骤B中使用基于薄板样条的TPS-RPM点配准模型，并采用确定性退火算法作为匹配的优化算法，分步交替求解对应度矩阵最优值和变换函数的最优变换参数。

进一步，所述步骤C中计算变形场Δ是根据一维B样条基函数以及控制网格点φ_i，j，k通过迭代方式求得，具体包括以下步骤：

C1、初始化均匀的B样条控制网格φ_init、残差δ以及最大迭代次数N₂；

C2、每一次迭代过程中，根据残差计算控制网格点增量并更新控制网格点，计算变换函数然后计算新的浮动点集更新残差；

C3、判断如果残差δ大于预设阈值ε₁或迭代小于最大迭代次数，继续步骤C2的迭代过程，否则迭代终止。

进一步，所述步骤D中采用基于光流场模型的Demons变形配准算法。

进一步，所述步骤D具体包括有以下子步骤：

D1、根据上一次迭代变形的浮动图像、参考图像上一次迭代变形的浮动图像的梯度场和参考图像的梯度场计算配准变形场增量；

D2、将本次迭代的配准变形场增量与上一次迭代变形场相加，得到本次的变形场；

D3、根据本次的变形场对上一次迭代得到的浮动图像进行变形，所得到的变形图像作为下一次迭代的浮动图像；

D4、重复执行上述步骤D1-D3得到配准变形场u。

进一步，所述步骤D4中的迭代终止条件为：终止参数小于等于1.0×10^-14。

本发明的有益效果是：本发明方法通过将基于特征的点配准方法和基于强度的图像配准方法结合在一起，尤其是通过B样条模型将点配准的结果转化成变形场，并以该变形场作为后续基于强度图像配准的初始变形场，从而将两种不同方式的配准方法结合在一起，巧妙地利用了点配准配准精度较高的特点，提高了感兴趣区域的配准精度，从而提高受照累积剂量的精度，实现了对图像中某些局部区域进行配准的精度干涉及调控。

附图说明

图1为本发明的主步骤流程图；

图2为本发明实施例中TPS-RPM点配准之前施源器区域的轮廓线点集重叠图；

图3为本发明实施例中TPS-RPM点配准之后施源器区域的轮廓线点集重叠图。

具体实施方式

下面结合附图对本发明的具体实施方式作进一步说明：

参照图1，一种联合点配准的局部区域精确变形配准算法，包括有以下步骤：

A、在浮动图像I_m和参考图像I_s中分别记录同一个感兴趣区域的轮廓线的三维浮动点集和三维参考点集其中K和N是点集中点的数目；

B、对三维浮动点集和三维参考点集利用基于特征点配准算法进行非刚性点配准；

进一步作为优选的实施方式，所述步骤B中使用基于薄板样条的TPS-RPM点配准模型，并采用确定性退火算法作为匹配的优化算法，分步交替求解对应度矩阵最优值和变换函数的最优变换参数。

其具体的优化过程需要交替迭代进行，在每次迭代的过程都要更新对应度矩阵和变换函数，其具体算法流程如下：

首先进行初始化：m_ai、f分别表示对应矩阵与薄板样条变换函数；T表示当次迭代的退火温度，T₀和T_final表示初始退火温度和最终退火温度，r表示退火速度；初始化薄板样条变化参数w＝0，d＝I；λ表示正则化参数；ε表示算法终止阈值，N表示算法最大迭代次数。

以下交替迭代进行的计算过程可参照以下伪代码段：

C、根据步骤B中计算得到的非刚性点配准后的浮动点集及配准前的浮动点集求解得到变形场(即变形场Δ) $\mathrm{\Δ}=$ $\{\stackrel{\→}{\mathrm{\Δ}}=({\mathrm{\Δ}}^x,{\mathrm{\Δ}}^y,{\mathrm{\Δ}}^z)|0\≤x\≤X,0\≤y\≤Y,0\≤z\≤Z\},$ 该过程可通过估算在图像域Ω＝(x，y，z)|0≤x≤X，0≤y≤Y，0≤z≤Z上一个大小为n_x×n_y×n_z的领域，且控制网格点φ_i，j，k间距为d_x，d_y，d_z的B样条网格来实现；

所述步骤C中根据一维B样条基函数以及控制网格点φ_i，j，k求得变形场Δ的计算公式为：

$\stackrel{\→}{\mathrm{\Δ}}\≡\mathrm{\ψ}(x,y,z)={\mathrm{\Σ}}_{l=0}^3{\mathrm{\Σ}}_{m=0}^3{\mathrm{\Σ}}_{n=0}^3{B}_l\left(u\right)B_m\left(v\right)B_n\left(w\right){\mathrm{\φ}}_{i+l,j+m,k+n}$

(式10)

其中：

(式11)

B_l，B_m和B_n为B样条基函数：

B₀(u)＝(1-u)³/6，

B₁(u)＝(3u³-6u²+4)/6，

B₂(u)＝(-3u³-3u²+3u+1)/6，

B₃(u)＝u³/6    (式12)

定义和的残差为 $\mathrm{\δ}=\{{\stackrel{\→}{\mathrm{\δ}}}_i=({\mathrm{\δ}}_i^x,{\mathrm{\δ}}_i^y,{\mathrm{\δ}}_i^z)|i=$ $\mathrm{1,2},...,K\},$ 并定义δ的一个子集δ′使得δ′中的元素位于控制网格点φ_i，j，k的4×4×4邻域中：

(式13)

假设只有上述残差δ′中的邻域点影响φ_i，j，k的大小，那么最小化残差δ的最小二乘解可由下式近似得到：

${\mathrm{\φ}}_{i,j,k}=\frac{{\mathrm{\Σ}}_c{w}_c{\mathrm{\φ}}_c}{{\mathrm{\Σ}}_c{w}_c^2}$     (式14)

其中：

${\mathrm{\φ}}_c=\frac{w_c\stackrel{\→}{{\mathrm{\δ}}^{\′}}}{{\mathrm{\Σ}}_{a=0}^3{\mathrm{\Σ}}_{b=0}^3{\mathrm{\Σ}}_{c=0}^3{w}_{\mathrm{abc}}},$     (式15)

w_c＝w_lmn＝B_l(v)B_m(u)B_n(w)，   (式16)

(式17)

进一步作为优选的实施方式，所述步骤C中计算变形场Δ是根据一维B样条基函数以及控制网格点φ_i，j，k通过迭代方式求得，具体包括以下步骤：

C1、初始化均匀的B样条控制网格φ_init、残差δ以及最大迭代次数N₂；

其中初始化残差 $\mathrm{\δ}=f\left({\stackrel{\→}{v}}_a\right)-{\stackrel{\→}{v}}_a;$

C2、每一次迭代过程中，根据残差计算控制网格点增量并更新控制网格点，计算变换函数然后计算新的浮动点集更新残差；

C3、判断如果残差δ大于预设阈值ε₁或迭代小于最大迭代次数，继续步骤C2的迭代过程，否则迭代终止。

上述步骤C2和C3采用迭代的方式来计算φ_i，j，k，即通过每次迭代更新后的φ_i，j，k计算变换函数ψ，再将变换至然后计算和的残差，如果残差||δ||₂≤ε₁，ε₁＝1.0×10^-3或迭代次数大于最大迭代次数N₂，迭代终止。

迭代过程可参照如下伪代码段：

D、根据变形场Δ对浮动图像I_m进行形变，得到初步变形后的浮动图像I′_m，然后基于图像强度的配准算法对初步变形后的浮动图像I′_m和参考图像I_s进行变形配准得到配准变形场u；其中浮动图像的计算公式为：I′_m＝I_m ^oΔ。

进一步作为优选的实施方式，所述步骤D中采用基于光流场模型的Demons变形配准算法。

进一步作为优选的实施方式，所述步骤D具体包括有以下子步骤：

D1、根据上一次迭代变形的浮动图像、参考图像上一次迭代变形的浮动图像的梯度场和参考图像的梯度场计算配准变形场增量；

具体计算公式为：

${\mathrm{du}}^{(k+1)}=\frac{(I_m^{\′\left(k\right)}-I_s){\▿I}_s}{{(I_m^{\′\left(k\right)}-I_s)}^2+{\left|{\▿I}_s\right|}^2}+\frac{(I_m^{\′\left(k\right)}-I_s){\▿I}_m^{\′\left(k\right)}}{{(I_m^{\′\left(k\right)}-I_s)}^2+{\left|{\▿I}_m^{\′\left(k\right)}\right|}^2}$     (式18)

上式中，du^(k+1)为配准变形场增量，为上次迭代变形的浮动图像，I_s为参考图像，为上次迭代变形的浮动图像的梯度场，为参考图像的梯度场，k为迭代循环次数的序号。

D2、将本次迭代的配准变形场增量与上一次迭代变形场相加，得到本次的变形场：

u^(k+1)＝u^(k)+du^(k+1)

上式中，u^(k+1)为配准变形场，du^(k+1)为配准变形场增量，k为迭代循环次数的序号。

D3、根据本次的变形场u^(k+1)对上一次迭代得到的浮动图像进行变形，所得到的变形图像作为下一次迭代的浮动图像；

D4、重复执行上述步骤D1-D3得到配准变形场u。

进一步作为优选的实施方式，所述步骤D4中的迭代终止条件为：终止参数ε_r小于等于1.0×10^-14。其中，ε_r表示为：

ε_r＝∑|du^(k-10)|/∑|u^(k-10)|-∑|du^(k)|/∑|u^(k)|

E、计算最终变形场u+Δ，并得到最终变形的浮动图像I_m ^o(Δ+u)。

以下结合图2-3说明本发明的一具体实施例，

读入CT₁图像和CT₂图像，CT₂为浮动图像，CT₂为参考图像。该实例是宫颈癌病人实施高剂量率近距离放疗期间获取的两套CT图像，在病人接受高剂量率近距离放疗时，需在其阴道内放置施源器，CT₁图像和CT₂图像是通过人工手动分割去除施源器后的CT图像，并对该区域赋予了空气的CT值(-1000HU)，目的是要通过变形配准，使得两幅CT图像形状匹配，并且使得施源器临近区域(临近区域是高剂量区域)得到较高的配准精度。两幅图像的分辨率均为256×256×54。CT₁图像和CT₂图像在分割出的施源器区域的轮廓线三维点集分别为 ${\stackrel{\→}{v}}_a,(a=\mathrm{1,2},...,K)$ 和 ${\stackrel{\→}{x}}_i,(i=\mathrm{1,2},...,N),$ 其中K＝2424，N＝2016，TPS-RPM点配准之前施源器区域的轮廓线如图2所示，轮廓线三维点集和的匹配程度不高。

采用基于薄板样条的TPS-RPM点配准算法对点集v_a和x_i进行非刚性点配准。首先初始化参数：初始退火温度T₀＝0.0002和最终退火温度T_final＝0.00002；退火速度r＝0.93；初始化薄板样条变化参数w＝0，d＝I；正则化参数λ＝1；最大迭代次数N＝5；算法终止阈值ε＝1.0×10^-3；浮动点集和参考点集的点数目分别为K＝2424，N＝2016。按前述第一段伪代码段中的算法迭代完成浮动点集和参考点集的非刚性配准。经过TPS-RPM点配准后，得到配准后的浮动点集TPS-RPM点配准之后施源器区域的轮廓线如图3所示，轮廓线三维点集和的匹配程度与图2明显不同。

初始化残差最大迭代次数N₂＝4，B样条控制网格的间距为d_x＝d_y＝d_z＝16，大小为n_x＝256，n_yn_x×n_y×n_z＝256，n_z＝54，按照前述第二段伪代码段中的算法计算变形场Δ。

利用步骤D计算的变形场Δ对CT₁进行形变，得到初步变形的浮动图像CT′₁＝CT₁ ^oΔ，然后采用于光流场模型的Demons变形配准算法对CT′₁和CT₃进行变形配准，得到配准变形场u，最终变形场等于u+Δ，并得到最终变形的浮动图像为CT₁ ^o(Δ+u)。

以上是对本发明的较佳实施进行了具体说明，但本发明创造并不限于所述实施例，熟悉本领域的技术人员在不违背本发明精神的前提下还可以作出种种的等同变换或替换，这些等同的变形或替换均包含在本申请权利要求所限定的范围内。

Claims (6)

1.一种联合点配准的局部区域精确变形配准算法，其特征在于：包括有以下步骤：

A、在浮动图像和参考图像中分别记录同一个感兴趣区域的轮廓线的三维浮动点集和三维参考点集；

B、对三维浮动点集和三维参考点集利用基于特征点配准算法进行非刚性点配准；

C、根据步骤B中计算得到的非刚性点配准后的浮动点集及配准前的浮动点集，通过B样条模型估算出变形场；

D、根据变形场对浮动图像进行形变，得到初步变形后的浮动图像，然后基于图像强度的配准算法对初步变形后的浮动图像和参考图像进行变形配准得到配准变形场；

E、计算最终变形场，并得到最终变形的浮动图像。

2.根据权利要求1所述的一种联合点配准的局部区域精确变形配准算法，其特征在于：所述步骤B中使用基于薄板样条的TPS-RPM点配准模型，并采用确定性退火算法作为匹配的优化算法，分步交替求解对应度矩阵最优值和变换函数的最优变换参数。

3.根据权利要求1所述的一种联合点配准的局部区域精确变形配准算法，其特征在于：所述步骤C中计算变形场是根据一维B样条基函数以及控制网格点通过迭代方式求得，具体包括以下步骤：

C1、初始化均匀的B样条控制网格、残差以及最大迭代次数；

C2、每一次迭代过程中，根据残差计算控制网格点增量并更新控制网格点，计算变换函数然后计算新的浮动点集，更新残差；

C3、判断如果残差大于预设阈值或迭代小于最大迭代次数，继续步骤C2的迭代过程，否则迭代终止。

4.根据权利要求1所述的一种联合点配准的局部区域精确变形配准算法，其特征在于：所述步骤D中采用基于光流场模型的Demons变形配准算法。

5.根据权利要求4所述的一种联合点配准的局部区域精确变形配准算法，其特征在于：所述步骤D具体包括有以下子步骤：

D1、根据上一次迭代变形的浮动图像、参考图像上一次迭代变形的浮动图像的梯度场和参考图像的梯度场计算配准变形场增量；

D2、将本次迭代的配准变形场增量与上一次迭代变形场相加，得到本次的变形场；

D3、根据本次的变形场对上一次迭代得到的浮动图像进行变形，所得到的变形图像作为下一次迭代的浮动图像；

D4、重复执行上述步骤D1-D3得到配准变形场。

6.根据权利要求5所述的一种联合点配准的局部区域精确变形配准算法，其特征在于：所述步骤D4中的迭代终止条件为：终止参数小于等于1.0×10^-14。