Technical Note: On the spatial correlation between robust CT-ventilation methods and SPECT ventilation.

PURPOSE: The computed tomography (CT)-derived ventilation imaging methodology employs deformable image registration (DIR) to recover respiratory motion-induced volume changes from an inhale/exhale CT image pair, as a surrogate for ventilation. The Integrated Jacobian Formulation (IJF) and Mass Conserving Volume Change (MCVC) numerical methods for volume change estimation represent two classes of ventilation methods, namely transformation based and intensity (Hounsfield Unit) based, respectively. Both the IJF and MCVC methods utilize subregional volume change measurements that satisfy a specified uncertainty tolerance. In previous publications, the ventilation images resulting from this numerical strategy demonstrated robustness to DIR variations. However, the reduced measurement uncertainty comes at the expense of measurement resolution. The purpose of this study was to examine the spatial correlation between robust CT-ventilation images and single photon emission CT-ventilation (SPECT-V).
METHODS: Previously described implementations of IJF and MCVC require the solution of a large scale, constrained linear least squares problem defined by a series of robust subregional volume change measurements. We introduce a simpler parameterized implementation that reduces the number of unknowns while increasing the number of data points in the resulting least squares problem. A parameter sweep of the measurement uncertainty tolerance, τ , was conducted using the 4DCT and SPECT-V images acquired for 15 non-small cell lung cancer patients prior to radiotherapy. For each test case, MCVC and IJF CT-ventilation images were created for 30 different uncertainty parameter values, uniformly sampled from the range 0.01 , 0.25 . Voxel-wise Spearman correlation between the SPECT-V and the resulting CT-ventilation images was computed.
RESULTS: The median correlations between MCVC and SPECT-V ranged from 0.20 to 0.48 across the parameter sweep, while the median correlations for IJF and SPECT-V ranged between 0.79 and 0.82. For the optimal IJF tolerance τ = 0.07 , the IJF and SPECT-V correlations across all 15 test cases ranged between 0.12 and 0.90. For the optimal MCVC tolerance τ = 0.03 , the MCVC and SPECT-V correlations across all 15 test cases ranged between -0.06 and 0.84.
CONCLUSION: The reported correlations indicate that robust methods generate ventilation images that are spatially consistent with SPECT-V, with the transformation-based IJF method yielding higher correlations than those previously reported in the literature. For both methods, overall correlations were found to marginally vary for τ ∈ [ 0.03 , 0.15 ] , indicating that the clinical utility of both methods is robust to both uncertainty tolerance and DIR solution.

INTRODUCTION

Computed tomography (CT)‐derived ventilation is an image processing modality that quantifies the apparent voxel volume changes within an inhale/exhale CT image pair as a surrogate for pulmonary ventilation. , With the exception of newer deep learning methods and methods based purely on image segmentation, there are two primary classes of CT‐ventilation methods: transformation based and intensity based. Both classes require, as preprocessing steps, lung volume delineation and deformable image registration (DIR). The DIR solution approximates respiratory motion and provides a spatial mapping between the inhale and exhale lung geometries. Intensity‐based methods estimate volume change using the Hounsfield Unit values of spatially corresponding inhale/exhale lung voxels. , Transformation‐based methods, in contrast, compute volume changes directly from the DIR‐defined spatial transformation via numerical approximation of the Jacobian Factor. While the images produced by transformation‐based methods have been shown to be highly variable with respect to numerical implementation, intensity‐based methods require key heuristic steps, including smoothing and vasculature segmentation. Taken together, these factors contribute to the suboptimal reproducibility of previously proposed CT‐ventilation methods. ,

We recently described two new CT‐ventilation algorithms designed to address reproducibility issues previously described in the literature , , : the intensity‐based mass conserving volume change (MCVC) method and the transformation‐based Integrated Jacobian Formulation (IJF) method. As opposed to traditional approaches, MCVC and IJF estimate the Jacobian factor of the DIR transformation from a series of spatially corresponding lung subregions. The numerical uncertainty in the subregional volume change measurements is modeled with Gaussian statistics, which allows for their uncertainty to be characterized and, consequently, controlled through the definition of an uncertainty tolerance parameter. In this context, the uncertainty tolerance is defined with respect to the uncertainties associated with the DIR solution and does not address uncertainties associated with CT acquisition. This strategy provides robustness to DIR variability at the expense of measurement resolution. While IJF and MCVC images were both shown to be robust to DIR variations, no comparisons to clinically accepted modalities have been reported.

The purpose of this study is to (a) introduce a parameterized numerical implementation of the MCVC and IJF methods and (b) characterize the spatial correlation between the resulting CT‐ventilation images and single photon emission computed tomography ventilation (SPECT‐V). Using a straightforward parameter sweep of the uncertainty tolerance, we examine the effect of the tolerance on the correlation between the resulting CT‐ventilation images and the SPECT ventilation images.

MATERIALS AND METHODS

Subregional volume change estimation

Computed tomography‐ventilation methods employ DIR to compute a spatial transformation which maps a reference image lung volume onto a target image lung volume . For a general subregion , the volume scaling factor under is described by the Jacobian factor:assuming is diffeomorphic. The transformation‐based IJF ventilation method numerically approximates Eq. (1) directly from the DIR solution as the sample mean of the membership oracle functionwhereand . Equation (2) is a “hit‐or‐miss” volume approximation and represents the number of “hits” that land within . The MCVC method, in contrast, represents an intensity‐based approach and requires the inhale/exhale CT image pair to be converted into two corresponding HU‐defined density functions, which we denote here as the reference image and the target image Assuming mass is conserved, the MCVC estimate is defined aswhere and represent the sample mean densities within and , respectively. As previously described in detail, the uncertainties in the Eqs. (2) and (4) measurements can be characterized using the standard errors of the corresponding sample means. , For IJF, it was shown that with 95% probabilitywhere the minimum hit count, , depends on and is a specified tolerance. For MCVC, it was shown that with 95% probabilityandfor . Equations (5), (6), (7) essentially state that the amount of uncertainty in the volume change estimates goes down as the size of goes up. Thus, the improved robustness of MCVC and IJF comes at the expense of measurement resolution.

Moving least squares ventilation image

A CT‐ventilation image, , requires computing the discretized variables.where represents the volume change of the undeformed unit volume voxel centered on . Because both IJF and MCVC obtain volume change estimates for a series of subregions their resulting mathematical formulations are equivalent. Specifically, for , the subdomain data acquisition process results in a linear system of equations that relates the unknown to the subdomain data:whereand the elements of contain the corresponding subregional estimates [Eq. (2) for IJF and Eq. (4) for MCVC]. Previous IJF and MCVC implementations solve a computationally intensive, large‐scale linear least squares problems defined by Eq. (9). In order to facilitate broader accessibility of the established methods and improve overall method robustness, we instead parameterize with moving least squares (MLS) to reduce the number of unknowns and, consequently, the overall computational and memory storage requirements. This approach also allows for the definition of an overdetermined data fitting problem.

Given a set of knot locations with corresponding scalar parameter values , the Shepard's class of moving least squares approximation is defined as:where the proximal weighting function is of the form

The Eq. (11) parameterization reduces the number of unknowns required to generate the volume change (ventilation) image from (the total number of voxels in the reference lung region of interest) down to (the number of knots used for the discretization).

Image Computation

A ventilation image, as defined by Eq. (11), can be recovered from the Eq. (9) subregional volume change estimates by solving the following constrained least squares problem:whereand

We point out that the data matrix can be constructed without explicitly constructing thereby reducing the overall memory requirements.

The Eq. (13) inequality constraints prevent physically impermissible negative volumes. Requiring is enough to guarantee that the resulting (Jacobian estimates) are strictly positive since the Eq. (11) MLS parameterization is simply a moving average operator applied to the . The equality constraint represents consistency with the global volume change, as measured on the full and lung volumes, which is often used as a validation metric for proposed ventilation methods. , , , , However, the constraint constant, , depends on the choice of method. For transformation‐based IJF,and for intensity‐based MCVC

The Eq. (13) problem structure is equivalent to the one first presented in Ref. [11] with two exceptions. First, the number of subregional estimates can be chosen such that the problem is highly overdetermined while the reduced number of total unknowns allows for the application of standard solvers (as opposed to the customized augmented Lagrangian method presented in Refs. [10, 11]), including the MATLAB lsqlin interior‐point method (release R2019a, The Mathworks Inc, Natick, Massachusetts, USA). Second, whereas previous implementations regularized the problem by penalizing the norm of the spatial gradient of , the proposed regularization model [second term in the Eq. (13) objective function] penalizes the variance of the moving least squares parameters with the mean defined as the average voxel volume change enforced by the equality constraint. This choice of regularization is derived from the recommendations of the AAPM Radiation Therapy Task Groupe No. 132, which state that large variations in the Jacobian are potentially indicative of error.

The original approaches described in Refs. [10, 11] both defined an underdetermined data fitting matrix and relied on spatially smoothing regularization to ensure a well‐posed problem. This approach, when applied to the MLS parameterization, had the potential to induce over‐smoothing in the resulting images. In contrast, the Eq. (13) formulation defines an overdetermined data fitting matrix, , and penalizes overall variance in the MLS recovered function, as opposed to local smoothness which is already inherent to MLS functions.

Numerical implementation

The maximum inhale and exhale phases from 4DCTs were used for this study. Lung masks were generated using a semi‐automated histogram segmentation (as done in Ref. [2]). A dart‐throwing algorithm was applied to in order to generate MLS knot locations, with approximately 30 mm uniform spacing. Similar to the number of cubic spline knots used for lung CT DIR, this procedure results in approximately 250–300 knots. An additional point cloud was similarly acquired with approximately 7 mm uniform spacing to serve as the subdomain locations for Eq. (9). This resulted in approximately 20 000 to 30 000 subdomain measurement points. Each initial subdomain was defined as a single voxel and then morphologically dilated with a 7 × 7 × 3 voxel structuring element until the tolerance criteria [Eq. (5) for IJF and Eqs. (6) and (7) for MCVC] were satisfied.

The spatial transformation is computed by applying the Quadratic Penalty DIR (QPDIR) algorithm to the inhale and exhale images. Briefly, QPDIR is an intensity‐based algorithm designed around a gradient‐free block coordinate descent strategy that essentially iterates between simple block matching and linear least squares operations to minimize the structural similarity index between an image pair. The implementation follows the description in Ref. [18], with the exception that an additional sum‐of‐squared difference term was included in the QPDIR objective function to improve lung mask alignment (as done in Ref. [19]).

Equation (13) is solved using the interior‐point method implemented in the MATLAB (release R2019a, The Mathworks Inc, Natick, Massachusetts, USA) optimization routine lsqlin.

Image data

A validation of the MLS IJF and MCVC methods was conducted using the simulation (treatment planning) 4DCT and single photon emission tomography (SPECT) ventilation images for 15 NSCLC patients who received radiotherapy at our institution. Data were retrospectively evaluated according to an IRB approved study (IRB 2016‐037, clinicalTrials.gov #NCT02528942). Patients received definitive radiotherapy (defined as prescription doses of 45 to 75 Gy) and a planned concurrent chemotherapy regimen. The majority of enrolled patients had stage III disease. Single photon emission tomography imaging was acquired prior to delivery of the first radiotherapy fraction, ensuring no clinical lung function changes occurred due to treatment between the simulation 4DCT and SPECT acquisitions.

A Philips Brilliance Big Bore CT (version 3.6.7) with a bellows system was used for respiratory correlated imaging. The 4DCT images were acquired with x‐ray tube settings of 120 kVp and 599 mAs, and reconstructed using phase binning to produce an average CT image and 10 phase indexed CT images. The phase images were indexed from 0% to 90% in steps of 10% where 0% indicates full inhalation and 50% indicates full exhalation on the breathing curve. Final images were then exported into the DICOM (Digital Imaging and Communications in Medicine) standard (512 x 512 pixels per 2D slice image, voxel dimensions of 1.27 mm × 1.27 mm × 3 mm).

SPECT ventilation images were acquired on a dual head Siemens Symbia SPECT/CT scanner (Siemens Medical Solutions, USA), using a parallel hole, high‐resolution collimator and an energy window of 15% at a centerline of 140 keV. 1 mCi Tc99m diethylenetriaminepentaacetic acid aerosol inhalation was used for the ventilation scan. SPECT acquisition was performed in steps of 6° for the entire 360° of rotation, with a 25 s collection time for each step. Total scan acquisition time, for most patients, was under 30 min. Free breathing, attenuation‐corrected CT images were subsequently recorded with 130 kVp, and 75–100 mAs (weight dependent) during continuous tidal respiration. The final reconstructed SPECT ventilation images were then exported into DICOM (64 × 64 pixel per 2D slice image, voxel dimensions 6.00 mm × 6.00 mm × 2.00 mm).

Uncertainty tolerance parameter sweep

A parameter sweep was conducted for both the MCVC and IJF methods in order to assess the effect of the uncertainty tolerance on the spatial correlation between SPECT ventilation and the resulting robust CT‐derived ventilation images. For each of the 15 test cases, a series of 30 MCVC and IJF CT‐ventilation images were computed using a uniformly sampled set of uncertainty tolerances ranging between and converted to absolute volume difference :

Each resulting image was spatially aligned with the SPECT ventilation image by first using affine registration to align the exhale 4DCT phase (on which the CT‐ventilation is computed) and the SPECT attenuation correction CT. The resulting affine transformation was then applied to the CT‐ventilation image and the voxel‐wise Spearman correlation was computed at the resolution of the SPECT ventilation, after applying a median filter with structuring element to the SPECT image (as done in Ref. [20]).

RESULTS

The median Spearman correlations taken across all 15 test cases for each uncertainty parameter value are presented in Figs. 1 and 2. The median SPECT‐V and MCVC correlation values ranged between 0.20 and 0.48 across the parameter sweep, with the higher median values being achieved for The median SPECT‐V and IJF correlation values ranged between 0.79 and 0.82 across the parameter sweep and were relatively consistent, though the highest values were achieved for

Fig. 1

The maximum, 75th percentile, median, 25th percentile, and minimum Spearman correlations across all 15 cases as a function of the uncertainty parameter for the IJF method. The correlation values remain relatively constant across the parameter sweep. [Color figure can be viewed at wileyonlinelibrary.com]

Fig. 2

The maximum, 75th percentile, median, 25th percentile, and minimum Spearman correlations across all 15 cases as a function of the uncertainty parameter for the MCVC method. After variations for smaller , the correlation values remain relatively constant between 0.05 and 0.15. [Color figure can be viewed at wileyonlinelibrary.com]

As summarized in Table I, the IJF and SPECT‐V Spearman correlations across the 15 test cases for the optimal ranged between 0.12 and 0.90. For MCVC, an optimal generated Spearman correlations that ranged between −0.06 and 0.84 across the 15 test cases. Wilcoxon signed‐rank test applied to the optimal IJF and MCVC correlations in Table I indicate that the IJF method outperforms MCVC with high statistical significance (P = 0.0062). Figure 3 illustrates the case with overall highest correlation between the two methods (Case 9).

Table I

Spearman correlations for optimal .

Case #	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	Median
IJF τ = 0.07	0.12	0.86	0.73	0.90	0.87	0.64	0.82	0.79	0.84	0.44	0.83	0.90	0.37	0.25	0.82	0.82
MCVC τ = 0.03	0.05	0.48	0.62	0.10	0.75	0.50	0.52	0.08	0.77	0.84	−0.06	0.33	0.20	0.06	0.54	0.48

The Spearman correlations between SPECT ventilation and both IJF and MCVC across all 15 test cases are listed, using the optimal uncertainty tolerances, as determined by the parameter sweep (Figs. 1 and 2), for each method. In all cases, the correlation is significantly different from zero (P < 0.001).

Fig. 3

The MCVC (top row) and IJF (bottom row) ventilation images computed using (left column) and (right column) for the test case with the overall highest correlation (Case 9, Table I) with the corresponding SPECT‐V (far right). The images were converted to percentile images for direct visual comparison. All images denote decreased function in the left upper lung. The Spearman correlation between SPECT‐V and MCVC for and are 0.79 and 0.77, respectively. For IJF, the correlations are both 0.85.

DISCUSSION

Integrated Jacobian Method and MCVC represent a new class of robust CT‐ventilation methods that seek to reduce the uncertainty associated with DIR variability. Mass conserving volume change has the added benefit of not requiring the advanced pulmonary vasculature mask needed by previous intensity‐based methods. However, the previous numerical implementations of robust methods required the solution of a constrained, large scale, linear least squares problem defined by an underdetermined data fidelity term and a spatial smoothness inducing regularization. In this work, we introduce a parameterized variant of the original IJF and MCVC methods that allows for overdetermined data fitting and the use of common optimization solvers. In addition, robust CT‐ventilation methods require the definition of an uncertainty tolerance parameter. The uncertainty tolerance reflects the methodology's general property that robustness is gained at the expense of measurement resolution. As such, the parameter's selection represents the trade‐off between data fidelity and numerical stability. Previous work provided general bounds on parameters values that maintained DIR robustness. In this study, we characterize the effect the parameter has on the resulting correlation between robust CT‐ventilation and SPECT‐V.

Across the parameter sweep, the IJF method demonstrated less variability than MCVC. This reflects the inherent characteristics of intensity‐ and transformation‐based methods. Whereas, intensity‐based methods operate on nonhomogeneous image data, transformation‐based methods operate on the DIR displacement field, which is often optimized for spatial smoothness. This is the case for the QPDIR method utilized in this study. As indicated by the Fig. 1 results, accurate estimates for the subregional IJF sample means can, in practice, be acquired with fewer data points than required by the Eq. (5) bound. As such, the IJF median correlation appears to be more stable with respect to when compared to those of the MCVC method. Moreover, IJF demonstrated a significantly higher correlation with SPECT‐V than MCVC (P = 0.0062). This result is not altogether surprising considering the underlying mathematical formulations on which the two methods are based. Integrated Jacobian formulation recovers the DIR Jacobian values from the geometric information captured by the DIR solution. The MCVC method also approximates the Jacobian. But similar to other HU‐based methods, MCVC is derived under the assumption that the density variations observed between inhale and exhale CT are caused solely by changes in air content. This assumption is known to be invalid due to variations in pulmonary blood mass that occur throughout the breath cycle and is a potential source of error for MCVC.

The 0.48 median correlation generated by the optimal tolerance value for MCVC (Table I) across the 15 test cases is consistent with the results from the Ventilation and Medical Pulmonary Image Registration Evaluation (VAMPIRE) study, where the highest performing method achieved a median of 0.49. However, the 0.82 median correlation generated by IJF represents a significant improvement over previously reported results, though the comparison with VAMPIRE is not directly applicable since its analysis included Galligas PET and Xenon CT. Validation based on nuclear medicine imaging is in general challenging due to the lower spatial resolution and aerosol deposition artifacts common to SPECT ventilation. Future work includes conducting further validation with other functional imaging modalities, including PET Galligas, hyperpolarized gas MRI (as done in Refs. [22, 23]), and data from the VAMPIRE study.

The mathematical framework governing the IJF and MCVC methods guarantee that similar DIR solutions will generate similar ventilation images. This property was derived and demonstrated numerically in our previous studies. , Considering the correlations reported in Table I, robust ventilation algorithms provide a numerical framework that (a) is reproducible with respect to image processing pipeline and (b) generates results consistent with SPECT ventilation. However, the framework does not guard against 4DCT acquisition artifacts, which can lead to substantial differences between DIR solutions generated from competing methods. As illustrated in Fig. 4, phase‐binning 4DCT artifacts introduce erroneous geometric and intensity information into the Eq. (13) calculation, which can corrupt the resulting IJF and MCVC ventilation images. CT‐ventilation computed from breathhold CT acquisitions, which do not require phase binning, should therefore be expected to produce higher quality results. However, for 4DCT, an area of future work is in extending the IJF mathematical framework to be incorporated into both the DIR and segmentation algorithm, with the goal of reducing the effect of phase‐binning artifacts.

Fig. 4

Top Row: The 4DCT‐Inhale phase (left) and SPECT (right) ventilation images for the case with the lowest over correlation in Table I (Case 1). Bottom Row: The IJF (left) and MCVC (right) superimposed on the 4DCT‐Exhale phase. Ventilation images were converted to percentile images for direct visual comparison. The inhale image possesses a phase‐bin artifact (blue arrows), which erroneously elevates the volume change signal in the lower right lung while erroneously suppressing it in the lower left. This leads to poor correlation with SPECT‐V.

CONCLUSION

We have assessed the spatial correlation between SPECT‐V and two robust CT‐ventilation methods that are dependent on the definition of an uncertainty parameter . Results demonstrate that the proposed numerical implementation yields consistent results for uncertainty parameter values ranging between , with the highest correlations being achieved with for IJF and for MCVC. While MCVC demonstrated correlations consistent with those previously reported in the literature, the IJF method's SPECT‐V correlations represent a significant improvement over previous results, despite known issues with SPECT validation. Immediate future work involves comparisons with other ventilation modalities, such as PET Galligas and Hyperpolarized gas MRI, that are generally considered to be superior to SPECT.

CONFLICT OF INTEREST

The authors have no conflict to disclose.