|Home | About | Journals | Submit | Contact Us | Français|
Intensity modulated radiation therapy (IMRT) for cancers in the lung remains challenging due to the complicated respiratory dynamics. We propose a shape-navigated dense image deformation model to estimate the patient-specific breathing motion using 4D respiratory correlated CT (RCCT) images. The idea is to use the shape change of the lungs, the major motion feature in the thorax image, as a surrogate to predict the corresponding dense image deformation from training.
To build the statistical model, dense diffeomorphic deformations between images of all other time points to the image at end expiration are calculated, and the shapes of the lungs are automatically extracted. By correlating the shape variation with the temporally corresponding image deformation variation, a linear mapping function that maps a shape change to its corresponding image deformation is calculated from the training sample. Finally, given an extracted shape from the image at an arbitrary time point, its dense image deformation can be predicted from the pre-computed statistics.
The method is carried out on two patients and evaluated in terms of the tumor and lung estimation accuracies. The result shows robustness of the model and suggests its potential for 4D lung radiation treatment planning.
Currently, computed tomography is the major imaging modality for comprehensive evaluation of the lung in respiratory motion due to significant advances in its both temporal and spatial resolution . Multi-slice CT enables acquisition of a sequence of temporal 3D images in one breathing cycle. Developments in computerized cancer treatment planning techniques such as intensity-modulated radiation therapy (IMRT) require precise and efficient calculation of dose distributions to maximize the tumor treatment while avoiding endangering the surrounding normal tissues. However, the complicated respiration dynamics limit the application of IMRT for lung cancer treatment and diagnosis . Clinical 4D modeling of the thorax region in respiratory motion remains to be a challenging task.
Various models have been proposed to simulate and predict the complicated free-breathing system. Recognizing the hysteresis and irregular breathing patterns, auxiliary devices, such as the pencil-beam navigators  and the spirometer , have been used to obtain extra modeling parameters as the surrogate for the motion estimation. Pure image based models [5, 6, 7, 8, 9, 10] have been adopted to characterize the respiratory motion in the whole imaging space or within the major organ region, mostly through image registration techniques. In  a mean motion model for lung is proposed based on image registration both intra- and inter-patients, where tidal volume estimated by the air content in the lung is used to account for individual differences. However, the model is more or less limited in the applications to regular and repeated breathing patterns due to the lack of surrogate signals; In  a patient-specific model is trained in the dense image deformations space, where the diaphragm position was manually labeled as the surrogate. It can be observed that the regions around the diaphragm contain the motions of large magnitudes, which makes the diaphragm position an efficient navigator of the deformation dynamics. Nevertheless a more globally involved surrogate is needed to account for the motions in more superior regions of lungs. Considering the large coverage and the high-contrast intensity feature of the air-filled lung fields in the image, the shape changes of the lung is a promising surrogate signal.
In this paper we propose a pure image based modeling method of estimating patient-specific respiratory model for lung cancer radiation therapy, using the shape of the lung as the navigator of the dense image deformation space.
Our method begins from calculating dense diffeomorphic deformations from the end of expiration (EE) phase image to all other phases (Sec.2.1). In the same images, shapes of the anatomical structures are automatically extracted (Sec.2.2). By correlating the shape changes with the temporally corresponding image deformation fields in their separate reparameterized spaces, we compute a linear mapping function which maps a shape model temporal change to a dense image deformation field (Sec.2.3). This relative computationally expensive training process can be done completely off-line. Given a target image at an arbitrary phase, we only need the extracted shape to estimate the corresponding image deformation using the trained statistics.
The state-of-the-art symmetric diffeomorphic image registration method developed by Avants et al.  is used to obtain the dense deformation fields of all phases. This method provides smooth deformation in a multi-resolution framework, resulting in both forward and inverse mappings regardless of which image is fixed and which image is moving.
In a respiratory sequence, each image is denoted by It where t [0, 1] with 0 as the starting phase and 1 as the end phase. We use the end of expiration (EE) phase IEE as the base image due to its high reproducibility. For each image It with dimension of I × J × K, the corresponding dense deformation field is denoted by Ht = [ut(0, 0, 0), …, ut(i, j, k), …, ut(I, J, K)], where is the displacement vector field as the result of registration with the base image IEE (HEE itself is zero everywhere).
The respiratory motions can be observed from the shape changes of the lungs (left lung and right lung) in a 4D CT clearly, especially near the diaphragm regions. As the major source of motion in the imaging field, the geometrical variations of lungs are used as the surrogate for the temporally corresponding dense image deformation. Furthermore the air filled regions have high contrast boundaries, which makes the automatic segmentation feasible.
Binary lungs are extracted using an automated segmentation pipeline, similar to the approach in . First, a binary threshold is chosen carefully to exclude airway regions and to separate the two lungs. A 3D math-morphological “ball rolling” operator fills the holes after thresholding and remove unnecessary details in connecting regions of bronchi structure, airway and lungs. An example of resulting binary images are illustrated in Fig.1a.
The surface point distribution model (PDM) representing the shape of lungs is denoted by , where is the ith point on the surface and N is the number of surface points for each lung.
In order to train statistics on PDMs of all phases, correspondences among the surface points across all phases must be satisfied. A direct way of producing such group-wise correspondence is to generate all the models via deforming a common initial model. We start from a semi-manual triangle mesh model of the lungs. A diffusion process deforms the triangle mesh along surface normal directions to obtain a good fit of the binary image contours. Correspondence among the mesh vertices is provided by the use of an underlying medial representation (m-rep)  as the initial shape model, as shown in Fig.1b. The m-rep provides a fixed number of sparse surface points (the spoke ends) with surface normals (unit spoke vectors) for the diffusion process. The vertices of each diffused mesh is a PDM; see Fig.1c. In our experiments each lung has 450 surface points.
For 3D CT lung images with typical resolution of 512 × 512 × 100 the image deformation fields are storage demanding. dimension of 3D shape representation for lungs is also in the order of thousands. Multiplying by the time dimension, the calculation is computationally prohibitive. Dimension reduction becomes the natural solution. We use principal component analysis (PCA) on all phases to reparameterize the shape and dense image deformation space respectively. The reparameterization of the image deformation space is described as
where μH is the sample mean, is the ith eigenvector in the dense image deformation space, EH is the matrix of eigen-vectors, is the corresponding coefficient or the projection score of Ht in the ith eigendirection, and kH is the number of eigenmodes. Similarly, we have the equation for the point sets forming the surrogate:
In our experiments with the training sample consisted of 10 phases, kH = 2 and kP = 5 cover more than 90% of the total variance in their spaces respectively. The rest of the variations are discarded as sampling noise.
Linear correlation between the surrogate and the image motion has been hypothesized in [3, 9] and was shown to be effective in modeling the respiratory dynamics. Based on the same hypothesis, linear correlation between the shape variation space and the image variation space within the breathing cycle is assumed. Given and , we first center the data by , . We then form the training data matrix X = [X0, …, Xt, …, X1] and Y = [Y0, …, Yt, …, Y1]. The linear mapping matrix M, which gives Y = XM, can be calculated via the standard multivariate regression as follows:
Given a target image at arbitrary phase t* with an automatically extracted shape of the lungs, we directly project the shape onto the shape prior to get , and the corresponding center is recovered by .
In our experiments, the respiratory correlated CT (RCCT) data sets are provided by a 4-slice scanner (lightSpeed GX/i, GE Medical System), acquiring repeat CT images for a complete respiratory cycle at each couch position while recording patient respiration (Real-time Position Management System, Varian Medical Systems). The CT images are retrospectively sorted (GE Advantage 4D) to produce a series of 3D images at 10 respiratory time points. The time resolution is 0.5 s. The CT slice thickness is 2.5 mm. For each time point, tumor contours by experts’ delineation are used as the clinician truth.
Sequences of RCCT for two patients ( patient A and B in the plots) are tested. We quantified the estimation errors by applying the predicted deformations to the gross tumor volume (GTV) of the EE phase and measured the results with the clinicians’ delineation. Fig.3 shows an example of GTV delineation and the predicted 3D tumor colored by the surface distance to the clinician truth.
Three measurements are adopted to investigate the results using standard boxplots (median, interquantile range, largest observation and outliers are plotted): The Dice coefficient measures the volume overlap between the prediction and the ground truth GTV; The average surface distance measures the overall geometry fitness; The 90% quantile surface distance accounts for large discrepancies. In comparison, we plot the estimation errors without deformation estimation (no registration), in which the tumors are assumed static and EE phase GTV are used for all time points. This static situation illustrates the mobility of the tumors. We also plot the estimation errors when directly applying the image deformation to the GTV via the same symmetric diffeomorphic registration used in training. The results are summarized in Fig.3. The maximum error of 0.25cm surface distance is in fact within one voxel size, 0.1cm × 0.1cm × 0.25cm. In both patients, our shape-navigated statistical model largely reduces the estimation error from the static (no registration) method and is more robust than the registration method as a result of the noise removal during the reparameterization.
We also evaluated the estimation error of the shape of lungs. Automatic binary segmentation results were used as the ground truth. Since the results on left and right lungs are quite consistent with each other, we use the average numbers of the two lungs in the boxplots of Fig.4. Results without motion modeling ( no registration) are not in comparison merely because the motion of the lungs are obvious to observe. Consistent with the GTV evaluations, our method yields good estimations on lungs, slightly better than the registration method in terms of the average surface distance.
In the evaluated cases the proposed model is able to estimate the dense image deformation within the breathing sequence robustly. The dimension reduction in both image and shape deformation spaces not only retains anatomically relevant variation but also reduces noise. The hypothesis of linear correlation between the shape and image deformation space is shown to be effective. Given 4D image sequences both at planning time and treatment time, our next goal includes estimating the respiratory motion at the treatment time based on the statics trained on the sequence of planning time. Also, it is possible to extend the methodology to inter-patent studies to include more training samples for more robust statistics.