Robust CT ventilation from the integral formulation of the Jacobian.

Computed tomography (CT) derived ventilation algorithms estimate the apparent voxel volume changes within an inhale/exhale CT image pair. Transformation-based methods compute these estimates solely from the spatial transformation acquired by applying a deformable image registration (DIR) algorithm to the image pair. However, approaches based on finite difference approximations of the transformation's Jacobian have been shown to be numerically unstable. As a result, transformation-based CT ventilation is poorly reproducible with respect to both DIR algorithm and CT acquisition method.
PURPOSE: We introduce a novel Integrated Jacobian Formulation (IJF) method for estimating voxel volume changes under a DIR-recovered spatial transformation. The method is based on computing volume estimates of DIR mapped subregions using the hit-or-miss sampling algorithm for integral approximation. The novel approach allows for regional volume change estimates that (a) respect the resolution of the digital grid and (b) are based on approximations with quantitatively characterized and controllable levels of uncertainty. As such, the IJF method is designed to be robust to variations in DIR solutions and thus overall more reproducible.
METHODS: Numerically, Jacobian estimates are recovered by solving a simple constrained linear least squares problem that guarantees the recovered global volume change is equal to the global volume change obtained from the inhale and exhale lung segmentation masks. Reproducibility of the IJF method with respect to DIR solution was assessed using the expert-determined landmark point pairs and inhale/exhale phases from 10 four-dimensional computed tomographies (4DCTs) available on www.dir-lab.com. Reproducibility with respect to CT acquisition was assessed on the 4DCT and 4D cone beam CT (4DCBCT) images acquired for five lung cancer patients prior to radiotherapy.
RESULTS: The ten Dir-Lab 4DCT cases were registered twice with the same DIR algorithm, but with different smoothing parameter. Finite difference Jacobian (FDJ) and IFJ images were computed for both solutions. The average spatial errors (300 landmarks per case) for the two DIR solution methods were 0.98 (1.10) and 1.02 (1.11). The average Pearson correlation between the FDJ images computed from the two DIR solutions was 0.83 (0.03), while for the IJF images it was 1.00 (0.00). For intermodality assessment, the IJF and FDJ images were computed from the 4DCT and 4DCBCT of five patients. The average Pearson correlation of the spatially aligned FDJ images was 0.27 (0.11), while it was 0.77 (0.13) for the IFJ method.
CONCLUSION: The mathematical theory underpinning the IJF method allows for the generation of ventilation images that are (a) computed with respect to DIR spatial accuracy on the digital voxel grid and (b) based on DIR-measured subregional volume change estimates acquired with quantifiable and controllable levels of uncertainty. Analyses of the experiments are consistent with the mathematical theory and indicate that IJF ventilation imaging has a higher reproducibility with respect to both DIR algorithm and CT acquisition method, in comparison to the standard finite difference approach.

Introduction

Computed tomography (CT) derived ventilation imaging is an image processing‐based modality, where image segmentation and deformable image registration (DIR) methods are applied to the temporally resolved phases of a four‐dimensional CT (4DCT) image set or breath‐hold CT pair in order to infer the lung voxel volume changes induced by respiratory motion.1 Since its inception, CT ventilation has been the subject of numerous validation studies, which primarily focused on comparison with established modalities such as 99mTc‐DTPA (diethylenetriamine pentaacetate) SPECT‐ventilation,2 99mTc‐MAA SPECT‐perfusion,3 68Ga‐MAA PET perfusion, 68Ga‐nanoparticles PET ventilation,4 and3He‐hyperpolarized MRI.5 The relative success of these studies has led to the development of applications, mainly in radiation oncology, such as functional avoidance radiotherapy planning6 and treatment response modeling.7

There are two approaches for computing CT ventilation: intensity based and transformation based.2 Intensity‐based methods estimate volume changes directly from the Hounsfield units of spatially corresponding inhale and exhale voxels according to the formulation proposed by Simon.8 However, with transformation‐based methods, estimates are derived from the Jacobian factor of the DIR recovered spatial transformation.9 Variations and combinations of the two approaches have also been proposed.10

Regardless of approach, CT ventilation is computed with image processing methods. As such, its physiological accuracy and clinical utility highly depend on the performance of the employed numerical algorithms. CT ventilation methods therefore face the same challenges arising in the field of scientific computing, where the goal is to develop numerical methods based on mathematical models of a physical phenomenon. Software implementations require both verification and validation of the numerical method and mathematical models before being confidently employed in practice.11 Despite the large number of validation studies, surprisingly little research has focused on the verification of CT ventilation, which in this case, refers to the characterization of the numerical approximation errors associated with the computation.

Recent work has shown that current methods for transformation‐based ventilation, which depend on finite difference approximations applied to the DIR‐recovered displacement field, are inherently nonstable numerically, meaning that small perturbations to the DIR solution can cause large relative changes in the resulting ventilation image.12 This property partly explains the reported sensitivity of CT ventilation with respect to DIR method13 and 4DCT acquisition artifacts14; in addition to its overall poor reproducibility.15 The purpose of this study is to introduce the Integrated Jacobian Formulation (IJF) method for computing DIR‐measured volume changes. The IJF method is based on robust estimates of regional volume change, which are computed by applying a sampling method to numerically integrate the regional Jacobian formulation. The resulting volume change estimates (a) respect the restrictions posed by the resolution of the digital grid and (b) have a quantifiable and controllable level of uncertainty. A Jacobian image, which depicts the volume change for each voxel within the reference image, is computed from the regional estimates by solving a constrained linear least squares problem.

This current work describes the motivation, mathematical framework and numerical implementation for the novel IJF method. Two sets of experiments are used to assess the performance of the IJF method with respect to (a) variable DIR solutions and (b) the correlation between ventilation images generated from 4DCT and 4D cone beam CT (4DCBCT) acquisitions.

Materials and methods

The Jacobian and finite difference approximations

The goal of CT ventilation methods is to estimate the volume changes, induced by the changes in air content within the lungs (breathing), for each voxel within a defined lung region of interest. As opposed to the HU‐based approach, transformation‐based CT ventilation approximates volume changes solely from a spatial transformation, ϕ, that describes the apparent lung motion within an inhale/exhale CT image pair. The spatial transformation is often described in terms of the voxel displacement field:where represent the regions of interest within the reference and target images, respectively. A standard assumption when computing ventilation images is the ability to generate both inhale and exhale lung masks (segmentations), which implies that Ω(, Ω( are both known and contain the same lung tissue.2 However, the true spatial transformation is not known and is approximated numerically by applying a DIR algorithm. The resulting DIR solution provides displacement vector estimates for the voxels in Ω(.

The Jacobian, J(x), is the first full derivative of ϕ:and it describes the volume scaling factor under:where Ω ⊆ Ω( and ϕ(Ω) ⊆ Ω( is the image of Ω under ϕ. For medical imaging applications, ϕ is assumed to be diffeomorphic, which implies:

Computing a volumetric transformation‐based CT ventilation image requires numerical approximation of the integral in Eq. (2) for a series of different subvolumes within Ω (. Current methods take Ω to be the unit voxel volume centered on .2 This allows for a volumetric midpoint rule approximation of the Eq. (2) integral:

Since voxels are spaced uniformly on the digital grid, the partial derivatives given by Eq. (1) can be approximated with a finite difference scheme. We refer to this method as the finite difference Jacobian (FDJ):where the entries in are obtained with forward difference approximations:and e is the standard basis vector. While the FDJ method has proven useful in qualitative applications, its poor reproducibility, with respect to both DIR method and CT acquisition, is well documented13, 15 and remains a limiting factor for applications dependent on precise quantitative analysis.

As demonstrated in Ref. 12, small perturbations to the DIR displacement field can cause large relative changes in the FDJ volume change estimate. The result follows from the fact that the general finite difference approximation error is O(h), where h is equal to the grid spacing. In the case of the digital voxel grid, h = 1, which implies that the uncertainty in the finite difference approximation in Eq. (6) is ±1. Intuitively, this uncertainty stems from the fact that DIR spatial accuracy resolution is limited to the digital grid. In other words, using expert‐determined landmarks, it is possible to determine if a reference voxel is mapped into the correct target voxel volume. However, it is not possible to determine where within the target voxel the correctly mapped position lies. Consequently, finite difference approximations can vary by as much as 1.0 for displacement estimates with the exact same spatial accuracy, as illustrated by Fig. 1.

Figure 1

Two transformations map the reference positions x and x + e (left) into the voxel volumes centered on (i,j) and (i,j + 1), respectively (right). While the transformations have the same spatial accuracy (i.e., map the reference positions into the same target voxel volumes), the forward difference approximations for ∂ ϕ/ ∂ x 1 and are 1.0 and 2.0, respectively. This reflects the O(1) error associated with finite differences on the voxel grid. [Color figure can be viewed at http://wileyonlinelibrary.com]

This result is problematic since inhale‐to‐exhale lung motion is primarily a contraction, meaning that the relative voxel volume change values are in general less than or equal to 1.0, and therefore, well within the error margin of the FDJ approximation.

With these issues in mind, our goal is to devise a numerical method for computing voxel volume changes, derived from a DIR spatial transformation, that is (a) computed with respect to the resolution of the digital grid and is (b) based on approximations with quantitatively characterized levels of uncertainty.

Robust regional volume change estimation

For a general subdomain Ω ∈ Ω(, Eqs. (2), (3), (4) imply:where the discretized variablesreflect that volume change is strictly positive.

The FDJ method directly estimates each v with finite differences and O(1) uncertainty. Rather than depending on finite differences to estimate volume changes for individual voxels, we propose robust regional volume change measurements by approximating vol(ϕ(Ω)) using the well‐known “hit‐or‐miss” algorithm.16

Sampling‐based numerical integration

The hit‐or‐miss algorithm is a sampling method that can be viewed theoretically as a simpler variant of the Monte Carlo integration method.16 Moreover, the method does not require an explicit representation of ϕ(Ω) to estimate vol(ϕ(Ω)). Since ϕ is assumed to be diffeomorphic, ϕ −1 exists and can be used to define a voxel membership oracle.17 In this case, the oracle answers the question: “Given x ∈ Ω(, is x ∈ ϕ(Ω)?” Mathematically, the membership oracle, ∀x ∈ Ω(, is defined as:

Since Ω represents a region of voxel volumes within Ω(, the oracle essentially operates on the components of ϕ −1 rounded to the nearest integer values. Or more specifically,where [i] represents the component‐wise rounding operation. As such, the oracle evaluations are equivalent for two DIR solutions with the same spatial accuracy.

The hit‐or‐miss method first computes the sampled average of the oracle function over an enclosing region with known volume. Since Ω( is known and ϕ(Ω) ⊆ Ω(, the estimated function average, , is given as:where N = |Ω(| and H is the number of voxels that “hit” ϕ(Ω). The corresponding volume estimation is then computed as the ratio of hits to samples multiplied by the volume of the enclosing region:

Combining Eqs. (7) and (12) results in a simple linear equation that relates individual voxel volume changes within Ω to the hit‐or‐miss volume approximation:

Uncertainty in regional volume change estimation

Since N is assumed to be large, similar to Monte Carlo error analysis,18 adopting Gaussian statistics allows for the uncertainty in the Eq. (11) average function value, and consequently the Eq. (12) volume estimation, to be characterized by the standard error of the mean:where s is the sample standard deviation. Applying Eq. (14) to the estimated mean in Eq. (11) yields:

Thus, the uncertainty in the Eq. (12) hit‐or‐miss volume estimate can be expressed as:

Equation (16) provides a probabilistic assessment of the volume approximation in terms of a confidence interval defined by the parameter β. For instance, setting β = 1.96 implies that with 95% probability:

General Monte Carlo methods reduce the uncertainty by raising the number of point samples, which implies the oracle function is precise for any sampling resolution. In this case, the problem is inherently discretized by the digital grid. Consequently, both DIR spatial accuracy and the Eq. (9) oracle are limited to the grid resolution. Therefore, the number of samples N = |Ω(| is fixed and the Eq. (17) uncertainty is determined solely by the number of hits H. Figure 2 shows that the relative uncertainty function for fixed N:is monotonically decreasing with respect to H. For H = 1, the uncertainty is approximately 100%, which mimics the O(1) error of forward finite differences. Conversely, the most trusted volume measurement is that of the complete target segmentation mask: N = |Ω(|.

Figure 2

Relative uncertainty in hit‐or‐miss volume change estimate as a function of the number of hits (x‐axis) for fixed N = 500 and β = 1.96. [Color figure can be viewed at http://wileyonlinelibrary.com]

Given a β ∗ corresponding to a specified confidence level and a desired relative error τ, there exists a unique H* ≤ N such that:This result follows from the fact that E is monotonically decreasing and E(N, β ∗) = 1.

Integrated Jacobian formulation

A full CT ventilation image, V(x), requires computing the discretized variables

Equation (13) represents a regional DIR‐based volume change measurement, the uncertainty of which is governed by Eq. (17), and provides one equation for M = |Ω(| total unknowns. Intuitively, a recovery algorithm can be built on the idea of computing enough regional estimates to define a full rank or overdetermined linear system of equations for the M unknowns. Applying Eq. (13) to a series of subregions results in the following linear system:whereand the elements of b are computed according to Eq. (12):

If the subregions are allowed to overlap, Eq. (21) is similar in structure to the standard image deblurring problem.19 Although in this case, the volume change image recovery problem is defined on an irregular domain, that is, the lung region of interest, while deblurring problems are typically defined on the full image domain.

The Eq. (21) subregions can be constructed to satisfy the uncertainty bound given by Eq. (19), but the resulting data values are still approximations. As such, Eq. (21) should not be treated as a hard constraint. However, global volume change, as measured by the volumes of the inhale/exhale masks, is often used as a validation metric for proposed ventilation methods.1, 2, 3, 20, 21 Thus, the global mass change:can be enforced as a hard constraint. Moreover, Eqs. (21) and (24) are not guaranteed to provide enough information to uniquely determine v. Therefore, a two‐norm penalization on the spatial gradient of V(x) is employed to both induce smoothness in the reconstructed image and guarantee the existence of a unique solution (i.e., regularize the problem). Combining these assumptions results in the following constrained linear least squares problem, the solution of which is the recovered volume change image:

The penalty parameter α dictates the amount of smoothness in the solution v* and the lower bound ϵ > 0 represents the minimum possible voxel volume. We refer to Eq. (25) as the IJF method for computing DIR‐induced volume changes and transformation‐based CT ventilation.

Numerical implementation

The IJF method requires constructing the subregion volume change estimates described in Eq. (21). To do this, we applied a dart‐throwing algorithm to x  ∈ Ω( in order to generate a point cloud of voxel locations uniformly spaced approximately 30 mm apart.22, 23 The Voronoi diagram generated from the point cloud and Ω( is used to define the initial subregions (Fig. 3). For each initial subregion, the Eq. (12) hit‐or‐miss estimate and the Eq. (18) uncertainty are computed. The initial subregion mask is morphologically dilated with a 5 × 5 × 5 structuring element until

Figure 3

An example of an initial subdomain decomposition (left) is defined as the Voronoi diagram generated from a voxel point cloud. The initial subregions are randomly colored for display purposes. The integrated Jacobian formulation method takes each individual subregion (middle image shows an example) and performs a morphological dilation until Eq. (26) is satisfied (right). [Color figure can be viewed at http://wileyonlinelibrary.com]

Equation (26) implies that the relative hit‐or‐miss uncertainty is <2.5% with 95% confidence. For all experiments, the Eq. (25) regularization parameter α = 0.5. Numerical experimentation and visual inspection indicated that solutions did not vary greatly for α between 0.1 and 0.5. For smaller values, the problem could potentially become ill‐conditioned, thereby resulting in possible numerical errors in the least squares solutions to Eq. (25). For larger values, the resulting solutions have the potential to become dominated by the regularization term, thereby resulting in oversmoothed images.

For the peak inhale and exhale 4DCT phases, lung masks, Ω(, Ω( were generated using a semi‐automated histogram segmentation (as done in Ref. 24). The 4DCBCT inhale/exhale phases were segmented manually by an imaging expert. The spatial transformation ϕ −1 is computed by applying the MILO DIR algorithm (see Ref. 24 for details) to the target and reference image. Eq. (25) is solved using the Augmented Lagrangian Method (see for 25 for details), implemented in MATLAB (release R2017b, The Mathworks Inc, Natick, MA, USA).

Image data

Two sets of experiments were conducted to assess the performance of the proposed IJF method. For the first experiment, ten publicly available 4DCTs from the http://www.dir-lab.com website (see Ref. 26, 27, 28 for acquisition details) were used to assess the reproducibility of the IJF method with respect to DIR solution.

For the second experiment, the simulation (treatment planning) 4DCT and setup verification 4DCBCT images for five NSCLC patients receiving radiotherapy at our institution were used to assess the robustness of the IJF method. Data were retrospectively evaluated according to an IRB‐approved study (IRB 2016‐037). 4DCBCTs were acquired 6 to 8 days after simulation and immediately prior to the delivery of the first radiotherapy fraction, ensuring no lung function changes occurred due to treatment between the simulation 4DCT and 4DCBCT acquisitions. The rationale for choosing 4DCT and 4DCBCT to demonstrate the robustness of the proposed method is that, due to the scatter artifact in CBCTs, CBCT images provide the most challenging scenario for comparing the generated CT ventilation images. Kipritidis et al. reported Spearman correlations between normalized intrafraction ventilation images generated from 4DCBCTs on the order of 0.6.29

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 × 512 pixels per two‐dimensional slice image, voxel dimensions of 1.27 mm × 1.27 mm × 3 mm).

A 4D VolumeView™ imaging sequence using Symmetry™ by Elekta (Elekta, Stockholm, Sweden) XVI (release version 5.0.2) was used to acquire 4DCBCT for treatment setup verification. For each 4DCBCT image, approximately 975 projection images were acquired, over a gantry rotational range of 200° with rotational speed of 67°/min. The S20 collimator cassette was used to define the field of view and imaging parameters were set to 120 kVp, 20 mA, 16 ms. The final 4DCBCT reconstructed images were then exported from XVI as a set of 10 phase indexed CT images in the Dicom format (135 × 135 × 42 voxels, voxel dimensions of 2 mm × 2 mm × 6 mm).

Results

The variability in the IJF method with respect to DIR solution was assessed using ten peak inhale and exhale 4DCT phases from the DIR‐LAB repository. Each image pair was registered twice using the same DIR algorithm (MILO24), but with different spatial smoothness parameters (MILO parameter γ = 0.5 and 0.6). The spatial accuracies of the two DIR solutions for each case are detailed in Table 1. The Pearson correlation coefficient between the pair of FDJ images generated from the two DIR solutions was computed for all cases. The same process was applied to the IJF generated images. No normalizations or postprocessing procedures were applied to the generated ventilation images. As detailed in Table 2, the average correlation values for the FDJ method was 0.83 (0.03), while for the IJF method the average correlation was 1.00 (0.00). These results are consistent with the FDJ numerical analysis findings detailed in Ref. 12, as well as in Sections 2 and 4. IJF operates on the rounded DIR displacement values [see Eq. (10)], which implies that the IJF Jacobian values are equivalent for DIR solution with the same spatial accuracy. The images for 4DCT 6, which provide a representative example of the Table 2 results, are displayed in Figure 4.

Table 1

Integrated Jacobian formulation (IJF) Subregion Summaries. Ten four‐dimensional computed tomography (4DCT) inhale/exhale phase test cases provided by http://www.dir-lab.com were used to assess the reproducibility of IJF and finite difference Jacobian methods with respect to registration accuracy (Table 2). A summary of the subregional domains used within the IJF method for deformable image registration (DIR) 1 (Table 2), including: the number of voxels contained in the target (inhale) image domain , the minimum number of hits, H , required to satisfy Eq. (26), the number of subregions, K, as defined by Eq. (21), and the minimum/maximum number of voxels within the K subregions

	|Ω^(T)|	H*	K	min|Ω_k|	max|Ω_k|
4DCT 1	972 538	6108	642	5160	13 370
4DCT 2	1 537 594	6123	976	5341	13 546
4DCT 3	1 191 488	6116	774	5321	14 201
4DCT 4	881 357	6104	558	4510	13 188
4DCT 5	1 063 434	6112	696	5112	12 855
4DCT 6	1 194 354	6116	666	4491	13 012
4DCT 7	1 458 702	6121	869	4452	13 210
4DCT 8	2 000 263	6128	1175	4401	13 145
4DCT 9	653 069	6090	445	4350	13 112
4DCT10	1 155 485	6115	727	4258	14 466

Table 2

Dataset four‐dimensional comouted tomography (4DCT): ten 4DCT maximum inhale/exhale phase pairs

	DIR 1 Avg. mm error	DIR 2 Avg. mm error	Mean absolute difference		Pearson correlation
			FDJ	IJF	FDJ	IJF
4DCT 1	0.72 (0.91)	0.71 (0.90)	0.03 (0.03)	0.00	0.88	1.00
4DCT 2	0.70 (0.90)	0.69 (0.90)	0.02 (0.03)	0.00	0.87	1.00
4DCT 3	0.91 (1.07)	0.91 (1.07)	0.05 (0.05)	0.00	0.80	1.00
4DCT 4	1.26 (1.28)	1.26 (1.23)	0.06 (0.06)	0.00	0.86	1.00
4DCT 5	1.18 (1.49)	1.29 (1.49)	0.07 (0.06)	0.00	0.83	1.00
4DCT 6	0.97 (0.96)	1.09 (1.00)	0.06 (0.05	0.00	0.80	1.00
4DCT 7	0.90 (0.98)	1.03 (0.99)	0.06 (0.05)	0.00	0.79	1.00
4DCT 8	1.10 (1.24)	1.14 (1.23)	0.05 (0.05)	0.00	0.82	1.00
4DCT 9	1.02 (0.97)	1.02 (0.96)	0.06 (0.06)	0.00	0.83	1.00
4DCT10	0.98 (0.98)	1.02 (1.03)	0.06 (0.07)	0.00	0.80	1.00
Avg. (std):					0.83 (0.03)	1.00 (0.00)

Ten 4DCT inhale/exhale phase test cases provided by http://www.dir-lab.com were registered by two versions of the MILO deformable image registration (DIR) algorithm. The spatial accuracies achieved by the two DIR approaches are given in average (SD) millimeter error with respect to 300 landmark point pairs. The Pearson correlation and mean absolute difference between the FD ventilation images computed from DIR 1 and DIR 2 are given, as are the correlations between the IJF ventilation images computed from DIR 1 and DIR 2.

Figure 4

The finite difference Jacobian (FDJ) images (top) and the integrated Jacobian formulation (IJF) images (bottom) computed for the same case with two deformable image registration (DIR) solutions (Table 1, four‐dimensional computed tomography 6). While there is a large visual difference between the FDJ images, the IJF images are nearly identical. The estimated Jacobian values indicate that the inhale‐to‐exhale lung motion recovered by the DIR algorithm is a contractive mapping.

Variability across modalities was assessed using 4DCT and 4DCBCT images for five lung cancer patients being treated with radiotherapy. Both the 4DCT and 4DCBCT were acquired prior to radiotherapy. For each patient, the FDJ and IJF images were computed from the 4DCT and 4DCBCT. No normalizations or postprocessing procedures were applied to the generated ventilation images. Affine registration was employed to spatially align the ventilation images, using the maximum inhale phase from the lower resolution 4DCBCT as the reference coordinate system. Pearson correlation coefficients between the spatially aligned IJF images generated from the 4DCB and 4DCT were computed, as they were for the FDJ images as well (Table 3). The IJF demonstrated good reproducibility between the 4DCT and 4DCBCT acquisitions with an average Pearson correlation of 0.77 (0.13). In contrast, the average correlation for the FDJ images was 0.27 (0.11). The highest IJF correlation of 0.95 was achieved in Case 1, which is illustrated in Fig. 5. The poorest correlation value occurred on Case 2, where a pronounced phase binning artifact is present near the diaphragm on the 4DCT image. Such an artifact is not present in the corresponding 4DCBCT (see Fig. 6).

Table 3

Four‐dimensional computed tomography (4DCT) vs 4D cone beam CT Jacobian values

	Mean absolute difference		Pearson correlation
	FDJ	IJF	FDJ	IJF
Case 1	0.11 (0.09)	0.03 (0.02)	0.41	0.95
Case 2	0.13 (0.10)	0.06 (0.04)	0.12	0.61
Case 3	0.17 (0.15)	0.05 (0.04)	0.21	0.68
Case 4	0.10 (0.08)	0.04 (0.03)	0.32	0.82
Case 5	0.19 (0.17)	0.05 (0.03)	0.29	0.80
Avg. (Std):			0.27 (0.11)	0.77 (0.13)

Finite difference Jacobian (FDJ) and integrated Jacobian formulation (IJF) images were created from the treatment planning (simulation) 4DCTs and the four‐dimensional cone beam CT images of five lung cancer patients. All images were acquired prior to radiotherapy. The high correlation values between the CB and CT IJF values indicate that the IJF is more reproducible with respect to acquisition modality than FDJ.

Figure 5

Finite difference Jacobian (FDJ) (top) and integrated Jacobian formulation (IJF) (bottom) images computed from the four‐dimensional cone beam computed tomograhy (4DCBCT) and 4DCT images for the same patient. (Table 3, Case 1), superimposed on the inhale CBCT phase. The 4DCT FDJ and IJF images were mapped onto the inhale CBCT phase via affine registration. While there is a large difference between the FDJ images (Pearson correlation 0.41), the IJF images are very similar (Pearson correlation 0.95). The estimated Jacobian values indicate that the inhale‐to‐exhale lung motion recovered by the deformable image registration algorithm is a contractive mapping.

Figure 6

The Integrated Jacobian Formulation image created from the four‐dimensional cone beam computed tomograhy (Left) and radiotherapy planning 4DCT (right) for a lung cancer patient (Table 2, Case 2). The middle image depicts the CT‐integrated Jacobian formulation image mapped onto the inhale CBCT phase. The 4DCT (right) contains a large phase binning artifact at the diaphragm, resulting in significant variation between the CB and CT IJF images. Table 2, Case 3 possess a similar phase binning artifact. The estimated Jacobian values indicate that the inhale‐to‐exhale lung motion recovered by the deformable image registration algorithm is a contractive mapping.

Discussion

The purpose of this work is to introduce a robust method for computing transformation‐based CT ventilation. The proposed methodology addresses the CT ventilation problem from an engineering/applied mathematics perspective, which notes that a mathematical model and its corresponding numerical implementation must undergo verification and validation before the methodology can be confidently employed in practice.11 Having defined a mathematical model, verification is the process of developing a numerical implementation that accurately solves the mathematical model. Validation is the process of evaluating whether the mathematical model and its verified numerical implementation accurately describe a physical phenomenon of interest. In this case, the physical quantity we wish to describe is pulmonary ventilation and the mathematical model is the Jacobian factor. The focus of this study is on the verification portion of the transformation‐based CT ventilation model. While several physiological factors contribute to lung volume change, the DIR recovered mapping only describes the apparent motion. Therefore, the volume changes recovered with IJF, or any other Jacobian‐based method, describe the concerted effect of all factors. As such, the proposed method's clinical utility cannot be directly inferred from this work on verification, and must be assessed separately with a rigorous validation study. This is an immediate area of our future research.

Current methods for transformation‐based volume change estimation rely on Jacobian approximations taken on individual voxels. As reported in the literature, the resulting ventilation images are not well reproducible with respect to either DIR algorithm or CT acquisition method.

The fundamental problem for FDJ ventilation is that accurately estimating the Jacobian for an individual voxel with finite differences requires a numerical precision that is higher than what DIR spatial accuracy assessments can measure. As a result, two DIR solutions with similar spatial accuracies can produce substantially different ventilation images, which is consistent with the findings in Ref. 12, 30. The standard error analysis for hit‐or‐miss integral estimation is also consistent with these findings. Consequently, the utility of CT ventilation has been limited to applications that do not rely on voxel‐by‐voxel accuracy of the generated CT ventilation images. For example, functional avoidance radiotherapy, which is one of the main clinical applications of CT ventilation, generally relies on regional/lobar accuracy rather than accuracy on a voxel level.31, 32, 33, 34, 35, 36, 37 Therefore, the development of robust and reproducible numerical CT ventilation methods is essential for broadening the scope of potential applications and decreasing the uncertainty associated with evaluating CT ventilation on a voxel level.

Considering that frameworks based on measuring individual voxel volume changes induced by a DIR solution are prone to numerical instability, our proposed methodology is based on computing a series of more reliable subregional volume change estimates. Thus, the IJF framework is premised on the ability to compute the warped subregion volume obtained by applying the DIR spatial transformation to a reference subregion. Estimating volumes of general multidimensional shapes is known to be a difficult problem.17 However, lung volume masks are not difficult to obtain and are often assumed to exist prior to DIR and CT ventilation computations in the radiotherapy setting. Moreover, the full lung volume is inherently discretized by the CT digital voxel grid. These characteristics allow for the warped subregion volume estimation problem to be cast as a simple sampling‐based integral approximation. In addition to providing a mathematical framework for characterizing uncertainty, the sampling‐based integral approach is defined in part by a membership oracle [Eq. (9)] that respects the spatial accuracy of the DIR solution (i.e., the oracle operates on the DIR solution rounded to the nearest integer). As a result, the oracle evaluations are equivalent for two DIR solutions with the same spatial accuracy. In this instance, the hit‐or‐miss sampling method and standard Gaussian statistics are used to characterize the uncertainty in the warped subregion volume estimates. However, it is certainly possible to incorporate more sophisticated statistical modeling and Monte Carlo machinery into the general IJF framework, or investigate the effect of adopting a Bernoulli distribution model for the hit‐or‐miss uncertainties. Although, considering that (a) N is assumed to be large and that (b) the probability of the subregional oracle function returning a “hit” (value of 1) will not be close to either zero or one, the use of Gaussian statistics is sound.38 Moreover, it may be possible to devise an alternative to the sampling approach based on applying level set or active contour methods to directly measure the warped subregional volumes.39

The numerical implementation of the IJF method involves making a series of subregional volume change estimations that satisfy an acceptable level of uncertainty in the relative approximation error [Eq. (26)]. As opposed to the FDJ method, individual voxel volume values are coupled together via the subregional estimates. Thus, the individual values are described as the solution to a constrained least squares problem [Eq. (25)]. The resulting IFJ ventilation images are robust to small changes in DIR solution, as evidenced by the numerical results described in Table 2, whereas the FDJ method's results are consistent with the numerical analysis of finite difference approximations (see Fig. 1). In particular, two DIR solutions obtained for 4DCT test case 9 had an equivalent average millimeter landmark error of 1.02, while the Pearson correlation between the corresponding FDJ images was 0.83. Variation with respect to the DIR solution occurred consistently within the FDJ images across all ten test cases. Conversely, the IJF method maintained a Pearson correlation of 1.00 between the images generated from the two DIR solutions on all ten cases.

The IJF method represents a robust numerical method for computing CT ventilation from DIR‐measured volume changes. However, the method is still susceptible to the uncertainties induced by patient breathing variations and the quality of the 4DCT acquisition. The effect of these uncertainties becomes apparent when comparing the IJF ventilation images generated from treatment planning 4DCT and 4DCBCT acquired for the same patient prior to radiotherapy. Across five test cases, the Pearson correlations detailed in Table 3 ranged between 0.61 and 0.95, with two cases scoring at or below 0.69, and three scoring at or above 0.80. As illustrated in Fig. 6, a large phase bin artifact is present near the diaphragm of the 4DCTs for the two low scoring cases. Recovered Jacobian values are near 1.0 in the regions around the artifact and close to the diaphragm, indicating little to no volume change in those areas. The corresponding CBCT IJF image does not possess the acquisition artifact and consequently does not depict a ventilation defect near the diaphragm. Considering that the high scoring cases did not contain apparent binning artifacts, the results suggest that differences between the IJF 4DCT and 4DCBCT ventilation images for the high scoring cases were primarily induced by variations in the patients’ breathing across the two imaging sessions, whereas for the lower scoring cases, 4DCT acquisition artifacts are a significant confounding factor.

In recent work,30 a mean Spearman correlation of 0.60 (0.23) between 4DCBCT and 4DCT ventilation was reported, whereas for intrafraction 4DCBCT‐generated ventilation images, a similar mean of 0.67 (0.20) was reported in Ref. 29. Moreover, the mean absolute differences demonstrated by the FD method (Table 3) ranged between 0.10 and 0.19, which is consistent with previously reported mean absolute differences measured between intrafraction 4DCBCT ventilation images.29 The 0.77 (0.13) average Pearson correlation measured between IJF 4DCBCT and 4DCT ventilation, as well as corresponding mean absolute differences ≤0.06, represent progress in comparing ventilation images generated from different modalities and at different time points. This improvement reflects both the inherent numerical instability properties of FDJ ventilation and the mechanisms for controlling DIR uncertainty built into the IJF method. We believe the robust numerical stability of the IJF method can aid in providing improved results to previously conducted validation studies and in future CT ventilation applications. As such, an immediate area of our future research is to investigate the correlation of IJF with previously acquired validation data (such as those in Ref. 3, 5, 40), as well as the potential dosimetric consequences of IJF on radiotherapy functional avoidance planning studies (such as those in Ref. 41, 42).

Since IJF is based on approximating subregional volume changes, another area of future work will focus on alternatives to the hit‐or‐miss sampling method for volume change estimation, including level set method and more sophisticated Monte Carlo approaches. Moreover, the current numerical implementation treats the Jacobian value for each voxel as an unknown in a least squares optimization problem. Finite element methods are well suited for integrated formulations and could substantially reduce the computational workload of the IJF method. Finally, the IJF optimization problem requires a regularization model for the recovered ventilation image. Models that allow discontinuities within the recovered solution, such as total variation, might be more appropriate considering that functional defects are likely to be contained within lung lobes.

Conclusions

The IJF method for estimating volume changes from a DIR spatial transformation is presented. The IJF method recovers individual voxel volume changes from a series of subregional volume change measurements. Uncertainty in the subregional estimates is characterized using Gaussian statistics and standard error analysis. Numerical experiments indicate that the IJF method is robust with respect to variations in DIR solution and demonstrates good reproducibility between 4DCT and 4CBCT acquisitions. An immediate area of future research is to apply a similar verification process to HU‐based ventilation methods.

Conflict of interest

The authors have no conflict to disclose.