A contour-guided deformable image registration algorithm for adaptive radiotherapy

The goal of demons registration is to have the intensity matching between static and moving image through deforming moving image. The demons algorithm can be cast into a generalized optimization framework (Vercauteren et al 2009):

$$\begin{eqnarray} &&E( {\bf u} ) = \frac{1}{2}\Vert I_0 \circ {\bf u} - T_0\Vert^2 + \frac{\alpha }{2}\Vert \nabla {\bf u}\Vert ^2 , \nonumber\\ &&{\bf u} = {\rm arg\, min}\;E( {\bf u} ) \end{eqnarray} \tag{ 1 }$$

Here, an optimal deformation vector field u is estimated by minimizing an energy function between a moving image I₀ and a target image T₀. α is a regularization parameter for controlling the smoothness of u over the image space.

In adaptive radiotherapy, I₀ often represents a planning CT (pCT) image and T₀ refers to a treatment CT (tCT) image. With physician delineated contours in pCT and edited contours in tCT, a constraint term defined by contours can be incorporated into equation (1) as

$$\begin{eqnarray} &&E( {\bf u} ) = \frac{1}{2}\Vert I_0 \circ {\bf u} - T_0\Vert^2 + \mathop \sum \limits_{i = 1}^k \frac{{\lambda _i }}{2}\Vert I_i \circ {\bf u} - T_i\Vert^2 + \frac{\alpha }{2}\Vert \nabla {\bf u}\Vert^2, \nonumber\\ &&{\bf u} = {\rm arg\, min}\;E( {\bf u} ) \end{eqnarray} \tag{ 2 }$$

Here, I_i and T_i denote modified images constructed by incorporating the ith contour pair set on pCT and tCT into original images I₀ and T₀, respectively; and k is the total number of modified contour pairs. The image modification is achieved with $I_i = \big( {1 + \frac{{M_i ( {\bf x} )}}{2}} \big)I_0$ and $T_i = \big( {1 + \frac{{N_i ( {\bf x} )}}{2}} \big)T_0$, where $M_i ( {\bf x} ) = \big\{ {\begin{array}{@{}cc@{}} \ssty{ \pm 1} & \ssty{{\bf x} \in O_i^p } [-4pt] \ssty 0 & \ssty{else} \end{array}}$ and $N_i ( {\bf x} ) = \big\{ {\begin{array}{@{}cc@{}} \ssty{ \pm 1} & \ssty{{\bf x} \in O_i^t } \ssty0 & \ssty{else} \end{array}}$. O^p_i (or O^t_i) represents a spatial domain that manifests a volume enclosed by the ith contour on pCT (or tCT). We refer O^p_i and O^t_i as regions of interest (ROIs) and I_i and T_i as modified images (MIs). M_i(x) and N_i(x) are positive when ROIs' intensity is higher than its surrounding and negative when lower. For simplification, λ_i is a regularization parameter for the ith contour pair set, which controls the influence of ROIs registration on the entire images registration.

Similar to the alternative optimization strategy proposed in Vercauteren et al (2009), an auxiliary variable c is introduced to equation (2) to enable a relaxation minimization

$$\begin{eqnarray} &&E( {{\bf c},{\bf u}} ) = \frac{1}{2}\Vert I_0 \circ {\bf c} - T_0\Vert^2 + \mathop \sum \limits_{i = 1}^k \frac{{\lambda _i }}{2}\Vert I_i \circ {\bf c} - T_i\Vert^2 + \frac{\beta }{2}\Vert {\bf c} - {\bf u}\Vert^2 + \frac{\alpha }{2}\Vert \nabla {\bf u}\Vert^2 , \nonumber\\ &&({\bf c},{\bf u}) = {\rm arg\, min}\;E( {{\bf c},{\bf u}} ) \end{eqnarray} \tag{ 3 }$$

The algorithm is minimizing candu in alternative steps. The first step starts with a given u and optimizes ${\rm arg\, min}_{\bf c}\; E_1 \left( {\bf c} \right) = \frac{1}{2}\Vert I_0 \circ {\bf c} - T_0\Vert^2 + \mathop \sum \nolimits_{i = 1}^k \frac{{\lambda _i }}{2}\Vert I_i \circ {\bf c} - T_i\Vert^2 + \frac{\beta }{2}\Vert {\bf c} - {\bf u}\Vert^2$, where optimization make steps from c = u. The second step solves the regularization by optimizing ${\rm arg\, min}_{\bf u}\; E_2 ( {\bf u} ) = \frac{\alpha }{2}\Vert\nabla {\bf u}\Vert^2 + \frac{\beta }{2}\Vert{\bf c} - {\bf u}\Vert^2$ with given c. In the second step, variety of regularizations have been proposed (Cachier and Ayache 2004). Here, we choose an optimal regularization reached by the convolution of the deformation field u with a Gaussian smoothing kernel G(α) (Vercauteren et al 2009).

In the first step, given the current deformation u, an updated (or displacement) deformation vector du can be calculated by minimizing

$$\begin{eqnarray} &&\fl {\rm arg}\mathop {\min }\limits_{{\rm d}{\bf u}} E_1 ( {{\rm d}{\bf u}} ) = \frac{1}{2}\Vert I_0 \circ {\bf u} \circ ( {{\bf Id} + {\rm d}{\bf u}} ) - T_0\Vert^2 \nonumber\\ &&+\, \mathop \sum \limits_{i = 1}^k \frac{{\lambda _i }}{2}\Vert I_i \circ {\bf u} \circ ( {{\bf Id} + {\rm d}{\bf u}} ) - T_i\Vert^2 + \frac{\beta }{2}\Vert {\rm d}{\bf u}\Vert^2 , \end{eqnarray} \tag{ 4 }$$

where I₀○u and I_i○u are deformed images with a given deformation u. For simplification, I₀○u is denoted as $\widetilde{I_0 }$ and I_i○u as $\widetilde{I_i }$. With a small displacement du, we can approximate the moving image term with a Taylor expansion: $I_0 \circ {\bf u} \circ ( {{\bf Id} + {\rm d}{\bf u}} ) = \widetilde{I_0 } \circ ( {{\bf Id} + {\rm d}{\bf u}} ) \approx \widetilde{I_0 } + \nabla \widetilde{I_0 }\cdot {\rm d}{\bf u}$ and $I_i \circ {\bf u} \circ ( {{\bf Id} + {\rm d}{\bf u}} ) = \widetilde{I_i } \circ ( {{\bf Id} + {\rm d}{\bf u}} ) \approx \widetilde{I_i } + \nabla \widetilde{I_i }\cdot {\rm d}{\bf u}.$ Then, equation (4) becomes:

$$\begin{equation} \fl {\rm arg}\mathop {\min }\limits_{{\rm d}{\bf u}} E_1 ( {{\rm d}{\bf u}} ) \approx \frac{1}{2}\Vert\widetilde{I_0 } + \nabla \widetilde{I_0 }\cdot {\rm d}{\bf u} - T_0\Vert^{2} + \mathop \sum \limits_{i = 1}^k \frac{{\lambda _i }}{2}\Vert\widetilde{I_i } + \nabla \widetilde{I_i }\cdot \,{\rm d}{\bf u} - T_i\Vert^{2} + \frac{\beta }{2}\Vert{\rm d}{\bf u}\Vert^2 \end{equation} \tag{ 5 }$$

In equation (5), we rewrite du = vn, where, v and n are the amplitude and the moving direction of the displacement vector du, respectively. Assuming that the displacement is determined on the projection du onto ${\bf n} = \frac{{\nabla \widetilde{I_0 }}}{{| {\nabla \widetilde{I_0 }} |}}$ (Cachier et al 1999), the optimization problem defined in equation (5) can be simplified from du = arg min E(du) to v = arg min E(v). Now the equation (5) becomes

$$\begin{equation} \fl {\rm arg}\mathop {\min }\limits_v E_1 ( v ) \approx \frac{1}{2}\Vert( {\nabla \widetilde{I_0 }\cdot {\bf n}} )v - ( {T_0 - \widetilde{I_0 }} )\Vert^{2} + \mathop \sum \limits_{i = 1}^k \frac{{\lambda _i }}{2}\Vert( {\nabla \widetilde{I_i }\cdot {\bf n}} )v - ( {T_i - \widetilde{I_i }} )\Vert^{2} + \frac{{\beta v^2 }}{2}\end{equation} \tag{ 6 }$$

Applying the first-order optimality condition $\frac{{\partial E_1 ( v )}}{{\partial v}} = 0,$we obtain a solution of equation (6):

$$\begin{equation} v = \frac{{( {\nabla \widetilde{I_0 }\cdot {\bf n}} )( {T_0 - \widetilde{I_0 }} ) + \mathop \sum \nolimits_{i = 1}^k \lambda _i ( {\nabla \widetilde{I_i }\cdot {\bf n}} )( {T_i - \widetilde{I_i }} )}}{{( {\nabla \widetilde{I_0 }\cdot {\bf n}} )^2 + \mathop \sum \nolimits_{i = 1}^k \lambda _i ( {\nabla \widetilde{I_i }\cdot {\bf n}} )^2 + \beta }};\end{equation} \tag{ 7 }$$

where β can be defined a spatial variable such as $\beta = ( {T_0 - \widetilde{I_0 }} )^2$ (Thirion 1998).

As described in section 2.1, CG-DIR is derived using a step of Taylor expansion and truncation, which implies that we can only allow a small and local displacement in registration. In order to reduce the magnitude of the displacement with respect to discretized image voxel size, a multi-scale strategy is adopted in CG-DIR algorithm. Registration iterations start with the lowest resolution, then the moving vectors obtained at a coarser level are up-sampled to serve as an initial solution at the next finer level. In this study, the entire registration is accomplished with two scales. The original images are only down-sampled once by a factor of 2 in each dimension. It was found that further down-sampling does not help improving either registration accuracy or efficiency. At coarse resolution levels, demons registration is adopted to generate an approximate DVF. At the finest level, the moving vectors are updated by the proposed CG-DIR algorithm. One major reason of adopting demons at coarse resolution level is to reduce computational burden. Also, the moving vectors obtained with demons at coarse resolution level are accurate enough to serve as an initial guess for the fine level CG-DIR registration.

Since the ROIs delineated by physicians are considered as the focus of our registration, the stopping criterion in CG-DIR is designed to place emphasis on the registration accuracy of these ROIs. Each ROI is expanded with 1 cm margin. The union of the expanded ROIs, referred to as UROI in this paper, is the designed region to estimate the stopping metric. We defined a root mean square error of intensity difference (RMS-ID) at the kth iteration in UROI as: $\delta ^k = \sqrt {1/N\mathop \sum \nolimits ( {I^k ({\bf x}) - T({\bf x})} )^2 ,} \;{\rm with}\;{\bf x} \in {\rm UROI}$. Then, a stopping metric is chosen as the change of RMS-ID in successive iterations: l^(k) = δ^k − δ^{k − 1}, and use l^{(k − 10)} − l^(k) ⩽ ε, where ε = 1.0 × 10⁻⁴, as our stopping criterion.

In the CG-DIR algorithm, there are three regularization parameters involved: α, β, and λ_i. These three regularization parameters have different functionality in the proposed CG-DIR algorithm. α, corresponding to the standard deviation of a Gaussian filter in our implementation, ensures a spatial smoothness of the deformation vector u. In our study, we empirically choose a unit standard deviation for the Gaussian filter. β, related to image noise, is estimated locally as: $\beta = ( {T_0 - \widetilde{I_0 }} )^2$ (Thirion 1998, Vercauteren et al 2009). λ_i is a parameter designed to balance the registration efforts between the modified image pair and the original image pair. The value of λ_i depends on the intensity ratio between modified image ROIs and original image ROIs. In this work, the value of λ_i is found by a trial-and-error strategy and λ_i ∼ 1.0 × 10⁻³ gives good deformation results for all cases studied.

The proposed CG-DIR algorithm was derived based on the original GPU-based demons algorithm (Gu et al 2010) and implemented on the NVIDIA CUDA enabled GPU platform. We adopt the following data parallel GPU kernels including: kernel 1. a low-pass filter kernel to smooth images and deformation vectors; kernel 2. a gradient kernel to calculate the gradient of images; kernel 3. a moving vector kernel to calculate displacement vectors and update deformation vectors; kernel 4. an interpolation kernel to deform images with moving vectors and update moving vectors; kernel 5. a comparison kernel to calculate the stopping criteria to terminate iterations properly. Here, the low-pass filter kernel is accomplished by convolving the image with a three-dimensional Gaussian filter of a unit standard deviation. Taking advantage of the separability property of the Gaussian kernel, the convolution is implemented in x,y, and z directions sequentially.

The CG-DIR algorithm is summarized as follows:

Algorithm A:

Initialization:

Set initial DVF u₀: =0;

Down-sample images (I₀, T₀) to the coarsest resolution;

A. Coarse-level registration:

Perform demons registration between down-sampled I₀ and T₀ images;

Up-sample moving vector u to a finer resolution level;

B. Finest-level registration:

Create modified images I_i and T_i;

 loop {over n until convergence}

 1. Given the DVF at the n-th iteration uⁿ, compute the displacement vector: ${\rm d}u^{n + 1} = \qquad\, v^{n + 1} {\bf n}^{n + 1}$, where v^{n + 1} is from equation (7) and ${\bf n}^{n + 1} = \frac{{\nabla I_0^n }}{{\left| {\nabla I_0^n } \right|}}$ (kernels 2 and 3);

 2. Smooth the displacement vector du^{n + 1} ← du^{n + 1}⊗G, here G is a Gaussian smoothing

operator (kernel 1);

 3. Update DVFs with a composite addition (Vercauteren et al 2009) ${\bf u}^{{\bf n} + 1} = {\bf u}^{\bf n} \circ \qquad\,( {{\bf Id} + {\rm d}{\bf u}^{{\bf n} + 1} } )$ (kernel 4), here I is an identity matrix;

 4. Smooth the updated DVF u^{n + 1} = u^{n + 1}⊗G (kernel 1);

 5. Warp images I₀○u^{n + 1} and I_i○u^{n + 1} (kernel 4);

 6. Calculate l⁽ⁿ⁾ to assess the stopping metric l^{(n − 10)} − l⁽ⁿ⁾ ⩽ ε (kernel 5).

end loop

The performance of CG-DIR is evaluated with clinical pCT and tCT image pairs of five head-and-neck cancer patients and one pelvis cancer patient. The tCT images were acquired 3–4 weeks after the first fraction of treatment due to the significant changes in anatomy. In this study, both pCT and tCT images are segmented to exclude couch, masks, and other accessories. All DIR analysis is conducted on the segmented images. The original pCT and tCT images have dimensions around ∼512 × 512 × 150. Due to the limitation of GPU memory, both tCT and pCT are down-sampled in axial planes and cropped in three dimensions. The images in this study have a resolution of 1.97 × 1.97 × 2.5 mm³. Representative organs with low contrast image intensity to background such as GTV, parotids, rectum, and bladder were selected for testing and evaluating the CG-DIR.

Direct evaluation of DVF accuracy is difficult in clinical cases due to the lack of ground-truth DVFs. However, other information, such as image intensity, structure volume, and contours, is often available. Hence, indirect evaluation approaches through the usage of available ground-truth information is commonly adopted. More specifically, the calculated DVF is used to deform the image and ROI volumes, or to propagate the contours. Then, the accuracy of DVF can be estimated in term of difference between the deformed and the target image intensities, overlapping rate of the deformed and the target ROI volumes, or distance between the propagated contours and the target contours.

In this study, the performance of the CG-DIR algorithm is first evaluated using image intensity difference. Two image intensity based assessment metrics are adopted: root mean square of intensity difference (RMS-ID) and correlation coefficient (CC) between the deformed pCT and the tCT. In the study, we calculate RMS-ID and CC over the entire image as well as the confined UROI.

The accuracy of the resulting DVF from the CG-DIR algorithm is further assessed using dice similarity coefficient (DSC). DSC assessment of volume overlapping, is a widely used similarity metric for evaluating deformable image registration results (Zhang et al 2007, Wang et al 2008, Xie et al 2008). DSC is defined as: ${\rm DSC} = \frac{{2\left( {V_D\mathop \cap \nolimits^ V_T } \right)}}{{\left( {V_D + V_T } \right)}}$, where V_D is the deformed volume of an ROI and V_T is the target volume of the same ROI. In this study V_T is the ROI volume enclosed by clinician manually edited contours. DCS measures the volume overlap of V_D and V_T, zero indicating no overlapping and one indicating an exact volume match. We apply this metric to solid organs like GTV and parotids.

The third type of comparison is contour matching, where the evaluated contours were propagated from pCT to tCT using the calculated DVF. Since DVFs generated by the CG-DIR algorithm only establish a 3D voxel-to-voxel correspondence map between pCT and tCT, when 2D axial-plane defined contour points on pCT are mapped onto rCT, points can easily cross the plane. Such a point-to-point mapping will destroy points order in contour definition and generate non-smooth contours. Extra efforts are needed to reorganize points and smooth contours. To avoid those extra tasks, we adopt a volume deformation strategy (Wang et al 2008): (1) deforming ROI volumes; (2) using morphological 'opening' operation to remove small islands; (3) extracting boundary to obtain the propagated contours.

Quantitatively, the contours comparison is assessed with surface-to-surface distance (SSD), which is defined as the shortest Euclidean distance of the deformed points on contours to the target surface. In this study, the deformed contours of bladder and rectum are evaluated with this metric. The main reason is that the accurate deformation of organ walls rather than the fillers for these two organs are more critical in therapy treatment. In particular, we choose several basic statistical quantities of SSD to quantify the bladder and rectum deformation in a pelvic cancer case: (1) mean value of SSD (M-SSD) measuring how close the deformed surface matched the target surface; (2) standard deviation of SSD values (SD-SSD) quantifying how deformed points scattered around M-SSD; and (3) 95 percentile SSD value gauging peak SSD values.