Deformable image registration in radiation therapy.

The number of imaging data sets has significantly increased during radiation treatment after introducing a diverse range of advanced techniques into the field of radiation oncology. As a consequence, there have been many studies proposing meaningful applications of imaging data set use. These applications commonly require a method to align the data sets at a reference. Deformable image registration (DIR) is a process which satisfies this requirement by locally registering image data sets into a reference image set. DIR identifies the spatial correspondence in order to minimize the differences between two or among multiple sets of images. This article describes clinical applications, validation, and algorithms of DIR techniques. Applications of DIR in radiation treatment include dose accumulation, mathematical modeling, automatic segmentation, and functional imaging. Validation methods discussed are based on anatomical landmarks, physical phantoms, digital phantoms, and per application purpose. DIR algorithms are also briefly reviewed with respect to two algorithmic components: similarity index and deformation models.

Introduction

The aim of radiation therapy (RT) is to deliver clinically effective dose to the target while minimizing dose to the surrounding normal tissues. Diverse advanced techniques have been introduced into the field of radiation oncology during the last two decades in order to get closer to this goal. One of the most impactful advancements is a sophisticated beam delivery method with intensity modulation (i.e., intensity-modulated radiation therapy and volumetric modulated arc therapy). Another instrumental advancement is combining imaging into treatment, such as on-line imaging and integrating multi-modality images. Modulated beams generate more conformal dose distribution to a target while sparing nearby normal structures even when those normal structures are located at the concave region of the target [1]. However, the conformal dose distribution is only achieved on the planning image set which is a snapshot of patient anatomy during a long treatment period. Inter- and/or intra-factional anatomy changes from the planning image set may deteriorate the dose conformality in actual delivered dose distribution. The use of more frequent on-line imaging, such as cone beam computerized tomography (CBCT) [2], CT on the rail [3], and integrated magnetic resonance imaging (MRI) [4], can detect these anatomic changes and aid in correcting or minimizing the effect of such changes. Multimodality images like positron emission tomography (PET) and MRI are used to define the target more accurately and can help reduce the clinical target volume margin [5,6]. As a consequence, the number of imaging data sets taken during the treatment has increased significantly in advanced radiation therapy.

The increased amount of imaging data allows for analysis of the relationship between a target response and delivered dose and/or a target shape/position change depending on the condition of nearby normal organs (e.g., whether the bladder and/or rectum is filled or empty). However, this is only possible when proper image registration is achievable. The imaging data may be taken within different timeframes (i.e., multi-temporal), from different imaging devices (i.e., multi-modality), and/or of certain populations of patients (i.e., multi-subject). By way of explanation, while two image sets are aligned on a reference coordinate, a pixel/voxel on the same coordinate may not necessarily represent the same anatomical structure. Methods used in image registration aim to resolve this issue.

Image registration can simply be described as a way to find the spatial correspondence between two or among multiple sets of images [7,8]. Consider two sets of images, a fixed one, F(x) and a moving one, M(x'). Image registration tries to find the best transform, T(x') that minimizes the difference between F(x) and M(T(x')). Transform is the sum of the local position vector in the moving image, x', and the displacement vector, u(x'). Thus, an ideal image registration between two image sets can be expressed as:

Image registration can be categorized into two groups: rigid and non-rigid. Non-rigid image registration is also known as deformable image registration (DIR). In rigid image registration (RIR), all pixels move and/or rotate uniformly so that every pixel-to-pixel relationship remains the same before and after transformation. In DIR, however, those pixel-to-pixel relationships change (Fig. 1).

Fig. 1.

Two image sets, a fixed image (A) and a moving image (B), were registered using image registration algorithms. Rigid image registration could not register the four sharp corners of the rectangle in the moving image into the rounded boundary in the fixed image. Deformable image registration locally deformed the four sharp corners with a different amount of deformation (or displacement). (C) The 2D deformation vector field (DVF) was displayed as blue arrows with the edge of the moving image object. The size and direction of the arrows represent the magnitude and direction of DVF. The magnitude of deformation is the largest at the corners and gradually decreases. (D) The deformed results with DVF.

RIR has long been used in radiation therapy (e.g., CT-to-MR information fusion). It is widely recognized that RIR is very effective in cases when no anatomic change is expected. Sometimes, however, patients do experience anatomical structure changes due to weight loss, tumor shrinkage, and/or physiological organ shape variation, which often cannot be handled by RIR at all. In comparison to RIR, DIR has a significantly greater number of degrees of freedom (DOF). Therefore, DIR can manage local distortion between two image sets (i.e, anatomical structure changes). A DIR algorithm consists of 3 core parts: objective function, optimization method, and deformation model. Thus, it can be categorized depending on the method used in each core part.

DIR has been studied since the early 1980s and its early clinical applications can be found in the field of brain surgery and neuroscience [9]. DIR had not been widely studied in RT until the late 1990s mainly due to limited computing power and the lack of interest. The size of a typical volumetric image in RT was much larger than that in brain surgery or neuroscience thus, required higher computing power to handle images in larger sizes. Up until the mid-1990s, the most common treatment techniques were 2D planning or 3D conformal radiation therapy. The dose distribution inside a field was assumed to be homogeneous so no complicated algorithm was needed to calculate the delivered dose. Patients often had only one or two volumetric image sets.

A new RT concept, adaptive radiation therapy (ART), was proposed in the late 1990s [10]. A key feature of ART was to update/modify a treatment plan based on the systematic feedback of measurements [11]. The systematic feedback could be an accumulated dose distribution delivered through a certain number of fractions during a multi-fraction treatment course [12] and/or the predicted patient anatomy at a certain point during the course of treatment [13]. Yan et al. [12] introduced a dose accumulation method based on a DIR. Shortly thereafter, the concept of DIR began to gradually adapt into RT. It has since been widely studied by virtue of advances in computing power, modulated beam delivery techniques, and on-line imaging.

In this review article, DIR techniques are discussed in terms of clinical application, validation and algorithm. The majority of articles reviewed were from major radiation therapy journals, while a few critical articles published in other fields (e.g., computer science and medical imaging) were also included. Because this article is mainly intended for clinicians, technical descriptions were simplified and the section of algorithm, the most technical part, was intentionally placed close to the end of the article not to discourage readers reading this article. Note a few of inevitable mathematical expressions were included for better understanding in explaining the types of DIR algorithms.

Application

Although there is no boundary to implement DIR in the field of RT, the application of DIR could be categorized into four major fields: dose accumulation, mathematical modeling, automatic segmentation, and functional imaging. As stated above, the main role of DIR is to define spatial correspondence between two considered image sets. The only difference among these four applications is how the generated correspondence information is used. Dose accumulation uses the information to map a dose distribution from updated images to a reference image. Mathematical modeling uses the information itself to find the statistic of motion or deformation in considered organs. Automatic segmentation uses the information to map a contour from a reference image to updated images. Finally, functional imaging uses the information to evaluate the functional capability of a certain organ (e.g., ventilation ability in a lung).

1. Dose accumulation

Dose accumulation [12,14-16] was the first application of DIR in radiation therapy, and was used in order to better estimate the actual delivered dose. The planning dose on a planning image is not the real delivered dose during treatment due to inter- and intra-fractional motion/anatomy changes. Yan et al. [12] considered the dose accumulation during a course of treatment and attempted to assess inter-fractional deviations. The displacement vector field (DVF) was calculated from the planning image set (i.e., moving image) to each daily image set (i.e., fixed image) using a biomechanical model based on the finite element method (FEM). The dose distribution of a region of interest (ROI) was estimated using the DVF and the planning dose distribution. Christensen et al. [14] proposed a method of dose accumulation between two different treatment modalities, external radiation therapy and intracavitary brachytherapy for cervix cancer patients. The DVF to be used for combining dose between external beam therapy and brachytherapy was calculated from the planning image set of the external beam therapy to each brachytherapy image set using the volumetric fluid registration (or viscous fluid flow). Schaly et al. [15] presented a dose accumulation procedure during external radiation therapy for prostate using the thin-plate splines (TPS). The dose difference was up to 29% between the planned dose and the accumulated one. Velec et al. [16] investigated the effect of breathing motion (intrafractional deviations) on dose accumulation by using an FEM-based DIR algorithm. They calculated the DVF from the end-of-exhale (EOE) phase to the end-of-inhale (EOI) phase. The DVF was then utilized to estimate inter-phase dose distributions based on weighted interpolation between the EOE and EOI dose distributions.

Commonly, a dose accumulation method deforms dose grid points based on the calculated DVF from DIR and then simply adds up the deformed doses of a number of fractions on reference grids. One of challenging problems of the current dose accumulation process is the lack of a clear description about mass change during treatment periods. The cumulative deposited energy estimated on the reference image set may not be consistent with the sum of actual deposited energy on each image set. Note that ‘energy’ is different from ‘dose’ and it is ‘dose’ times ‘mass’ (i.e., E = D × M).

2. Mathematical modeling

Mathematical modeling is similar to a process of obtaining a probability density function for an object in 3D space. By defining inter- and intra-patient correspondence from multiple imaging data sets using DIR, a mathematical model of organ motion and deformation [13,17-19] could be generated and used for ART. Sohn et al. [13] proposed a method of modelling inter-fractional organ motion and deformation using repeated CT images of a prostate cancer patient. They defined the correspondences among CT images using the FEM based DIR. The organ motion was simulated using the principal component analysis (PCA). To eliminate the dependency of manual segmentation on segmenting subjects and improve segmentation efficiency, Nguyen et al. [17] introduced a population liver motion model using the FEM based DIR. They considered the summation of the involved patients’ liver structures as a population liver model. The population liver model was rigidly registered to the liver structure of each individual patient based on PCA in the EOE first, then non-rigidly registered using the FEM based DIR program to obtain the population liver model at EOE. The population liver model at EOE was rigidly registered to the liver structure of the patient at EOI based on center-of-mass. Then it was deformed to finally obtain the population liver motion model. The authors claimed that the model could be used in online motion assessment. Budiarto et al. [18] built a population prostate motion and deformation model using CT scans of prostate patients in order to describe the patient-specific local probability distributions of motions. The correspondence between inter- and intra-patient prostates was found using DIR based on thin-plate-spline robust-point-matching (TPS-RPM). The population model was described using a limited number of eigenmodes obtained from PCA. Oh et al. [19] proposed a method to model spatial and temporal geometric changes using weekly MR images of cervix canter patients. They defined the inter- and intra-patient correspondences through a common coordinate using a 3D surface DIR algorithm based on a parametric active contour. A population based anisotropic margin model was generated based on the calculated correspondences.

3. Automatic segmentation

Automatic segmentation [20-24] is also proposed as an application of DIR and has used pre-defined reference segment data and correspondence data based on DIR. Shekhar et al. [20] proposed the automatic segmentation methods using the B-spline based free-form deformation (i.e., DIR). They tested their method using 4D CT images of lung and abdomen cancer patients. In the automatic segmentation procedure discussed by another group [21-24], contours in a reference image were deformed to CT image sets obtained during treatment course by DVF defined by the ‘Demons’ DIR algorithm or the fact free-form DIR algorithm via calculus of variations. They evaluated their method on head and neck cancer patients [21,22]; head and neck, prostate, and lung cancer patients [23]; and breast cancer patients [24].

4. Functional imaging

The functional imaging technique had been proposed using the volume changes of each voxel in 4D images [25-27]. Guerrero et al. [25] proposed a method to simulate ventilation image using 4D CT images and DIR algorithm based on optical flow. The volume change of each voxel between the exhale image and inhale image was calculated based on the correspondence information from a DIR algorithm. Yaremko et al. [26] expanded the idea to 4D CT images and proposed intensity-modulated radiation therapy optimization using a lung functional image set generated based on 4D CT and DIR. A similar approach was tested by Yamamoto et al. [27].

In these approaches, DIR was applied to multiple anatomic image sets and the outcome of such process was functional information. Thus, it is termed with ‘functional imaging’ but different from typical functional imaging based on a functional imaging modality (e.g., PET or MRI unit).

Validation of DIR Algorithm

The result of deformable image registration is a deformation vector field which describes the correspondence between the fixed and moving images. The known correspondence information between the two image sets is essential for validating DIR algorithms. To define the correspondence information as a ground truth, there are three commonly used methods: (1) based on anatomical landmarks, (2) physical phantoms, and (3) digital phantoms.

1. Anatomical landmarks

Anatomical landmarks based validation is the most intuitive thus, considered the best approach. The performance of DIR algorithms can be directly estimated and compared by using anatomical landmarks in actual patient images. However, the anatomical landmarks are not enough to cover whole image areas. DIR is a rather local registration than global one in contrast to RIR. It is impossible to measure the accuracy of registration where the landmarks do not exist. It is difficult to define anatomical landmarks in certain organs, such as the rectum, small bowel, uterus, bladder, etc. Some researchers injected lipiodol markers [28] and used implanted stents or seeds [29]. The deformation results may be different with and without these markers, especially for DIR algorithms using image intensity similarity metrics, because these markers have different image intensity from surrounding tissues. Some researchers made attempts to collect patient image sets with manually identified anatomical landmarks information and then shared that information publically [30].

Brock and Deformable Registration Accuracy Consortium [31] published results of a multi-institution DIR accuracy comparison using lung, liver, and prostate patient images based on anatomical structures localized by radiation oncologists. In this paper, a total of 21 institutions sent deformation results using their own DIR algorithms and the authors analyzed the accuracy of deformation results based on the anatomical structures. They considered slice thickness of 4D CT data (2.5 mm) as a performance criterion. For lung 4D CT image sets, all replied 21 algorithms satisfied the criterion along all directions by the virtue of the consistent contrast between the two image sets. However, the maximum errors were up to 12 mm. For liver 4D CT, 7 out of 20 algorithms satisfied the criterion along three directions. The maximum errors were up to 13 mm. For the liver MRI-CT image sets, only 3 institutions sent the deformation results and the mean absolute error ranged from 1.1 to 5.0 mm. The prostate MRI image sets were found to be very challenging due to artifacts caused by gas in the rectum along with substantial deformation of the prostate. There were 3 investigators who did submit their prostate results.

Kadoya et al. [32] published a multi-institutional comparison study using 10 thoracic cancer patients with manually defined anatomical landmarks. Each institution reported DIR results of 10 patients’ data using commercial DIR softwares—RayStation (RaySearch Laboratories, Stockholm, Sweden), MIM Software (Cleveland, OH, USA), and Velocity (Varian Medical Systems, Palo Alto, CA, USA). The institution number which reported the results using RayStation, MIM software, and Velocity was 4, 5, and 3, respectively. In RayStation group, the institution used a hybrid DIR with ‘focus ROI’ option showed the highest accuracy. The results from 5 institutions of MIM group showed similar accuracies from each other within the group, and it was likely because there were few parameter settings that could be changed by the users. One institution manually applied rigid registration before applying deformable registration and obtained the best accuracy among the 5 MIM institutions. Velocity group exhibited relatively lower accuracy than RayStation and MIM at least with the given tested cases. The longest calculation time was reported by one of Velocity group but, it was still in the order of 3 minutes.

2. Physical phantom

The use of physical phantoms was proposed to evaluate DIR algorithms by generating known deformation [33]. Using a phantom, two or more sets of images were obtained with different statuses of deformation. The DVF was calculated based on the phantom correspondence information and was considered the ground truth. Next, a DIR algorithm was evaluated using the phantom image sets and the DVF from the DIR algorithm was compared with the ground truth. The distribution and density of the known correspondence information could be modified and optimized by changing the phantom design. The drawbacks of the phantom were its insufficient anatomical structures and unrealistic deformation.

Kashani et al. [34] compared the accuracy of DIR algorithms using image sets of a physical phantom proposed in [33]. They obtained inhale and exhale image sets using their anthropomorphic thoracic phantom which included a chest wall, skeleton, and compressible section. The phantom had 48 small, manually localized plastic markers. They tested 8 DIR algorithms from 6 institutions. All tested algorithms used an image intensity based similarity metric. One interesting finding was that while some methods showed relatively uniform accuracy throughout the whole phantom, others exhibited local dependency.

3. Digital phantom

The approach of digital phantoms is similar to that of physical phantoms. One image set is considered as a reference image set. Another image set is synthesized using a pre-defined DVF and considered a deformed image set. Then, a DIR algorithm is used to deform the reference image set to the deformed image set. The digital phantom gives the densest DVF information. Unlike the anatomical landmark or physical phantom approaches, however, the pre-defined DVF may be physically inappropriate. Pukala et al. [35] compared five commercialized DIR software programs using the digital phantom approach with 10 head and neck patient cases.

4. Application purpose

Another validation method is based on the purpose of the DIR application, like contour propagation and dose accumulation. Hardcastle et al. [36] used the physician-drawn contours as a validation metric to evaluate DIR algorithms. They propagated the gross tumor volume and OAR ROIs from a planning CT to a treatment CT for 22 patients based on the DVF from DIR algorithms. The quality of the propagated ROI was evaluated based on the Dice volume overlap score. Yeo et al. [37] developed a tissue-equivalent, deformable gel phantom. They obtained an initial CT image along with 3 deformed CT images and calculated the accumulated dose based on the DVF from 11 DIR algorithms. The accuracy of the dose accumulation was evaluated with the measured dose distribution based on the gamma index.

Types of DIR Algorithms

The process of image registration was concisely explained before reviewing the types of DIR algorithms. The goal of image registration is to identify the best transformation vector, T(x') where x' is a position vector (i.e., pixels) in a moving image, M, which maximizes the similarity between M(T(x')) and F(x) where F is a fixed image. Before starting the DIR process, the moving image is globally aligned with the fixed image based on a rigid body registration (6 DOFs: translation and rotation along x-, y-, and z-axis) or affine translation (12 DOFs: translation, rotation, scale, and shear along x-, y-, and z-axis). Next, the local area of the moving image is registered to the fixed image using DIR algorithms. The DIR process starts iteratively with respect to the optimization algorithms. In an iteration process, the DVF is generated according to a transformation model and the deformed moving image is generated. The objective function is also updated with the deformed moving image. The optimization algorithms attempt to locate local/global minima. Fig. 2 shows the flow chart of the process. DIR algorithms can be categorized based on three main parts; objective function, transformation model, and optimization algorithm.

Fig. 2.

Flow chart of deformable image registration process. The similarity index is calculated with given a fixed image and a moving image. The optimization algorithm tries to maximize the similarity index by changing deformation vector field (DVF) and the moving image is deformed based on the DVF. The similarity index is recalculated with the deformed moving image and the fixed image. This optimization process is done iteratively until the improvement of the similarity index reaches its target.

As mentioned earlier, the medical application of image registration had been studied in the field of brain surgery and neurosciences since the early 1980s [9]. Since then, a vast amount of image registration methods have been proposed. Review articles have been published in the following different fields; computing science [38], medical image [8], and radiation oncology [7,39]. Zitova and Flusser [38] reviewed the image registration algorithms in terms of objective functions (the feature- and area-based method). Sotiras et al. [8] meticulously reviewed many DIR algorithms based on transformation models describing deformation, matching criteria (objective functions), and optimization methods. Kaus and Brock [39] also categorized the DIR algorithms in terms of a similarity metric (objective functions) and transformation models. In this section, the DIR algorithms are briefly reviewed with respect to objective functions and transformation models.

1. Objective functions

An objective function in DIR methods depends on how to define the similarity between two image sets. Objective functions are commonly characterized into 2–3 categories: intensity-based (area-based or iconic), feature-based (contour-based or geometric), and combing of two objective functions.

1) Intensity-based objective function

The intensity-based objective function calculates the degree of similarity using image intensity. The assumption of this similarity index is that the pixel values (image intensity) of the same anatomical area are similar among image sets. This is a sound assumption with single-modality images. Several intensity-based objective functions have been proposed as a similarity index [7]; sum of squared difference (SSD) of image intensity, correlation coefficient (CC), and mutual information (MI).

Consider the image intensity of a fixed image and a deformed moving image, I(x) and I'(x'), respectively. The simplest objective function using image intensity is the SSD of image intensity between two images;

where Ω is the overlapped area of a fixed image and a deformed moving image and N is the number of voxels in Ω. SSD is sensitive to image intensity and works well with noiseless images from the same modality. It is possible to implement the optimization with fast convergence by virtue of the least-square form of the SSD.

CC between two images is less sensitive to noise in images than SSD and can be defined as such:

where and is the mean intensity in the fixed image and the deformed moving image, respectively.

CC can be implemented in an efficient computation because of its least-square form.

The concept of MI is from the information theory. The image registration is considered as a process to maximize the amount of shared information in two images. Hill et al. [7] thoroughly explained the concept of MI from the information theory. The MI between the fixed image, F and the deformed moving image, M’ can be calculated as such:

where H(F) and H(M') are defined as entropy of the fixed and the deformed moving images, respectively. H(F, M') is joint entropy between two images. P and P are the probabilities of intensity values at voxel x and x’ in two images, respectively. P is the joint probability of intensity values at voxel x and x’ occurring together.

MI is considered the standard intensity-based objective function for DIR algorithms with multimodality images [39].

2) Feature-based objective function

One of the advantages of the feature-based objective function is that the feature does not depend on the image intensity. This objective function is suitable for multimodality image registration. However, a significant drawback is that defining features in images can be quite difficult. Manually defined pairs of points of interest (POIs) or contours of ROIs are considered as the features to calculate objective functions. The manual definition of image features introduces inter- and intra-observer dependency.

The easiest way to implement the feature-based objective function is to calculate the squares of the distances between paired POIs in the fixed image and the moving images. Yan et al. [12] used the difference between squares of the distances of adjacent elements; is the distance from i-th volume element of the meshed organ to its one of adjacent j-th volume element in the moving and fixed image, respectively. For calculating this, they assumed that the meshed organ in the moving image set had the same volume as the corresponding organ in the fixed image set.

Brock et al. [40] used surface projection from the surface of a moving ROI to the surface of a fixed ROI. The projection was primarily perpendicular to the moving surface and the projection data was used as a boundary condition of their DIR algorithm. Oh et al. [19] applied a gradient vector flow field as an external force field. The surface meshes of ROIs had a minimum energy at the surface of target ROIs (contours).

3) Hybrid objective function

Some researchers proposed a hybrid objective function in order to solve the limitations of image-intensity based and feature-based objective functions. Weistrand and Svensson [41] proposed a DIR algorithm of which objective function was a linear combination of four terms: image similarity, grid regularization, shape based regularization, and penalty. Among these four terms, the first and last are similarity terms and the second and third are regulating terms to prevent physically unreasonable deformation. The image similarity terms are based on CC between a fixed image and a deformed image. The penalty term is the distance between predefined features in the moving and fixed image sets. The predefined features can be POIs and ROIs. The term is the sum of two distances. The first is the Euclidean distance between a vertex in the deformed triangular mesh and a signed distance from a fixed ROI contour. The second is the distance between corresponding POIs. The grid regularization term aimed to prevent physically unreasonable deformation by minimizing Dirichlet energy. The shape based regularization term is based on the triangular meshes of predefined contours. Because of the shape based regularization term, their algorithm worked well even with noisy images such as CBCT.

2. Transformation model

The correspondences of POIs or ROIs between two images are updated based on the transformation model to maximize the similarity index. If the transformation model can consider only several motion parameters, such as translation and rotation, the updated POIs or ROIs can only be changed globally. A deformation model requires significantly large number of motion parameters to achieve local transformation. With such large number of parameters, a DIR algorithm needs huge amount of memory and computation time for POIs/ROIs update.

The model can be categorized into two; a parametric model and a non-parametric model. A parametric model generates fine DVF as a linear combination of their basis functions. Typical parametric transformation models are a group of spline models. One of the advantages of the parametric model is that the local change of a point can be made by nearby points within a certain distance. This property significantly reduces computation time and memory required. In contrast, non-parametric models calculate transformation vectors of all points thus, require longer computation time and larger memory than parametric models do.

The deformation vector of a position vector x' in a moving image, u can be expressed as u(x') after k-th iterations. The deformation vector transforms x' into x in the fixed image.

The B-splines methods model the deformation as a linear combination of basis-splines [20]. The moving image with the size of N×N×N is deformed by interpolating the displacement vector (φ) of control points on a rectangular grid, C×C×C The deformation vector at a pixel on the moving image can be expressed as such:

where u, v, and w are the distances from x' to the nearest control points, (i, j, k), respectively. B represents the l-th basis of B-spline.

A multiresolution approach is possible by changing the grid spacing of control points.

TPS modeled the bending of a thin metal sheet by orthogonal directional forces. The control points are represented by the forcing area and the force represents the displacements from control points in the moving image and their corresponding control points in the fixed image. The control points are manually or automatically selected in the fixed and moving images [34]. The deformation vector at a pixel on the moving image can be expressed as such:

where coefficient matrix A and B, and the weighting factor, are calculated based on the n-pairs of matched control points, and U is a basis function to measure the distance from x' to i-th control points in p in the moving image. The deformation quality of B-spline and TPS depends on the accuracy of the control point’s correspondence.

The original ‘Demons’ algorithm used gradient information of a static fixed image (∇I(x)) in order to generate the Demons force that deformed the moving images [42]. Unlike the spline models which interpolated or smoothed the displacement vectors of a certain amount of control points, the Demons algorithm generated the local Demons force at each voxel and generated DVF at each voxel. The DVF at a position in the moving image is expressed as such:

where is the gradient of the fixed image at a position xk which is k-th corresponding position in the fixed image with x' in the moving image and is the intensity difference as external force. A limitation of the original Demons algorithm is that the force may not be strong to deform the moving images in the low gradient region.

The viscous fluid flow model considers the image deformation as viscous fluid [8]. This model allows the ROI in the moving image to be deformed to the ROI in the fixed image while maintaining the topology of the deforming structures, even during large nonlinear deformations such as the insertion of tandem and ovoid [14]. The deformations of viscous fluid flow are continuous, one-to-one, and differentiable. The model is governed by a Navier-Stokes equation and the displacement vector at a time, t and location, x' in the moving image can be described as such:

where δ is the time increment and v(x, t) is the deformation velocity. In the viscous fluid flow model, the force that deforms the moving image is proportional to the deformation velocity.

The force in the linear elastic model is proportional to the deformation itself, much like a spring. The linear elastic model is governed by Hooke’s law of elasticity, f = -kx , where f is the restoring force, k is the spring constant, and x is the change of spring length. If one voxel/pixel is deformed in the moving image, that deformation introduces the force and the other voxels are also deformed until forces archive equilibrium. This force and displacement relationship among all voxels in the moving image can be expressed as such:

where K is a stiffness matrix, u is displacement matrix of all considered element, and F is an elastic force matrix [43].

Challenges and Future Directions

A brief survey of the up-to-date DIR technology was provided in this article. The main focus was on its clinical applications; however, validation methods and major algorithms were also reviewed. In addition, a concise summary was provided in Table 1.

Table 1.

Summary of reviewed DIR application

Study	Transformation model	Application	Site
Yan et al. [12]	FEM-based linear elastic	Dose accumulation	Prostate cancer
Christensen et al. [14]	Viscous fluid flow	Dose accumulation	Cervix cancer
Schaly et al. [15]	Thin-plate splines	Dose accumulation	Prostate cancer
Velec et al. [16]	FEM-based linear elastic	Dose accumulation	Lung cancer
Sohn et al. [13]	FEM-based linear elastic	Mathematical modeling	Prostate cancer
Nguyen et al. [17]	FEM-based linear elastic	Mathematical modeling	Liver cancer
Budiarto et al. [18]	Thin-plate splines	Mathematical modeling	Prostate cancer
Oh et al. [19]	Parametric active contour	Mathematical modeling	Cervix cancer
Shekhar et al. [20]	B-splines	Automatic segmentation	Lung cancer and abdomen cancer
Chao et al. [21]	Demons algorithm	Automatic segmentation	Head and neck cancer
Lee et al. [22]	Calculus of variance	Automatic segmentation	Head and neck cancer
Wang et al. [23]	Commercial algorithm (Pinnacle)	Automatic segmentation	Head and neck, prostate, and lung cancer
Reed et al. [24]	Demons algorithm	Automatic segmentation	Breast cancer
Guerrero et al. [25]	Optical flow	Functional imaging	Breath hold CT of lung
Yaremko et al. [26]	Optical flow	Functional imaging	4D CT lung
Yamamoto et al. [27]	Calculus of variance	Functional imaging	4D CT lung

DIR, deformable image registration; FEM, finite element method; CT, computed tomography.

It is generally accepted that DIR has the potential to greatly improve the current radiation therapy process and accelerate the realization of personalized medicine in radiotherapy. To reach this desired outcome, however, there are still several aspects that pose challenges to overcome: (1) automatic image registration; (2) more rigorous validation methodology; (3) robust planning taking DIR uncertainty into account; and (4) ultimate automatization of critical processes in treatment workflow

DIR must deal with tasks such as complex non-linear and local distortion management, multi-modality image registration and multi-dimensional image registration, all of which make automatic image registration particularly challenging. Thus, algorithm development for automatic image registration has been a critical topic for many researchers, especially imaging scientists and computer software engineers. Another interesting issue in medical imaging is that organ deformation can be due to not only deformation itself but also actual mass change. Mass variation may not be a significant problem in the deformation process itself; however, it can cause extreme difficulty in mapping voxel-to-voxel radiation doses. The authors believe this is the most challenging obstacle to overcome in DIR application for radiation therapy.

Although significant attention has been given to DIR validation, establishing more rigorous validation methods still remains an open problem. One solution may be to construct an ideal benchmark phantom that can be voxel-by-voxel deformed, mas-variated and imaged in multi-modalities (e.g., CT, MR, and PET). To our best knowledge, there has not been such a phantom reported.

While DIR can bring opportunities of response evaluation and cumulative dose estimation, even during the treatment course, non-negligible uncertainties still exist. Considering that some treatment parameters could change based on such information in ART, it is necessary to analyze the reliability of the ART approach based on DIR. Robust planning mechanisms that do account for such DIR uncertainties can increase the confidence of ART practice.

Regardless of how impressive a new technology may be, the excessive time and effort put forth into implementing it will often prevents from being commonly adopted. Optimization of workflow is as critical as the effectiveness of the technology itself. Therefore, the authors believe that properly automatizing/integrating DIR applications into an overall workflow is of upmost importance.