PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Med Phys. Author manuscript; available in PMC 2010 April 8.
Published in final edited form as:
PMCID: PMC2851165
NIHMSID: NIHMS110779

Iterative sorting for four-dimensional CT images based on internal anatomy motion

Abstract

Current four-dimensional (4D) computed tomography (CT) imaging techniques using multi-slice CT scanners require retrospective sorting of the reconstructed two-dimensional (2D) CT images. Most existing sorting methods depend on externally monitored breathing signals recorded by extra instruments. External signals may not always accurately capture the breathing status and may lead to severe discontinuity artifacts in the sorted CT volumes. This paper describes a method to find the temporal correspondences for the free-breathing multi-slice CT images acquired at different table positions based on internal anatomy movement. The algorithm iteratively sorts the CT images using estimated internal motion indices. It starts from two imperfect reference volumes obtained from the unsorted CT images; then, in each iteration, thorax motion is estimated from the reference volumes and the free-breathing CT images; Based on the estimated motion, the breathing indices as well as the reference volumes are refined and fed into the next iteration. The algorithm terminates when two successive iterations attain the same sorted reference volumes. In 3 out of 5 patient studies, our method attained comparable image quality with that using external breathing signals. For the other two patient studies, where the external signals poorly reflected the internal motion, the proposed method significantly improved the sorted 4D CT volumes, albeit with greater computation time.

Keywords: 4DCT, internal motion indices, motion estimation, image registration

I. INTRODUCTION

There has been a significant trend in recent years toward development and utilization of computed tomography (CT) techniques that generate multiple volumetric reconstructions, each one descriptive of a state of breathing. To reconstruct such time-resolved CT volumes of patients that breathe freely during a scan, different scanning protocols using multi-slice CT scanners are being investigated [14]. Although subtle differences exist among those scanning protocols, their general ideas are the same, and can be described by an oversampling-sorting process. Oversampling here means that at each table position the X-ray gantry continuously rotates for one to two breathing cycles. Multiple CT slices are reconstructed from the acquired projection data at each position. Temporally coherent CT slices are then sorted and stacked to form “4D” CT volumes, as illustrated in Fig. 1. The sorting process usually depends on external breathing signals that are recorded synchronously with the scan by some motion monitoring system. The recorded breathing signals may reflect the skin motion [2, 5], the skin tension [6] or the tidal volume measured orally[1]. Such external breathing indices may not always accurately represent the internal motion status [7, 8]. Using an external breathing signal that poorly correlates with the actual thoracic motion may cause severe tissue mismatches in the retrospectively sorted CT volumes. Note that 4D CT imaging can also be based on cone-beam CT scanners. Those scanners have large flat-panel detectors, leading to more consistent volume sampling, but they usually rotate very slowly (1 min per rotation) because they are mounted on linac systems. Methods have been proposed to use such scanners to image the thorax under free breathing conditions [912]. In this paper we mainly focus on 4D multi-slice CT imaging.

FIG 1
Illustration of 4D CT imaging methods.

The sorting requires a good respiratory-motion-correlated signal. Sometimes external breathing signals are not available or poorly correlate with the actual internal anatomy motion. Existing methods for extracting breathing indices directly from the image include tracking the center of mass (COM) [13, 14], correlating of a region of interest (ROI) between adjacent slices in consecutive table positions [3], calculating the internal air content [15] or estimating the diaphragm’s superior-inferior (SI) position in the cone-beam projection views [9, 16]. In 4D multi-slice CT images, tracking the diaphragm transition is inapplicable because the structure is present only in slices near the bottom of the thorax. The COM or correlation metrics may help identify the phase in one breathing cycle of the acquired images. However, respiratory motion usually varies in amplitude, duration and shape from cycle to cycle, so the reproducibility of the motion with respect to phase may be poor, unlike the more reproducible cardiac motion. It has been reported that phase sorting often results in more artifacts in the stacked thorax CT volumes than amplitude sorting does due to weak reproducibility of breathing motion [1719]. However, the amplitudes of COM or correlation signals are incomparable between different volumetric structures. Therefore, they may be insufficient to facilitate amplitude sorting for 4D multi-slice CT imaging.

We believe a need exists to have a system to improve on external or simple internal sorting for cases in which these methods yield unacceptable artifacts. In this paper, we develop a retrospective method that finds breathing indices for all slices based on internal anatomy motion, without using any externally recorded surrogates of breathing motion. The estimated motion indices are not real-time. However, the purpose of this work is to build 4D CT volumes with fewer artifacts to facilitate more accurate treatment planning. Details of this method are described next followed by the experiment results on five patient studies.

II. MATERIAL AND METHODS

A. Data acquisition

A General Electric (GE) 8-slice Lightspeed CT scanner was used to acquire CT data. The gantry speed was 0.5 second per rotation. The scanner was operated in an axial cine mode. In this mode, the scanner continuously scans the patient at one position for a certain time interval, then the x-ray beam is automatically turned off and the table moves to the next position, where the CT scanner resumes another round of continuous scan. This process repeats to cover the whole predetermined length of the body. Usually the cine duration (the scan time at each table position) is set to the maximum observed breathing period of the patient plus the overhead time of a full gantry rotation, to ensure that the scan covers an entire respiratory cycle. Patients’ breathing cycle may vary from 3 to 6 seconds. During the scan, a respiratory signal was synchronously recorded by a Varian real-time position management (RPM) system, which tracked the motion of a marker placed on the abdominal surface. Note that the RPM signals are not necessary for our internal motion based sorting method. However, RPM signals were recorded for comparison with our estimated internal motion indices.

B. Iterative sorting method

For each table position, multiple 8-slice CT volumes are reconstructed that temporally sample the moving anatomy for at least one full breathing cycle. We call the 8-slice CT volumes in each table scan period a group of free-breathing sub-volumes. Assume N table positions are scanned, then the whole data set contains N groups of free-breathing sub-volumes, denoted fn,k(x), x [set membership] R3, n = 1, (...), N, k = 1, (...), K, where n is the index for scan position, k is the temporal index and K is the total number of temporal samples in each table scan period. The nth group of subvolumes have an axial range of znd/2 ≤ zzn + d/2, where d is the axial aperture of the scanner at isocenter, and zn is the central axial coordinate of the nth group sub-volumes. To obtain 4D CT volumes from the unsorted images, we estimate internal motion that can provide sorting indices for all the images.

Here is an overview of the iterative internal motion-based sorting method.

Step 0. Extract simple preliminary breathing indices based on centroids of the unsorted images.

Step 1. Sort two reference volumes at relatively deep exhale and inhale states, denoted fex(x) and fin(x), x [set membership] R3, using the simple breathing indices.

Step 2. Find the full deformation during inhalation, Dfull, by registering the two reference volumes of the patient, fex(x) and fin(x).

Step 3. For each table position, estimate internal motion indices by iteratively updating a motion model to best match the deformed reference volume fex(x) to each group of moving sub-volumes.

Step 4. Sort two CT volumes at near end-of-exhale and end-of-inhale states based on the normalized internal motion indices.

If they differ from the previous reference volumes fex(x) and fin(x), replace them and go back to Step 2; otherwise, go to the final step 5

Step 5. Do a final amplitude-based sorting to form 4D CT volumes.

A registration step is necessary before Step 3 because the motion model we will describe later is built upon the extreme deformation during inhalation. To sum up, this method starts from two imperfect reference volumes, refines the internal motion-based indices to sort two reference volumes that contain fewer artifacts, and perform another round of motion estimation until two successive iterations attain same inhale and exhale reference volumes. Typically two to three rounds are sufficient for this process to converge. We now describe each step in detail.

1. Step 0 and 1: extract a simple breathing index and sort two initial reference volumes

We need inhale and exhale reference volumes of the patient to characterize the overall patient deformation during inhalation. However, such volumes are not available directly from the acquired CT images. In Step 0 we extract a simple breathing index for each sub-volume and in Step 1 we sort out two reference volumes using this simple breathing index. We treat the y-axis (anterior-posterior (AP) direction) centroid of the 8th slice of each sub-volume as the initial breathing index for that sub-volume. It is calculated as follows,

cn,k=jiyifn,k(xi,yi,8))jifn,k(xi,yi,8),
(1)

where cn,k denotes the initial breathing index of fn,k(·). We then stack the sub-volumes having largest or smallest index values at all positions to form two CT volumes. Other slices may also be used for calculating the initial breathing index. However, we need to determine whether the index with a peak value corresponds end-inhale or end-exhale states and using a boundary slice can save us effort in resolving this relationship, which is automatically decided based on the following property. The Chest always expands during inhale. Assuming n = 1 denotes the most superior sub-volume (close to the neck) and the AP coordinates are labeled from anterior to posterior, as shown in Fig. 2, then the peak of cn occurs near the end-exhale state and the valley occurs near the end-inhale state. However, the abdomen may expand or contract during inhale. To determine the states of the peak indices for the inferior sub-volumes, we examine the correlation coefficient (CC) of the y-axis centroids of the 1st and 8th slice of the moving sub-volumes. If the two sequences are poorly correlated (e.g., CC < 0.6), we will treat the peak of cn at the same state of the valley of cn−1, assuming the 1st slice of the nth sub-volumes is adjacent to the 8th slice of the (n − 1)th sub-volumes. So starting from the 1st sub-volume, the breathing states of the peak and valley of the initial breathing indices can be decided and propagated through subsequent sections.

FIG 2
A slice of the acquired CT image.

2. Step 2. Registration

In step 2 we register the two reference CT volumes to estimate the extreme deformation during inhale. This extreme deformation is needed in the next step for building the motion model. Various image registration methods have been developed in recent years [2023]. We use a B-spline based image registration method [24], but any other method that has been successfully applied to medical image registration can be used here.

The B-spline based deformation model is represented as follows,

Dfull(x;θ)=iθiβ(xxiΔx),
(2)

where θ are the B-spline knot coefficients of the deformation field, β(x) denotes the tensor product of cubic B-spline functions, xi is the ith spatial knot location, and Δx controls the width of the B-spline functions. Wide B-spline functions tend to capture more global deformations, but poorly represent local deformations. Narrow B-spline functions better describe local deformations, but they have more unknown parameters, complicating the optimization. We recommend placing two B-spline knots along the axial axis in each sub-volume region to capture more local deformations. For example, for an 8-slice CT with a slice thickness of 2.5 mm, we set the knot spacing in the axial direction to be 4 pixels, corresponding to 1 cm. We used 1.6 cm knot spacings in the left-right (LR) and AP directions

During registration, we deform the exhale reference volume fex(x) to match the inhale reference volume fin(x). We use a sum of squared differences (SSD) similarity term. Because the estimation problem is ill-posed, it is necessary to include some regularization on the estimated deformations. Since the cubic B-spline based deformations are intrinsically smooth, a smoothness regularization here is superfluous. We also would like to have the deformation be invertible, or at least locally invertible, so we adopt the regularization used in [25, 26] that encourages local invertibility by penalizing non-positive values of the Jacobian determinants of the transformation x + Dfull(x; θ). The cost function thus contains a sum of squared differences (SSD) similarity term and a Jacobian penalty term:

θ^=argminθ(x(fin(x)fex(x+Dfull(x;θ)))2+γR(J(;θ))),
(3)

where γ controls the tradeoff between the fidelity term and the penalty term and J(x; θ) denotes the Jacobian determinant of x + Dfull(x; θ). The penalty function R(.) is calculated as follows:

R(J(·;θ))=xg(J(x;θ)),g(J)={12(J0.05)2,J<0.050,     J0.05.
(4)

This penalty term helps pull the search directions of the estimator away from deformations having negative Jacobians.

We use the gradient descent algorithm to search for the parameter values that minimize the cost function. The multi-resolution technique is also applied in the optimization to avoid local minimum problems [27]. We start the registration from the downsampled-by-2 images and then use the results to initialize the registration of images with finer resolutions.

3. Step 3. Estimate the internal motion indices

Step 3 is the key task of this iterative sorting process, in which we find the internal motion-based breathing indices for all sub-volumes, denoted αn,k, n = 1, (...), N, k = 1, (...), K. We estimate each vector αn = (αn,1, (...), αn,K) from the exhale reference volume fex(x) and the free-breathing sub-volumes fn,k(x). The estimators iteratively update the parameters of a motion model to best match the reference volume to a sequence of moving sub-volumes. The motion parameters essentially represent the internal motion indices.

Now we describe the estimator. We start from the motion model and then explain the cost function that the estimator minimizes.

The motion model we use is based on a motion proportionality assumption [28]. We assume the displacement of each voxel at any time is proportional to its full movement from end-exhalation to end-inhalation. Based on this assumption, we express the motion model during the scan time of the nth sub-volume as follows:

Tn(x;t)=x+αn(t)Dfull(x),
(5)

where αn(t) denotes the scalar proportionality sequence. We allow αn(t) to be negative or greater than 1 because the estimated deformation Dfull(x) may not be the extreme displacement of each voxel during inhale. This proportionality motion model is imperfect. However, the ultimate goal for this work is not to find the precise motion of every voxel. We use the proportionality sequence αn(t) as a breathing index for the internal motion.

For simplicity, we parameterize the continuous proportionality function αn(t) using rect basis functions as follows:

αn(t)=k=1Kαn,krect(tk).
(6)

Based on the motion model (5) and (6), we optimize each vector of proportionality parameters αn by dynamically deforming the reference volume fex(x) according to the motion model to best match the nth group of free-breathing sub-volumes fn,k(x). We again use the SSD as the image matching criteria because the images are within the same modality. Respiratory motion usually changes slowly during free breathing, except in unusual like sneezing or coughing. Because the rect basis function does not ensure smoothness itself, we include temporal regularization into the estimator. To estimate αn, we use a cost function Ψn that contains two terms as follows:

α^n=argminαnΨn(αn),    Ψn(αn)=k=1KLn(fn,k(),fn,k(;αn,k))+λR(αn),
(7)

where the deformed reference volume at time tk is

fn,k(x;αn,k)=fex(x+αn,kDfull(x)),

and Ln(·) denotes the data fidelity term measured by the SSD over the field of view (FOV) of the nth sub-volume:

Ln(fn,k(),fn,k(;αn,k))=xFOVn(fn,k(x)fn,k(x;αn,k))2.

The roughness penalty term R(·) discourages rapidly changing motion estimates, and is measured by the differences of neighboring points as follows:

R(αn)=k=1K1(αn,k+1αn,k)2.

The scalar λ controls the tradeoff between the data fidelity and regularization terms.

Because the similarity is calculated only within the FOV of the sub-volume, it is not necessary to deform the whole reference volume when evaluating the cost function (7). To save computation time, we warp only the volume in the region that is slightly larger than the FOV of the sub-volumes. For example, if the extreme deformation of the nth sub-volume has a maximum absolute value of r along the axial direction, we deform the reference volume only within the axial range of znd/2 − rzzn + d/2 + r, because we assume the changes outside that range have negligible effect on calculating Ln. We use the conjugate gradient method [29] to minimize the cost function (7).

4. Step 4: Update the two reference volumes

In this step, we use the estimated {αn, n = 1, (...) , N} to form two CT volumes at near end-of-exhale and end-of-inhale states, which will be used in the next iteration of motion estimation. We first normalize each sequence αn as follows,

sn=αnαnminαnmaxαnmin×100%,n=1,,N,
(8)

where αnmax and αnomin denote the minimum and maximum values of αn over k. We then simply stack those sub-volumes having breathing indices that are closest to 80% as an inhale reference volume, and closest to 0% as an exhale volume. We use the CT volume at 80% inhale state because it appeared more consistent than the “full” inhalation state for different breathing cycles.

We normalize each breathing index sequence αn before sorting the exhale and inhale CT volumes to compensate for mismatch in the imperfect reference volumes. This normalization helps stabilize the whole process, i.e., it will find similar final motion indices even if started from different initial reference volumes. Fig. 3 illustrates an object that expands along one direction. Case 1 and case 2 select time t3 and t2 as the reference inhale state respectively. The estimated motion indices for these two cases will be quite different. However, the normalization equalizes them. Thus the estimated final breathing indices can be robust to imperfect initial reference volumes. Due to these normalizations, two rounds of motion estimations are generally sufficient to arrive at the final motion indices. Although this process is robust to imperfect reference volumes, extremely “bad” initial reference volumes will still fail with this algorithm. The bottom line is that each sub-volume stacked in the inhale reference volume should be at a deeper inhale state than that in the exhale reference volume.

FIG 3
Illustration that the normalization (8) of each estimated breathing index sequence can improve robustness to imperfect reference volumes.

5. Step 5: Final Sorting

Using the final internal motion indices associated with all the free-breathing sub-volumes at all table positions, we divide the sub-volumes into several breathing state bins to form 4D CT volumes. Various binning methods have been proposed [1, 2, 9, 18, 30], but they mainly fall into two categories: amplitude-based sorting and phase-based sorting. Amplitude-based sorting methods bin the data using the values and directions (inhale or exhale) of the breathing indices. Phase-based sorting methods usually bin the data according to phase-angles determined by some transformations on the breathing signal from the time domain to a phase domain. Phase-based sorting results in larger mismatch because of the insufficient motion reproducibility with respect to phase [18]. We chose amplitude sorting for this application. We first normalize the breathing signals using the following formula,

sn,k=αn,kαminαmaxαmin×100%,
(9)

where αmin and αmax denote the minimum and maximum values of α over both n and k. This normalization still preserves the relative amplitude variations among different breathing cycles. We then define a set of breathing state bins, each with an amplitude value in [0%, 100%] and a direction called ascending or descending. The amplitude value represents how deep the breathing is, and the direction indicates whether it is during inhalation or exhalation. We associate to each state bin the sub-volume whose corresponding breathing index is the closest to the bin amplitude value and that has the correct bin direction. We determine the direction of each motion index αn,k by comparing its value with its neighboring points αn,k−1 and αn,k+1. For the experimental results that follow, we specified 11 breathing states with these bin values:

bin.value=[1,0.8,0.6,0.4,0.2,0,0.2,0.4,0.6,0.8,1]×ratioratio=median(α1max,,αΝmax),
(10)

and these bin directions:

bin.direction=[1,1,1,1,1,1,0,0,0,0,0,],

where 1 and 0 represent “descending” and “ascending” respectively. Under this assignment, bin[1], bin[6] and bin[11] correspond to start of exhale, end of exhale and end of inhale respectively. Because the bin values are adjusted by the ratio defined in (10), offsets between the actual states using different kinds of breathing indices are minimized, hence the comparisons we report later between the internal motion indices and the external indices on sorting the images can be fair.

III. RESULTS AND DISCUSSION

A. Experiment results

We applied the iterative sorting method to five sets of 4D CT patient data. The CT images had a resolution of 0.98 cm × 0.98 cm and a slice thickness of 2.5 mm. The data was obtained from an existing ”standard” 4D CT imaging protocol. The X-ray tube current was 100 mA. The cine durations were about 6 seconds. Since the gantry speed was 0.5 second per rotation, a new set of 180 degree projections became available every 0.25 second for image reconstruction. This yielded about 20 distinct sub-volumes per table position.

Fig. 4 shows examples of the extracted centroid-based simple breathing indices for several positions of one patient. The ranges of the curves of the superior positions are smaller because the motion of the upper thorax during breathing is less. However, clear ascending and descending trends exist. The smallest and largest values in these breathing signals occur either near the end-inhalation or end-exhalation state. We stacked those sub-volumes associated with peak or base indices to form two initial reference CT volumes.

FIG 4
The extracted simple breathing signals of patient 1 based on centroid tracking. Positions 4 is near the neck and position 15 is near the abdomen.

Using these two reference volumes, the rest of the procedure described in Sect II was implemented. For Step 2, we set the regularization parameter γ in (3) to be 8×105 . This parameter was set based on our previous registration experiments [25, 26]. For Step 3, we set the regularization parameter λ to be 10. It was manually tuned using the first data set. We started λ from 1 and scaled it up by 5 until there were no abrupt bumps in the estimated proportionality sequence. For both registration and motion estimation, the reference volumes were downsampled by 2 in the transaxial plane to save computation time. The downsampled images provide enough information because of the small number of parameters in the proportionality motion model (5). Using a finer resolution would require significantly longer computation time but with very minor improvements if any.

We experimented with five patient data sets. In two of them, the RPM signals correlated with the patients’ breathing motion sufficiently to have minimized sorting artifacts and our method attained comparable image quality. The other three data sets (labeled patient 1,2 and 3) are further described in this section.

For patients 1 and 2, the RPM signals resulted in severe mismatches in the 4D CT volumes, as seen in Fig. 5b and Fig. 6b. We compared the estimated internal motion indices and the RPM breathing indices for those two patients in Fig. 7 and Fig. 8 respectively. In both plots these two sequences show similar transitions of inhalation and exhalation, but there also exist considerable discrepancies between them. Sorting the slices using the internal motion indices resulted in significantly better tissue consistency at the sub-volume boundaries. As shown in Fig. 5 and Fig. 6, most of these artifacts such as “flatness” and discontinuities in the RPM-sorted CT volumes were corrected.

FIG 5
Sorted CT volumes of patient 1 using recorded RPM indices (a) and internal motion indices (b). From upper left to lower right, the patient exhale and then inhale. Severe tissue mismatches are marked by arrows.
FIG 6
Sorted CT volumes of patient 2 using recorded RPM indices (a) and internal motion indices (b). From upper left to lower right, the patient exhale and then inhale. Severe tissue mismatches are marked by arrows.
FIG 7
The estimated internal motion breathing signals (’+’) and the recorded external RPM signals (’.’). Both signals were normalized according to (9). Each piece of curve represents the breathing signal for one scan position. ...
FIG 8
The estimated internal motion breathing signals (’+’) and the recorded external RPM signals (’.’) for patient 2. Both signals were normalized according to (9). Each piece of curve represents the breathing signal for one ...

For patient 3, neither the RPM indices nor the internal motion indices could find acceptable deep-inhale CT volumes. The mismatch occurred at the boundary between the 7th and 8th groups of sub-volumes, as shown in Fig. 9. We did a manual stack of the sub-volumes between the 7th and 8th groups. We found that no data in the 8th group could match the relatively deep-inhale sub-volumes in the 7th group because the patient breathed much more shallowly when acquiring the scans of the 8th group. So this failure is mainly due to data insufficiency rather than the proposed sorting algorithm.

FIG 9
End-of-inhale CT volumes of patient 3 using recorded RPM indices (a) and internal motion indices (b). Both contained mismatch at the bronchi.

B. Discussion

This paper describes an algorithm that forms 4D CT volumes by sorting using internal motion estimates. This 4D CT method does not require external breathing signals that may be inaccurate in measuring the actual overall tissue movement. As shown in the two “bad” examples (patient 1 and 2) presented above, the inaccuracy of the external RPM signals led to significant mismatches in the sorted CT volumes. The proposed internal motion based sorting method led to better consistency for both examples.

Our current implementation of this algorithm spent about 40 min on Step 2 (registration) and 20 min on Step 3 (motion estimation) on a single Pentium 3 GHz processor respectively for the first iteration. For the later iterations, registration and motion estimation can start from the results of the previous iteration hence requiring less computation time. When we strictly followed the termination condition, i.e., two successive iterations find the same reference volumes, four experiments needed three iterations and one needed two iterations, indicating that the second iteration already arrives at the final motion indices for most cases. This agrees with the conclusion we made in Step 4 that two iterations are sufficient for this iterative sorting process. The computation time is still long for routine clinical usage. However, in cases when the recorded breathing signals fail to produce acceptable 4D CT volumes, it may be worthwhile to use this method to sort the CT slices rather than rescan the patient, which would involve more X-ray exposure and still may not guarantee a successful reconstruction. From this point of view, the longer computation of the proposed method may be clinically acceptable. Certainly more effort will be put in reducing the computation time such as using a faster image registration algorithm, and parallelization could accelerate it significantly.

Our internal motion estimation is based on a proportionality motion model (5). This one dimensional, spatially-linear motion model cannot describe perfectly the actual trajectory of any point during breathing. However, the proportionality parameter can summarize the “average” deformation of all points, and hence can be a useful index. Although there are missing data at the boundary slices for CT image registration, we expect its effect on the accuracy of the estimated motion indices to be small because the motion estimator uses very few parameters. A higher-dimensionality motion model could better describe the actual 3D thorax motion. However, for sorting purposes, using lower dimensionality can simplify establishing the motion correspondences. From this point of view, we may also think that the proportionality motion represents the projection of the actual higher-dimensional motion onto a lower one-dimensional linear space to facilitate sorting. Similarly, after Xu et al. [19] obtained the deformation of each voxel by registering each slice to a high-resolution breath-hold reference volume, they still needed extra signal processing to generate a 1D breathing signal from the high-dimensional deformation field to enable a convenient motion phase synchronization among the CT images. The processing included an averaging of the deformations of all voxels and a principle component analysis on the 3D vector obtained from averaging. Because their 1D breathing motion signal could not facilitate an amplitude sorting, their final 4D CT volumes are generated by deforming the reference volume according to deformation fields that were smoothed at the sub-volume boundaries. So for sorting-based 4D CT imaging techniques, simple but reasonably descriptive motion models like (5) seem sufficient.

Our method improves the sorting of 4D CT images obtained from an existing ”standard” imaging protocol. It would be an interesting future study to investigate how this new method performs for alternative scan protocols, including helical scan and lower dose studies. Although ultra-low-dose protocols may suffice for motion estimation, may not be relevant clinically because the 4D CT scans are used for treatment planning and the image quality should be sufficient for accurate delineation of critical structures.

IV. CONCLUSION

We developed a method to sort free-breathing multi-slice CT images according to internal anatomy motion. Our method completely eliminates the reliance on any externally recorded surrogate of breathing motion. Patient studies showed that the internal motion indices estimated by our method resulted in better or at least comparable consistency in the stacked thorax CT volumes. Longer computation is required for this method because of image registration and motion estimation, hence we suggest it as a backup solution when the external surrogates of breathing motion result in poor sorting.

ACKNOWLEDGMENTS

This work was supported in part by NIH grant P01 CA59827-06A1.

Contributor Information

Rongping Zeng, EECS Dept., University of Michigan, Ann Arbor, MI 48109-2122.

Jeffrey A. Fessler, EECS Dept., University of Michigan, Ann Arbor, MI 48109-2122.

James M. Balter, RadOnc Dept., University of Michigan, Ann Arbor, MI 48109-2122.

Peter A. Balter, The University of Texas M.D. Anderson Cancer Center, Houston, TX 77030.

References

1. Low DA, Nystrom M, Kalinin E, Parikh P, Dempsey JF, Bradley JD, Mutic S, Wahab SH, Islam T, Christensen G, Politte DG, Whiting BR. A method for the reconstruction of four-dimensional synchronized CT scans acquired during free breathing. Med. Phys. 2003 June;30(6):1254–1263. [PubMed]
2. Rietzel E, Pan T, Chen GTY. Four-dimensional computed tomography: Image formation and clinical protocol. Med. Phys. 2005 April;32(4):874–889. [PubMed]
3. Pan T, Lee T-Y, Rietzel E, Chen GTY. 4D-CT imaging of a volume influenced by respiratory motion on multi-slice CT. Med. Phys. 2004 February;31(2):333–340. [PubMed]
4. El Naqa IM, Low DA, Christensen GE, Parikh PJ, Song JH, Nystrom MM, Lu W, Deasy JO, Hubenschmidt JP, Wahab SH, Mutic S, Singh AK, Bradley JD. Automated 4D lung computed tomography reconstruction during free breathing for conformal radiation therapy. Proc. SPIE 5369, Medical Imaging 2004: Physiology, Function, and Structure from Medical Images. 2004:100–106.
5. Vedam SS, Keall PJ, Kini VR, Mostafavi H, Shukla HP, Mohan R. Acquiring a four-dimensional computed tomography dataset using an external respiratory signal. Phys. Med. Biol. 2003 January;48(1):45–62. [PubMed]
6. Kleshneva T, Muzik J, Alber M. An algorithm for automatic determination of the respiratory phases in four-dimensional computed tomography. Phys. Med. Biol. 2006 August;51(16):N269–N276. [PubMed]
7. Yan H, Yin F-F, Zhu G-P, Ajlouni M, Kim JH. The correlation evaluation of a tumor tracking system using multiple external markers. Med. Phys. 2006 November;33(11):4073–4084. [PubMed]
8. Tang J, Dieterich S, Cleary KR. Respiratory motion tracking of skin and liver in swine for Cyberknife motion compensation. spie-5367. 2004:729–734.
9. Sonke J-J, Zijp L, Remeijer P, Herk M. Respiratory correlated cone beam CT. Med. Phys. 2005 April;32(4):1176–1186. [PubMed]
10. Zeng R, Fessler JA, Balter JM. Estimating 3-D respiratory motion from orbiting views by tomographic image registration. IEEE Trans. Med. Imag. 2007 February;26(2):153–163. [PMC free article] [PubMed]
11. Zeng R, Fessler JA, Balter JM. Respiratory motion estimation from slowly rotating X-ray projections: Theory and simulation. Med. Phys. 2005 April;32(4):984–991. [PubMed]
12. Li T, Schreibmann E, Yang Y, Xing L. Motion correction for improved target localization with on-board cone-beam computed tomography. Phys. Med. Biol. 2006 January;51(2):253–268. [PubMed]
13. Kachelriess M, Sennst D-A, Maxlmoser W, Kalender WA. Kymogram detection and kymogram-correlated image reconstruction from subsecond spiral computed tomography scans of the heart. Med. Phys. 2002 July;29(7):1489–1503. [PubMed]
14. Larson AC, White RD, Laub G, McVeigh ER, Li D, Simonetti OP. Self-gated cardiac cine MRI. Mag. Res. Med. 2004 January;51(1):93–102. [PMC free article] [PubMed]
15. Lu W, Parikh PJ, Naqa IME, Nystrom MM, Hubenschmidt JP, Wahab SH, Mutic S, Singh AK, Christensen GE, Bradley JD, Low DA. Quantitation of the reconstruction quality of a four-dimensional computed tomography process for lung cancer patients. Med. Phys. 2005 April;32(4):890–901. [PubMed]
16. Zijp L, Sonke J-J, Herk M. Extraction of the respiratory signal from sequential thorax cone-beam X-ray images. Proc. Intl. Congr. Comp. in Radiotherapy. 2004
17. Vedam SS, Keall PJ, Kini VR, Mohan R. Determining parameters for respiration-gated radiotherapy. Med. Phys. 2001 October;28(10):2139–2146. [PubMed]
18. Lu W, Parikh PJ, Hubenschmidt JP, Bradley JD, Low DA. A comparison between amplitude sorting and phase-angle sorting using external respiratory measurement for 4D CT. Med. Phys. 2006 August;33(8):2964–2974. [PubMed]
19. Xu S, Taylor RH, Fichtinger G, Cleary K. Lung deformation estimation and four-dimensional CT lung reconstruction. Acad Radiol. 2006 September;13(9):1082–1092. [PubMed]
20. Toga AW, Banerjee PK. Registration revisted. J. Neurosci. Meth. 1993 June;48(1):1–13. [PubMed]
21. Maintz JBA, Viergever MA. A survey of medical image registration. Med. Im. Anal. 1998 March;2(1):1–36. [PubMed]
22. Hill DLG, Batchelor PG, Holden M, Hawkes DJ. Medical image registration. Phys. Med. Biol. 2001 March;46(3):R1–R47. [PubMed]
23. Kybic J, Unser M. Fast parametric elastic image registration. IEEE Trans. Im. Proc. 2003 November;12(11):1427–1442. [PubMed]
24. Unser M. Splines: A perfect fit for signal and image processing. spmag. 1999 November;16(6):22–38.
25. Kim J, Fessler JA. Nonrigid image registration using constrained optimization. SIAM Conf. Imaging Sci., Abstract Book. 2004
26. Kim J. PhD thesis. Ann Arbor, MI: Univ. of Michigan; 2004. Intensity based image registration using robust similarity measure and constrained optimization: applications for radiation therapy. 48109-2122, Ann Arbor, MI.
27. Thevenaz P, Unser M. Optimization of mutual information for multiresolution image registration. IEEE Trans. Im. Proc. 2000 December;9(12):2083–2099. [PubMed]
28. Zeng R, Fessler JA, Balter JM. A simplified motion model for estimating respiratory motion from orbiting views. Proc. SPIE 6512, Medical Imaging 2007: Image Proc. 2007:1–8. pages 651240. [PMC free article] [PubMed]
29. Press WH, Flannery BP, Teukolsky SA, Vetterling WT. Numerical recipes in C. New York: Cambridge Univ. Press; 1988.
30. Keall PJ, Starkschall G, Shukla H, Forster KM, Ortiz V, Stevens CW, Vedam SS, George R, Guerrero T, Mohan R. Acquiring 4D thoracic CT scans using a multislice helical method. Phys. Med. Biol. 2004 May;49(10):2053–2068. [PubMed]