A prospective gating method to acquire a diverse set of free-breathing CT images for model-based 4DCT

Four dimensional computed tomography (4DCT), also referred to as respiratory-correlated CT, utilizes multiple projections of the same tissues acquired at different points in the breathing cycle to reconstruct a series of breathing-gated volumetric images (Lu et al 2006). 4DCT was developed to improve the accuracy of therapeutic radiation delivery to tumors in the lung and abdomen, which can be displaced by as much as 3 cm during treatment due to breathing motion (Barnes et al 2001, Keall et al 2006).

4DCT has seen widespread adoption in the field of radiation oncology as a tool to manage breathing motion in treatment planning; implementations are offered by the major scanner manufacturers and 4DCT datasets are recognized by many treatment planning systems. It has also shown promise as a functional imaging technique for measuring pulmonary ventilation (Guerrero et al 2006, Vinogradskiy et al 2013, Brennen et al 2015), and more recently as a means to estimate the elastic properties of lung tissue for bio-mechanical modeling (Hasse et al 2017). Much of the literature on ventilation imaging and lung elastography relates to the context of radiation therapy—both techniques are generally proposed as a means to identify healthy lung tissue for preferential sparing in treatment planning.

Outside of radiation therapy, 4DCT could potentially be employed to characterize lung function for the phenotyping and response monitoring of chronic obstructive pulmonary disease (COPD) (Yamamoto et al 2011, Murphy et al 2012), a progressive lung disease which is the third leading cause of death in the United States (Hoyert and Xu 2012). COPD is typically diagnosed by measuring the ratio of forced expiratory volume in one second to the full vital capacity (FEV1/FVC). FEV1/FVC is unable to provide information about the physiology causing impaired respiratory dynamics (Jögi et al 2011). Single photon emission computed tomography (SPECT) ventilation perfusion scintigraphy (V/Q) is the current gold standard for lung functional imaging, but 4DCT is less expensive, more widely available, and has superior spatial resolution.

In most cases, 4DCT effectively characterizes respiratory induced motion, however, commercial implementations suffer from two major limitations: susceptibility to image artifacts if the patient breathes irregularly during acquisition (Yamamoto et al 2008), and a high imaging dose—up to four times higher than standard thoracic CT (Matsuzaki et al 2013). Multiple strategies have been proposed to reduce 4DCT artifacts and imaging dose. Castillo et al demonstrated that oversampling led to the fewer artifacts in cine 4DCT, but imaging dose increased (Castillo et al 2015). In Tian et al (2011), the authors developed a temporal non-local means approach to reconstruct gated images using fewer projections and acquired with low mAs, however their method does not treat inter-cycle respiratory variability. Pan et al developed a prospective technique to mitigate artifacts in cine 4DCT by stopping irradiation when irregular breathing occurs and selectively re-acquiring images only at the necessary couch positions (Pan et al 2017). Although their method achieved reduced imaging dose compared to repeating an entire scan if severe artifacts occur, the current standard clinical practice, it does not provide a means to lower imaging dose below that of a typical 4DCT.

Respiratory motion model-based 4DCT techniques have been developed that mitigate or avoid the artifacts commonly seen in commercial 4DCT protocols (McClelland et al 2013). In contrast to commercial protocols, which generally function by retrospectively sorting projections into bins for independent reconstructions, model-based methods aim to describe the correspondence between internal tissue displacements, typically found using deformable image registration, and an externally measurable breathing surrogate signal using a motion model. Examples of previously published models include polynomial functions of one or more surrogate signals (Low et al 2005), B-splines (McClelland et al 2006), and principal components of deformation vector fields (Li et al 2011, Xu et al 2015). More recent developments in respiratory motion modeling include work by McClelland et al (2017), where the authors present a unified framework for image registration and model parameter estimation, and the kernel ridge regression approach proposed by Geimer et al (2017) that used abdominal surface motion as a multidimensional breathing surrogate signal.

Our group previously introduced a model-based technique, termed '5DCT', which accurately characterizes tumor motion in the presence of irregular breathing (Low et al 2013, Thomas et al 2014). Our technique has been demonstrated to be accurate to within 2 mm using a landmark analysis (O'Connell et al 2015), and to produce sorting artifact-free images which accurately describe tissue motion to within 2 mm on average in a 16 patient cohort (Dou et al 2015).

Model-based 4DCT techniques require images that capture different respiratory states. Unlike other techniques described in McClelland et al (2013) that use commercial 4DCT or cine protocols to obtain image data, 5DCT uses repeated fast helical scanning as the patient breathes freely to sample the respiratory cycle. Although the free breathing images do not correspond to a single respiratory phase as a commercial 4DCT reconstruction would, each individual axial slice can be assigned a breathing amplitude and rate measurement, given by the bellows signal and its derivative at that slice's acquisition time. Blurring and doubling artifacts are mitigated by using a fast helical protocol. Repeated scans are performed to collect observations throughout the breathing cycle. The previously published research protocol acquires 25 scans in succession. Each scan uses a low mA and fast rotation time technique such that the total dose of all 25 scans is similar to that of a commercial 4DCT scan. Although the current protocol ensures a high probability that the breathing cycle will be adequately sampled, redundant image data is acquired because the scans are not selected prospectively. In earlier work, we determined that using unsynchronized scanning required at least ten scans to guarantee acquisition of images with sufficient variation in the breathing phase to reliably generate accurate motion models (Thomas et al 2016). Lowering the number of free-breathing scans would allow for a reduction in imaging dose beyond what is currently achievable with commercial 4DCT. While Thomas et al demonstrated that the number of scans could be reduced without substantially increasing the model residual if the scans were acquired at sufficiently diverse respiratory phases, no method was proposed to select such a set of scans. Ruan et al introduced an objective function to quantify the phase coverage of a set of breathing surrogate signals measured during each scan, which will subsequently be referred to as respiratory trajectories (Ruan et al 2015). In this work, we present a novel adaptive method to prospectively select a desirable set of 5 to 10 trajectories by triggering image acquisition based on real-time analysis of the breathing surrogate signal using that objective function. We hypothesize that using the proposed method, the number of scans can be reliably reduced while maintaining acceptable model accuracy. The method was tested by simulating image acquisition using the measured breathing traces of five patients and evaluating the phase coverage of the sampled respiratory trajectories. Although the analysis in this work focuses on 5DCT, the proposed technique could be used in conjunction with any respiratory motion model that can utilize free-breathing CT image data.

We previously found that the phase coverage of a set of respiratory trajectories correlated with model accuracy only for five or more scans, and that ten scans yielded adequate sampling of the breathing cycle. Therefore, the goal of a prospective scanning approach would be to select between five and nine scans with a desirable spread of respiratory trajectories. We had two primary objectives for this study: firstly, to determine the minimum phase coverage for sets of N scans, where 5  ⩽  N  ⩽  9, necessary to achieve mean model accuracy within 0.5 mm of the corresponding 25 scan model; and secondly, to test a prospective protocol and determine the number of scans necessary to reach the previously determined minimum phase coverage. Sections 2.1–2.3 describes the patient datasets and the process of finding the phase coverage thresholds. Sections 2.4 and 2.5 describe the prospective method and the simulation used to validate it.

2.1. Image acquisition and motion modeling

Five lung cancer patients were imaged under the published 5DCT protocol, described in detail in Thomas et al (2014), using 64-slice CT scanners (Biograph TruePoint, Definition AS; Siemens Healthcare, Forchhiem, Germany). A low dose technique (40 effective mAs), short rotation time (0.33 or 0.285 s) and a pitch of 1.5 or 1.2 were used to scan each patient 25 times while breathing was simultaneously monitored with an abdominal bellows (Model 76513NM12G; Lafayette Instrument, Lafayette, IN), and a pressure to voltage transducer sampled with LabVIEW (National Instruments, Austin, TX) software. The acquisition protocol was designed such that tissues were imaged for a sufficiently short amount of time (approximately one third of a second) to treat the axial slices of the free-breathing images as measurements of instantaneous tissue positions. Each axial slice was assigned a breathing amplitude and rate according to the bellows signal at the time of acquisition. Importantly, because scans were acquired during free-breathing, the respiratory phase of the images varied by slice—rather than an image depicting a single point in the breathing cycle, as in a gated 4DCT reconstruction, the respiratory phase of each scan was described by a curve as shown in figure 1(b).

Zoom In Zoom Out Reset image size

The first scan, selected as the reference geometry, was deformably registered to the remaining 24 scans using the 'deeds' algorithm, proposed in Heinrich et al (2013) and demonstrated to achieve an average registration error of 1.43 mm over 10 cases using expert selected landmarks. The algorithm recovers large tissue displacements well, and models the sliding motion between the lungs and the rib cage. Additionally, it does not require a lung segmentation. The deformation vector fields (DVF) output by the DIR provided measurements of tissue displacement from the reference geometry to the remaining 24 observed breathing states. The DVFs and corresponding amplitude and rate values were used to solve for the voxel-specific parameters relating tissue motion to breathing amplitude, v, and rate, f, in the 5D motion model (Low et al 2005):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \vec{X}=\overrightarrow{{{X}_{0}}}+\vec{\alpha }v+\vec{\beta }f.\nonumber \end{align} \tag{ 1 }$$

Here, $\vec{X}$ denotes the tissue displacement, $\overrightarrow{{{X}_{0}}}$ its position in the reference geometry, $\vec{\alpha }$ relates inhalation depth to surrogate amplitude, and the parameter $\vec{\beta }$ relates hysteresis motion to the derivative of the surrogate signal.

For each patient, model parameters were solved for using all 25 DVFs to build a 'full' motion model. A cuBLAS batched QR factorization routine was used to find least-squares solutions for model parameters on a graphics processing unit (Tesla K40; NVIDIA, Santa Clara, CA). After fitting $\vec{\alpha }$ and $\vec{\beta }$, equation (1) allows for the computation of a DVF to transform the reference image to arbitrary amplitude v and rate f, producing a breathing-gated image at any user-selected breathing phase without sorting artifacts; no binning of projections is performed because the model parameters are fit using free-breathing fast helical images.

Model residual was calculated by comparing the tissue displacements computed using equation (1) to the 25 displacements given by DIR, and taking the root mean square of the difference in displacements for each voxel. An overall figure of merit was calculated by averaging the residual values for lung voxels where the measured motion was  ⩾5 mm. The motion criteria was included to prevent voxels which underwent little to no motion, such as regions of the upper lung, from biasing the accuracy measure. Importantly, the residual is representative of how well the motion model agrees with the DIR and does not account for uncertainty in the registration, ie differences between the measured displacements and the actual physical motion.

2.2. Quantifying the phase coverage of sampled respiratory trajectories

To faithfully characterize tissue motion with a reduced number of scans it is necessary to image each slice location's tissues at breathing states that are well-distributed throughout the respiratory cycle. We adopt the previously published phase coverage metric described in Ruan et al (2015):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M=\underset{z}{\mathop \sum }\,\underset{i=1}{\overset{N-1}{\mathop \sum }}\,{{\left( {{{\tilde{v}}}_{i+1,z}}-{{{\tilde{v}}}_{i,z}} \right)}^{0.5}}\nonumber \end{align} \tag{ 2 }$$

where $\tilde{v}$ is the vector of the N measured breathing amplitudes for a slice location z, shown as values lying on the dashed line in figure 1(b), sorted in ascending order.

2.3. Relationship between phase coverage and model accuracy

2.4. Prospective gating method

The proposed method consisted of three elements: an overall triggering scheme, real-time breathing phase estimation, and a decision on whether or not the patient is currently breathing in a reproducible manner. From previous work, it was observed that sets of trajectories exhibiting a high M value generally resembled a pattern of phase-shifted waveforms and therefore it was hypothesized that adequate phase coverage could be achieved by staggering the phase at which scans are started. A respiratory profiling method (Ruan et al 2009) was used to estimate the breathing phase at each timepoint. A 'representative breath' (White et al 2013) was formed based on training data and subsequently updated each new period. In order to quantify the reproducibility of the breathing pattern, the normalized cross-correlation between the representative breath and the most recent respiratory surrogate measurements was continually computed. An overview of the method is shown in figure 2.

Zoom In Zoom Out Reset image size

2.4.1. Breathing phase estimation

A detailed description of the technique used for real-time breathing phase estimation is given in Ruan et al (2009). In brief, in order to characterize hysteresis motion, an augmented state vector for each time-point t_i was formed by taking the current breathing amplitude, v(t_i) and a previous observation v(t_i  −  kd_s), where k is the lag expressed in number of samples and d_s is the sampling interval. An ellipse was continually fit to points in the augmented state space within a sliding window encompassing (v(t_i), v(t_i  −  kd_s)) and approximately one period of previous observations. Transition events marking the end of a period were detected using Poincaré sectioning. An estimated phase was assigned to each time-point by linearly extrapolating from the most recent marker event.

2.4.2. Reproducibility evaluation

For each patient, the first 30 s of the measured breathing trace was used as training data. The trace was segmented into individual breaths and the mean period was taken as the period of the representative breath. The amplitude was found by linearly interpolating the amplitude values at each time-point in the segmenting breaths and computing the weighted average, where the weights emphasized the most recently observed breaths. The resulting waveform was smoothed using Savitzky–Golay filtering (Savitzky and Golay 1964) to produce the representative breath. After each marker event, the representative breath was updated to include the most recent period. In order to avoid bias due to outliers or incorrectly detected periods, breaths which exceeded two standard deviations in amplitude or period from the mean were not used to update the representative breath.

The normalized cross-correlation (Lewis 1995) between the representative breath and the most recent respiratory surrogate signal values, used as the template, was computed. The length of the template was set to half of the breathing period. The correlation coefficient provided an indicator of breathing reproducibility at the current time point. If the coefficient exceeded a preset threshold a scan was eligible to be triggered if the breathing phase criterion was satisfied.

2.4.3. Triggering scheme

Sets of respiratory traces with high phase coverage were observed to generally resemble phase-shifted waveforms. A triggering scheme was designed to stagger the starting phase of each free-breathing scan with the aim of collecting a desirable set. The breathing cycle was divided uniformly into ten contiguous phase intervals, with each interval i representing 10% of the breathing period. Initially, all ten phase intervals are marked as 'unsampled'. No restrictions are placed on the phase of the first scan—acquisition is triggered as soon as the breathing reproducibility criterion is met. The interval containing the phase point at which the first scan was triggered is marked as 'sampled' and removed from consideration for triggering subsequent scans. Of the remaining nine intervals, a subset of candidate intervals is chosen to maximize the minimum distance (defined as the number of intervals), allowing for wrap-around because the representative breath is periodic. There can be more than one candidate interval—if there are multiple unsampled regions with the same minimum distance to a sampled region, then scan k can be triggered in any of them. An example case is shown in figure 3, where the four candidate intervals for the next acquisition are darkly shaded, the sampled intervals are unshaded, and the ineligible intervals are lightly shaded.

Zoom In Zoom Out Reset image size

Once the set of candidate intervals is determined, acquisition is triggered only when the reproducibility criterion is satisfied and the current phase is within a candidate interval. After five scans are acquired, the phase coverage metric M is computed. If M is satisfactory, acquisition is ended. Otherwise, scanning continues until a satisfactory M is obtained or ten scans acquired.

2.5. Simulation study

3.1. Phase coverage and model accuracy

The average full model residual over the five patients in this study was 1.06  ±  0.2 mm. The mean minimum M values for subsets of N scans which achieved a model residual within tolerance were 0.368  ±  0.020, 0.339  ±  0.020, 0.310  ±  0.020, 0.290  ±  0.020, 0.268  ±  0.022 for N  =  5, 6, 7, 8, and 9, respectively. These M values corresponded to the 91st, 65th, 39th, 8th and 1st percentiles on average among the M values of sets of N non-prospectively acquired scans.

Figure 4 displays the distribution of residual for reduced models for one patient. Note that for N  =  5 scans, few of the subsets are able to achieve a residual within tolerance. The proportion of subsets that would yield an acceptable model residual increases as the number of scans increases. When N  =  9, almost all subsets result in tissue motion being acceptably modeled.

Zoom In Zoom Out Reset image size

3.2. Simulation study

The average minimum M values for N  =  5, 6, 7, 8 and 9 were used as stopping criteria in the simulation study. For four out of the five patients, the prospective method was able to select a set of respiratory trajectories with an acceptable M value in five scans. For the remaining patient, acceptable phase coverage was achieved after six scans. Percentile scores describing the relative M value with respect to subsets of actually acquired scans for the five patients were 90.1, 100.0, 100.0, 99.1, and 98.6.

Figure 5 shows the M values for different numbers of prospectively selected scans for one patient, along with M values of subsets of the non-prospectively selected scans, and the minimum M value necessary for sufficient accuracy to provide context. When N  =  5,6,7, very few subsets of the non-prospective scans have an M value greater than the threshold. For these cases, the difference in phase coverage achieved by the prospectively selected scans compared to non-prospectively selected is most apparent.

Zoom In Zoom Out Reset image size

We have described a prospective method to acquire a set of free-breathing CT scans such that their associated respiratory surrogate signals are well-spread throughout the breathing cycle in order to facilitate the fitting of a correspondence model between the respiratory surrogate signal and tissue displacement for model-based 4DCT. The results presented here suggest that using the proposed approach, the number of scans required for accurately describing tissue motion with 5DCT can be reduced from the previously published 25 to as few as six with less than a 0.5 mm increase in mean residual. The proposed method can be employed to select free-breathing scans for any respiratory motion model.

The primary motivation behind the proposed method is to reduce the imaging dose associated with our model-based 4DCT technique, 5DCT, by eliminating acquisition of redundant data which does not improve modeling accuracy. In this study modeling accuracy was assessed by computing the residual, defined as the root mean square difference in tissue positions given by the model and image registration. While the residual is indicative of model performance, it does not translate directly to a measurement of error. A thorough assessment of error, including the effects of uncertainty in image registration and the relationship between the respiratory surrogate signal and tissue motion, would require a evaluation such as a landmark analysis.

Uncertainty in the correlation between the respiratory surrogate signal and internal tissue motion could potentially reduce the effectiveness of the prospective scan selection method and degrade the accuracy of the resulting motion model. A major complicating factor in estimating this correlation is respiratory hysteresis. In the real-time phase estimation technique employed in the proposed method, hysteresis is treated by increasing the dimensionality of the surrogate signal and defining each surrogate measurement by its current amplitude as well as the amplitude observed some fixed amount of time ago (Ruan et al 2008). Similarly, the 5D model accounts for hysteresis by describing tissue motion as a function of surrogate amplitude as well as its derivative.

In figure 1(b), the set of five scans with the highest M value for patient 5 are shown. The 5D model, equation (1), constrains each voxel's motion to lie within a plane with respect to amplitude and rate. The two planes in figure 1(c) demonstrate the difference in the full model parameters, fit using all 25, versus the reduced model, fit using the high M subset of 5. Although overall agreement is good, note that the inclusion of deeper breaths, which drive the M value to be higher, biases the reduced model toward estimating greater displacement. The redundant data incorporated into the full model fitting damps this effect.

We incorporated a delay of two seconds between a trigger signal and the beginning of data acquisition into our simulation. With the current scanner software at our institution (Syngo, Siemens Healthcare; Forchheim, Germany), there is a time lag of 2 s between the operator initiating a scan and acquisition. Future work will be centered upon implementing the method in real-time on a workstation separate from the CT console, and triggering will be accomplished by visually signaling to the operator when to manually initialize a scan. Ideally the proposed method could be incorporated into the scanner software to eliminate the delay, similar to commercially available cardiac gating schemes, but we chose to model a more feasible implementation in this study.

The reproducibility evaluation described in section 2.4.2 aims to avoid the effects of irregular breathing by preventing imaging under such conditions, but it has obvious limitations. An abrupt change in respiration after triggering but during or shortly before image acquisition could not be protected against. Such variance is unpredictable even for human operators. Although the model-based 4DCT technique used in this work is not susceptible to image artifacts due to irregular breathing, the proposed prospective selection method may require additional scans before the stopping criteria is met compared to a more regularly breathing case, and consequently the potential imaging dose savings would be limited.

Typically, binning-based 4DCT reconstruction methods require an irradiation period equal to or longer than the breathing period in order to reconstruct gated images. In contrast, model-based techniques construct a continuous description of tissue motion from a limited set of observations acquired at different respiratory states. For 5DCT, the model-based approach employed in this work, fast helical free-breathing scans are used as observations. An approximate volume computed tomography dose index (CTDI_vol) value for a single free-breathing scan for 5DCT is 2.7 mGy (McCollough et al 2008). Simply reducing the number of scans from 25 to six allows for the acquisition of 5DCT with less than 25% of the current dose. Further reduction could be achieved by combining the proposed technique with low-dose iterative CT reconstruction (Baumueller et al 2012). The substantial savings afforded by the proposed technique have the potential to reduce the total imaging dose of model-based 4DCT into the range where it is clinically acceptable to use as a tool for non-radiation therapy applications.

A method to prospectively select a set of free-breathing fast helical CT scans which are well-spread throughout the respiratory cycle has been described and tested in a five patient simulation study. Using the proposed technique, respiratory phase coverage which suggests model accuracy within 0.5 mm on average of a 25-scan model could be achieved in six or fewer scans.

The Tesla K40 used for this research was donated by the NVIDIA Corporation. This work was supported in part by funding from Siemens Healthcare.

None.