Accelerated gradient-based free form deformable registration for online adaptive radiotherapy

The intensity-based registration algorithms would cause registration error because of inconsistent intensities between the two modalities. We will illustrate this with the famous intensity-based registration algorithm, the Demons algorithm (Wang et al 2005, Rogelj and Kovacic 2006). The Demons algorithm generates displacement field by $\text{d}\vec{u}{{(\vec{x})}^{(k)}}$ through iteration:

$$\begin{eqnarray}&&\vec{u}{{(\vec{x})}^{(k+1)}}={{G}_{\sigma}}*(\vec{u}{{(\vec{x})}^{(k)}}+\text{d}\vec{u}{{(\vec{x})}^{(k)}})\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\text{d}\vec{u}{{(\vec{x})}^{(k)}}=-\left(\frac{\left[B\left(\vec{x}^{\prime} \right)-A(\vec{x})\right]\vec{\nabla}A(\vec{x})}{|\vec{\nabla}A(\vec{x}){{|}^{2}}+{{\gamma}^{2}}{{\left[B\left(\vec{x}^{\prime} \right)-A(\vec{x})\right]}^{2}}}+\frac{\left[B\left(\vec{x}^{\prime} \right)-A(\vec{x})\right]\vec{\nabla}B\left(\vec{x}^{\prime} \right)}{|\vec{\nabla}B\left(\vec{x}^{\prime} \right){{|}^{2}}+{{\gamma}^{2}}{{\left[B\left(\vec{x}^{\prime} \right)-A(\vec{x})\right]}^{2}}}\right)\end{eqnarray} \tag{ 2 }$$

where k is the current number of iterations, G_σ is the Gaussian filter, which is used to smooth the displacement field in order to maintain the geometric continuity of the deformation. $A\left(\vec{x}^{\prime} \right)$ is the static image, $B\left(\vec{x}^{\prime} \right)$ is the moving image at the k-th iteration and γ is used to control the displacement step length.

Figure 1 shows the FBCT to CBCT registration error due to the intensity inaccuracy for a lung example. The regions manually drawn by the red line in figure 1(b) represent some typical areas where CBCT intensity is inaccurate. It is just a sketch, rather than a precise measure. As can be seen from figure 1(c), the misalignment is evident before non-rigid registration, which is highlighted by the red arrows (a)–(e). After registration by the Demons algorithm, most of the regions are well matched and they are highlighted by the green arrows (a)–(c) in figure 1(d). But the Demons algorithm unrealistically distorts the tissues with severe scatter artifacts. In particular, before the Demons registration, the triceps borders of the FBCT image were already on the left of the CBCT triceps borders (see red arrow (e) in figure 1(c)). However, the triceps borders further dilate and seriously mismatch the triceps borders in the CBCT image after the Demons registration. This phenomenon can be easily observed by comparing the same region indicated by the red arrow (e) in figures 1(c) and (d).

Zoom In Zoom Out Reset image size

In figure 2, the deformation force acting on the borders of the triceps is shown to explain how the mismatch occurs. For simplicity, we take the green dashed line in figure 1(a) as an example. In this example, the triceps borders in FBCT were already on the left of the CBCT triceps borders (see region indicated by the red arrow (e) in figures 1(c) and 2). However, because the intensity of the triceps contours on CBCT is greater than that of FBCT, both B(x₁)–A(x₁) and $\vec{\nabla}A(\vec{x})$ are greater than 0, the passive force at x₁ is along −x direction (left direction). So the passive force drives the triceps borders to dilate. The active force acting on the CBCT triceps borders has a similar effect. So the triceps borders further dilated until the anti-directional deformation force produced by the low-intensity FBCT pixels dragged out of the triceps region stopped it. Although the intensity is inconsistent between the two modalities, the gradient is approximately consistent, particularly on the border of tissue or organs (Wiegert et al 2007, Zhen et al 2012). Because large gradient usually only exists on the border of tissue or organs, registration gradient field can effectively decrease the registration error due to CBCT intensity inaccuracy. In order to decrease the aforementioned effect by the intensity inaccuracy, a more accurate gradient-based deformable registration algorithm was developed.

Zoom In Zoom Out Reset image size

2.2.1. Gradient-based free form registration algorithm construction.

The registration problem can be mathematically modeled as minimizing the cost function $\varepsilon ({\overset \rightharpoonup {u}})$ (Troy and Kris 2013)

$$\begin{eqnarray}&&\varepsilon \left({\overset \rightharpoonup{u}}\right)={{\alpha}_{i}}{{D}_{i}}\left(A\left({\overset \rightharpoonup{x}}\right),B\left({\overset \rightharpoonup{x}}\right)*{\overset \rightharpoonup{u}}\left({\overset \rightharpoonup{x}}\right)\right)+{{\alpha}_{i+1}}{{R}_{i+1}}\left({\overset \rightharpoonup{x}}\right)+\ldots\end{eqnarray} \tag{ 3 }$$

where D_i denotes a data similarity term, R_{i + 1} denotes a regularization term, α_i is a weight on the ith term. There are several regularization terms, of which one of the most desired is smoothness as it directly affects the smoothness of the image to which it is applied. It is involved in our model. In this paper, the gradient-based deformable registration problem is preliminarily described as finding the displacement fields ${{\vec{u}}_{0}}(\vec{x})$ that minimize the energy functional $\varepsilon \left(\vec{u}\right)$:

$$\begin{eqnarray}&&\varepsilon \left(\vec{u}\right)=\underset{\vec{x}\in {{R}^{3}}}{\int} \left[\left\|\nabla B\left(\vec{x}^{\prime} \right)-\nabla A(\vec{x})\right\|^{2}+\lambda \underset{i=1}{\overset{3}{\sum}} \underset{j=1}{\overset{3}{\sum}} \frac{\partial {{u}_{i}}}{\partial {{x}_{j}}}\right]\text{d}\vec{x}\end{eqnarray} \tag{ 4 }$$

where $\nabla$ is the gradient operator, $~\left\|\cdot \right\|$ is the Euclidean norm, $\partial {{u}_{i}}/\partial {{x}_{j}}$ indicates the partial derivative of the component of the displacement $\vec{u}$ along the ith coordinate axis with respect to the component of $\vec{x}$ along the jth coordinate axis.

From the analysis in section 2.1, we know that the gradient vector can be used as the similarity measure to avoid the effect of CBCT intensity inaccuracy. So in this model the minimization of gradient distance $\left\|\nabla B\left(\vec{x}^{\prime} \right)-\nabla A(\vec{x})\right\|$ matches the two images in the gradient. For the smooth item, which is usually a local neighborhood filter, just like a Gaussian filter in the Demons algorithm, it directly affects the smoothness of the image. In this model, $\lambda \underset{i=1}{\overset{3}{\sum}} \underset{j=1}{\overset{3}{\sum}} \partial {{u}_{i}}/\partial {{x}_{j}}$ plays this role implicitly in an average way. However, it is likely that some pixels used in this average process do not situate on the same gray level as the processed one. As illustrated in figure 3(a), their participation will blur the image edge and degrade the registration performance. In order to alleviate this effect, we exclude these pixels using the LPA-ICI algorithm when the image proceeds by the multi-resolution strategy (Sun and Han 2009). The anisotropic neighborhood, which belongs to its center pixel, is shown in figure 3(b). Since the pixels in the anisotropic neighborhood belong to the same intensity level as its center one, these pixels are more suitable for smoothing the intensity of the center pixel. The minimization of the smooth item ${{\left(\lambda \underset{i=1}{\overset{3}{\sum}} \underset{j=1}{\overset{3}{\sum}} \partial {{u}_{i}}/\partial {{x}_{j}}\right)}_{\text{LPA-ICI}}}$ ensures a smooth and continuous displacement field. So, we modify (4) as:

$$\begin{eqnarray}&&\varepsilon \left(\vec{u}\right)=\underset{\vec{x}\in {{R}^{3}}}{\int} \left[{{\left(\left\|\nabla B\left(\vec{x}^{\prime} \right)-\nabla A(\vec{x})\right\|^{2}+\lambda \underset{i=1}{\overset{3}{\sum}} \underset{j=1}{\overset{3}{\sum}} \frac{\partial {{u}_{i}}}{\partial {{x}_{j}}}\right)}_{\text{LPA}-\text{ICI}}}\right]\text{d}\vec{x}\end{eqnarray} \tag{ 5 }$$

Zoom In Zoom Out Reset image size

Through the calculus of variations (Xu 2000), we have:

$$\begin{eqnarray}&&\lambda {{\nabla}^{2}}\vec{u}(\vec{x})-\left\|\nabla B\left(\vec{x}^{\prime} \right)-\nabla A(\vec{x})\right\|{{\vec{g}}_{s}}(\vec{x})=0\end{eqnarray} \tag{ 6 }$$

where ${{\nabla}^{2}}$ is the Laplacian operator and ${{\vec{g}}_{s}}(\vec{x})=\partial \left\|\nabla B\left(\vec{x}^{\prime} \right)-\nabla A(\vec{x})\right\|/\partial \vec{x}.$ Then a finite difference scheme (Press et al 1992, Lu et al 2004) was used to solve (6) in which the iterative formula of deformable registration is given by:

$$\begin{eqnarray}&&\text{d}{{\left(\vec{u}(\vec{x})\right)}_{n}}^{(k)}=\frac{\lambda {{\left({{\nabla}^{2}}\vec{u}(\vec{x})\right)}_{n}}-\left\|\nabla B\left(\vec{x}^{\prime} \right)-\nabla A(\vec{x})\right\|{{\left({{{\vec{g}}}_{s}}(\vec{x})\right)}_{n}}}{\lambda +{{\left({{\left({{{\vec{g}}}_{s}}(\vec{x})\right)}_{n}}\right)}^{2}}}\end{eqnarray} \tag{ 7 }$$

2.2.2. The bi-directional force to diffusion process.

As can be seen from (6) that ${{\vec{g}}_{s}}(\vec{x})$ and $\left\|\nabla B\left(\vec{x}^{\prime} \right)-\nabla A(\vec{x})\right\|$ have to be re-calculated in each iteration, which greatly increases the computational burden of the GFFD algorithm compared to the intensity-based algorithms. So it is inefficient to utilize gradient information (${{\vec{g}}_{s}}(\vec{x})$) to drive the deformation alone. Similar to the Demons algorithm, we assume that the diffusion is bi-directional. That is to say 'GFFD' will produce two forces at any point in the image. One force allows the 'moving' object to diffuse into the corresponding 'static' image and another allows the 'static' object to diffuse into the corresponding 'moving' object. Different from the Demons algorithm, we will not introduce the 'active' force, but combine ${{\vec{g}}_{s}}(\vec{x})$ with ${{\vec{g}}_{k}}\left(\vec{x}^{\prime} \right),$ which allows the 'static' object to diffuse into the corresponding 'moving' object to make a 'bi-directional' force. Similar to ${{\vec{g}}_{s}}(\vec{x}),{{\vec{g}}_{k}}\left(\vec{x}^{\prime} \right)$ is defined as ${{\vec{g}}_{k}}\left(\vec{x}^{\prime} \right)=\partial \left\|\nabla B\left(\vec{x}^{\prime} \right)-\nabla A(\vec{x})\right\|/\partial \vec{x}^{\prime} .$ To make both ${{\vec{g}}_{s}}(\vec{x})$ and ${{\vec{g}}_{k}}\left(\vec{x}^{\prime} \right)$ information work and accelerate the convergence, we modified (7) as:

$$\begin{eqnarray}&&\text{d}\left(\vec{u}(\vec{x})\right)_{n}^{(k)}=\alpha \frac{\lambda {{\left({{\nabla}^{2}}\vec{u}(\vec{x})\right)}_{n}}-\left\|\nabla B\left(\vec{x}^{\prime} \right)-\nabla A(\vec{x})\right\|{{\left({{{\vec{g}}}_{s}}(\vec{x})-{{{\vec{g}}}_{k}}\left(\vec{x}^{\prime} \right)\right)}_{n}}}{\lambda +{{\left({{\left({{{\vec{g}}}_{s}}(\vec{x})-{{{\vec{g}}}_{k}}\left(\vec{x}^{\prime} \right)\right)}_{n}}\right)}^{2}}}\end{eqnarray} \tag{ 8 }$$

where the parameter α was used to adjust the step size. For digital images, the 3D integer vector was used to represent the image pixel coordinates. The approximation of ${{\nabla}^{2}}\vec{u}(\vec{x})$ can be expressed as:

$$\begin{eqnarray}&&{{\nabla}^{2}}\vec{u}(\vec{x})=\frac{\sum\nolimits_{t=-1}^{t=1}{{\left(-1\right)}^{t+l}}u\left(i+t,j,k\right)+\sum\nolimits_{t=-1}^{t=1}{{\left(-1\right)}^{t+l}}u\left(i,j+t,k\right)+\sum\nolimits_{t=-1}^{t=1}{{\left(-1\right)}^{t+l}}u\left(i,j,k+t\right)}{6}\end{eqnarray} \tag{ 9 }$$

The equation used to calculate static image gradient was:

$$\begin{eqnarray}&&\nabla A\left(i,j,k\right)=\frac{1}{2}{{\left(\begin{array}{c} A\left(i+1,j,k\right)-A\left(i-1,j,k\right) \\ A\left(i,j+1,k\right)-A\left(i,j-1,k\right) \\ A\left(i,j,k+1\right)-A\left(i,j,k-1\right) \\ \end{array}\right)}^{T}}\left(\begin{array}{c} {\vec{I}} \\ {\vec{m}} \\ {\vec{n}} \\ \end{array}\right)\end{eqnarray} \tag{ 10 }$$

where $\vec{l},\ \vec{m}\ \text{and}\ \vec{n}$ denote the unit vectors along three coordinate axes. The moving image gradient $\nabla B\left(i,j,k\right)$ was calculated by (10) and the moving image gradient $\nabla B\left(\vec{x}^{\prime} \right)$ at any point was obtained by tri-linear interpolation. The formulae to calculate ${{\vec{g}}_{s}}\left(i,j,k\right)\ \text{and}\ {{\vec{g}}_{k}}\left(i,j,k\right)$ are respectively:

$$\begin{eqnarray}&&{{\vec{g}}_{s}}\left(i,j,k\right)=\frac{1}{2}\left(\begin{array}{*{35}{l}} \|\nabla B\left(\vec{x}^{\prime} \right)-\nabla A\left(i+1,j,k\right)\|-\|\nabla B\left(\vec{x}^{\prime} \right)-\nabla A\left(i-1,j,k\right)\| \\ \|\nabla B\left(\vec{x}^{\prime} \right)-\nabla A\left(i,j+1,k\right)\|-\|\nabla B\left(\vec{x}^{\prime} \right)-\nabla A\left(i,j-1,k\right)\| \\ \|\nabla B\left(\vec{x}^{\prime} \right)-\nabla A\left(i,j,k+1\right)\|-\|\nabla B\left(\vec{x}^{\prime} \right)-\nabla A\left(i,j,k-1\right)\| \\ \end{array}\right)\left(\begin{array}{*{35}{l}} {\vec{l}} \\ {\vec{m}} \\ {\vec{n}} \\ \end{array}\right)\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{{\vec{g}}_{k}}\left(i,j,k\right)=\frac{1}{2}\left(\begin{array}{*{35}{l}} \|\nabla B\left(i+1,j,k\right)-\nabla A(\vec{x})\|-\|\nabla B\left(i-1,j,k\right)-\nabla A(\vec{x})\| \\ \|\nabla B\left(i,j+1,k\right)-\nabla A(\vec{x})\|-\|\nabla B\left(i,j-1,k\right)-\nabla A(\vec{x})\| \\ \|\nabla B\left(i,j,k+1\right)-\nabla A(\vec{x})\|-\|\nabla B\left(i,j,k-1\right)-\nabla A(\vec{x})\| \\ \end{array}\right)\left(\begin{array}{*{35}{l}} {\vec{l}} \\ {\vec{m}} \\ {\vec{n}} \\ \end{array}\right)\end{eqnarray} \tag{ 12 }$$

Similarly, ${{\vec{g}}_{k}}\left(\vec{x}^{\prime} \right)$ was also obtained by tri-linear interpolation. In order to further improve the computation efficiency, the proposed algorithm was implemented by the GPU parallel computing. This is the objective of the next subsection.

The multi-resolution strategy and the bi-directional force can improve computation efficiency, but still cannot meet the clinical requirement needs. One way to further speed up the algorithm efficiency is to port them from the central processing units to the GPU. The gradient-based algorithm has thus been implemented on the GPU in OpenCL framework for parallel computing. The GPU used in this paper has 240 processors. The parallel-executed function is called a 'kernel', which was programmed by an extended C language.

The gradient-based algorithm comprises five GPU kernels and a CPU program. The CPU program controls the main program flow and invokes the GPU kernels for parallel computing. The five GPU kernels are InitializeMap kernel, DownSample kernel, GetGradient kernel, Registration kernel and UpSampleMap kernel, respectively. The InitializeMap kernel was used to initialize the displacement field in the video memory to 0. The DownSample kernel was used to down-sample the images. The pixel value of the down-sampled images is the mean value of all pixels in the down-sampling region. The GetGradient kernel was used to calculate the gradient field of the down-sampled images, and the Registration kernel implemented the process of updating the displacement vector by one-time iteration on a down-sampled pixel. The procedure of the gradient-based algorithm is described as follows:

• (a)  
  First, the FBCT and CBCT images were transferred from the host memory to the video memory;
• (b)  
  All the displacement field values were set to 0 through the InitializeMap kernel;
• (c)  
  The image down-sampling was done by the DownSample kernel;
• (d)  
  The gradient field was calculated by the GetGradient kernel;
• (e)  
  At each down-sampling ratio, the displacement field was updated for 16 times via the Registration kernel;
• (f)  
  After 16 iterations were completed, we went to step 3 to acquire the next down-sampling ratio image by the DownSample kernel and started the image registration at the new down-sampling ratio;
• (g)  
  After the registration was finished, the displacement field was transferred to the host memory.

There are 16 KB registers on each multi-processor in the GTX 285. The number of threads set in each block needs to be determined in combination with the number of registers occupied by each Kernel. The total registers should not exceed the 16 KB limit; otherwise it will affect the calculation performance.

Six CT data sets (three for head (cases 1–3) and three for lung (cases 4–6)) from 6 patients and lung phantom data were collected. Each CT data set contains a planning FBCT image sequence and an on-line CBCT image sequence. The planning FBCT images were acquired using a brilliance big bore 16 slice CT (Philips Medical System). The CBCT images were acquired using an on-board imager (OBI, Varian Medical Systems, Palo Alto, CA). The size of the image and the voxel in FBCT and CBCT are summarized in table 1. Prior to image registration, FBCT and CBCT images were re-sampled to the same voxel size. FBCT and CBCT images were re-sampled at 0.80 × 0.80 × 2.5 mm³ in case 1. All the other CBCT images were re-sampled in accordance with corresponding FBCT pixel size. The pixel intensity was acquired by the tri-linear interpolation in the re-sampling step.

All experiments were performed on a desktop computer with a 1.87 GHz Core2 6320 (Intel Corporation) CPU and 2G RAM. To reduce the computation time, both the Demons algorithm and the GFFD algorithm were implemented by OpenCL parallel programming. An ENGTX285 graphics card (ASUS Corporation) was used. It had 1G video memory and a GTX285 graphics processing unit (Nvidia Corporation) with 240 stream-processors. In the GFFD algorithm, the parameter λ was set to 80 000 and the number of iterations was set to 16 for each down-sampling ratio. The Demons method was applied to these cases for comparison purposes. A Gaussian smoothing kernel with standard deviation σ was set to 2.0 and the number of iterations was set to 16 for each down-sampling ratio.

In this section, the accuracy of the GFFD algorithm is demonstrated with phantom data. The CT and CBCT slices from the same lung phantom case are shown in figure 4. To get the known deformation field, the evaluation framework of the non-rigid registration algorithm proposed by Urschler et al (2007) was introduced. In this deformation model, two parameters, the grid size α and the possible extent of the random deformation β that is applied to each grid point, are adjustable. Three sets of distorted CBCT images with different degrees of synthetic deformations were generated. The deformations are transformed based on the grid size α = 64 and with maximum deviations (β) of 4, 6 and 8, respectively. To evaluate the registration algorithm quantitatively, the mean absolute differences (MAD) and standard deviations (SD) of displacement field between the synthetic ground truth and the computed deformation field were computed.

$$\begin{eqnarray}&&\text{MAD}=\frac{1}{N}\underset{i=1}{\overset{N}{\sum}} |{{W}_{i}}-{{K}_{i}}|\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\text{SD}=\sqrt{\frac{1}{N-1}\underset{i=1}{\overset{N}{\sum}} {{\left({{W}_{i}}-\text{MAD}\right)}^{2}}}\end{eqnarray} \tag{ 14 }$$

where N represents the total number of pixels, W_i denotes the ith vector of the calculated deformation field and K_i defines the ith vector of the known displacement field. The matching extent of the edges between the deformed image and the target image was assessed with the mean absolute error (MAE), which is defined as follows:

$$\begin{eqnarray}&&\text{MA}{{\text{E}}_{\text{edge}}}=\frac{1}{N}\underset{i=1}{\overset{N}{\sum}} \left|I_{1}^{\text{edge}}(i)-I_{2}^{\text{edge}}(i)\right|\end{eqnarray} \tag{ 15 }$$

where $I_{1}^{\text{edge}}(i)\ \text{and}\ I_{2}^{\text{edge}}(i)$ are the binary Canny edge images of image I₁ and I₂. All of the pixels used in this experiment are on the image edge (right–left and anterior–posterior directions) because (1) the accurate deformation of their organ walls rather than the fillers of these two organs is more critical in therapy treatment (Gu et al 2013); (2) In the lung phantom, the largest region is smooth (see figure 4) which directly influences the veracity calculation results on both the Demons and GFFD algorithms.

Zoom In Zoom Out Reset image size

Tables 2 and 3 show the MAD, SD and MAE_edge of the registration errors obtained with the two algorithms. From these two tables, we can see that both the Demons and GFFD algorithms reduce MAD, SD and MAE_edge between CT and CBCT images and the GFFD method outperforms the Demons algorithm in terms of MAD, SD and MAE_edge.

Different from the Demons algorithm, the similarity is assessed with the 3D gradient vector in the GFFD method. In order to discuss the impact of voxel size on the GFFD algorithm, the lung phantom data set was taken and reconstructed at different voxel size. In fact, the thickness can only be selected as 1.5, 3 and 6 mm and the resolution of CT scans on the transverse section cannot be changed in the clinical brilliance big bore 16 slice CT. Image size and voxel size are summarized in table 4.

As we know, both CBCT and CT should have the same spatial resolution in order to proceed to the registration process. Therefore, resampling is necessary if two images have different resolution. In this paper, the CBCT images were resampled in accordance with corresponding FBCT pixel size. The registration results with different voxel size are shown in figure 5. As can be seen, when the CT voxel size is much smaller than CBCT (case 1), the deformed CT image edges mismatch the CBCT edges. When the CT voxel size is much bigger than CBCT, the registration result is good, but the intensity of the registered image does not represent the true intensity of the current region. The reason is that when the FBCT voxel size is much bigger or smaller, it may have different structures included in the same region between CT scans and CBCT scans. More specifically, we can illustrate it via the similarity measure-3D gradient vector. From equation (10), it can be seen that six points (three slices) are needed if we want to compute the gradient at a point. In Case 1 as an example, if we want to compute the gradient of the image shown in figure 6(b), the other two slices shown in figures 6(a) and (c) are needed. Although the CT image and CBCT image in figure 6(b) represent the same region of the lung phantom, the gradient is very different. Thus, the gradient mismatch and registration error occur. But it is easier to solve compared with the impact of inconsistent intensity between two modalities of the intensity algorithm. If the voxel size of the CT image is similar to that of the CBCT image, there is no difficulty. Fortunately, this is often the case in practice, as listed in table 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The displacement field generated by the registration algorithm was applied to the markers manually labeled on the FBCT images. The principle of marker selection is that we can clearly and accurately find it in both CT and CBCT images. An effort was made to distribute the landmarks as uniformly as possible throughout the images. So, the marker is mainly located at the vertices of the structure (e.g. cervical vertebrae, oropharynx, ramus of mandible) in the head and vascular as well as bronchial bifurcations in the lung. The CT and corresponding CBCT image were imported into the in-house-developed software (see section 4) and initially aligned through the iso-center to facilitate manual selection of landmark pairs. After the marker was defined in the CT images, two experts independently labeled the corresponding markers. Two groups of labeled results were compared by the experts and the landmarks were removed if the same corresponding landmark was labeled on different slices. The final corresponding landmarks were decided according to the consensus result. In order to evaluate the precision of the manual landmark selection, inter-observer variability, characterized by the MAD (13) and SD (14) between the target points labeled by two experts in CBCT were computed. For three head cases, the MAD and SD were 0.74 and 0.65, respectively. The maximum landmark displacement between two experts was 1.583 mm. For three lung cases, the MAD and SD were 0.83 and 0.85, respectively. The maximum landmark displacement between two experts was 1.612 mm. The MAD (13) and SD (14) between the reference markers in the CT and the corresponding markers in the CBCT images were used as the registration error to assess the accuracy of the registration algorithm. The difference between the two groups was determined by the t-test with p < 0.05, which is considered as a statistically significant difference and a 95% confidence interval (CI) of MAD was compared (Chen et al 2011). As illustrated in table 5, for the head, the MAD obtained by the GFFD algorithm (95% CI: 0.87 ~ 1.15) is always smaller than the Demons method (95% CI: 1.22 ~ 1.58). For the lung, we can get a similar result. Meanwhile, a statistically significant difference between the Demons algorithm and the GFFD algorithm (p < 0.05) is observed for both head and lung cases. The experimental results demonstrate that the GFFD is more accurate than the Demons algorithm. Furthermore, from table 5 we can conclude that the GFFD is more robust than the Demons method in terms of standard deviations.

In section 2.1, the registration error of the triceps (case 5) due to intensity inaccuracy in the CBCT images was shown in figure 1(d) and the reason was discussed by analyzing the deformation force on the borders of the triceps. It is indicated that the inconsistent intensities between two modalities can affect the accuracy of the registration. But, after registration by the GFFD algorithm, the regions with scatter artifacts indicated by the green arrows (a) and (b) in figure 7 are well-registered. The result indicates that the GFFD algorithm can avoid the effect of the CBCT intensity inaccuracy and get more accurate results.

Zoom In Zoom Out Reset image size

An important parameter in the GFFD algorithm is the smooth factor λ. A larger λ ensures the continuity of the displacement field, but limits the deformability of the displacement field while a smaller λ allows for dramatic deformation but may result in the loss of geometric continuity in the deformed image. Figure 8 shows the convergence and final results of the MAD and the SD in the iterative process when different values of λ were adopted in the GFFD algorithm for all lung cases. Figure 8 shows when λ is above 80 000, the smaller λ is, the faster convergence of the mean error will be; yet when λ is smaller than 40 000, the registration algorithm becomes unstable due to the weak constraints. Therefore, in order to compromise between the continuity and deformability of the registration, the value of λ was set to 80 000 in the GFFD algorithm.

Zoom In Zoom Out Reset image size

On the other hand, regarding computational expense, the GFFD requires the same number of iterations (16 iterations) as that of the Demons algorithm to achieve an accurate registration result. But due to the current state-of-the-art implementations of the Demons algorithm and re-calculation of the gradient information in the iteration of the GFFD algorithm, the computational cost of the GFFD algorithm is higher than that of the Demons algorithm. But compared with the GFFD algorithm without the bi-directional force (the displacement field is generated by (7) through iteration), the number of iterations could be approximately reduced by half. With the use of the GPU, the computation time of the GFFD algorithm for each group of clinical data ranged from 8 to 13 s.

After deformable registration of planning CT and CBCT, some applications with on-line ART by using the deformation fields can be executed. Automatic organ propagation and direct aperture optimization-based plans for improving the treatment process are such applications. An in-house-developed software experiment platform was used for on-line ART. It is composed of five modules: Dicom input and output, image registration, dose calculation, plan optimization and evaluation. The dose verification results show that the error was within 3% between calculation results and actual measurement results. In order to evaluate the plan, three plans for the same head and neck patient were compared: (1) the original plan (the 'CT-plan'); (2) the the original plan delivered to the CBCT (the 'unadapted-plan'); and (3) the re-optimized plan.

According to the deformation field, the target and organ contours were automatically propagated from the planning CT to CBCT in order to save tedious time of organ delineation. Figure 9 illustrates the automatic propagation of the head and neck patient. As can be seen from figure 9, the contours of the planning CT directly overlaid on the CBCT, mismatch the tissues and organs particularly in the left parotid and body indicated by the red arrow. While the contours propagated by the GFFD algorithm is good. The volume overlapping was assessed by the Dice similarity coefficient (DSC). The DSC value between the contours propagated by the GFFD algorithm and the contours edited by the oncologist are 96 and 94% for the parotid and GTV separately. The radiation oncologists concluded that the results are helpful for clinical re-planning.

Zoom In Zoom Out Reset image size

A multi-leaf collimator (MLC) adjustments method based on the deformation fields (Ahunbay et al 2010) was applied to the on-line modification of the radiotherapy plan, so as to reflect the changes in the anatomical structure. The result is illustrated in figure 10. As can be seen, there is a small amount of planning target volume outside of the leaf edge before aperture-optimization (middle case), while it provided significantly complete target coverage after optimization (right case).

Zoom In Zoom Out Reset image size

After completing the re-optimization process, three DVH were analyzed for radiation therapy, in order to determine whether the adapted plan is even necessary and whether the adapted plan can achieve dose coverage at a level comparable to the original CT-plan (Wu et al 2008). As shown in figure 11, it would cause significant variation in dose coverage to the plan target volumes (PTV), if the original CT-plan delivered to the anatomy-of-the-day directly. For the un-adapted plan, the dose to volume of PTV above the prescribed dose (45 Gy) is 92%. This is unacceptable. For the adapted plan, the dose to volume of PTV above the prescribed dose (45 Gy) is 96%. It can satisfy the dose coverage criterion. For the organ at risk, the results were complicated. The DVH of the right parotid (parotid-R) for the un-adapted plan was better compared to the other two plans because the majority of the parotid-R volume is far away from the beam-coverage region. On the other hand, the un-adapted plan delivers a much higher dose to the spinal cord and left parotid (parotid-L). For the spinal cord, the maximum dose was more than the prescribed dose (42 Gy), while the dose is below this limiting value in both the CT-plan and adapted plan.

Zoom In Zoom Out Reset image size

The adapted plan based on the GFFD registration was compared with that based on the Demons. A new re-optimized plan was acted on the contours propagated by the Demons and the GFFD algorithm to produce these two plans. For the Demons algorithm, the values of the DSC for GTV between the propagated contours and the contours edited by the oncologist are 90%. This value increases to 94% when the contours are propagated by the GFFD algorithm. Compared with the Demons algorithm, the propagation accuracy of the GFFD increases significantly, while there is very little improvement in the DVH (figure 12). This is mainly because almost all the propagated contours are located within the high dose region (PTV). The adapted plan based on the Demons delivers more dosage to the parotid-L, while there is no difference for the parotid-R. The results indicate that the adapted plan based on the GFFD registration algorithm is better than or equivalent to the adapted plan based on the Demons registration.

Zoom In Zoom Out Reset image size