Rapid estimation of 4DCT motion-artifact severity based on 1D breathing-surrogate periodicity

Purpose: Motion artifacts are common in patient four-dimensional computed tomography (4DCT) images, leading to an ill-defined tumor volume with large variations for radiotherapy treatment and a poor foundation with low imaging fidelity for studying respiratory motion. The authors developed a method to estimate 4DCT image quality by establishing a correlation between the severity of motion artifacts in 4DCT images and the periodicity of the corresponding 1D respiratory waveform (1DRW) used for phase binning in 4DCT reconstruction. Methods: Discrete Fourier transformation (DFT) was applied to analyze 1DRW periodicity. The breathing periodicity index (BPI) was defined as the sum of the largest five Fourier coe ffi cients, ranging from 0 to 1. Distortional motion artifacts (excluding blurring) of cine-scan 4DCT at the junctions of adjacent couch positions around the diaphragm were classified in three categories: incomplete, overlapping, and duplicate anatomies. To quantify these artifacts, discontinuity of the diaphragm at the junctions was measured in distance and averaged along six directions in three orthogonal views. Artifacts per junction (APJ) across the entire diaphragm were calculated in each breathing phase and phase-averaged APJ, defined as motion-artifact severity (MAS), was obtained for each patient. To make MAS independent of patient-specific motion amplitude, two new MAS quantities were defined: MAS D is normalized to the maximum diaphragmatic displacement and MAS V is normalized to the mean diaphragmatic velocity (the breathing period was obtained from DFT analysis of 1DRW). Twenty-six patients’ free-breathing 4DCT images and corresponding 1DRW data were studied. Results: Higher APJ values were found around midventilation and full inhalation while the lowest APJ values were around full exhalation. The distribution of MAS is close to Poisson distribution with a mean of 2.2 mm. The BPI among the 26 patients was calculated with a value ranging from 0.25 to 0.93. The DFT calculation was within 3 s per 1DRW. Correlations were found between 1DRW periodicity and 4DCT artifact severity: − 0 . 71 for MAS D and − 0 . 73 for MAS V . A BPI greater than 0.85 in a 1DRW suggests minimal motion artifacts in the corresponding 4DCT images. Conclusions: The breathing periodicity index and motion-artifact severity index are introduced to assess the relationship between 1DRW and 4DCT. A correlation between 1DRW periodicity and 4DCT artifact severity has been established. The 1DRW periodicity provides a rapid means to estimate 4DCT image quality. The rapid 1DRW analysis and the correlative relationship can be applied prospectively to identify irregular breathers as candidates for breath coaching prior to 4DCT scan and retrospectively to select high-quality 4DCT images for clinical motion-management research.


INTRODUCTION
Four-dimensional computed tomography (4DCT) based on respiratory correlation is widely applied in radiation oncology to define the internal tumor volume (ITV) for thoracic and abdominal lesions. 1 The ITV approach is simple and effective in treating a moving target but often exposes normal tissue along the tumor motion path to the full prescribed dose. Tumor motion tracking or 4D radiotherapy (4DRT) 2,3 is a solution that spares normal tissue by focusing the radiation beam on the target volume. However, due to the presence of motion artifacts in 4DCT, the delineated gross tumor volumes in different breathing phases can vary by up to 90% 4 -110%, 5,6 especially for small lesions with large motion ranges. This variation may affect the defined ITV, particularly if the artifact occurs at the extreme respiratory phases, such as the full inhalation phase that is known for irreproducibility. 7 When considering 4DRT treatment planning, 6,8 scan artifacts in all respiratory phases are expected to have similar impact on the planning target volume, leading to greater treatment uncertainties.
Motion artifacts in clinical 4DCT images were previously reported as common phenomena since patients rarely breathe reproducibly. 7,[9][10][11][12] These artifacts cause anatomical distortion, resulting from breathing irregularities during 4DCT scans. In a cine scan, breathing irregularity occurring within one bed position breathing may cause blurring artifacts, since image projections from a narrow range of time points are binned into phases based on a respiratory motion surrogate. From one bed scan to the next, breathing irregularities cause discontinuity of the 4DCT image at the junction, leading to such artifacts as incomplete, overlapping, and duplicate anatomies, which are much pronounced than motion blurring. 13 A patient study consisting of 50 clinical 4DCT images of the thorax and abdomen showed that 90% of these images contained at least one motion artifact with a mean magnitude of 12 mm. 13 Because breathing irregularities produce most of motion artifacts using the binning method, the motion artifacts in respiratory-correlated 4DCT are often referred as binning artifacts. 14,15 As a result of recognizing the severity and frequency of motion artifacts and realizing the importance of artifact reduction for both current and future clinical needs, effort has been devoted to correcting or minimizing binning-associated artifacts. Such studies include the use of internal, feature-based, or multipoint motion surrogates for binning [16][17][18] or of using direction-dependent displacement or the motion velocity for binning. [19][20][21] Other studies correct motion artifacts by modeling respiratory motion with deformable image registration [22][23][24][25][26] or probability density functions. 15 More recently, advanced CT scanners, such as 64-slice fan-beam or 320-slice cone-beam CT scanners, 15,22 have been developed; these circumvent some binning artifacts with a shorter rotation time and a longitudinally widened 4D scan volume. However, these methods are not yet implemented in most clinics.
In this study, we report a correlative relationship between motion artifacts and breathing irregularities. Two new indexes are introduced to quantify the motion-artifact severity (MAS) in 4DCT images and the breathing periodicity index (BPI) of the 1D respiratory waveform (1DRW), which is the respiratory surrogate used for retrospective binning. Twentysix free-breathing 4DCT images and corresponding 1DRW curves were studied by manual measurement of the MAS and by using discrete Fourier transformation (DFT) method to calculate the BPI, respectively. The established correlation could be used as to anticipate image quality from a breathing trace acquired prior to 4DCT scanning or to further develop a gated 4DCT scanning to suspend image acquisition during periods of irregular breathing.

METHODS AND MATERIALS
Twenty-six patients who had received a 4DCT during quiet unprompted breathing at simulation for radiotherapy of thoracic and abdominal lesions were selected for this study. The images were acquired in cine mode using an 8-slice CT scanner (Discovery, GE Healthcare, Milwaukee, WI) and using real-time position management system (RPM, Varian Oncology Systems, Palo Alto, CA) as an 1DRW surrogate. These images were binned into ten phases using proprietary software (GE Advantage PET/CT Simulator). In this study, only the motion artifacts around the diaphragm are reported. Figure 1 shows the workflow diagram used to quantify motion-artifact severity in the phase-binned images, the periodicity of the RPM breathing trace, and their correlation.

2.A. Fourier analysis for 1DRW (RPM) periodicity
The time range of the RPM (1DRW) curve was selected from the entire breathing trace to correspond to the bed positions, in which the diaphragm was scanned in all breathing phases. The superior-inferior locations of the diaphragm within 4DCT field of view were used to select the corresponding RPM time range. Fourier analysis on the selected breathing trace was used to obtain the BPI. Mathematically, a 1DRW was considered as a 1D n-element array v, which was processed to remove its constant component by subtracting its mean: v = v − mean(v). After Fourier transformation, the Fourier coefficients of v were sorted by absolute magnitude in a descending order as F(i), i = 1, 2, ..., m. The BPI is defined as namely, it is defined as the ratio of the energy (F(i) 2 ) of the largest five Fourier components (BPI-5) to the total Fourier energy. In this study, k was defaulted to five and BPI was defaulted to BPI-5. BPI value ranges from 0 to 1. Based on Parseval's theorem, the zero-mean array v is conserved by Fourier transformation. Thus, BPI-5 quantifies how closely the first five Fourier components represent the original array v (1DRW). For a highly periodic function, the corresponding BPI is close to one. Conversely, the energies of a nonperiodic function are scattered over a far wider spectral domains and the associated BPI-5 should be small. Figure 2 illustrates the difference between a nonperiodic function (BPI-5 = 0.18) and a periodic sinusoidal function (BPI-5 = 0.95) with minor noise.
In addition to BPI-5, other BPI-k (k = 2 to 8) were evaluated, including Pearson correlation among themselves and the trend of BPI-k values as a function of k. BPI-5 was used to study the correlation with the MAS from the phase-binned 4DCT.
From the Fourier analysis, the mean breathing period (T) was calculated as below.
where the total number of samples (N), sampling rate ( f s ), and the index of the dominant Fourier coefficient (c F ) in Fourier frequency spectrum were used. The period was used to calculate the mean velocity of the diaphragm during respiration. The Fourier analysis program, which can process all 1DRW data files in batch mode, was implemented in  on a desktop computer (Windows XP, Core duo CPU E8500 at 3.1 GHz with 2 GB of RAM).

2.A.1. Estimation of mean diaphragmatic displacement and velocity
The mean diaphragm location was determined by using six pivot points on the diaphragm in each of the ten phases. 27 Three points on the right and left sides of the diaphragm were measured; these points were chosen to be the most inferior points at the anterior and posterior sides of the diaphragm, and the most superior point at apex point of the diaphragmatic dome (Fig. 3). The maximum diaphragm displacement between the full exhalation (FE) and full inhalation (FI) phases was calculated (D Max Diaphragm ) and the mean diaphragmatic velocity (v Diaphragm ) was the D Max Diaphragm divided by the mean period (T) where Z represents the superior-inferior position of the six points, FE is 50% phase and FI is averaged between 0% and 90% phases. These patient-specific motion parameters were used to "normalize" the MAS, which minimized patient specificity and allowed motion artifact comparability among patients.

2.B. Measurement of motion artifacts in 4DCT images
Motion artifacts in 4DCT were visually identified and manually quantified using an in-house treatment planning system (TPS) following the procedure in Fig. 1. In cine 4DCT images, motion artifacts that distort the diaphragmatic appearance at F. 3. Illustration of three pivot points on one lateral side of the diaphragm. Six points in total on both sides of the diaphragms were measured to define the mean diaphragm position. This image shows a respiratory phase (yellow) and the full exhalation (red) as the reference phase. Starting with respiratory phase 0% (full inhalation), all motion artifacts were measured in three orthogonal planar images. At a 2-bed junction in the axial view, a crosshair was placed at the center of the diaphragm. Within the sagittal and coronal views, corresponding to the crosshair cuts, discontinuities due to the artifacts were clearly seen, as shown in Fig. 4. Within this fixed location, all artifact magnitudes were measured in the superior-inferior (SI), anterior-posterior (AP), and medial-lateral (ML) directions. Because of asymmetry of the diaphragm and its motion (Figs. 3 and 4), discontinuities measured may not be the maximum amplitude of the artifacts. This procedure was applied consistently to every junction in the diaphragmatic range on both sides of the diaphragm. If no artifact was observed at a particular scan junction, the artifact amplitude was set to zero in calculating the average artifact per junction (APJ). Upon completion of the 0% phase image, the process was repeated for the remaining phase (10%-90%) images.
The magnitudes of each artifact in the SI, AP, and ML directions were recorded as follows: For INC and DUP artifacts, the SI magnitude, or the distance between the actual edge of the artifact and visually estimated edge without artifact, was measured with the "ruler" tool of the TPS and recorded as distance (cm)

2.C. Quantification of MAS
The magnitudes of artifacts in the SI, AP, and ML directions in coronal and sagittal images were averaged at each scan junction on both lateral sides of the diaphragm. Within the longitudinal scan range of the diaphragm, APJ were calculated for each respiratory phase. The MAS was defined as the averaged APJ (in cm) among all respiratory phases (n = 10).
The MAS depends not only on breathing irregularity but also on patient-specific diaphragm motion amplitude. Therefore, to make MAS depend solely upon the breathing irregularity, two normalized MAS indexes were introduced to the maximum diaphragmatic displacement and the diaphragmatic velocity, defined in Eqs. (6) and (7).
(b) MAS V is MAS normalized to mean diaphragmatic velocity [Eq. (4)] These three motion-artifact severities (MAS, MAS D , and MAS V ) were analyzed and correlated to the periodicity (BPI-5) of the corresponding 1DRW.

2.D. Correlation analysis between BPI in 1DRW and MAS in 4DCT
The Pearson correlation coefficient (r), defined as the covariance (cov) of two variables divided by individual standard deviation (σ), was used to quantify the linear relationship between BPI and MAS values, namely, where MAS can also be MAS D or MAS V . This equation was also used to check the correlations among BPI-5, BPI-4, and BPI-3, and among MAS, MAS D , and MAS V . For linear fitting of two variables (BPI, MAS), in which both variables contain non-negligible uncertainties, singular value decomposition (SVD) was applied to treat BPI and MAS equally in the total least-square linear regression.

3.A. Phase-dependent motion artifacts
The amplitude of the motion artifacts is phase dependent, as shown in Fig. 5(A). The smallest artifact magnitude occurs around the full exhalation phase, while the largest magnitude with greatest variation happens at or near the full inhalation phase. The average motion artifact over all respiratory phases and all 26 patients is 2.2 ± 1.1 mm. The histogram of motion artifacts including all phases and all patients is shown in Fig. 5(B), illustrating the distribution of the amplitude of artifacts, which is similar to a Poisson process with a mean of 2.2 mm [ Fig. 5(B)].

3.B. Periodicity and number of Fourier components
Among the 26 patients studied, the average BPI-5 is 0.75 ± 0.18, ranging from 0.25 to 0.93, as shown in Table I. Figure 6 shows BPI-k (k = 2-8) as a function of the number of largest Fourier coefficients in 26 patients. At k = 5, the slopes of the curves become steady, suggesting that BPI-5 is appropriate and effective. Figure 7(A) shows that the linear correlation increases slightly from (BPI-5 and BPI-3) to (BPI-5 and BPI-4), suggesting that k = 5 should be used even though their difference is small. For MAS, Fig. 7(B) depicts that MAS D correlates with MAS V (r = 0.92) better than with MAS (r = 0.82). The patient-averaged breathing period is 4.1 ± 1.3 s, ranging from 1.7 to 6.5 s. Table I shows the values of 26 patients, together with mean diaphragm displacement (0.8 ± 0.5 cm, ranging from 0.4 to 2.1 cm) and velocity (0.2 ± 0.1 cm/s, ranging from 0.1 to 0.4 cm/s). After patient-specific motion factor is removed from the MAS, the MAS D and MAS V should be primarily a result of breathing irregularities. Figure 8 shows a negative correlation of −0.53 between the MAS and BPI-5. However, MAS D and MAS V , which were normalized to the displacement and velocity of the diaphragm, respectively, become comparable among patients, resulting in improvement in correlation with BPI-5 substantially to −0.71 and −0.73, respectively. One patient outlier (not shown) was removed from this study because of a possible mismatch between 4DCT images and 1DRW dataset.

3.D. Timing requirements for the 1DRW and 4DCT analyses
The averaged time for Fourier analysis was about 3 s per 1DRW, including the time for loading the data, displaying the 1DRW curve, and saving the result. However, it would take a trained researcher at least 2-3 h to process the ten phases for each patient's 4DCT image set and record the data into an existing template. 111717-6 T I. Three motion-artifact severities (MAS = APJ, MAS D , and MAS V ) for the 4DCT images of 26 patients, whose maximum diaphragm displacement, mean motion period, and mean velocity are shown.

4.A. Breathing irregularity as the primary source of motion artifacts
Breathing irregularities have long been recognized as the primary source of motion artifacts in respiratory-correlated 4DCT. 2,13 For example, using a mobile phantom with a spherical ball and motion driven by real patient respiratory waveforms, the true volume of the sphere was estimated from F. 6. Breathing periodicity index (BPI-k) changes as a function of number of largest Fourier coefficient (k) for 26 patients. Generally, after k = 5, the slopes become steady.
4DCT imaging and significant deviation from its true value was shown. 7,11 However, in contrast to phantom studies, delineating the absolute volume variation in patients is nearly impossible since the initial tumor volume within patients is unknown; only relative variations of the tumor volume within a 4DCT image can be assessed. 6,7 The impact of breathing irregularities on 4DCT image quality is not entirely predictable and varies depending upon the differences of breathing between two adjacent bed positions for cine scans. This suggested that statistical means would be effective in quantifying the severity of motion artifacts in 4DCT as well as its relationship with breathing irregularities observed in the 1DRW. Thus, we studied the 1DRW-4DCT relationship using the linear correlation approach. We found that MAS is lowest at full exhalation, while higher at other phases. In addition, the standard deviation around exhalation is the smallest since the passive relaxation process tends to be reproducible [ Fig. 5(A)]. This is consistent with the general knowledge on reproducibility within the breathing cycle. 2 However, as the artifact severity is also a function of other motion factors, such as displacement or velocity: generally the greater the motion, the larger the artifact. In other words, if an irregular breather has little diaphragm motion, the magnitude of artifacts will be limited by the motion amplitude. To remove the motion amplitude factor and make the irregularity the sole source of the artifacts, we introduced MAS D and MAS V . Our results (Fig. 8 and Table I) clearly depict a stronger correlation between BPI-5 (the cause) and MAS D /MAS V (the consequences) than between BPI-5 and MAS. The correlation is negative because periodicity is opposite to irregularity. This is the first time that a quantitative relationship is established, which not only links the two random phenomena but also quantifies the 1D-to-4D correlative relationship.
The artifact distribution from this 26-patient study shows artifacts are a common event (Figs. 5 and 7). All patients have nonzero APJ, ranging from 1 to 5 mm. This is consistent with the previous 50-patient study, in which 90% of patient 4DCTs contained motion artifacts. 13 Note that APJ value is small because it averages with zero where motion artifact F. 8. Correlation between breathing periodicity index (BPI-5) and motionartifact severity MAS (A), MAS D (normalized to displacement) (B), and MAS V (normalized to velocity) (C). As the uncertainties in BPI and MAS are similar, a linear regression treating BPI and MAS equally is applied using SVD (singular value decomposition). Nevertheless, the linear fitting does not affect the correlation between BPI and MAS.
is not observable. We note that the RPM system contains a predictive filter, which is supposed to provide an assessment of breathing periodicity. However, this is a proprietary software tool that provides little information about its use of the RPM data. Therefore, this filter could not be utilized in a quantitative study.

4.B. Potential clinical and research applications
With the knowledge of the relatively high correlation between patient breathing irregularities and resulting motion artifacts, it becomes possible to predict the severity of motion artifacts in clinical 4DCT images. Two potential applications could be implemented in the clinic to improve 4DCT image quality with reduced motion artifacts. The BPI can be applied to predict artifact severity (1) prior to 4DCT scan or (2) during 4DCT scan. In both applications, appropriate actions can be taken when a patient's BPI score is or becomes poor.
In the first scenario, a patient's breathing regularity could be quantified after a short (<60 s) session of RPM monitoring and the BPI score can be calculated on-the-fly or at the end of the test. If BPI is lower than 0.85, then breath coaching may be an option. Breath coaching has been used clinically to reduce breathing irregularities. 28 Three general breathing irregularities can affect 4DCT image quality: (A) frequency, (B) amplitude, and (C) breathing pattern [variation in chest (costal) and abdominal (diaphragmatic) muscles involvements]. Audio coaching regulates breathing frequency, 29 but it often makes patients breathe deeper, resulting in greater tumor motion. Video coaching regulates both frequency and amplitude 28 but has no control in maintaining the same breathing pattern. So far, no coaching method is known to regulate breathing pattern variations. In addition, breath coaching requires a patient's voluntary effort, which introduces a new source of uncertainty including variations (learning curve) between 4DCT simulation and the course of treatment. Due to medical conditions, some patients may have difficulties in breathing regularly and/or following the coaching instructions. 28 Therefore, although coaching may improve breathing regularity, it adds another layer of workload with associated uncertainty.
In the second scenario, on-the-fly DFT analysis of 1DRW can be performed to monitor the BPI variation using graphics processing unit (GPU)-based parallel computing. BPI as a function of time can potentially monitor patient breathing behavior in near-real time and identify sudden irregularities during the scan. A sudden change in the BPI curve between two adjacent bed positions could be used as a trigger to temporarily suspend 4D scanning and reassume the scan after BPI becomes normal. This could minimize breathing irregularity at the image acquisition level to gain improved 4DCT image quality.
From the clinical research perspective, 4DCT provides a rich dataset about respiratory motion for retrospective patient studies. However, image quality is of paramount importance to set a solid foundation for studying respiratory motion. Manual selection of patient 4DCT images is subjective, tedious, and labor-intensive. Our finding of a correlative relationship could serve as a rapid screening method to identify a large number of 4DCT images with minimal motion artifacts.

4.C. Limitations and future directions
The correlation between 1DRW breathing periodicity and 4DCT motion-artifact severity is established for cine-mode phase-binned 4DCT images and is useful for clinics that use similar 4DCT scanning protocols. For helical scan or amplitude-binned 4DCT, this result may not be directly applicable. Although the methodology is generally valid it has to be "recalibrated" for different types of 4DCT image sets.
In this study, we focused on the distortional artifacts but ignored motion blurring, which also contributes to overall motion artifacts. Motion blurring occurs within the same bed scan and is difficult to quantify without knowledge of the actual anatomy although it is feasible to scan the patient with breath holding as a control for some patients, or to create a theoretical model to account for the residual motion artifacts in 4DCT. 11 A global Fourier analysis of 4DCT was also previously applied to assess overall motion artifacts, including the motion blurring, but only provide binary classification of image quality. 30 The DFT-based periodicity definition uses sinusoidal functions to represent a 1DRW curve, and it is well-known that Fourier transformation produces multiple sine/cosine waves to mimic a curve with possible discontinuities. In rare circumstances, such as a periodic but discontinuous 1DRW, BPI may not fully represent periodicity, although it can catch principal components. In addition, in case a severe baseline drift in 1DRW occurs, which could result in large motion artifacts, the corresponding Fourier component (with a low frequency) should be excluded from BPI-5. In this study, we have not seen any case with severe baseline drift in the 26 patients. Nevertheless, it is worthwhile to perform a further investigation in the future.
In both clinical practice and clinical research, future application of this 1D-4D method could improve patient care and facilitate screening for high-quality images with minimal artifacts. The time saving in estimating 4DCT image quality is tremendous and the ability to differentiate regular breathers from irregular breathers prior to 4DCT scans would assist in making individualized decisions to fit patient respiration. We can take advantage of the periodicity of those with steady breathing cycles while providing help to those who need it. Moreover, the speed performance of Fourier analysis can be further improved using GPU for parallel computing to provide real-time results. This performance is essential in monitoring BPI on the fly during 4DCT scanning for a potential "gating" of the 4DCT scan using the BPI to set an irregularity threshold.

CONCLUSION
We have established a correlation between 1DRW breathing periodicity and 4DCT motion-artifact severity, by defining and quantifying these two quantities (BPI and MAS). We found that the correlation is about −0.70, which can be applied to predict the quality of 4DCT images based on 1DRW alone, before or during 4DCT scanning. This method can be useful in rapidly selecting high-quality 4DCT images for clinical research with a tremendous time saving factor. Further investigation is necessary to implement this method in the radiotherapy clinic to make an adaptive decision at the time of 4DCT simulation, and subsequent planning and treatment.