1 Introduction

Lung cancer is the leading cause of cancer death (\(25.9\%\) of all cancer deaths) in the United States [11]. Approximately \(50\%\) of the patients receive radiation therapy (RT) during the course of lung cancer [4]. The goal of RT is to deliver high-energy radiation beams to the tumors to kill cancer cells. However, irradiation of the surrounding healthy lung tissue during RT causes lung toxicity to about \(5\% \) to \(20\%\) of lung cancer patients [10]. To reduce this radiation-induced lung injury, functional avoidance RT has been proposed to minimize irradiation of high-function lung regions [3, 8, 9, 13,14,15,16,17, 19].

In functional avoidance RT planning, the high-function lung regions can be identified by imaging of lung function. In this paper, we use ventilation as a synonym for lung function since the main function of lung is for gas exchange. Clinical standard ventilation modalities such as single photon emission computed tomography (SPECT) and positron emission tomography (PET) have been used in functional RT planning [3, 9, 13, 14]. Although studies have shown the possibility of reducing dose to high-function tissue in RT using SPECT and PET, these techniques are often limited by low spatial resolution, high cost, long scan time, and/or low accessibility to patients [18]. Recently, ventilation images derived from four-dimensional computed tomography (4DCT) data have been used in functional avoidance RT [5, 8, 15, 17, 19]. One advantage of CT ventilation imaging is that it only requires processing of 4DCT data and acquisition of a 4DCT scan is often included in the treatment planning for lung cancer. Therefore, CT ventilation imaging is less expensive and more accessible to patients than clinical standard ventilation imaging techniques. Moreover, 4DCT has a shorter scan time and can be used to generate CT ventilation images with a higher spatial resolution than clinical standard lung function modalities.

CT ventilation images may be estimated by local expansion ratio (LER) of the lung [6, 12]. We define the LER at each voxel as the ratio of the maximum to the minimum local lung volume in a breathing cycle. Most CT ventilation imaging algorithms use pairwise image registration to find a one-to-one correspondence map between the end exhale (0EX) CT image and the end inhale (100IN) CT image. Pulmonary ventilation at each voxel can then be estimated by the Jacobian determinant of the correspondence map at that voxel [12], which we refer to as LER3D.

However, the lung may have out-of-phase ventilation, i.e., local lung volume change is out of phase with respect to global lung expansion and contraction. In such a situation, the LER3D measure which only uses the 0EX and 100IN phases may underestimate the LER quantity. In this paper, we proposed a new measure of LER by all phases of 4DCT, which we refer to as LER4D. Note that both LER3D and LER4D provide voxel-wise measurement of lung function. The purpose of this paper is to quantify the amount of out-of-phase ventilation of the lung.

2 Methods

2.1 Image Acquisition

This study evaluated fourteen human subjects undergoing radiation therapy and was approved by the University of Wisconsin-Madison institutional review board. Two 4DCT scans were acquired for each subject before radiation treatment, with a 5-min break between two scans. The 4DCT scans were acquired on a Siemens EDGE CT scanner using 120 kV, 100 mAs per rotation, 0.5 s tube rotation period, 0.09 pitch, 76.8 mm beam collimation, 128 detector rows, and reconstructed slice thicknesses ranging between 0.6 and 3 mm. Musical melody and voice instruction guidance were played throughout the scan to improve the repeatability of the respiratory pattern. Each 4DCT data set was reconstructed into 10 breathing phases, with 20\(\%\) (20IN), 40\(\%\) (40IN), 60\(\%\) (60IN), 80\(\%\) (80IN) and 100\(\%\) (100IN) inspiration phases and 80\(\%\) (80EX), 60\(\%\) (60EX), 40\(\%\) (40EX), 20\(\%\) (20EX) and 0\(\%\) (0EX) expiration phases.

2.2 Tissue-Volume Preserving Deformable Image Registration

Given a fixed image and a moving image, the goal of image registration is to find a one-to-one correspondence map between the two images. Image registration is an optimization problem whose objective function is a combination of the difference between the fixed image and the deformed image, and the smoothness of the correspondence map. In this paper, we focus on pairwise registration of volumetric CT images.

We denote the fixed and moving CT lung images by \(I_0:\varOmega \rightarrow \mathbb {R}\) and \(I_1:\varOmega \rightarrow \mathbb {R}\), respectively, where the closed and bounded set \(\varOmega \subset \mathbb {R}^3\) is the image domain. The CT images in Housfield unit (HU) can be converted into tissue ratio images by

$$\begin{aligned} \frac{HU - HU_{air}}{HU_{tissue} - HU_{air}} = \frac{HU + 1000}{1000} \end{aligned}$$
(1)

where the HUs of the tissue and the air are approximately \(HU_{tissue} = 0\) and \(HU_{air} = -1000\).

We denote the tissue ratio images associated with \(I_0\) and \(I_1\) by \(R_0\) and \(R_1\), respectively, i.e., \(R_0=\frac{I_0 + 1000}{1000}\) and \(R_1=\frac{I_1 + 1000}{1000}\). The HU of CT lung image varies with tissue density change during breathing (see Eq. 1). To take into account this variation of CT intensity, we use the sum of squared tissue volume difference (SSTVD) similarity metric [1, 7, 20].

The input to the registration algorithm are the tissue ratio images \(R_0\) and \(R_1\). The moving image \(R_1\) is deformed by the transformation \(\phi : \varOmega \rightarrow \varOmega \) by the operation of \(\phi \) acting on \(R_1\) denoted by \(\phi \cdot R_1\). This action is defined as follows: \(\phi \cdot R_1 \triangleq |J_{\phi }|\times R_1\circ \phi \), where \(|J_{\phi }|\) is the Jacobian determinant of \(\phi \). It can be shown that this definition is a group action. The SSTVD similarity metric is then given by

$$\begin{aligned} C_{SSTVD} = \int _{\varOmega }\Big (R_0(x) - |J_{\phi }|(x)\times R_1(\phi (x))\Big )^2 dx \end{aligned}$$
(2)

The regularization constraint is given by

$$\begin{aligned} Reg(\phi ) = \int _{\varOmega }\Big ( c_1(\nabla \cdot \nabla )u(x) + c_2\nabla (\nabla \cdot u(x)) \Big ) dx \end{aligned}$$
(3)

where \(\nabla = [\frac{\partial }{\partial x_1}, \frac{\partial }{\partial x_2}, \frac{\partial }{\partial x_3}]^{T}\), \(\nabla \cdot \) is the divergence operator, \(u = \phi - \text {Id}\) is the associated displacement vector field, and \(c_1\) and \(c_2\) are constants.

We choose \(c_1 = 0.75\) and \(c_2 = 0.25\) in this paper. The objective function is given by

$$\begin{aligned} C_{total} = C_{SSTVD} + \lambda \cdot Reg(\phi ) \end{aligned}$$
(4)

where the variable \(\lambda \) is used to balance the weights put on similarity cost and regularization cost.

The image registration algorithm used in this paper has been shown to have sub-voxel accuracy [2].

2.3 Local Expansion Ratio by Two 4DCT Phases (LER3D)

The SSTVD deformable image registration algorithm is used to find a plausible correspondence map (transformation) \(\phi \) from the 0EX CT image to the 100IN CT image. The Jacobian matrix \(J_{\phi }\) of the transformation \(\phi \) is given by

$$\begin{aligned} J_{\phi } \triangleq \begin{bmatrix} \frac{\partial \phi _1}{\partial x_1}&\frac{\partial \phi _1}{\partial x_2}&\frac{\partial \phi _1}{\partial x_3} \\ \frac{\partial \phi _2}{\partial x_1}&\frac{\partial \phi _2}{\partial x_2}&\frac{\partial \phi _2}{\partial x_3} \\ \frac{\partial \phi _3}{\partial x_1}&\frac{\partial \phi _3}{\partial x_2}&\frac{\partial \phi _3}{\partial x_3} \end{bmatrix} = \mathbb {I}_3 + J_u. \end{aligned}$$
(5)

where \(\mathbb {I}_3\) is the \(3 \times 3\) identity matrix and \(J_u\) is the Jacobian of the displacement field \(u = \phi - \text {Id}\).

Reinhardt et al. [12] proposed to estimate local expansion ratio of the lung at each voxel x by the Jacobian determinant of \(\phi \) (LER3D), i.e.,

$$\begin{aligned} LER3D(x)\triangleq |J_{\phi }(x)| \end{aligned}$$
(6)

2.4 Local Expansion Ratio by All 4DCT Phases (LER4D)

The LER3D measure may underestimate LER when out-of-phase ventilation happens. To solve this problem, we propose to estimate LER by all 4DCT phases, which we refer to as LER4D. We perform a nonrigid registration from each breathing phase to the 0EX phase and denote the resulting Jacobian determinant images by \(\mathbb {J}_1,\cdots ,\mathbb {J}_N\), where N is number of breathing phases of a 4DCT scan. The LER4D measure at each point \(x\in \varOmega \) is given by the ratio of the maximum to the minimum local lung volume:

$$\begin{aligned} LER4D(x) = \max _{\begin{array}{c} i\in \{1,\cdots ,N\} \end{array}} \mathbb {J}_i(x) \Big / \min _{\begin{array}{c} i\in \{1,\cdots ,N\} \end{array}} \mathbb {J}_i(x) \end{aligned}$$
(7)

By definition, the LER4D measure is always greater or equal to the LER3D measure, and we hypothesize that it may provide a more accurate estimate of local lung expansion ratio.

2.5 Out-of-Phase Ventilation

The lung has out-of-phase ventilation when local lung volume change is out of phase with respect to global lung volume change. By definition, out-of-phase ventilation occurs whenever the LER4D measure is greater than the LER3D measure. We define the lung to have out-of-phase ventilation if \(LER4D\ge T\times LER3D\), where \(T>1\) a threshold value. The threshold T is used to reduce the effect of noise. We choose the out-of-phase threshold T to be 1.05, i.e., the LER4D measure is greater than or equal to \(5\%\) of the LER3D measure. We chose 1.05 as the threshold value since the standard deviation of the ratio of two LER3Ds for repeated 4DCT scans is about 0.05.

3 Results

For each of the 14 human subjects, we compute the LER3D measure and the LER4D measure by Eqs. 6 and 7, respectively. Figure 1 shows the cumulative 2D histogram of LER3D image versus LER4D image for all 14 subjects. A logarithmic scale is used for visualization. Overlaid on this histogram are the functions \(y = x\) and \(y = 1.05x\). Notice all points lie above the \(y = x\) solid line since \(LER4D\ge LER3D\) by our definition. Points that lie above the \(y = 1.05x\) dashed line are defined to be out-of-phase ventilation and points that lie between the two lines are defined to be in-phase ventilation.

Fig. 1.
figure 1

Cumulative 2D histogram of LER3D images versus LER4D images for 14 subjects. A logarithmic scale is used for visualization. Low-function (high-function) regions are defined as regions that have less (greater) than \(10\%\) expansion and are denoted by the dashed lines at 1.1. Region A is where the lung has in-phase ventilation, region B, C and D are regions where the lung has out-of-phase ventilation. The lung is considered as low-function by both LER3D and LER4D in region B, the lung is considered as low-function by LER3D while high-function by LER4D in region C, the lung is considered as high-function by both LER3D and LER4D in region D. We use P(A), P(B), P(C), and P(D) to denote the percentages of the voxels in regions A, B, C and D, respectively.

Full size image

We divide the 2D plane in Fig. 1 into four regions: A, B, C and D, where region A corresponds to in-phase ventilation, and regions B, C and D correspond to out-of-phase ventilation. Low-function (high-function) regions are defined as regions that have less (greater) than \(10\%\) expansion and are denoted by the dashed lines at 1.1. On average for all 14 subjects, \(80.7\%\) of all voxels are in region A, i.e., \(80.7\% \) of the lung has in-phase ventilation. Conversely, \(19.3\%\) of the lung has out-of-phase ventilation. \(3.8\%\) of all voxels are in region B, i.e., on average \(3.8\%\) of the lung volume is labeled as low-function by both of the LER3D and LER4D measures and at the same time has out-of-phase ventilation. \(9.6\%\) of all voxels are in region C, i.e., on average \(9.6\%\) of the lung volume is labeled as low-function by the LER3D measure while high-function by the LER4D measure. \(5.9\%\) of all voxels are in region D, i.e., on average \(5.9\%\) of the lung volume is labeled as high-function by both of the LER3D and LER4D measures and at the same time has out-of-phase ventilation.

Table 1 summarizes the percentages of the lung volume for regions A, B, C and D for each of the 14 subjects. This table shows that every subject has out-of-phase ventilation in more than \(10\%\) of the lung volume. The last row of Table 1 shows the mean (± standard deviation) over all 14 subjects for all regions.

Figure 2 shows the location of the tumor and the spatial distribution of regions A, B, C and D of two subjects. Both subject 1 and subject 8 have less out-of-phase ventilation than the average. Note that for subject 1, large regions near the tumor have out-of-phase ventilation. This means it is possible that some of the tissue that was classified as low-functioning by LER3D should have been classified as high-functioning when using LER4D due to out-of-phase ventilation. For this subject, using the LER4D measure instead of the LER3D measure will make a difference to the functional avoidance treatment plan. However, for subject 8, out-of-phase regions are far from the tumor and using the LER4D measure will not make much difference to functional avoidance plan derived from LER3D.

Table 1. Percentages of the lung volume for Regions A, B, C and D.
Full size table
Fig. 2.
figure 2

Spatial distribution of regions A, B, C and D of two subjects. Region A is where the lung has in-phase ventilation; region B is where the lung has out-of-phase ventilation and is labeled as low-function by both LER3D and LER4D; region C is where the lung is labeled as high-function by LER4D while low-function by LER3D; region D is where the lung has out-of-phase ventilation and is labeled as high-function by both LER3D and LER4D.

Full size image

4 Discussion and Conclusions

We are conducting a clinical trial at UW-Madison (NCT02843568) that uses both conventional and functional avoidance treatment plans to treat patients with lung cancer. In this trial, the LER3D method is used in functional avoidance RT to minimize radiation delivered to the regions of high-function. The goal of this trial is to show functional avoidance RT can reduce damage to healthy lung tissue compared to conventional treatment. In total, 120 patients will be recruited in this clinical trial; half will be treated with conventional RT and the other half will be treated with functional avoidance RT. In this paper, we have shown that out-of-phase ventilation was common in all 14 human subjects that we investigated. When large difference between LER4D and LER3D happens near the tumor (e.g. subject 1), then using the LER4D measure may preserve more healthy lung tissue than the LER3D measure in functional avoidance RT.

In conclusion, we investigated the out-of-phase ventilation of 14 human subjects. Our study shows that all subjects had out-of-phase ventilation in more than \(10\%\) of the lung volume. Our results show that on average \(19.3\%\) of the lung has out-of-phase ventilation; \(3.8\%\) of the lung has out-of-phase ventilation and is labeled as low-function by both of LER3D and LER4D; \(9.6\%\) of the lung is labeled as low-function by LER3D while high-function by LER4D; and \(5.9\%\) of the lung has out-of-phase ventilation and is labeled as high-function by both LER3D and LER4D.