Markerless Four-Dimensional-Cone Beam Computed Tomography Projection-Phase Sorting Using Prior Knowledge and Patient Motion Modeling: A Feasibility Study.

AIM: During cancer radiotherapy treatment, on-board four-dimensional-cone beam computed tomography (4D-CBCT) provides important patient 4D volumetric information for tumor target verification. Reconstruction of 4D-CBCT images requires sorting of acquired projections into different respiratory phases. Traditional phase sorting methods are either based on external surrogates, which might miscorrelate with internal structures; or on 2D internal structures, which require specific organ presence or slow gantry rotations. The aim of this study is to investigate the feasibility of a 3D motion modeling-based method for markerless 4D-CBCT projection-phase sorting.
METHODS: Patient 4D-CT images acquired during simulation are used as prior images. Principal component analysis (PCA) is used to extract three major respiratory deformation patterns. On-board patient image volume is considered as a deformation of the prior CT at the end-expiration phase. Coefficients of the principal deformation patterns are solved for each on-board projection by matching it with the digitally reconstructed radiograph (DRR) of the deformed prior CT. The primary PCA coefficients are used for the projection-phase sorting.
RESULTS: PCA coefficients solved in nine digital phantoms (XCATs) showed the same pattern as the breathing motions in both the anteroposterior and superoinferior directions. The mean phase sorting differences were below 2% and percentages of phase difference < 10% were 100% for all the nine XCAT phantoms. Five lung cancer patient results showed mean phase difference ranging from 1.62% to 2.23%. The percentage of projections within 10% phase difference ranged from 98.4% to 100% and those within 5% phase difference ranged from 88.9% to 99.8%.
CONCLUSION: The study demonstrated the feasibility of using PCA coefficients for 4D-CBCT projection-phase sorting. High sorting accuracy in both digital phantoms and patient cases was achieved. This method provides an accurate and robust tool for automatic 4D-CBCT projection sorting using 3D motion modeling without the need of external surrogate or internal markers.

INTRODUCTION

Target control and normal tissue complication probabilities are highly correlated to the target localization accuracy in radiation therapy.[1,2] Accurate localization of moving targets is challenging due to the respiratory motion of the normal tissue and target.[3] Conventionally, three-dimensional-cone-beam computed tomography (3D-CBCT) has been used for target localization before the treatment. However, 3D-CBCT does not capture the motion information and only provides an average location of the moving target, which may not be sufficient for high-precision radiotherapy treatments. More recently, 4D-CBCT has been developed for 4D localization of moving targets. 4D-CBCT reconstructs patient on-board 4D images by first sorting the acquired cone beam projections into different respiratory phases.[4-9] The accuracy of the projection sorting determines the accuracy of the reconstructed 4D-CBCT images, which in turn affects the tumor localization accuracy using 4D-CBCT.

Previous projection sorting algorithms can be classified into the following categories: (1) External surrogate-based sorting. This method uses the motion of external surrogates such as abdominal displacement or lung air volume to determine the respiratory phase. Commercially available products include the Varian real-time position management (RPM) system (Varian Medical Systems, Palo Alto, CA, USA)[10,11] and strain-gauged Anzai belt (Anzai Medical Systems, Tokyo, Japan).[12,13] Spirometry was also implemented to detect changes in lung air volume.[14] All of these external surrogates were assumed to correlate with the displacement of internal anatomy, but studies have demonstrated that uncertainties exist in this relationship,[15-18] which can potentially lead to sorting errors. (2) Internal surrogate-based sorting. To mitigate the limitations of external surrogates, several internal “markerless” techniques were proposed. One example is known as the Amsterdam shroud (AS) technique,[19] in which the diaphragm motion is enhanced by converting 2D projections into the so-called AS image and then breathing signal is extracted by temporal derivative. However, this technique requires a visible oscillating structure present across all projection angles, which may not always be available. Berbeco et al.[20] proposed a technique based on the radiological path-length change with lung volume expansion/retraction. The breathing-phase information is extracted by analyzing the image intensity fluctuation of fluoroscopic images in the middle of lung[20] or region of interest in CBCT projections.[21] However, the angular dependence of image intensities may affect the accuracy of the sorting. Siochi[22] developed a motion tracking method by converting the mega-voltage CBCT rotating coordinate system into a patient coordinate system. The method requires the diaphragm to be present in the projections images, which may not always be available for upper lobe lesions. Lewis et al.[23] proposed a template-based tumor tracking algorithm that registers each phase of 4D-CT to CBCT and generates digitally reconstructed radiographs (DRRs) to match with the projections. This method could be affected by the “out of range” problem when the on-board motion range exceeds that of the template. Vergalasova et al.[24] proposed a Fourier domain-based method to extract breathing signal from CBCT projections. The accuracy of this method is degraded for 4D-CBCT scans with faster gantry rotation speeds due to the angular dependence of the image intensities.

The aforementioned internal markerless phase sorting techniques either require specific anatomical features,[19,22,23] or are affected by gantry rotation angle or speed.[20,21,24] Recently, our group and others have proposed using prior information and motion modeling for CBCT volume reconstruction with limited angle cone beam projections.[9,25] The CBCT volume can be estimated by iteratively matching DRRs of deformed prior CT to on-board projection images. Our previous results demonstrated that a single projection is not adequate for CBCT reconstruction when there is substantial breathing pattern change from simulation to on-board treatment.[8,9] However, we hypothesized in this study that the motion model coefficients solved from a single projection are adequate for projection-phase sorting for 4D-CBCT reconstruction.[26] Since the motion model coefficients are correlated with the overall 3D motion information of the patient, they can potentially provide more robust information for projection sorting without the need of specific anatomical features or gantry rotation speeds. In our study, the method was developed and evaluated using XCAT simulation with different scenarios and lung cancer patient CBCT data with different gantry rotation speeds and scanning modes.

METHODS

Motion modeling from prior four-dimensional-computed tomography

As illustrated in Figure 1, patient motion model is generated from 10-phase planning 4D-CT images. One phase of the planning 4D-CT (e.g., end of exhalation phase) is defined as the “prior.” The other nine phases are deformed to the prior using Velocity AI (Varian Medical Systems, Palo Alto, USA) to generate nine 3D deformation vector fields (DVFs) which describe patient respiratory deformations. Principal component analysis (PCA) is then used to extract the principal components {PC} of the deformation fields, which represent the major deformation patterns of the patient. Patient respiratory deformation at any instant can be represented as a linear combination of the principal components. In this study, we used the first three principal components corresponding to the three largest eigenvalues as they were proven to be sufficient in depicting lung motion:[27-29]

Figure 1

Flowchart of patient motion modeling and on-board principal component analysis coefficient calculation

Where, D is the deformation field vector at any instant, Davg is the average of the nine deformation fields initially extracted from 4D-CT, PC is the i principal component, and wi represents the corresponding coefficients for the principal components.

Solving the principal component analysis coefficients for each single cone beam computed tomography projection

The on-board CBCT volume at any instant is approximately considered as a deformation from the CT “prior.” The deformation field map D is represented as a weighted linear combination of the principal motion modes {PCi} as shown in Eq (1).

The weighting coefficients {wi} are solved by minimizing the objective function f (w), which is defined based on the data fidelity constraint:

Where, NCC is normalized cross correlation between two 2D images, M is projection operator calculating DRR from new CBCT volume, CBCT is the new volume after applying estimated deformation map onto prior volume Iprior, and OBI is the single CBCT on-board projection image.

The PCA coefficients {wi} are solved iteratively using the gradient descent method with a maximum of ten iterations and ten backtracking line searches for each iteration. The PCA coefficients describe the degree of 3D patient deformation caused by the respiratory motion, and therefore can be used as surrogates for 4D projection sorting. In our study, the primary PCA coefficients along the superoinferior (SI) and/or anteroposterior (AP) directions were used for CBCT projection-phase sorting due to the more substantial respiratory motions along these two directions.

Phase sorting based on principal component analysis coefficients

The PCA coefficients {wi} are solved for all projections and the primary coefficient w1 in the SI or AP direction was used for phase sorting. The coefficients are plotted for all projections and the valley projections are defined as the peak-inspiration phase (0% or 100%). The phase for each projection is defined as the percentage of the respiratory cycle that has passed since peak inspiration.

XCAT digital phantom study

The 4D digital extended-cardiac-torso (XCAT) phantom[30] was used to test the algorithm. Based on the human anatomical database from the National Library of Medicine, XCAT uses nonuniform rational B-spline surfaces to model highly realistic and detailed human anatomical structures and images. It can be utilized to generate 4D images according to respiratory profiles and anatomical parameters specified by the user. XCAT has been widely used in multiple studies as an algorithm and as an evaluation and verification tool.[9,24,27,31-33]

The respiratory motion of the body volume is controlled by two breathing curves: diaphragm and chest wall curves. The diaphragm curve controls both diaphragm motion and motions of the lung, liver, stomach, and other downstream organs. It generally determines motion in the SI direction. The chest wall curve controls AP expansions of the body, the ribcage, and lungs. In this study, both curves were simplified as sinusoidal curves with 5 s respiratory period. The peak-to-peak amplitude of the two curves was set to 3 and 2 cm, respectively. A spherical lesion of 30 mm diameter was inserted in the middle of the lung.

A ten-phase 4D-CT was then simulated as the prior knowledge. We used 40 keV monochromatic CT source to approach the effective energy of 120 kVp polychromatic spectrum used in clinical CT scanners. CT volume of each phase was composed of 256 × 256 × 150 voxels, with 1.67 mm isotropic resolution. The end-expiration phase was selected as the prior volume as the body structure is most stable in this phase.

As shown in Table 1, nine XCAT scenarios were generated to simulate 4D-CBCT volumes. Scenario 0 represents no breathing pattern or anatomical changes from prior 4D-CT, while scenarios 1–8 include changes in tumor motion amplitude (S1: 3 cm to 2 cm), tumor diameter (S2: 30 mm to 20 mm, S3: 30 mm to 40 mm), tumor spatial shifts (S4: 8 mm SI, S5: 8 mm AP, and S6: 5 mm in SI, AP, and LR), and motion removal in SI or both SI and AP directions (S7, S8).

Table 1

Four-dimensional digital extended-cardiac-torso scenarios simulating changes in respiratory motion amplitude, tumor dimension and position, and respiratory motion directions

XCAT scenario	Tumor diameter (mm)	Diaphragm motion amplitude SI (mm)	Chest wall motion amplitude AP (mm)	Tumor position SI (mm)	Tumor position AP (mm)	Tumor position LR (mm)
S0	30	30	20	0	0	0
S1	30	20	12	0	0	0
S2	20	20	12	0	0	0
S3	40	20	12	0	0	0
S4	30	20	12	8	0	0
S5	30	20	12	0	8	0
S6	30	20	12	5	5	5
S7	30	20	0	0	0	0
S8	30	0	0	0	0	0

XCAT: 4D digital extended-cardiac-torso; AP: Anteroposterior; SI: Superoinferior; LR: Left-right; 4D: Four dimensional

These 4D-CBCT volumes were set as “ground truth” on-board images to simulate on-board projections to be sorted by the proposed method. Totally around 200 projections within 20 respiratory cycles were simulated over the 200° scan angle. The source to isocenter distance was set to 100 cm, and isocenter of detector distance was set to 50 cm. Each projection contains 512 × 384 pixels with each pixel of 0.78 mm in both dimensions.

Patient study

In addition to the XCAT digital phantom studies, the proposed method was evaluated on five sets of patient data that were previously acquired under an investigational 4D-CBCT study.[34,35] The 4D-CBCT projections were acquired with an X-ray flat panel detector mounted orthogonally to the MV treatment gantry on the Varian Trilogy system (Varian Medical Systems, Palo Alto, CA, USA). Two of the patients were scanned with full-fan projections over a 200° arc at slow-gantry speeds of 0.71 and 0.6°/s, with total scan lengths of 4.5 and 5.7 min and frame rates of 7 and 5 fps, respectively. Three of the patients were scanned with half-fan projections over 360° arc at normal gantry speed of 6°/s, with total scan length of 1 min and frame rate of 10, 15, and 15 fps, respectively.

Phase sorting evaluation

In order to evaluate the phase sorting accuracy of the proposed technique, respiratory phases were derived from the ground truth respiratory signal in the digital phantom study and visual inspection of the on-board projections in patient study. Projections that correspond to peak inspiration were identified and recorded as 0% or 100%. The projection phases in-between peak inspirations were linearly interpolated as before. This “manual phase sorting” served as the gold standard for comparison.

The projection-phase sorting accuracy was evaluated by three parameters. One is the average difference in phases across the entire dataset, which is defined as:

Where, n is the total number of projections, i is the index of projection, phase (imanual) is the phase value sorted by manual method for the ith projection, and phase (i) is the phase value sorted by PCA coefficient based method for the ith projection.

The other two parameters are the percentage of projections with a phase difference within 10% or 5%:

Where, n, phase (imanual), and phase (i) had the same definition.

The three quantitative evaluations were performed on all the XCAT simulation scenarios and five patient data sets. All image processing and data analyses were performed with in-house developed MATLAB programs (The Mathworks Inc., Natick, MA, USA).

RESULTS

XCAT study

To verify the effectiveness of the PCA coefficients in representing the respiratory motion in the original prior 4D-CT images, we calculated the coefficients along all the three directions for the nine DVFs extracted from the 4D-CT, as shown in Figure 2. Sinusoidal respiration motion patterns are manifested in the coefficients of the first principal component along all directions. Second and third PCA coefficients (w2, w3) did not show clear sinusoidal pattern since the 2nd and 3rd principal components capture much smaller portion of the variances of the respiratory deformations. The coefficients for the SI and AP directions showed clearer breathing patterns than those for the lateral direction, due to the more substantial respiratory motions along SI and AP directions in the XCAT phantom. In the following studies, we will use the PCA coefficient of the first principal component (w1) for projection sorting.

Figure 2

Principal component analysis coefficients of nine deformation vector fields of prior four-dimensional-computed tomography. From top to bottom, the principal component analysis coefficients (w1, w2, and w3) correspond to deformable field vectors along different directions (superoinferior, anteroposterior, and lateral). From left to right, the principal component analysis coefficients correspond to the first three primary principal components. The X-axis represents the phase number of the nine deformation vector fields used to calculate principal component analysis coefficients

Figure 3 shows the calculated primary PCA coefficients (w1) for each CBCT projection in XCAT scenario 0 in three orthogonal directions. Full-fan scanning with 200° rotation was simulated, with a rotation speed of 3.3°/s and 0.9° projection interval. The PCA coefficients calculated for all projections formed a sinusoidal cycle pattern, which was most clear in SI and AP directions and less apparent in lateral direction.

Figure 3

Primary principal component analysis coefficient (w1) for three orthogonal motion directions for XCAT scenario 0. Each dot represents principal component analysis coefficient calculated for one on-board projection. From top to bottom: superoinferior, anteroposterior, and lateral directions, respectively

Figure 4 shows the effect of changes in motion amplitude and tumor anatomy from prior 4D-CT to on-board 4D-CBCT scans on the calculated PCA coefficients. Figure 4a shows that, with lower motion amplitude, the oscillation pattern of PCA coefficient stayed the same, while the amplitude decreased (red dotted line). Figure 4b illustrates the sensitivity of PCA coefficient in detecting the changes in respiratory motion directions. The PCA coefficients showed no cyclic pattern in the AP direction when the respiratory motion in AP direction was removed (red), and there was no cyclic pattern in both AP and SI directions when the respiratory motions in both AP and SI directions were removed (green). Figure 4c demonstrates that changes of tumor diameter from 30 mm to 20 mm or 40 mm did not change the cyclic patterns or amplitudes of the PCA coefficients (green and red). Figure 4d shows that shifts of tumor central location in three orthogonal directions (blue, green, and red) did not change the cyclic pattern or amplitudes of the PCA coefficients.

Figure 4

Primary principal component analysis coefficients (w1) of XCAT scenarios 0–8 in superoinferior and anteroposterior directions. (a) Change in motion amplitude; (b) changes in motion directions (superoinferior and anteroposterior, superoinferior only, and no motion); (c) changes in tumor size (20 mm, 30 mm, and 40 mm); and (d) shifts in tumor central location in superoinferior, anteroposterior, and diagonal directions

Table 2 shows the phase sorting accuracy of the nine XCAT scenarios. For all the eight quantifiable scenarios, the average phase difference was 2% or less. All of the projections in the eight scenarios showed a < 10% phase difference from the ground truth. Similar accuracies were observed using PCA coefficients in SI and AP directions.

Table 2

Principal component analysis coefficient-based phase sorting accuracy for nine four-dimensional digital extended-cardiac-torso scenarios

XCAT	Feature	DVF coordinate
		SI			AP
		Average phase difference (%)	Percentage of phase difference < 10%	Percentage of phase difference < 5%	Average phase difference (%)	Percentage of phase difference < 10%	Percentage of phase difference < 5%
S0	OB motion same as prior	0	100	100	0	100	100
S1	Lower motion amplitude	0.5	100	95	2	100	80
S2	Tumor size decreased	1	100	90	1.5	100	85
S3	Tumor size increased	2	100	80	0.5	100	95
S4	Tumor shifted in SI	0.5	100	95	1.5	100	85
S5	Tumor shifted in AP	1	100	90	2	100	80
S6	Tumor shifted in SI, AP, and LR	1.5	100	85	0.5	100	95
S7	AP motion removed	0.5	100	95	N/A	N/A	N/A
S8	AP and SI motion removed	N/A	N/A	N/A	N/A	N/A	N/A

XCAT: 4D digital extended-cardiac-torso; AP: Anteroposterior; SI: Superoinferior; LR: Left-right; DVF: Deformation vector field; N/A: Not available; 4D: Four dimensional; OB: On-board

Figure 5 shows a representative OBI and final DRR of the deformed prior CT images in a patient case with half-fan mode and fast gantry rotation (~1 min). The anatomy showed in the final estimated DRR [Figure 5b] agrees well with the OBI [Figure 5a]. The normalized cross correlation between OBI and final DRRs optimized in Equation 2 was above 0.9 for most of the projections, as shown in Figure 5c.

Figure 5

Patient cone beam computed tomography on-board projection and final digitally reconstructed radiograph after optimization. (a) On-board projection; (b) final digitally reconstructed radiograph after optimization; and (c) normalized cross correlation between on-board projection and digitally reconstructed radiographs for all 885 projections

Figure 6 displays the solved PCA coefficients (SI direction, primary) for all the five CBCT patients. Respiratory motion patterns were clearly manifested by the changes in the coefficients in both full-fan slow gantry rotation cases (patients 1 and 2) and half-fan fast gantry rotation cases (patients 3–5). As expected, only the SI direction coefficients showed clear breathing patterns, possibly because SI is perpendicular to the OBI projection direction and is the major breathing motion direction.

Figure 6

Principal component analysis coefficients in the SI direction for all the five cone beam computed tomography patients. Patients 1–2, full-fan cone beam computed tomography slow gantry rotation (~5 min); patients 3–5, half-fan cone beam computed tomography fast gantry rotation (~1 min). The temporal resolution in the X-axis is not to the same scale because of different gantry rotation speeds and frame rate settings for each patient

Figure 7 shows the evaluation process for the PCA coefficient-based sorting for patient 1. The sorting result based on PCA coefficients is shown in blue, with 0/100% phases corresponding to the valleys in the coefficient plot in Figure 6. Manual phase sorting result is shown in red, with 0/100% phases corresponding to visually identified peak inspiration projections. The phase differences were < 10% for most of the projections.

Figure 7

Comparison of principal component analysis coefficient-based sorting with manual phase sorting. Top, phase sorting result using principal component analysis coefficient; middle, phase sorting result by visual inspection; and bottom, phase difference between the two methods

Table 3 shows the accuracy of PCA coefficient-based phase sorting for the five lung cancer patient 4D-CBCT datasets. The average phase differences ranged from 1.62% to 2.23%. The percentage of projections with phase sorting difference < 10% ranged from 98.4% to 100% and ranged from 88.9% to 99.8%

Table 3

Principal component analysis coefficient-based phase sorting accuracy for five lung cancer patients

Patient number	Average phase difference (%)	Projections with phase difference (%)
		< 10%	< 5%
P1 (full fan)	2.23	99.4	88.9
P2 (full fan)	1.81	98.7	93.4
P3 (half fan)	1.62	100.0	99.8
P4 (half fan)	1.82	100.0	96.2
P5 (half fan)	1.63	98.4	96.6

DISCUSSION

The results presented in this article have demonstrated the feasibility of using patient motion modeling-based PCA coefficients as an internal markerless surrogate for 4D-CBCT phase sorting. High phase sorting accuracy was achieved in both XCAT digital phantoms and real 4D-CBCT lung cancer patient data sets. The method avoids using any internal marker or external surrogates which may cause complications such as pneumothorax or errors due to potential mismatch between external surrogates and internal organ motions.

As demonstrated in the XCAT digital phantoms, this method showed robustness against patient breathing or anatomical changes from prior to on-board imaging, including changes in tumor size, location, and breathing motion amplitudes. By controlling the motion directions and changing the breathing amplitudes in the digital phantom, we further demonstrated that the primary PCA coefficient can also be indicative of the major motion directions and amplitudes.

The proposed motion modeling and PCA coefficient-based method utilizes prior 4D-CT, 3D deformation maps, dimensional reduction, on-board CBCT projection matching, and iterative optimization for phase sorting. Instead of using 1D- or 2D-based information as in previous methods,[10-12,14,20,21,24] 3D motion information is used in searching for the correct breathing phase. Using information in this higher dimension may have led to several advantages of this method.

First, different from the Fourier transform-based phase sorting method[24] or image intensity based methods,[21] this method does not have gantry rotation speed dependence. This advantage has been demonstrated in our 4D-CBCT patient datasets, in which a high accuracy was achieved with both fast and slow rotation speeds. This feature allows us to increase the gantry rotation speed to substantially reduce the scanning time and imaging dose of 4D-CBCT.[36]

Second, no postfiltering of the surrogate is required because, rather than using the summation of X-ray attenuation information in intensity-based methods, the PCA coefficient represents the overall 3D deformation based on motion modeling. Therefore, no correction for the slow changing signal due to projection path-length angular change is needed, which is sensitive to variations in external structure or internal structure deformations.[21,24]

Finally, for half-fan CBCT patient cases, lateral projections may have very limited 2D information. Using 3D motion modeling and matching, this method have achieved decent sorting results in these areas, which could be one of the contributing factors for its overall high accuracy.

This method was inspired by our earlier work in CBCT image reconstruction using the motion modeling method. Due to the suboptimal accuracy of reconstruction using motion modeling alone, Zhang et al.[9] used free-form deformation to further improve the reconstruction accuracy after motion modeling. While 3D motion modeling itself may not be sufficient for 3D volumetric reconstruction, our study demonstrated that 3D motion modeling is accurate for 1D breathing signal extraction and thereafter used for projection-phase sorting. The method can be used for 4D-CBCT reconstruction using either full-scan angles or limited scan angles.[7,37,38]

XCAT results for scenarios S0, S1, S7, and S8 showed potential correlation between the values of PCA coefficients and breathing motion amplitudes, which suggests the potential to use the PCA coefficient as surrogates for amplitude-based sorting.[13,39]

The limitation of this method may be the calculation speed. On a Dell Optiplex 7010 PC with 3.4 GHz CPU and 8 GB RAM, 3–6 min is required for each projection-phase calculation. However, the PCA coefficient calculation for the projections is completely independent and thus can be paralleled. In addition, graphics processing unit acceleration can potentially accelerate the calculation speed substantially (< 1 s).[32]

In conclusion, the feasibility of using the PCA coefficient of a patient respiratory motion model as an internal markerless surrogate for phase sorting of 4D-CBCT projections has been demonstrated through both simulation and patient studies. The percentages of projections sorted within 10% phase difference were 100% and 98% for XCAT and patient data, respectively. The average phase differences between sorted phases and the gold standard were below 3% for both XCAT and patient data. Overall, the proposed method is promising in providing robust markerless automatic projection sorting for 4D-CBCT scans.