Robust HU-based CT ventilation from an integrated mass conservation formulation.

Computed tomography (CT) ventilation algorithms estimate volume changes induced by respiratory motion. Existing Hounsfield Unit (HU) methods approximate volume change from the measured HU variations between spatially corresponding voxel locations within a temporally resolved CT image pair, assuming that volume changes are caused solely by changes in air content. Numerical implementations require a deformable image registration to determine the inhale/exhale spatial correspondence, a preprocessing lung volume segmentation, a preprocessing high-intensity vessel segmentation, and a post-processing smoothing applied to the raw volume change estimates obtained for each lung tissue voxel.
PURPOSE: We introduce the novel mass-conserving volume change (MCVC) method for estimating voxel volume changes from the HU values within an inhale/exhale CT image pair. MCVC is based on subregional volume change estimates that possess quantitatively characterized and controllable levels of uncertainty. MCVC is therefore robust to small variations in DIR solutions and the resulting ventilation images are overall more reproducible. In contrast to existing HU methods, MCVC does not require a preprocessing lung vessel segmentation or pre/post-processing Gaussian smoothing.
METHODS: Subregional volume change estimates are defined in terms of mean density ratios. As such, the corresponding uncertainty is characterized using Gaussian statistics and standard error analysis of the sample density means. A numerical solution is obtained from the MCVC formulation by solving a constrained linear least squares problem defined by a series of subregional volume change estimates. Reproducibility of the MCVC method with respect to DIR solution was assessed using expert-determined landmark point pairs and inhale/exhale phases from 10 four-dimensional CTs (4DCTs) available on www.dir-lab.com. MCVC was also compared to the robust Integrated Jacobian Formulation (IJF), a transformation-based ventilation method.
RESULTS: The ten Dir-Lab 4DCT cases were registered twice with the same DIR algorithm, but using different degrees of freedom (DIR 1 and DIR 2). Standard HU ventilation (HUV) and MCVC ventilation images were computed for both solutions. The average spatial errors (300 landmarks per case) for DIR 1 ranged between 0.74 and 1.50 mm, whereas for DIR 2 they ranged between 0.68 and 1.18 mm. Despite the differences in spatial errors between the two DIR solutions, the average Pearson correlation between the corresponding HUV images was 0.94 (0.03), while for the MCVC images it was 1.00 (0.00). The average correlation between MCVC and IJF ventilation over the ten test cases was 0.81 (0.14), whereas for HUV and IJF it was 0.56 (1.11).
CONCLUSION: While HUV is robust to DIR solution, its implementation depends on heuristic Gaussian smoothing and vessel segmentation. MCVC is based on subregional volume change measurements with quantifiable and controllable levels of uncertainty. The subregional approach eliminates the need for Gaussian smoothing and lung vasculature segmentation. Numerical experiments are consistent with the underlying mathematical theory and indicate that MCVC ventilation is more reproducible with respect to DIR algorithm than standard HU methods. MCVC results are also more consistent with the robust IJF method, which suggests that incorporating robustness leads to more consistent results across both DIRs and ventilation algorithms.

Introduction

Computed tomography (CT) derived ventilation is an image processing based modality where information provided by segmentation and deformable image registration (DIR) algorithms is leveraged to quantify the apparent volume change within a CT inhale/exhale image pair.1, 2 There are two primary classes of ventilation algorithms. Transformation‐based ventilation algorithms recover volume changes directly from the Jacobian factor of the DIR estimated spatial transformation.3 Intensity (or Hounsfield Unit)‐based ventilation methods are dependent on the CT values and estimate volume changes by analyzing the CT value variations between spatially corresponding voxels within the inhale/exhale CT pair.1 While validation of CT ventilation has been a major focus of research,4, 5, 6, 7 surprisingly little work has focused on numerical analysis of CT‐ventilation algorithms. Recently, we characterized the numerical properties of transformation‐based methods8 and presented a robust algorithm for estimating the Jacobian factor of a DIR spatial transformation.9 However, a similar analysis for intensity‐based ventilation methods has yet to be presented in the literature. As demonstrated in Ref. [8], a ventilation method can give substantially different results on the same image data when supplied with different DIR solutions. Addressing the current lack of reproducibility in CT ventilation has the potential to increase the significance of clinical validation studies and move the modality closer to widespread clinical use. Thus, as opposed to clinical validation, the focus of this work is on the development of a HU ventilation method that possesses robust numerical properties and is reproducible with respect to DIR solution.

Image formation in CT derives from the measurable attenuation of x rays traversing an imaged body within the confines of a well‐defined geometry and incident photon fluence. Reconstructed images reflect the spatial distribution of x‐ray attenuation, primarily due to photoelectric and Compton interactions for the clinically relevant peak kilovoltage potentials used to generate the incident x rays. Diagnostic radiology and radiotherapy simulation scanners are calibrated such that inherently quantitative image values assigned to discretized volume elements (voxels) convey explicit physical significance according towhere HU is the Hounsfield Unit,10 is the linear attenuation coefficient of water, and represents the average linear attenuation coefficient for the material contained within the voxel centered at . CT is primarily an anatomic imaging modality, although the use of hyperdense radiographic contrast agents and/or multispectral incident energy fluence have facilitated direct imaging of physiologic metrics such as pulmonary ventilation from Xenon‐CT11 or pulmonary perfusion by dual energy CT with iodinated contrast.12

Since 2000, there has been considerable interest in the use of image processing techniques leveraging physical characteristics of the image acquisition process to compute regional volume change as a surrogate metric for ventilation. The approach was first articulated in a report by Simon,13 which outlined a formalism for expressing specific lung compliance in terms of regionally correlated HU in images acquired at two different distending pressures. This approach is derived from the basic assumption that lung tissue can be physically approximated as a linear combination of air and tissue components. This allows for the volume of a region of interest, to be expressed in terms of the respective proportions of each component:where and represent the relative proportion of air and tissue contained in . Owing to the linearity of the HU scale relative to , the fractional proportion can be inferred under the two‐component model according to the linear form

Assuming an idealized scenario in which and , the relative proportion of air to the total volume of a given region can be approximated as:where is averaged within . Thus, a regional estimate for specific compliance can easily be computed by using Eq. (3) to approximate the change in air content between spatially corresponding regions within a pair of images, if the distending pressure is known. In addition to the Eq. (2) material approximation, this approach implicitly assumes that observed image intensity changes are caused solely by changes in air content.

Considering that pressure differences between image acquisitions are not readily or routinely controlled in clinical subjects, Guerrero et al. modified Simon's formulation to describe specific volume change in terms of spatially registered voxel‐level HU between inhale/exhale BH‐CT image pairs.1 Notably, a modification was required to account for an observable variation in lung tissue mass between respective inhale and exhale segmented lung volumes. The material density estimate used for this correction, which is derived from the Eq. (2) model, is defined (in units:) as:

In a latter study,14 Eq. (4) was used to demonstrate cyclic variation in apparent lung mass on four‐dimensional CT (4DCT), with magnitude variation on the order of 40–50  between peak phases in three lung cancer patients. The observed fluctuation was attributed to global variation in lung perfusion over the respiratory cycle. However, with respect to estimating volume changes, the assumption that variations in HU are attributed only to change in local air content implies that tissue mass must be conserved.14 As such, a global correction factor, which assumes a uniform perfusion effect, was applied to enforce mass conservation between segmented lung volumes.

Further examination of the Guerrero's original voxel model reveals an additional consideration that becomes applicable as the size of registered volumes approaches that of the underlying image discretization. The expansion in lung volume between exhale and inhale breathing states implies that ventilated tissue elements contained within a single voxel on exhale will necessarily expand to be contained within multiple voxels on inhale. Without correction, the model would result in overestimation of the spatial distribution of ventilation by treating each redundant pair of registered voxels as independently expanding tissue elements. Castillo et al.2 addressed this issue by reformulating the voxel model strictly in terms of exhale grid points and showed that multiple registered voxels are appropriately treated according to the average expanded intensity. From this perspective, the accuracy of the volume change estimation is dictated by the accuracy of the estimated intensity mean, taken over all the intensity values mapped into each exhale voxel.

Numerical implementations of HU‐ventilation require a deformable image registration algorithm to generate a map of spatially corresponding voxel locations within a CT inhale/exhale image pair. Since volume change, as estimated by the HU voxel model, is determined for each exhale voxel solely from the intensities of individually mapped inhale voxels, the method is susceptible to small errors in the DIR solution. As a result, high intensity blood vessels are often segmented out of the HU ventilation calculation to prevent erroneous DIR mappings from corrupting the volume change measurement.2 Moreover, heuristically applied pre‐ and post‐processing Gaussian smoothing is often required to counteract the effects of CT noise and further reduce the effect of DIR errors.2 There is currently no accepted or general formalism for applying the smoothing or for interpreting its effect on the validity of the resulting ventilation image. These issues coupled with the heuristic components have hindered the reproducibility of HU‐ventilation15, 16 and, as a consequence, limited its clinical applicability.

In contrast to HU methods, transformation‐based CT ventilation methods estimate volume changes solely from the DIR by numerically approximating the Jacobian factor of the recovered spatial transformation on individual voxel volumes. Empirical analyses of these Jacobian methods, which depend on finite difference approximations applied to the DIR recovered displacement field, reported poor overall reproducibility16 with respect to DIR method,17 and 4DCT acquisition artifacts.18 These findings are consistent with numerical analysis theory describing the approach’s instability with respect to DIR solution.8 In recent work, we derived a framework for incorporating robustness into the Jacobian calculation.9 As opposed to traditional single‐voxel approaches, the resulting Integrated Jacobian Formulation (IJF) method recovers the ventilation image from a series of robust subregional volume change estimates. The subregional measurements are formulated as the sampled average of a membership oracle function defined over the lung region of interest. As such, their uncertainty can be characterized, and subsequently controlled, by examining the standard error of the sample mean. IJF demonstrated substantially improved reproducibility with respect to DIR solution and improved correlations between ventilation images generated from CT data acquired during different imaging sessions.9

The purpose of this work is to introduce a robust numerical method for computing HU‐based ventilation using the density formulation in Eq. (4) and the integrated form of the mass conservation assumption. As with Simon's original model, the proposed mass‐conserving volume change (MCVC) formulation ignores the effects of perfusion and assumes that mass is conserved between the inhale and exhale lungs. Unlike current approaches based on single‐voxel measurements, MCVC estimates the Jacobian factor of the DIR transformation using the HU‐defined density values within a series of spatially corresponding lung subregions, similar to the robustness strategy employed within the IJF method. Moreover, similar to Castillo's approach,2 subregional volume change is estimated from: (a) the mean density value within the subregional volume and (b) the mean density value of the voxels mapped into the subregional volume by the DIR solution. The numerical uncertainty can therefore be quantified as the standard error of the sample density means and the subregional volumes can be constructed such that the data estimates possess a specified level of acceptable uncertainty. Thus, unlike any existing HU methods, our proposed methodology characterizes and controls the numerical uncertainty in the data measurements used to generate the ventilation image. A full ventilation image, which provides volume change values for each lung voxel, is recovered from the subregional estimates by solving an inequality constrained‐linear least squares problem that is similar in structure to those arising in image‐deblurring problems.19 Since MCVC is based on robust subregional estimates, the method does not require the heuristic Gaussian smoothing or lung vasculature segmentation needed by traditional intensity methods.

Materials and Methods

Subregional volume change estimates

Applying Eqs. (1) and (4) to an inhale/exhale CT image pair results in two corresponding density functions, which we denote here as the reference image and the target image . For our purposes, the distinction is irrelevant (i.e., inhale and exhale can correspond to reference and target or vice versa). A spatial transformation computed with a DIR algorithm maps the reference image lung volume onto the target image lung volume . As such, the target lung volume is the image of under :

Mathematically, for a general subregion , the volume scaling factor under is described by the Jacobian factor:where is the first full derivative of and is assumed to be diffeomorphic. In contrast to HU methods, transformation‐based CT‐ventilation methods compute volume change images by numerically approximating Eq. (6) directly from the DIR solution . Our goal here is to derive a robust numerical approximation for the Jacobian factor using CT‐density values.

Assuming the volume change between and is caused solely by a change in air content (and ignoring the effects of perfusion) implies that mass is conserved:

Standard Monte Carlo integration methods, under mild assumptions, express an integral in terms of the integrand average over a desired volume.20 Applying the same approach to Eq. (7) results inwhere

Thus, the unknown Jacobian factor can be expressed in terms of image density means by substituting Eq. (6) into Eq. (8):

The true density means are not known and must be approximated by sample means:

andwhere each represents a discrete voxel location.

In general, is not explicitly provided by DIR. However, the DIR solution is assumed to be diffeomorphic,21 which implies the existence of . Thus, voxels within can be identified with the help of a membership oracle function22:so thatand

Finally, the MCVC estimate is defined in terms of the unknown Jacobian factor and the sample mean estimates:

Uncertainty in the subregional density means

The uncertainties associated with the density sample means, assuming Gaussian statistics and large N, can be characterized by the standard error, which for Eq. (12) is defined as:

Equation (18) can be used to provide a probabilistic assessment of the density mean approximation in terms of a confidence interval defined by a parameter . For instance, setting implies that with 95% probability:where represents the true mean.

Considering that the Eq. (4) density formula is linearly scaled, with zero corresponding to the density of air and corresponding to the density of water, it is reasonable to assume that for lung subregions. This impliesand the standard error can be bounded with respect to the number of samples:

Substituting the bound into Eq. (19) gives an uncertainty estimate that depends solely on the number of samples:

A similar result can be obtained for the Eq. (13) approximation using the same analysis. Therefore, given a specified tolerance and , ifwherethen with greater than 95% probability:

This follows from Eq. (21) and the fact that when Eq. (23) is satisfied.

General sampling methods reduce measurement uncertainties by raising the number of point samples within the volume of interest. This approach assumes that the integrand function is precise for any sampling resolution. However, in this case, the integrand is the density image, which is inherently discretized by the digital voxel grid. Therefore, increasing the number of samples can only be accomplished by growing the size of the subregion . Figure 1 shows a plot of as a function of and illustrates the reduction in uncertainty as grows.

Figure 1

For , the uncertainty tolerance is given as a function of the corresponding voxel sample size, as defined by Eq. (24). [Color figure can be viewed at http://wileyonlinelibrary.com]

The reduction in uncertainty of the volume change measurement comes at the expense of measurement resolution. In other words, because volume change is assessed with respect to the subregion size, volume change variations occurring within the subregion are coalesced into a single measurement. Thus, individual voxel volume changes must be computed from the subregional volume change measurements using a numerical algorithm.

Mass‐conserving volume change formulation

A full CT ventilation image, , requires computing the discretized variables:where represents the volume change of the undeformed unit volume voxel centered on . The general subregional volume change estimate given by Eq. (17) can be used to provide one equation for the total unknowns. Thus, intuitively, a recovery algorithm can be built on the idea of computing enough subregional estimates to define a full‐rank or overdetermined linear system of equations for the unknowns. Applying Eq. (17) to a series of subregions results in the following linear system:

whereand the elements of contain the subregional density means:

Even though the Eq. (27) subregions can be constructed to satisfy the uncertainty constraints defined by Eq. (23), the resulting data values are still approximations. As such, Eq. (27) should not be treated as a hard constraint. Moreover, Eq. (27) is not guaranteed to provide enough information to uniquely determine . Therefore, a two‐norm penalization on the spatial gradient of is employed to both induce smoothness in the reconstructed image and guarantee the existence of a 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 and the lower bound represents the minimum possible voxel volume size. We refer to Eqs. (27)–(30) as the MCVC estimation method.

Numerical implementation

The MCVC method requires constructing the subregion volume change estimates described in Eq. (27). To do this, we applied a dart‐throwing algorithm to the to generate a point cloud of voxel locations uniformly spaced approximately 5 mm apart.23, 24 The initial reference subregions are defined as the 7 × 7 × 3 voxel neighborhood (approximately 7 × 7 × 7 mm for typical voxel sizes of thoracic CTs) centered on . Each initial subdomain is morphologically dilated with a 7 × 7 × 3 structuring element and the corresponding is recomputed until Eq. (23) is satisfied for and a user‐defined . This dilation process has the effect of adding more samples to the mean estimates, thereby reducing uncertainty in the measurements. When not otherwise specified, and the Eq. (30) regularization parameter . Numerical experimentation and visual inspection indicated that solutions did not vary greatly for between 0.05 and 2.0. As with all regularization strategies, the problem could potentially become ill‐conditioned for smaller regularization parameter values, thereby resulting in possible numerical errors. For larger values, the resulting solutions have the potential to become dominated by the regularization term, which would result in over smoothed images.

For the peak inhale and exhale 4DCT phases, lung masks, were generated using a semi‐automated histogram segmentation (as done in Ref. [2]). The MCVC derivation and formulation only require the lung ROI mask. Therefore, no lung vessel delineation is needed.

The spatial transformation is computed by applying the QPDIR algorithm (see Ref. [25] for details) to the target and reference image. Equation (30) is solved using the augmented Lagrangian method (see for Ref. [26] for details), implemented in MATLAB (release R2017b, The Mathworks Inc, Natick, MA, USA).

Image data

Ten publicly available radiotherapy treatment planning 4DCTs from the http://www.dir-lab.com website (see Ref. [27, 28, 29] for acquisition details) were used to assess the reproducibility of the MCVC method and the original HU ventilation method2 (HUV) with respect to DIR solution. Correlation with IJF was also assessed using the same test cases. In all cases, the reconstructed CT slice thickness is 2.5 mm, whereas the in‐plane dimension ranged from (0.97 × 0.97)–(1.16 × 1.16) mm2.

Image registration

Quadratic penalty DIR (QPDIR) is an intensity‐based DIR algorithm designed around a gradient‐free block coordinate descent strategy that essentially iterates between simple block matching operations and linear least squares solves to minimize the structural similarity index between an image pair.25 The algorithm was applied to the lung ROIs twice, first using a coarse displacement field parameterization with fewer degrees of freedom (2000), and then again using a fine displacement field parameterization with more degrees of freedom (20 000). Intuitively, for nonlinear lung deformation, a coarse displacement field parameterization is unable to capture the highly nonlinear nature of lung deformation, which in turn leads to a reduction in DIR spatial accuracy. The measured accuracy differences between the two DIR solutions are representative of the variance in DIR results detailed on http://www.dir-lab.com across over roughly 30 algorithms applied to the same test cases.

HUV and IFJ implementation

HUV was computed by first segmenting the lung vasculature with a semi‐automated seed growing method. A 9 × 9 × 3 voxel preprocessing smoothing filter was applied to the segmented CT data, and a 9 × 9 × 3 voxel post‐processing smoothing filter was applied to HU‐measured volume change estimates. This approach represents a common implementation strategy for traditional HUV methods.2, 15, 30 IJF was computed using the exact description in Ref. [9]

Results

In order to assess the reproducibility of the HUV and MCVC ventilation with respect to DIR solution, the ten test cases were each registered twice with the QPDIR method. The first DIR solution was computed with a smaller number of degrees of freedom (approximately 2000 for each case), while the second was computed with a larger number (approximately 20 000 for each case). The resulting mean spatial errors for the two solutions, as measured by 300 expert‐determined landmarks per case, are provided in Table 1 and are indicative of the differences between the two solutions.

Table 1

Data set four‐dimensional computed tomography (4DCT): Ten 4DCT maximum inhale/exhale phase pairs.

	DIR 1 avg. mm error	DIR 2 avg. mm error	Mean absolute difference			Pearson correlation
			HUV		MCVC	HUV	MCVC
4DCT 1	0.76 (0.96)	0.70 (0.91)	0.01 (0.02)	0.00 (0.00)		0.97	1.00
4DCT 2	0.74 (0.89)	0.68 (0.90)	0.01 (0.02)	0.00 (0.00)		0.94	0.99
4DCT 3	1.03 (1.09)	0.87 (1.05)	0.02 (0.02)	0.00 (0.00)		0.95	1.00
4DCT 4	1.50 (1.34)	1.18 (1.22)	0.04 (0.05)	0.00 (0.00)		0.95	1.00
4DCT 5	1.40 (1.54)	1.07 (1.40)	0.03 (0.04)	0.00 (0.00)		0.90	0.99
4DCT 6	1.21 (1.08)	0.90 (0.97)	0.06 (0.08)	0.00 (0.00)		0.94	1.00
4DCT 7	1.28 (1.11)	0.81 (0.91)	0.04 (0.05)	0.00 (0.00)		0.93	1.00
4DCT 8	1.26 (1.26)	1.03 (1.22)	0.04 (0.05)	0.00 (0.00)		0.96	1.00
4DCT 9	1.16 (0.98)	0.89 (0.95)	0.04 (0.05)	0.00 (0.00)		0.95	1.00
4DCT 10	1.15 (1.09)	0.83 (0.94)	0.05 (0.06)	0.00 (0.00)		0.88	1.00
Avg. (std):						0.94 (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 QPDIR algorithm. The spatial accuracies achieved by the two deformable image registration (DIR) approaches are given in average (std) millimeter error with respect to 300 landmark point pairs. The Pearson correlation and mean absolute difference between the Hounsfield unit ventilation (HUV) ventilation images computed from DIR 1 and DIR 2 are given, as are they are for the mass‐conserving volume change (MCVC) ventilation images computed from DIR 1 and 2.

For each case, the two DIR solutions were used to compute two HUV images. The Pearson correlation between the two HUV images was computed and the same process was applied to the MCVC method. All HUV correlations were computed on the lung volume excluding the vessel segmentation, which defines lung parenchyma regions. MCVC correlations were computed on the full lung ROI. Results, which are detailed in Table 1, indicate that both methods demonstrated good reproducibility across the two DIR solutions, with the HUV method generating an average correlation of 0.94 (0.03) and the MCVC method generating an average correlation of 1.00 (0.00). The HUV and MCVC images for 4DCT 6 are illustrated in Fig. 2.

Figure 2

Coronal slices of the mass‐conserving volume change (MCVC) and Hounsfield unit ventilation (HUV) for a patient with malignant airway stenosis (Table 1, four‐dimensional computed tomography 6). HUV appears noisier than the robust MCVC image, but the images are qualitatively similar.

The uncertainty tolerance for the MVCV method was , which resulted in a minimum subvolume size of [Eq. (24)]. However, the subregion construction procedure (see Section 2.4) only guarantees that the subregion will contain . Table 2 details the number of subregions generated and the range in sizes. In order to demonstrate the effect of on the MCVC subregional measurements in Eq (29), different values for were used to generate estimates for the representative case 4DCT 6. Box plots of the resulting subregional estimates are provided in Fig. 3 and corresponding coronal slices from the resulting MCVC images are illustrated in Fig. 4.

Table 2

Mass‐conserving volume change (MCVC) subregion summaries.

	K	minΩ_k	maxΩ_k	minϕΩ_k	maxϕΩ_k
4DCT 1	5323	6150	10571	6159	12954
4DCT 2	8422	6148	10906	6167	12700
4DCT 3	6403	6149	10906	6155	13342
4DCT 4	4500	6147	10894	6154	16068
4DCT 5	5856	6147	10571	6151	14303
4DCT 6	5705	6148	11058	6156	16787
4DCT 7	7208	6147	11160	6148	15746
4DCT 8	11386	6147	11075	6147	16969
4DCT 9	3956	6150	11077	6154	14870
4DCT 10	6875	6147	11706	6148	14919

A summary of the subregional domains used within the MCVC method for DIR 1 (Table 2), including: the number of total subregions, K, and the minimum and maximum sizes of the reference subregions, , and the corresponding mapped subregions,. For these experiments, N* = 6147, (as defined in Eq. (24)), which corresponds to

Figure 3

Box plots of the subregional volume change estimates [Eq. (29)] obtained for different levels of (uncertainty) for test cast four‐dimensional computed tomography 6. The plots were created using the MATLAB (release 2017b, The Mathworks Inc., Natick, Massachusetts, United States) boxplot function. On each box, the central mark indicates the median, and the bottom and top edges of the box indicate the 25th and 75th percentiles, respectively. The whiskers extend to the most extreme data points not considered outliers. While smaller values of reduce uncertainty in the measurements, the spatial resolution of measurements is also reduced. This is illustrated by the reduction in both outliers and data variation as decreases. [Color figure can be viewed at http://wileyonlinelibrary.com]

Figure 4

Corresponding coronal slices from mass conserving volume change images generated from the same computed tomography images and DIR but with two different uncertainty tolerances. The image created using the stricter tolerance () is less spatially varying than the one generated from the more relaxed .

IJF ventilation was computed using the DIR 2 solutions for all cases. The Pearson correlations between HUV and IJF, HU and MCVC, and MCVC and IJF (all using the DIR 2 solution), are given in Table 3. Correlations between IJF and MCVC were computed over the entire lung ROI, whereas correlations with HUV were computed on the lung volume, excluding the vessel segmentation, on which the HUV is defined. For all but one case, the correlation between MCVC & IJF was higher than those computed with HUV.

Table 3

Correlation between Hounsfield unit ventilation (HUV), Integrated Jacobian formulation (IJF), and mass‐conserving volume change (MCVC).

	HUV MCVC	HUV IJF	MCVC IJF
4DCT 1	0.75	0.62	0.88
4DCT 2	0.62	0.44	0.76
4DCT 3	0.65	0.53	0.88
4DCT 4	0.74	0.67	0.87
4DCT 5	0.53	0.37	0.65
4DCT 6	0.72	0.58	0.76
4DCT 7	0.72	0.64	0.92
4DCT 8	0.69	0.59	0.92
4DCT 9	0.77	0.72	0.96
4DCT 10	0.56	0.44	0.53
Avg. (STD)	0.67 (0.08)	0.56 (0.11)	0.81 (0.14)

The correlation between the HUV, MCVC, and IJF ventilation images generated from deformable image registration 2 for all ten test cases is given.

In order to the assess the sensitivity of the reported MCVC reproducibility results to the selected parameter values, for all ten test cases, 16 MCVC images were created using all value combinations between and for both DIR 1 and DIR 2. The average Pearson correlation between the DIR 1 and 2 MCVC images across all ten test cases for each parameter set combination are given in Table 4. These results indicate that the MVCV images maintain robustness for any combination of and .

Table 4

Mass‐conserving volume change (MCVC) reproducibility (Pearson correlation) with respect to and .

α	τ
	0.400	0.100	0.025	0.010
0.05	0.969 (0.014)	0.989 (0.005)	0.998 (0.001)	0.999 (0.001)
0.25	0.986 (0.009)	0.996 (0.002)	0.999 (0.001)	1.000 (0.000)
0.50	0.991 (0.007)	0.997 (0.002)	0.999 (0.001)	1.000 (0.000)
2.00	0.996 (0.004)	0.999 (0.001)	1.000 (0.000)	1.000 (0.000)

Each entry represents the average (std.) Pearson correlation between the MCVC images computed from deformable image registration (DIR) 1 and DIR 2, given the specified values for and , across all ten test cases. The results indicate that the MCVC method is robust to DIR solution for varying parameter value combinations.

Discussion

The proposed MCVC method is based on robust subregional measurements, as opposed to current approaches, which are based on single voxel measurements. Specifically, MCVC uses the basic mass conservation model to express the Jacobian factor, of the DIR transformation in terms of the HU‐defined density values within a series of spatially corresponding lung subregions. As such, the MCVC method estimates the same quantity as transformation‐based Jacobian ventilation methods. However, the MCVC volume change numerical approximations are computed from subregional mean density estimates. Assuming Gaussian statistics, the approximation uncertainty can be characterized by the standard error of the mean. As a consequence, robustness can be incorporated into the volume change approximation by constructing subregions that are guaranteed to achieve a specified level of uncertainty.

As evidenced by Table 1, the Gaussian pre‐ and post‐processing smoothing utilized by HUV methods necessarily introduces robustness with respect to DIR solution. However, there is no consensus or theoretical justification for choosing the size of the smoothing kernel. While the subregional density mean estimates used within the MCVC are conceptually similar to the HUV Gaussian smoothing procedure, the effect of volume size is well characterized by the standard error [Eq. (18)]. In fact, further analysis on the standard error yields the precise volume size, needed to guarantee a specified level of uncertainty [Eq. (24)]. Thus, represents the resolution of the MCVC subregional volume change approximation, given a specified uncertainty tolerance.

The uncertainty tolerance represents the tradeoff between the volume change measurement's resolution and its fidelity. As illustrated in Fig. 3, smaller reduces uncertainty, but the resulting measurements must be acquired over a larger subvolume size, thereby reducing the millimeter resolution of the volume change estimation and inducing smoothness in the resulting MCVC image. As shown in Fig. 4, the MCVC images generated using depicts more spatial variation than the ones computed with . However, the derivation of is independent of the CT image resolution. In other words, is constant with respect to the total number of CT image voxels. Therefore, the millimeter resolution of the MCVC volume change measurements is proportional but less than the resolution of the utilized CT images. Put another way, for a fixed , a higher spatial resolution in the CT images implies higher resolution MCVC volume change measurements.

Automated vessel segmentation is a notoriously difficult problem.31 Consequently, HUV typically employs a semi‐automated process to generate the required vessel mask. For MCVC, precise density mappings are not crucial since subregional measurements are inherently robust to small errors in the DIR solution, as demonstrated by the 0.00 (0.00) mean absolute differences between the MCVC images computed from two DIR solutions (Table 1). Furthermore, while the impact of subregional size is well characterized, the MCVC formulation [Eq. (30)] requires a regularization model and corresponding penalty parameter . Variability with respect to the uncertainty and regularization parameters was assessed using a parameter sweep. The average Pearson correlations computed across all ten test cases for 16 different combinations of and were computed and detailed in Table 4. The results indicate that the MVCV images maintain robustness for any combination of and , with the lowest measured average correlation being 0.97 ().

Our recently developed IJF method for transformation‐based ventilation is conceptually similar to the MCVC in that both methods incorporate robustness through subregional volume change estimation. From an engineering/applied mathematics perspective, incorporating robustness into the numerical method is in line with the definition of verification, which is the process of developing a numerical method that accurately solves a mathematical model. In this case, the mathematical model is the Jacobian of the DIR transformation. Validation, on the other hand, is the process of assessing whether the mathematical model accurately describes an observable physical phenomenon, which in this case, is ventilation. In our view, studies of verification and validation need to be (a) decoupled and (b) performed in the correct order (i.e., verification prior to validation). The work presented here is aimed at the verification component. As such, MCVC clinical utility cannot be directly inferred from this verification work, and must be addressed separately with a rigorous validation study. Moreover, since DIR robustness is maintained for a wide range of parameter values (Table 4), optimal values should be selected with respect to clinical data. This could potentially be accomplished by adopting a data science framework where the optimal parameter set is defined as the one which yields the best correlation with an accepted clinical modality, such as hyperpolarized He3 MRI. However, validation of CT ventilation is itself a challenging problem. Two examples of uncertainties within available current validation data sets and approaches are (a) lack of reproducibility and (b) lack of reliable ground‐truth. With respect to (a) available data sets32 do not provide a measure of reproducibility, therefore, it would be challenging to tell if the validation results are within the error‐bounds of the data set. With regard to (b), other lung function imaging modalities (SPECT, MRI, PET) suffer their own shortcomings (clumping artifact, poor spatial resolution, etc.) which make the designation of a true gold standard challenging.

As reported in Table 3, MCVC is more consistent with the robust IJF method than the HUV method. This result suggests that incorporating robustness into the formulation leads to more consistent results across algorithm class, as well as with respect to DIR solution. In fact, the average correlation between MCVC and IJF was 0.81 (0.14). Measured differences between the methods presumably reflect the known flaws in both the IJF and MCVC models. As discussed in Ref. [9], discrepancies in the lung geometry, such as those induced by 4DCT acquisition artifacts, are detrimental to IJF performance. As illustrated in Fig. 5, Case 10 (Table 3), which produced the lowest MCVC and IJF correlation, contains a phase‐bin 4DCT artifact that has the visual effect of cropping an apparent tumor. Despite describing a primarily expansive exhale‐to‐inhale lung motion, the resulting IJF volume change erroneously computes an area of contraction near the artifact. This reflects the apparent shrinkage of the tumor volume caused by the artifact. Accounting for motion and reconstruction artifacts is currently not possible by any numerical methodology and is typically accounted for by manually reviewing acquired 4DCTs.33 Even though geometric artifacts will intuitively have less influence on MCVC, due to its strict reliance on CT intensities, the effects of perfusion are potentially significant14, 34 and could also contribute to the discrepancies between IJF and MCVC.

Figure 5

Top Row: Spatially corresponding coronal slices from the Case 10 (Table 1) inhale/exhale images. A phase binning artifact has the apparent effect of artificially truncating the upper portion of a tumor (yellow arrows). Bottom Row: The IJF Jacobian estimation (left) and the mass conserving volume change Jacobian estimation (right) are dissimilar near the artifact. Near the artifact, the Integrated Jacobian Formulation image indicates an erroneous contractive motion induced by the tumor truncation.

Like IJF, MCVC approximates the Jacobian of the DIR transformation, which itself describes respiratory induced motion. Several physiological factors contribute to lung volume change within the respiratory cycle, but like all Jacobian ventilation methods, MCVC describes the effect of all physiological factors combined. Recent work has demonstrated that there are two factors contributing to the apparent density variations between inhale and exhale CT images: (a) changes in air content and (b) changes in mass, which presumably describe variations in the spatial distribution of pulmonary blood mass induced by respiratory motion.34 While the effect of blood mass variation was shown to be measurable, its effect on density variation was also shown to be substantially less than that of air content changes.34 Determining the degree to which blood mass variations effect the estimated Jacobian values is not known. This question and the mathematical theory needed to account for mass change variations are the focus of our immediate future work.

Conclusion

We introduced the MCVC method for estimating Jacobian factor volume changes from HU defined densities. The MCVC method recovers individual voxel volume changes from a series of subregional volume change measurements, which are defined in terms of mean density ratios. The uncertainty in the subregional estimates is characterized using Gaussian statistics and standard error analysis applied to the sample density means. As a result of the subregional approach, MCVC does not require the heuristic approach to Gaussian smoothing or lung vasculature segmentation needed by traditional HU methods. Numerical experiments indicate that the MCVC method is robust with respect to variations in DIR solution and that MCVC demonstrates good correlation with the robust IJF transformation‐based ventilation method. This indicates that incorporating robustness into the ventilation method leads to more consistent results across algorithm class as well as with respect to DIR solution.

Conflict of Interest

The authors have no conflict to disclose.