Evaluation of Demons- and FEM-Based Registration Algorithms for Lung Cancer

Abstract

We evaluated and compared the accuracy of 2 deformable image registration algorithms in 4-dimensional computed tomography images for patients with lung cancer. Ten patients with non-small cell lung cancer or small cell lung cancer were enrolled in this institutional review board-approved study. The displacement vector fields relative to a specific reference image were calculated by using the diffeomorphic demons (DD) algorithm and the finite element method (FEM)-based algorithm. The registration accuracy was evaluated by using normalized mutual information (NMI), the sum of squared intensity difference (SSD), modified Hausdorff distance (d_{H_M}), and ratio of gross tumor volume (r_GTV) difference between reference image and deformed phase image. We also compared the registration speed of the 2 algorithms. Of all patients, the FEM-based algorithm showed stronger ability in aligning 2 images than the DD algorithm. The means (±standard deviation) of NMI were 0.86 (±0.05) and 0.90 (±0.05) using the DD algorithm and the FEM-based algorithm, respectively. The means of SSD were 0.006 (±0.003) and 0.003 (±0.002) using the DD algorithm and the FEM-based algorithm, respectively. The means of d_{H_M} were 0.04 (±0.02) and 0.03 (±0.03) using the DD algorithm and the FEM-based algorithm, respectively. The means of r_GTV were 3.9% (±1.01%) and 2.9% (±1.1%) using the DD algorithm and the FEM-based algorithm, respectively. However, the FEM-based algorithm costs a longer time than the DD algorithm, with the average running time of 31.4 minutes compared to 21.9 minutes for all patients. The preliminary results showed that the FEM-based algorithm was more accurate than the DD algorithm while compromised with the registration speed.

Introduction

Four-dimensional computed tomography (4D CT) has been widely used in radiation therapy to acquire time-resolved volumetric images ^{1 –3} and access patient-specific breathing motion for determining individual safety margins. ^{4 –8} Since the thoracic organs expand and contract nonuniformly with the breathing, deformable image registration (DIR) is required to extract the motion and anatomical information from 4D CT images.

Several different algorithms have been proposed and validated for nonrigid registration of 4D CT images. ^{9 –17} For example, Wang et al ⁹ proposed an accelerated “demons” algorithm and validated it on patients with prostate, head-and-neck, and lung cancer. An “active” force along with an adaptive force strength adjustment during the iterative process was introduced to accelerate the registration algorithm. The results showed that the accelerated demons algorithm had a significant potential for delineating and tracking doses in target and critical structures during computed tomography (CT)-guided radiotherapy. Wu et al ¹⁰ developed a registration method that could be used to create discontinuous vector fields at the boundaries of anatomic subregions and overcame the difficulties in most registration algorithms, which assumed that the deformation field was smooth and continuous everywhere. The authors indicated that the lung was allowed to slide almost independently along the chest wall, resulting in a discontinuous deformation vector field, which could not be generated using standard deformable registration methods. In their study, registration was performed independently on each subregion, with a boundary-matching penalty used to prevent gaps. Both B-spline and demons registration algorithms were used to implement and test this method. The approach was validated on 4D CT data sets of 4 patients with lung cancer for registration of the end-inhalation and end-exhalation volumes. The results showed that the registration achieved consistently improved accuracy for both the B-spline and the demons algorithms in both moving and less-moving subregions. Retzel et al ¹¹ registered 4D CT data to a reference respiratory phase using a deformable registration based on B-splines. Anatomic landmarks were selected across 10 respiratory phases in 5 patients to assess registration performance. The average registration error for 5 landmarks for each of the patients was 2.1 mm. Velec et al ¹² used a biomechanical model-based DIR algorithm, previously developed by Brock et al ¹³ to retrospectively deform each exhale CT to inhale CT. The deformation map was combined with exhale and inhale dose grids from the treatment planning system to accumulate dose over the breathing cycle. The diffeomorphic demons (DD) algorithm proposed by Vercauteren et al ¹⁴ has been recognized as a robust and reliable method. ¹⁵ This algorithm is based on the classic demons, ¹⁶ which is a well-established technique for nonrigid registration and computes the pixel displacements driven by the image intensity difference. Instead of performing the optimization on the displacement vector field (DVF) as in the classic demons algorithm, the DD algorithm combines the Lie group framework on diffeomorphisms and optimizes for the Lie group. Thus, this method ensures to generate invertible output transformation. Furthermore, most nonrigid registration methods are established by solving an energy function minimization problem or a partial differential equation. The finite element method (FEM) ¹⁷ is a powerful technique to solve these problems, which commonly divides the problem domain into elements and further solves the equations element-by-element. The advantage of FEM-based method is the ability of dividing the domain into nonuniform elements. Thus, it can reduce the computational time by applying small elements where accurate registration is required and large elements in areas of less importance.

The goal of our study is to evaluate and compare the accuracy of the DD- and FEM-based algorithms. This is a novel study because no investigation has been presented to report the comparison between the 2 algorithms. Ten 4D CT data sets of patients with lung cancer were included in this study. The CT image at phase T = 0% (CT₀₀) was chosen as the reference image, while the CT images at other respiratory phases were treated as the floating images successively. To speedup the registration process, maintain robustness, and avoid local extremum, the reference and floating images were downsampled into 3 resolution levels, forming 2 multiresolution image pyramids, and the multiresolution registration was performed iteratively from the coarsest level to the finest level using the DD- and FEM-based registration algorithms. The DVFs relative to the reference image calculated by the 2 DIR algorithms were extracted for comparison. To quantitatively evaluate the 2 algorithms, 3 most commonly used similarity measures, normalized mutual information (NMI), the sum of squared-intensity differences (SSD), and modified Hausdorff distance (d_{H_M}), were used in this study. Furthermore, the computation time of the 2 algorithms was recorded and compared. The registration schemes used in this study were evaluated and compared quantitatively and qualitatively.

Materials and Methods

Patient Cohort and Imaging Study

Ten patients (5 male, 5 female, age 34~73) who had nonsmall cell lung cancer (NSCLC; n = 5) or small cell lung cancer (SCLC; n = 5) were enrolled in this prospective study, which has been approved by the Institutional Review Board of Shandong Cancer Hospital and Institute. Table 1 summarizes the clinical characteristics of the patients. All the patients underwent CT scans on the same day for treatment planning and were positioned head-first-supine with the arms raised above the head. Vacuum bags were used to immobilize the patients during image acquisition.

Table 1.

Summary of Patients’ Characteristics and Measurements.

Patient	Age	Gender	Cancer Site	Tumor Location	Staging
1	70	M	NSCLC	R-L Lobe	T1
2	68	F	NSCLC	L-L Lobe	T2
3	67	F	SCLC	L-U Lobe	T1
4	34	F	NSCLC	R-L Lobe	T3
5	68	F	SCLC	R-U Lobe	T1
6	67	M	NSCLC	L-L Lobe	T2
7	55	M	SCLC	L-U Lobe	T2
8	54	M	NSCLC	R-L Lobe	T3
9	60	F	SCLC	R-L Lobe	T2
10	73	M	NSCLC	L-L Lobe	T1

Abbreviations: M, male; F, female; NSCLC, non-small cell lung cancer; SCLC, small cell lung cancer; R-L, right-lower; L-L, left-lower; L-U, left-upper; R-U, right-upper.

For each patient, a 4D CT scan was performed under uncoached free breathing condition on a 16-slice CT scanner (Philips Brilliance Bores CT, Andover, Massachusetts) equipped with real-time position management (RPM) system (Varian Medical Systems, Inc, Palo Alto, California) and Advantage 4D software (GE Healthcare, Milwaukee, Wisconsin). The respiratory signal was recorded with the RPM gating system by tracking the trajectory of infrared markers placed on the patient’s abdomen. Each CT image from the scanner was labeled by the time tag according to the respiratory signal. The reconstructed 4D CT images were sorted into 10 respiratory phases based on tags by the Advantage 4D software, with 0% corresponding to end inhalation and 50% corresponding to end -exhalation. The imaging parameters were as follows: voltage/current: 120 kV/290 mA, slice thickness: 2.5 mm, gantry rotation: 0.5 s/cycle, reconstruction matrix: 512 × 512, and field of view: 450 to 500 mm.

Diffeomorphic Demons Algorithm

The complete registration framework can be treated as an optimization problem that aims at finding the deformation field to obtain a reasonable alignment between images. And it can also be summarized by a model with an energy function consisting of a similarity function and a regularized term. ¹⁸

For a fixed image (F), a moving image (M), and a transformation field (s), the energy function of the DD algorithm with respect to the update field (u) can be calculated using Equation 1:

$E\left(\mathit{u}\right)=F-M\circ {\left(s+\mathit{u}\right)}^2+\frac{{\mathit{\sigma }}_i^2}{{\mathit{\sigma }}_x^2}{\mathit{u}}^2,$

1

where ${\mathit{\sigma }}_i^2$ and ${\mathit{\sigma }}_x^2$ are image noise and the corresponding spatial uncertainty, respectively. The operator ○ denotes composition. The first term on the right side represents the similarity metric, which is the sum of absolute differences (SAD) between the fixed (F) and moving (M) images. The right-most term is used for regularization. As aforementioned, the registration work could be summarized by a model with an energy function consisting of a similarity function and a regularized term. The idea of this method is to consider the regularization criterion as a prior on the smoothness of the transformation (s). Instead of requiring the point correspondences between image pixels be exact realizations of the transformation, one allows some errors at image point.

The energy function calculated in Equation 1 can be linearized using the first-order Taylor expansion as follows:

$E\left(\mathit{u}\right)=F-M\circ s+J{\mathit{u}}^2+\frac{{\mathit{\sigma }}_i^2}{{\mathit{\sigma }}_x^2}{\mathit{u}}^2,$

2

$J=-\left(\mathrm{\nabla }F+\mathrm{\nabla }\left(M\circ s\right)\right)/2,$

3

where J is used for performing an efficient second-order minimization symmetric registration, with $\mathrm{\nabla }$ accounting for the gradient of the images. We can obtain the update field with the minimization of energy function, satisfying $\mathrm{\nabla }E\left(\mathit{u}\right)=0$ .

$\mathit{u}=-\frac{F-M\circ s}{J^2+\frac{{\mathit{\sigma }}_i^2}{{\mathit{\sigma }}_x^2}}$

4

Typically, the Gaussian smoothing is used as the regularization term for the displacement field ( u ). To generate an invertible output transform, a diffeomorphic extension was proposed by Vercauteren et al. ¹⁴ The extension adapted the optimization procedure to a space of diffeomorphic transformation. And the update of the current transformation s was done through the exponential map on the Lie group with an update step u:

$s=s\circ \exp \left(\mathit{u}\right).$

5

So the iterative registration algorithm can be organized as follows:

For the current transformation s, compute the corresponding field u by minimizing the energy function E( u ) with respect to u .

Typically, use the Gaussian kernel to regularize the update field u , satisfying u ← k_G × u .

Let $c=s\circ \exp \left(\mathit{u}\right)$ .

Typically, use the Gaussian kernel to regularize c, satisfying $s\leftarrow k_G\ast c$ .

In our study, 3 levels of resolution consisting of 20 iterations each were adopted. A Gaussian kernel with σ = 3 and σ = 2 were used for the update field regularization ( u ) and deformation field regularization (c), respectively. Those parameters were chosen through trials by considering balancing of the registration accuracy and the convergence speed of the DD algorithm.

FEM-Based Registration Algorithm

Given 2 images, a fixed image $F\left(\mathbf{x}\right),\mathbf{x}\in {\mathrm{\Omega }}_F$ and a moving image $M\left(\mathbf{x}\right)$ , $\mathit{x}\in {\mathrm{\Omega }}_M$ where ${\mathrm{\Omega }}_F$ and ${\mathrm{\Omega }}_M$ account for the fixed image domain and the moving image domain, respectively. The purpose of image registration is to find a deformation Φ(x), which best aligns the moving image with the fixed image according to a specific similarity measure. ¹⁹ And the energy functional used in FEM-based algorithm is as follows:

$\begin{array}{l}\mathrm{\Pi }={\int }_{{\mathrm{\Omega }}_F}(F(\mathit{x})-M\circ \mathrm{\Phi }(\mathit{x}){)}^{2}d\mathit{x}+\mathit{\alpha }\sum _i\frac{(\mathrm{\Phi }({\mathit{p}}_i^F)-{\mathit{p}}_i^M{)}^2}{2{\mathit{\varsigma }}_i^2}\\ +\frac{1}{2}{\int }_{{\mathrm{\Omega }}_F}{\mathit{\sigma }}^T\mathit{ϵ}d\mathit{x}\end{array},$

6

$\mathrm{\Pi }={\mathrm{\Pi }}^1+{\mathrm{\Pi }}^2+{\mathrm{\Pi }}^3,$

7

where the first term represents the similarity of 2 images and the second term aligns the predefined landmarks in the fixed image ( ${\mathit{p}}_i^F$ ), with the corresponding landmarks in the moving image ( ${\mathit{p}}_i^M$ ). Parameter α defines the importance of the landmark constraint and parameter ς _i represents the weighting of each corresponding set of landmarks. The last term is a linear regularizer, where ∊ is a strain tensor and σ is a stress tensor. In this study, we used the SAD instead of Gee probabilistic similarity metric for ${\mathrm{\Pi }}^1$ , whereas the 2 terms Π² and Π³ were the similar to the ones used by Gee. Our goal is to find arg min_ΦΠ.

In this study, anatomical landmark feature pairs were identified manually across 10 respiratory phase images. A Matlab-based software named Assisted Point Registration of Internal Landmarks (APRIL, The MathWorks Inc., Natick, Massachusetts) was employed to facilitate the selection of landmarks. ^20,21 The corresponding feature points were determined by an expert via manually designing the feature correspondence on the target image. For each patient, a reference set of lung landmark feature points was generated on the end-inhale phase image (CT₀₀) from the 4D CT data set. The vessels and bronchial bifurcations were typically included in the feature points, which were selected systematically on the 10-phase images by the expert. About 10 feature points per axial slice were selected and propagated throughout the respiratory phases as shown Figure 1, which ensured the collection of >300 feature points for all cases. For each case, the feature points were visually reviewed by the observer a second time and location adjusted according to the end-inhale phase image if necessary. The points were then used to verify the accuracy of the FEM-based algorithm in the study.

Figure 1.

Representative 4D CT images in coronal and sagittal views of patient 1. Dash lines are drawn to facilitate the visualization of respiratory motion. Arrows indict the tumor. 4D CT indicates 4-dimensional computed tomography.

A 3-level hierarchy same as the DD algorithm was used for the registration procedure with a decreasing Young modulus schedule of (1.00, 0.10, 0.01) and a constant Poisson ratio of (0.5, 0.5, 0.5). The FEM parameters were chosen according to the previous literature ¹⁹ by considering balancing the registration accuracy and the convergence speed of the algorithm.

Quantitative Evaluation of Registration

During the registration process, the DVFs relative to the reference image (CT₀₀) were calculated by the 2 algorithms. To demonstrate the ability of these 2 algorithms visually, each DVF was delineated and compared.

Four commonly used metrics, NMI, SSD, d_{M_H}, and registration speed were employed to compare the efficiency of the 2 registration methods. The NMI can be calculated as follows:

$\mathrm{N}\mathrm{M}\mathrm{I}\left(I_1,I_2\right)=\frac{H\left(I_1\right)+H\left(I_2\right)}{H\left(I_1,I_2\right)}\mathrm{,}$

8

where H(I₁), H(I₂), and H(I₁,I₂) denote the information entropy of 2 images and the joint entropy, respectively. The greater NMI between 2 images, the better alignment they have.

The SSD between 2 images is one of the simplest voxel similarity measures, and it is minimized during registration. For voxel location X _A in image A, within an overlap domain Ω _AB ^T containing N voxels:

$\mathrm{SSD}=\frac{1}{N}\sum _{X_A\in {\mathrm{\Omega }}_{A,B}^T}{\left|A\left({\mathit{X}}_A\right)-B^T\left({\mathit{X}}_B\right)\right|}^2\mathrm{.}$

9

The Hausdorff distance which can measure the distance of two images from each other was commonly used as the similarity metric to evaluate the registration results. Given 2 images I and J, it is defined as:

$d_H\left(I,J\right)=\max \left(h\left(I,J\right),h\left(J,I\right)\right),$

10

where h(I, J) and h(J, I) are the Hausdorff distance of I → J and J →I, respectively, and defined as:

$h\left(I,J\right)=\underset{i\in I}{max}h\left(i,J\right)=\underset{i_n\in I}{max}\underset{j_n\in J}{min}{i}_n-j_n,$

11

$h\left(J,I\right)=\underset{j\in J}{max}h\left(j,I\right)=\underset{j_n\in J}{max}\underset{i_n\in I}{min}{j}_n-i_n\mathrm{,}$

12

where i and j represent pixels in image I and image J, respectively, and $∥\cdot ∥$ represents the calculation of distance between 2 points calculated, using the Euclidean distance as metric in this study.

Since the traditional Hausdorff distance is sensitive to image noise, in this study, we used the modified Hausdorff distance proposed by Dubuisson et al, ²² which can avoid the deviation caused by the interference of some noisy pixels. The modified Hausdorff distance (d_{M_H}) was defined as:

$d_{H\mathrm{\_}M}\left(I,J\right)=\frac{1}{N_I}\sum _{i\in I}h\left(i,J\right)\mathrm{,}$

13

where N_I is the number of feature points in image set I.

The r_GTV between reference image and deformed phase images using 2 methods was also measured and compared. The GTVs on the reference image (GTV _r ) and deformed phase images (GTV _f ) were delineated by an oncologist. The r_GTV value was calculated as follows:

$r_{GTV}=\frac{\left|GT{V}_r-GT{V}_f\right|}{GT{V}_r}(\%).$

14

In our study, images at phase ${\mathrm{C}\mathrm{T}}_{00}$ were served as the reference, and the rest images at the other phases were treated as the phase images, respectively. To directly illustrate the efficacy of the registration methods, we calculated all the average measurements between the reference image (CT₀₀), and the phase images calculated by NMI, SSD, d_{H_M}, r_GTV and registration speed for all the patients.

Results

Figure 2 shows a representative 4D CT data set of patient 1. From the tumor locations in the lung at 10 phases, we can observe the respiratory motion over the breathing cycle.

Figure 2.

Examples of the DVFs relative to a specific reference phase (CT₀₀) calculated by the DD algorithm (a) and the FEM-based algorithm (b). DVF indicates displacement vector field; DD, diffeomorphic demons; FEM, finite element method.

For each patient, 9 phase images were registered to the reference (CT₀₀), and 9 sets of DVF recording the movement of each voxel in phase images to corresponding voxel in the reference were obtained. Figure 3 demonstrates the DVFs in the lung area obtained from the DD algorithm (Figure 3A) and the FEM-based algorithm (Figure 3B) in patient 1. It was visible that the corresponding displacement fields showed differences between the 2 methods, with greater displacements appearing in the FEM-based algorithm.

Figure 3.

Examples of the displacement trajectories (right) of a voxel in the lung tumor (indicted by arrows) spanning the phases from maximum inhalation (CT₀₀) to maximum exhalation (CT₅₀) determined by the 2 DIR algorithms compared with a landmark trajectory treated as a golden standard. The tracked feature point was indicted by the red points. DIR indicates deformable image registration (Color version of the figure available online).

Figure 1 shows an example of the displacement trajectory (right) of a pixel identified by a radiation oncologist in the tumor spanning the phases from maximum inhalation (CT₀₀)to maximum exhalation (CT₅₀) using the 2 DIR methods. The tracked points were indicted by red points. We can observe that the displacement trajectories show large differences.

Figure 4 shows an example of sampled DVFs between CT₅₀ and CT₀₀ confined in the lung area in coronal plane. As aforementioned, CT₀₀ (a) served as the reference image to obtain the DVFs calculated by the FEM-based (b) and DD (c) algorithms.

Figure 4.

The coronal DVFs confined in the lung area calculated by the FEM-based (b) and DD (c) algorithms for patient 1. The CT image at phase T = 0% (a) was used as the reference image to derive the DVFs. DVF indicates displacement vector field; FEM, finite element method; DD, diffeomorphic demons, CT, computed tomography.

To validate the registration efficacy, we calculated the subtractions of phase images before/after registration and reference image. Figure 5 gives an example of subtraction images in 3 orthogonal directions (superior–inferior, anterior–posterior, and medial–lateral) between ${\mathrm{C}\mathrm{T}}_{00}$ and CT₅₀ of patient 1. Figure 5A illustrates the fusion of phase CT₀₀ and phase CT₅₀, with the reference image in gray and the phase image in hot metal. From the coronal and sagittal views, we can find that the tumor and diaphragm locations in the 2 images show large differences mainly caused by the respiratory motion. Figure 5(B-D) shows the subtractions of reference and phase images before and after registration using the DD algorithm and the FEM-based algorithm, respectively. The results implied that the subtraction of reference and original phase images was large, especially in the tumor area (indicted by circles) and lower lobe in the lung. And this subtraction was greatly reduced after the registration using the DD algorithm, but the differences in the lower lobe of lung and tumor locations could not be ignored. The phase image and the reference image showed good consistency after using the FEM-based registration algorithm as shown in Figure 5(D), excluding some specific regions. It was hard to see any difference (indicted by circles) in tumor location between the 2 phases. So the preliminary results showed that the FEM-based algorithm showed greater ability in aligning 2-phase images of lung 4D CT than the DD algorithm.

Figure 5.

Representative subtractions of reference image and phase image before (b) and after using the DD algorithm (c) and the FEM-based (d) algorithm. The first line (a) shows the fusion images of reference image ( ${\mathrm{C}\mathrm{T}}_{00}$ ) and original phase image (CT₅₀). Tumors in the reference image and the phase image were contoured. DD indicates diffeomorphic demons; FEM, finite element method.

Table 2 summarizes the measurement results of alignment between CT₀₀ and CT₅₀ from all the patients. Of all patients, the FEM-based registration method showed more accuracy than the demons-based method. The means (±standard deviation [SD]) of NMI were 0.70 (±0.10), 0.82 (±0.06), and 0.85 (±0.07) before and after registration using the DD and FEM-based algorithms, respectively. The means of SSD were 0.016 (±0.010), 0.011 (±0.009), and 0.009 (±0.009) before and after registration using the DD and FEM-based algorithms, respectively. The means of d_H were 0.17 (±0.13), 0.06 (±0.03), and 0.04 (±0.03), before and after registration using the DD and FEM-based algorithms, respectively. The means of r_GTV were 4.5 (±1.2)% and 3.3 (±1.1)% after using the DD and FEM-based algorithms, respectively. However, the FEM-based algorithm costs a longer time on the registration compared to the DD algorithm (35.8 vs 25.3 minutes). As aforementioned, the FEM-based algorithm solved the problem on an element-by-element basis by partitioning the problem domain into elements controlled by nodes and assembles the contributions of each element into a global system of equations. So the computational complexity is large when using high-resolution FEMs.

Table 2.

Comparison of NMI, SSD, d_{H_M}, and Running Time on Registration Between $\mathrm{C}{\mathrm{T}}_{00}$ and $\mathrm{C}{\mathrm{T}}_{50}$ using the DD and FEM-Based Algorithms.

Patient	NMI			SSD			dH (pixel)			Running Time (min)		rGTV, %
	Before	After		Before	After		Before	After
		DD	FEM		DD	FEM		DD	FEM	DD	FEM	DD	FEM
1	0.57	0.79	0.82	0.012	0.007	0.005	0.09	0.04	0.03	25	40	3.4	2.0
2	0.76	0.88	0.92	0.008	0.006	0.004	0.12	0.04	0.02	22	35	4.2	3.5
3	0.89	0.92	0.95	0.014	0.010	0.008	0.19	0.06	0.04	20	30	3.2	2.4
4	0.74	0.85	0.88	0.003	0.003	0.001	0.08	0.06	0.03	21	32	6.8	5.0
5	0.67	0.76	0.80	0.027	0.008	0.006	0.04	0.03	0.01	28	39	3.5	2.3
6	0.72	0.86	0.90	0.010	0.007	0.006	0.47	0.05	0.03	26	36	5.6	3.8
7	0.68	0.79	0.81	0.003	0.002	0.001	0.04	0.02	0.01	30	38	5.5	4.0
8	0.54	0.70	0.72	0.016	0.013	0.010	0.27	0.11	0.08	32	42	5.0	4.2
9	0.72	0.82	0.80	0.028	0.020	0.018	0.20	0.10	0.09	24	35	4.0	3.7
10	0.70	0.85	0.86	0.034	0.032	0.030	0.15	0.08	0.06	25	31	3.5	2.0
Mean	0.70	0.82	0.85	0.016	0.011	0.009	0.17	0.06	0.04	25.3	35.8	4.5	3.3
SD	0.10	0.06	0.07	0.010	0.009	0.009	0.13	0.03	0.03	3.86	4.00	1.2	1.1

Abbreviations: NMI, normalized mutual information; SSD, sum of squared-intensity difference; d_{H_M,} modified Hausdorff distance; d_H, Hausdorff distance; DD, diffeomorphic demons; FEM, finite element method; r_GTV, ratio of gross tumor volume; SD, standard deviation.

Table 3 summarizes the averaged measurement results over 9 phase bins from all the patients. On average, the FEM-based registration method showed more accuracy than the demons-based method: The means of NMI were 0.75 (±0.08), 0.86 (±0.05), and 0.90 (±0.05) before and after registration using the DD and FEM-based algorithms, respectively. The means of SSD were 0.008 (±0.06), 0.006 (±0.003), and 0.003 (±0.002) before and after registration using the DD and FEM-based algorithms, respectively. The means of d_H were 0.91 (±2.5), 0.04 (±0.02), and 0.03 (±0.03) before and after registration using the DD and FEM-based algorithms, respectively. The means of r_GTV were 3.9 (±1.1)% and 2.9 (±1.1)% after using the DD and FEM-based algorithms, respectively. However, the FEM-based algorithm costs a longer time compared with the DD algorithm (32 vs 22 minutes).

Table 3.

Comparison of NMI, SSD, d_{H_M} and Running Time Measured by Averaging Over 9 Phases Using the DD and the FEM-Based Algorithms for All Patients.

Patient	NMI			SSD			dH (pixel)			Running Time (min)		rGTV (%)
	Before	After		Before	After		Before	After
		DD	FEM		DD	FEM		DD	FEM	DD	FEM	DD	FEM
1	0.65	0.82	0.85	0.008	0.006	0.003	0.05	0.03	0.02	23	35	3.0	1.8
2	0.80	0.90	0.96	0.006	0.005	0.003	008	0.03	0.02	20	30	3.5	3.0
3	0.91	0.92	0.96	0.010	0.008	0.006	0.12	0.04	0.03	18	28	2.6	2.0
4	0.78	0.86	0.94	0.002	0.001	0.001	0.06	0.04	0.02	18	30	5.2	4.6
5	0.72	0.88	0.92	0.012	0.009	0.005	0.02	0.01	0.01	25	35	3.0	2.0
6	0.76	0.89	0.93	0.008	0.006	0.002	0.35	0.03	0.02	20	32	5.0	3.5
7	0.74	0.84	0.86	0.002	0.001	0.001	0.02	0.01	0.01	27	32	5.2	3.5
8	0.60	0.76	0.80	0.012	0.010	0.006	0.20	0.08	0.02	25	36	4.5	4.0
9	0.76	0.83	0.85	0.014	0.006	0.002	0.15	0.07	0.06	20	30	3.6	3.2
10	0.74	0.86	0.92	0.020	0.010	0.005	0.10	0.06	0.04	23	26	3.0	1.5
Mean	0.75	0.86	0.90	0.008	0.006	0.003	0.91	0.04	0.03	21.9	31.4	3.9	2.9
SD	0.08	0.05	0.05	0.006	0.003	0.002	2.50	0.02	0.03	3.14	3.2	1.01	1.1

Abbreviations: NMI, normalized mutual information; SSD, sum of squared-intensity difference; d_{H_M,} modified Hausdorff distance; DD, diffeomorphic demons; FEM, finite element method; d_H, Hausdorff distance; r_GTV, ratio of gross tumor volume; SD, standard deviation.

Discussion and Conclusion

In this study, we evaluated the accuracy of the DD and FEM-based algorithms for 4D CT registration in 10 patients with lung cancer. The results demonstrated that the FEM-based algorithm showed greater ability in aligning 2 images than the DD algorithm. However, this algorithm requires a longer time mainly due to the ideal of nonuniform element partition. This work provided a fundamental research for DIR in radiation therapy.

This pilot study included a limited number of patients and assessed only patients with lung cancer. A larger pool of patients is needed in future studies to answer the following questions: (1) Is this comparison effective for cancers in other locations treated by radiation therapy, such as liver or esophageal cancers? (2) Are there any other methods to better evaluate the accuracy of DIRs? (3) How does the level of resolution affect the convergence of the cost function since we only employed three levels in registration procedure? There is no controversy about the fact that the deformation of lung is more easily affected by respiratory motion compared with liver or esophagus referred earlier. However, it could be estimated that the comparison results might be different if liver or esophagus was used to validate the performance of the 2 presented DIR algorithms due to the suboptimal contrast of liver 4D CT and the small deformation of esophagus 4D CT. In addition, poor soft-tissue contrast in liver 4D CT will definitely affect the selection of landmarks using the APRIL software.

Quantitative validation of a DIR algorithm is difficult due to the lack of a golden standard. Recently, there is no established method to validate the performance of registration algorithm. In our study, we qualitatively evaluated the DD algorithm and the FEM-based algorithm by comparing the magnitudes of motion vectors between the phase images and the reference image. In addition, the subtractions of the reference and deformed images calculated by the DD algorithm and the FEM-based algorithm were generated for directly visualized evaluation. Furthermore, 3 commonly used measurement metrics, NMI, SSD, and d_{H_M} were chosen to quantitatively validate the registration algorithms. All measurements were calculated before and after deformable registration between the reference image and the phase image, which could be used for determining the more accurate algorithm. A more robust method to evaluate the registration performance may be achieved by introducing a golden standard, such as anatomic landmarks, fiducial markers, or contour positions. ¹³ Furthermore, we can evaluate the efficiency of the registration according to different parts of the lung, since we observed that the lower lung had a greater motion compared with the upper lung.

The material properties for the lung were determined according to the previous literature. ^13,19 Brock et al ¹³ indicated that the optimal value of Poisson ratio (v) was difficult to determine for the lung, since trends varied between patients as well as between medial–lateral, anterior–posterior, and superior–inferior directions. The optimal value of v was assigned to 0.45 across all the patients in our study as indicated in Brock study. So the differences in optimal v between different directions and different patients could introduce potential errors. It is of great interest in our future work to investigate the potential benefit of different material properties to enhance fixed boundary conditions. Another drawback was the computational costs of the FEM-based registration method. In the current implementation, the in-plane 4D CT data set consisting of 512 × 512 pixels, with the initial element size of 16 × 16 requiring more than 1 hour using a standard PC (Inter(R) Core(TM) i7-3770 central processing unit @ 3.40 GHz and 8 G RAM). Therefore, multiresolution technique was used in this study to speedup the computation of the FEM-based algorithm. The computation time had been reduced to a half-time, which was comparable with the DD algorithm and meanwhile contained good accuracy. In addition, we found that, for both algorithms, the registration accuracy gradually improved and the convergence speed gradually slowed with the increase in the SD (σ) of Gaussian kernel. However, after σ was greater than 1, both the convergent speed and the registration accuracy were very similar with different $\mathit{\sigma }$ values. Furthermore, the FEM algorithm could remain balanced with good accuracy and high convergence speed with initial element size of 16 × 16 in this study.

The deformation and distortion of thoracic organs and tumors are mainly caused by respiratory motion. With the development of 4D CT techniques, we can acquire volumetric images combined with the breathing signals throughout the breathing cycle. In 3-dimensional (3D) treatment planning, the gross target volume is delineated manually by the oncologists in the planning system and then expanded to a clinical target volume (CTV) by adding a margin to cover expected microexpansion. Furthermore, the CTV is expanded to the planning target volume by adding an internal margin (IM) according to the clinical criterion by considering the target movement induced by the respiratory motion. The IMs are usually determined according to experiences or reported literatures. These margins are relative to cancer, category specific rather than patient specific. Thus, adding safety margins using this method is approximately unsafe. The advent of 4D CT technique has solved these problems and can be used for assessing patient-specific breathing motion and adding individual safety margins.

Furthermore, the DIR algorithms can be used to evaluate the dosimetric effect of respiratory motion on dose delivered to the tumor target and organs at risk (OARs). ^{23 –27} Calculation of 3D dose in radiotherapy usually utilizes 1 static phase (breath-hold 3D CT, end-exhale or end-inhale of 4D CT, average intensity projection, and maximum intensity projection). A major disadvantage is insufficient evaluation of the dose distribution to tumor target and OARs. Recently, we can perform 4D dose calculation by transferring the DVFs obtained from the DIR procedure to a 3D dose matrix calculated from 3D plans on the exhale or inhale CT of 4D CT images. Deformed phase dose on the reference image was calculated by deforming the 3D phase dose on each phase with the corresponding phase DVF, which can be calculated as follows:

$D_i^{x^{\prime },y^{\prime },z^{\prime }}=D_i^{x,y,z}\times DV{F}_i^{x,y,z},$

15

where the parameters $D_i^{x,y,z}$ and $D_i^{x^{\prime },y^{\prime },z^{\prime }}$ denote the 3D phase dose calculated on the i th phase image and the deformed dose matrix on the reference image in the ith phase. And the $DV{F}_i^{x,y,z}$ represents the DVF in the ith phase. The 4D dose is then calculated to accumulate deformed phase dose on the reference image with equal weighting (1/number of total phases), as follows:

$D_{4D}^{x^{\prime },y^{\prime },z^{\prime }}=\sum _i{w}_i\times D_i^{x^{\prime },y^{\prime },z^{\prime }},$

16

where $D_{4D}^{x^{\prime },y^{\prime },z^{\prime }}$ represents the accumulated 4D dose on the reference image, and the w_i. denotes the weighting factor of the ith phase. So we can perform 4D dose calculation on the basis of DIR.

So choosing a reliable and high-efficient DIR algorithm for the 4D images is crucial and significant for radiation therapy. It is of great interest to evaluate substantial registration algorithms and determine the most effective one for the 4D CT alignment and validate our point-to-point dose accumulation in our future study.