1. IntroductionChronic obstructive pulmonary disease (COPD) refers to a collection of inflammatory lung diseases including chronic bronchitis and emphysema [1]. COPD continues to be the third leading cause of death worldwide [1], with an ever-increasing disease burden [2]. People with COPD may experience high morbidity, high healthcare costs, poor quality of life, activity limitation, and exacerbations [3]. Although COPD is not curable, available treatments can help relieve symptoms, improve exercise capacity, improve the quality of life, reduce the risk of death, and reduce the cost of healthcare [4]. Spirometry—a common pulmonary function test—is currently the gold standard for diagnosing and staging COPD. Based on spirometry, the Global Initiative for Chronic Obstructive Lung Disease (GOLD) recommends grading the disease into four categories, ranging from GOLD 1 to GOLD 4 [5]. However, spirometry alone cannot capture the heterogeneous manifestations of COPD, which calls for better diagnostic methods. Towards advancing our understanding of the disease, several multi-center studies such as SPIROMICS [6] and the Genetic Epidemiology of COPD, COPDGene [7] are underway. Along with spirometry, these studies rely on computed tomography (CT) scans of the subjects, which have led to the development of several imaging biomarkers for characterizing COPD.

Biomechanical properties of the lung can be used to characterize lung function and are important because they provide information on whether or not the lung is functioning normally or abnormally. One way to indirectly measure the biomechanical properties of the lung at the local level is to analyze the pointwise correspondences between inspiration and expiration CT image volumes. This can be achieved using image registration. Image registration is used to find the point-to-point correspondence between two images in the form of a transformation or deformation vector field (DVF). Biomechanical measurements (also called biomechanical biomarkers) can then be extracted from the transformation to describe the local expansion or contraction of the lung during the breathing cycle.

A complementary approach to extracting biomarkers from CT images is to compute CT intensity-based disease biomarkers. Intensity-based disease biomarkers are computed from CT image intensities with or without using image registration. These methods extract biomarkers at each CT voxel based on the voxel intensity or based on the intensity pattern (texture) of a local region of voxels. Intensity-based biomarkers are widely used to infer disease patterns, but they are susceptible to image noise and do not provide biomechanical measurements. An example of this type of analysis includes defining emphysematous regions by thresholding the inspiratory scan below −950 Hounsfield units (HU), while the cut-off for air-trapping regions in an expiratory scan is −856 HU. This technique is often referred to as CT densitometry, and is prone to noise arising from changes in volume of acquisition, CT dosage, and scanning parameters [8,9,10]. To overcome these challenges, texture-based CT analysis methods have been developed to assess parenchymal integrity and its relation to disease progression. Sørensen et al. [11] developed a CT texture classification method to assess COPD. A similar method was used to assess pulmonary emphysema using local binary patterns [12]. These methods rely on globally annotated CTs, and lack local functional information related to parenchymal disintegration. A texture-based approach called CALIPER was used to quantify interstitial lung abnormalities and showed good correlations between the fraction of abnormal lung texture and lung functional measurements [13]. More recently, a deep convolutional neural network (CNN) was employed for staging COPD and predicting disease progression [14]. Although this approach was able to learn effective representations associated with COPD severity and progression, it failed to provide insights into the regional distribution of lung function abnormalities.Several image-registration-based biomechanical biomarkers have been derived from the DVF to understand, diagnose, and stage COPD. Galbàn et al. [15] used CT intensity and image registration to define the parametric response mapping (PRM), a two-dimensional histogram that measures the amounts of normal tissue, small airways disease, and emphysema within the lung. Amelon et al. [16] extracted three biomechanical indices from the DVF of pulmonary inspiration-to-expiration registration to characterize local and global lung deformation: the Jacobian determinant (J); the anisotropic deformation index (ADI), and the slab-rod index (SRI) [16]. Bodduluri et al. [17] evaluated the predictive performance of J, ADI, and strain extracted from image registration and compared their performance with conventional CT texture and densitometry [18,19] features. They demonstrated that for a complex task such as COPD severity prediction, the biomechanical features performed better than the conventional CT texture and density features. The current study differs from that of Bodduluri et al. in many ways. The work of Bodduluri et al. analyzed data from the Genetic Epidemiology of COPD (COPDGene) study whereas the current study analyzes data from the SPIROMICS study. The SPIROMICS study collected CT scans at total lung capacity (TLC) and residual volume (RV) whereas the COPDGene study collected CT scans at functional residual capacity (FRC), as opposed to RV. Registering TLC to RV CT images presents a more challenging problem than registering TLC to FRC because there is a comparatively larger shape change from TLC to RV than TLC to FRC. The current study extends the work of Bodduluri et al. from a single image-registration method to multiple image-registration algorithms. Finally, the current study analyzes nearly twice as many data sets compared to that of the previous study by Bodduluri et al.In another study, Bhatt et al. [20] showed an agreement between registration-based mechanics and spirometric measures of lung function and concluded that dual volume (inspiration–expiration) biomechanical measures are better indicators of declining lung function and emphysema than spirometry alone. Another important work by Bodduluri et al. [21] identified the mean of J to be significantly associated with different measures of lung function including forced expiratory volume in 1 s (FEV1), emphysema, and six-minute walk distance (6MWD). In a recent work by Pan et al. [22], correlations between air-trapping regions or emphysema regions with mean J and mean ADI were found. This work showed that the trends of mean J and ADI decreased monotonically as GOLD stages increased for one registration algorithm. While the aforementioned studies point towards clinical effectiveness of registration-based mechanics, they were limited to global features derived from these biomechanical measures. Recently, Chaudhary et al. [23] showed the predictive capabilities of global statistical features extracted from J, ADI, and SRI across the four registration methods used in this study. The work by Chaudhary et al. [23] differs from the current paper in that it was limited to the development of classification methods for predicting COPD GOLD stage.Several image registration algorithms have recently been proposed for analyzing CT lung images of COPD subjects, such as the total variation regularization method [24], key-points-based method [25], and Markov random field (MRF)-based discrete optimization [26]. In addition, deep learning image registration techniques bear the potential to perform image registration orders of magnitude faster than traditional iterative image registration techniques [27,28,29,30,31,32].The goal of the current study was to understand and evaluate the robustness of different lung biomechanical measures across different image registration methods over a varying disease severity. In this study, we compared and contrasted five diverse image registration algorithms with different similarity costs and transformation models to investigate the impact of different algorithms on the biomechanical measures extracted from the DVF. Image registration performance was evaluated using the lung lobe Dice coefficient (LDC) to measure the overlap of the five lobes after registration; the worst 10% surface error (W10SE) to measure how well the lung surfaces aligned after registration; and the vessel tree position error (VTPE) and symmetric closest skeleton error (SCSE) to measure how well the lung vessel trees aligned. The transformations were used to compute the J, ADI and SRI measures [16] in order to analyze and characterize the lung function of individuals with COPD at varying GOLD stages. 2. Materials and Methods 2.1. DataThe CT images used in this study were part of the SubPopulations and Intermediate Outcome Measures in COPD Study (SPIROMICS) with institutional review board (IRB) approval number 201003733 [6]. SPIROMICS is an ongoing prospective cohort study designed to identify novel clinical stratifications in subjects with COPD. CT images were collected at baseline, one-year, three-year, and five-year follow-ups from fourteen university-based clinical centers across the United States [6].In this study, we analyzed CT images from current and former smokers across varying degrees of disease severity, as defined by GOLD stage. The GOLD disease staging system for COPD ranges from GOLD 1 (mild), to GOLD 2 (moderate), GOLD 3 (severe) and GOLD 4 (very severe) [5]. GOLD 0 subjects are asymptomatic smokers without airflow obstruction but at risk for COPD due to their smoking history [6]. We analyzed the baseline scans of 245 subjects chosen randomly from the 14 clinical sites, with 49 subjects from GOLD 0, 50 subjects from GOLD 1, 49 subjects from GOLD 2, 50 from GOLD 3, 47 subjects from GOLD 4. At each visit, a pair of 3D breath-hold CT scans were acquired; one at total lung capacity (TLC), and the other at residual volume (RV) [33]. The SPIROMICS imaging protocol uses a CT dose index to standardize exposure across scanners in different sites. The slice collimation was set at 0.6 mm, rotation time 0.5 s and pitch to 1.0. Philips B, GE Standard, and Siemens B35 reconstruction kernels were used [33]. The resolution of the CT scans was approximately 0.6×0.6×0.5 mm3, and the image size was 512×512 per slice, with 500 to 600 slices per image. Full details on the scanning protocol are described by Sieren et al. [33]. 2.2. PreprocessingFigure 1 shows the preprocessing pipeline applied to each CT image volume.All CT image volumes were resampled to be isotropic with spacing 1×1×1 mm3. The image volumes were cropped based on the 3D bounding box containing the union of the lung regions of the inspiration and expiration scans to reduce computer memory requirements and computation time. A multi-resolution convolutional neural network [34] was used to generate a segmentation of the entire lung from the CT volume. The lobes were segmented from the CT volume using the FissureNet deep learning method [35,36]. We used the recommended parameters from the cited papers for the lung and lobe segmentations.The lobe segmentations were tessellated into triangles using the surface mesh algorithm [37] implemented in the Computational Geometry Algorithms Library (CGAL) (https://www.cgal.org (accessed on 15 September 2017)) that implements the Variational Shape Approximation (VSA) method [38]. The mesh size was set to 2 mm2 with smoothing to give accurate representation of the surfaces. The blood vessel trees were segmented from the CT images using the vesselness filter developed by Jerman et al. [39]. In this method, the eigenvalues of the Hessian matrix are computed from the intensity values of the CT image at each voxel location. Tubular structures are identified at voxel locations that have one near-zero eigenvalue and two non-zero eigenvalues with similar magnitudes. Vessels were detected at different scales by computing the Hessian matrix at different scales. For the vessel segmentation, we started from the parameters recommended in the cited papers, then slightly fine-tuned them. A binary vessel segmentation was computed from the vesselness probability map by thresholding. The threshold was computed using Otsu’s method. Next, the skeletons were extracted from the binary vessel segmentation using the binary 3D thinning algorithm implemented with Insight Toolkit (ITK) in C++ [40] with default parameters. The TLC and RV CT data sets used in this study were collected in register with each other. Therefore, no affine registration was performed in this study before applying the nonrigid image registration algorithms used in this study. 2.3. Image Registration AlgorithmsWe selected five image-registration algorithms to cover a wide range of image registration algorithms. In general, image registration algorithms consist of three major components: an overall cost function to be minimized, a transformation model, and an optimization method. In this study, we focus on the two former components. In general, the overall cost function is a linear combination of a similarity cost and a regularization cost of the transformation. There are various choices for these two components, and therefore different combinations. In terms of similarity cost functions, we can divide them into two major categories: intensity-based and feature-based. Commonly employed intensity-based similarity cost functions used for matching CT scans include, but are not limited to the sum of squared difference (SSD) [41,42], cross correlation (CC) [43], mutual information (MI) [44], and the sum of squared tissue volume difference (SSTVD) [45,46,47]. Examples of intensity-based featured-matching cost functions are the SSD vesselness measure [45], and SSD of lobar segmentations [48] of the lung. For shape-based feature matching, the shapes could be corresponding or non-corresponding points, curves, and surfaces. A popular example is the Iterative Closest Point (ICP) algorithm, which can be used to match point clouds [49], curves (contours), and surfaces [50]. Another group of shape-matching methods adopt the concept of currents and varifolds, initially proposed by Charon et al. [51], and extensively studied by Durrelman et al. [52,53,54,55,56,57]. The idea of applying currents to pulmonary registration was first proposed by Gorbunova et al. [58], and later Pan et al. extended the study by using varifolds representation [59,60]. The advantages of representing shapes for registration purpose with currents and varifolds are (1) point correspondence is not required to define the distance between landmarks, curves, surfaces, and intensity with currents and varifolds representations; and (2) landmarks, curves, surfaces, and intensity can be unified in one framework.The transformation model plays a critical role in an image registration algorithm since it determines the properties of the DVF. The aim of this study was to register TLC inspiration scans to RV expiration scans of the lungs. Non-rigid image registrations—also often referred to as deformable image registration (DIR)—can be classified based on physical models, interpolation theory models, knowledge-based models and task-specific models [61]. In this manuscript, we classify the transformations as either small or large deformation models. The main differences between small and large deformation models are that the large deformation model allows for more curved particle paths from the moving image to target image, and guarantees a one-to-one correspondence between the moving and target images. The large deformation diffeomorphic metric mapping (LDDMM) model [41] refers to the large deformation model throughout this manuscript. The DVFs that are parameterized by basis functions are categorized into small deformation model, such as Fourier series [42], Thin-Plate spline [62], B-spline [46,63,64,65], etc.To cover a wide category of image registration algorithms, we selected five algorithms as follows: the Sum-of-Squared-Tissue Volume-Difference (SSTVD) [45], Geodesic Density Regression (GDR) [66,67], Greedy Symmetric Normalization (GSyN) [43], Pulmonary blood Vessel and lobe Surface Varifold (PVSV) [51,57,60], and the Population Learning followed by One Shot Learning (PLOSL) [32]. 2.3.1. SSTVD

The SSTVD algorithm was selected as one of the algorithms used because it was developed particularly for registering pulmonary CT and won the Computed Tomography Ventilation Imaging Evaluation 2019 (CTVIE19) Grand Challenge at the American Association of Physicists in Medicine (AAPM) 2019 annual meeting. The SSTVD algorithm is currently being used in the Functional Avoidance Radiation Therapy Clinical trial NCT02843568 and has been used to treat more than 50 subjects in that trial. The SSTVD cost function models the local intensity changes seen in the CT images of the lung due to breathing and provides good correspondence information between inspiration and expiration CT scans.

The overall cost of the SSTVD algorithm is the sum of three terms: (1) the sum of squared tissue volume difference (SSTVD) similarity cost for matching pulmonary CT intensity [46,47], (2) the sum of squared vesselness measure difference (SSVMD) similarity cost for matching lung vessels, and (3) a regularization cost on the cubic B-Splines parameterized DVF to guarantee a smooth and plausible transformation [45]. The main idea behind SSTVD is that lung tissue volume remains relatively constant over the breathing cycle while the volume of air in the lung does not. This method can be considered as a lung-specific intensity-based and feature-based combined small deformation algorithm. The cost function of SSTVD algorithm [22,45,46,47,68,69] is given by:

C=γ11|Ω|∫Ω(If(x)−|Jφ−1(x)|Im(φ−1(x)))2dx       +γ21|Ω|∫Ω(Ifvm(x)−Imvm(φ−1(x)))2dx+γ3∫Ω||L(u(x))||2dx,

(1)

where If and Im are the fixed and moving tissue density images, respectively; Ifvm and Imvm are the fixed and moving vesselness images, respectively; L is the linear elasticity differential operator of the type L=(−αΔ+γ)βIn×n; Ω⊂R3 is the domain of If and Im; The |Jφ−1(x)| term is the Jacobian determinant of φ−1 and used to accommodate intensity changes in the CT images due to changing air content by ensuring the total tissue volume of the lung remains constant; φ−1(x) is the transformation from the fixed to the moving coordinate system; u(x)=x−φ−1(x) is the displacement field; and γ1, γ2 and γ3 are weights that control the relative importance of each term in the cost function. 2.3.2. GDR

The GDR algorithm was chosen since it is an example of a large deformation image registration algorithm that uses the SSTVD cost function. The SSTVD and GDR image registration algorithms were chosen to compare and contrast the effect of using small vs. large deformation transformation models. The SSTVD and GDR image registrations also differed in that the SSTVD algorithm included the SSVMD similarity method while the GDR did not.

The Geodesic Density Regression (GDR) algorithm [66,67] finds the geodesic path in image space to align pulmonary CT images. The overall cost of the GDR includes the SSTVD term and the regularization term on the velocity fields, which are used to parameterize the DVF as a flow of diffeomorphism. The differences between GDR in this study and the previous SSTVD method are that the GDR algorithm is based on the LDDMM large deformation model [41], and uses the SSTVD similarity cost but does not use the SSVMD similarity cost. GDR can be categorized as a lung-specific intensity-based large deformation model. The cost function of the GDR is given by

C=γ11|Ω|∫Ω(If(x)−|Jφ−1(x)|Im(φ−1(x)))2dx+γ2∫01||Lv(t)||2dt,

(2)

where the first term is the SSTVD cost function using the same definitions as in Equation (1). The time-varying velocity field v(t) parameterizes the transformation φ by the ordinary differential equation (O.D.E.):

∂∂tϕv(x,t)=vt(ϕtv(x)),

(3)

where t∈[0,1], ϕ0=Identity and φ=ϕ1v (see [41] for details). The second term in Equation (2) is the regularization cost and is used to find the geodesic path in image space between If and Im. The large deformation model estimates a diffeomorphic transformation φ that guarantees the transformation is smooth with a smooth inverse. A diffeomorphic transformation preserves the topological properties of objects in an image, i.e., connected structures remain connected, disjoint structures remain disjoint, and the smoothness of curves and surfaces are preserved. The relative importance of each term compared to each other are controlled by weights γ1 and γ2. 2.3.3. GSyNThe Greedy Symmetric Normalization (GSyN) algorithm [43] is part of the open source ANTs (Advanced Normalization Tools) toolbox [70]. GSyN aligns two images with the deformable diffeomorphic transformation (LDDMM) model [43] in a symmetric scheme. The GSyN algorithm differs from the SSTVD and GDR methods in that it estimates the transformation from moving image to fixed image by deforming both images to the mid-point between the two images. The GSyN algorithm fits into the intensity-based large-deformation model image registration algorithm class. The transformation model of GSyN differs from the GDR algorithm in that GDR uses a time-varying velocity field whereas the GSyN uses a stationary velocity field. The difference is that a stationary velocity field is constant for all time t in Equation (3) and therefore uses less parameters to represent the transformation compared to the non-stationary velocity field parameterization and has fewer degrees of freedom than a time-varying velocity field.

The GSyN algorithm used the normalized cross correlation (NCC) similarity cost function compared to the SSTVD similarity cost used by the GDR algorithm. The NCC cost function does not accommodate changes in CT intensity due to breathing in contrast to the SSTVD cost function. This means that the NCC cost function may not always provide the proper correspondences between the moving and target images. To understand this, assume two regions of the parenchyma have the same intensity in the moving image. Next, assume that one region expanded and the other region did not as a result of breathing. The intensity of the expanded region will become darker due to an increase in air filling the regions (i.e., reduced tissue density) and the other region will show no intensity change. The NCC cost is then forced to do its best to match one intensity in the moving image with two different intensities in the target image.

The overall GSyN cost function is defined as

C=γ1∫ΩNCC(Im(φm(x)),If(φf(x)))dx+γ2∫00.5(||Lvm(t)||2+||Lvf(t)||2)dt,

(4)

where x is a coordinate of the mid-point coordinate system of Im and If. Note that we use the fixed and moving notation to describe the GSyN method to be consistent with the other registration algorithms even though both images can be thought of as moving. The transformations φm and φf are from the mid-point coordinate system to the coordinate systems of Im and If, respectively. Equation (3) is used to parameterize φm and φf by vm and vf, respectively. NCC is the normalized cross-correlation cost function [43]. The relative importance of each term compared to each other are controlled by weights γ1 and γ2. 2.3.4. PVSVThe Pulmonary blood Vessel and lobe Surface Varifold (PVSV) algorithm [60] is a shape feature-based (LDDMM) large deformation registration approach that aligns varifold representations of the lung blood vessel skeletons and lobe surfaces. The process for extracting the blood vessel trees and lobe surfaces used for the PVSV algorithm is the same as described in Section 2.2. The skeleton of each branch of the vessel tree was parameterized by a quadratic line segment using a least squares fitting process. The skeletons of the vessel trees and the lobe surfaces were then represented with delta Dirac varifolds using the procedure in [60]. Varifolds were introduced to the field of computational anatomy to overcome the orientation problem of currents [51]. The advantages of matching shape with a varifolds representation are: it does not require pre-knowledge of the point correspondence between shapes, and it can match shapes with large geometric changes.

The PVSV algorithm was chosen to determine whether meaningful biomechanical properties could be extracted using only lung surfaces and vessel tree represented via varifolds. The advantage of using lung surfaces and vessel trees for registration is that it reduces the amount of information needed to represent the lung and that these features are invariant to CT intensity changes caused by breathing. Another advantage of this approach is that it can accurately match surface and vessel structures but it has the disadvantage of needing to interpolate the transformation between these features.

The vessel tree skeletons and lung surfaces are registered as a shape complex with varifolds representation [57], and the overall PVSV cost function is given by

C=γ1||φ*Cm−Cf||W′2+γ2||φ*Sm−Sf||W′2+γ3∫01||Lv(t)||2dt,

(5)

where Cm and Cf are the blood vessel tree skeletons represented with varifolds, while Sm and Sf are the lung surfaces with varifolds representation. φ*Cm is the push forward action of the transformation φ on the varifold Cm that transports Cm into the space of Cf, and same for φ*Sm. The relative importance of each term is controlled by the weights γ1, γ2 and γ3. 2.3.5. PLOSLThe PLOSL is a fast unsupervised-learning-based framework developed for 3D pulmonary CT based on population learning (PL) and one-shot learning (OSL) [32]. It uses the same tissue volume preserving and vesselness constraints similarity metrics as the SSTVD method. PLOSL has two training stages: (1) PL which serves as a base model; (2) OSL which generates an individual specific model. The PLOSL algorithm has been shown to produce state-of-the-art deep-learning-based image registration performance [32] with comparable accuracy to traditional iterative image registration algorithms and other deep-learning-based methods for lung CT registration.The cost function for the PLOSL network consists of three components: SSTVD, SSVMD, and a regularization cost as follows:

C=γ11|Ω|∫ΩIf(x)−|Jφ−1(x)|Im(φ−1(x))2dx       +γ21|Ω|∫ΩIfvm(x)−Imvm(φ−1(x))2dx+γ3∫Ω||∇u(x)||L22dx,

(6)

where the first two terms are the same as the two terms of the SSTVD method, and the third term is a regularization term and the operator ∇=[∂∂x1,∂∂x2,∂∂x3]. 2.4. Image Registration Parameters

The registration algorithms employed in the study were fundamentally different and therefore, the corresponding parameters and multi-resolution setups were different for each method. All the registration algorithms used a coarse-to-fine multi-resolution framework to prevent themselves from becoming stuck in local minima, speed up computation, and achieve better registration performance. Parameters were determined by optimizing image registration performance by hand using a few typical data sets and hundreds of different parameter combinations for each algorithm. We optimized the algorithms independently to achieve the best performance. The following algorithm-specific parameters were used in this work.

The multi-resolution framework used for the SSTVD algorithm consisted of six image resolutions at 1/8, 1/8, 1/4, 1/4, 1/2 and 1 of the original resolution, respectively. The corresponding B-Spline node spacings were 8, 4, 8, 4, 4 and 4 mm. The cost function described in Equation (1) used weights γ1=γ2=1. The differential operator used was L=−0.75∇2−0.25∇(∇·).

The multi-resolution framework employed for the GDR algorithm consisted of three image resolutions equal to 1/8, 1/4 and 1/2 of the original resolution, respectively. The corresponding standard deviation of the deformation Gaussian kernel size for the transformation were 60, 30, and 15 mm, respectively. A total of ten time points were used for the large deformation model, and the weight on the image cost γ1=0.06, and the regularization cost γ2=1 at each resolution.

The multi-resolution framework used for GSyN algorithm consisted of five image resolutions. The input images were smoothed with a Gaussian kernel with standard deviation 4, 3, 2, 1, and 0 mm and down sampled by a factor of 1/16, 1/8, 1/4, 1/2, and 1, respectively.

The multi-resolution framework used for PVSV algorithm consisted of two resolutions. At the coarse resolution, the shape kernel for the vessel tree skeletons and lung surfaces were set to 5 and 15 mm, respectively. At the fine resolution, the shape kernel for the vessel tree skeletons and lung surfaces were set to 3 and 10 mm, respectively. The corresponding deformation kernel sizes were 40 and 30 mm, respectively. The weights on the vessel tree and surface cost terms were set to γ1=10,000 and γ2=25 at coarse resolution, respectively, and to γ1=4 and γ2=0.015, respectively, at the fine resolution. The weight on the regularization term is γ3=1 for the two resolutions.

The population learning U-Net of PLOSL was trained with 2042 subjects randomly selected from the SPIROMICS database. There was no overlap between training data used in PLOSL and the data sets (247 subjects) analyzed in this study. The weights on the terms of the cost function Equation (6) were set to γ1=1, γ2=1, and γ3=0.01, respectively. The number of one-shot iterations was set to 50. The methods were implemented in TensorFlow version 2.3.0 and run on an NVIDIA GeForce GTX 2080Ti GPU with 11 GB of memory. The ADAM optimizer with a learning rate of 10−4 was used for training. 2.5. Image Registration Performance

A total of four measures were used to evaluate image registration performance: (1) lung Lobe Dice Coefficient (LDC); (2) Worst 10% Surface Error (W10SE); (3) Vessel Tree Position Error (VTPE); and (4) Symmetric Closest Skeleton Error (SCSE).

LDC measures the overlap of the moving and fixed lung lobes. The Dice Coefficient (DC) of two overlapping regions L1 and L2 is given by DC=2|L1∩L2||L1|+|L2|. The human lung has five lobes: right upper lobe, right middle lobe, right lower lobe, left upper lobe, and left lower lobe. The LDC is the average DC of the five lobes.

W10SE is an error measurement of the alignment of two lung surfaces. The W10SE captures the average surface error where two surfaces disagree the most. The process for extracting the lung lobe surface triangulation used for this evaluation is described in Section 2.2. For each vertex p of the triangles that parametrize the surface S1, its closet point q on S2 is found, and the distances between p and q are stored. It is worth noting that q is not necessarily a vertex of S2 and may instead be on the face of a triangle. As a next step, the distances for all the vertices of S1 are sorted from large to small, and the mean of the largest 10% distance is computed, which is denoted as W10S1→S2. Similarly, we compute the mean of worst 10% distance in the opposite direction, i.e., we rank the closest-point distances computed from each vertex of S2 to S1, then compute the mean of worst 10% again, and denote it as W10S2→S1. The W10SE is defined as: W10SE=12(W10S1→S2+W10S2→S1).VTPE measures the error of the lung blood vessel alignment. The process for extracting the blood vessel tree skeleton used for this evaluation is described in Section 2.2. The closest-point distances from the fixed to the moving segmentation were computed for each voxel in the fixed segmentation and averaged. This process was repeated to compute the average closest-point distance from the moving to the fixed segmentation. The VTPE is defined as the average of these two closest-point distances.

SCSE measures the alignment of the skeletons of the vessel tree. SCSE is similar to VTPE, except the symmetric closest-point distance is computed using the vessel tree skeletons instead of binary vessel segmentations.

2.6. Biomechanical MeasuresThis section briefly summarizes the calculation of J, ADI, and SRI [16]. Let F(x) be the deformation gradient tensor at the point x in the lung, i.e., the gradient of the transformation from the inspiration-to-expiration CT images at x. For each point x in the lung, the quantity FTF is the right Cauchy–Green deformation tensor with corresponding eigenvalues λi2 such that λ1≥λ2≥λ3.

The Jacobian determinant J(x) measures the volume change at x and can be computed as J=λ1λ2λ3. A value of J(x)>1 corresponds to local volume expansion at x and a value of J<1 corresponds to local volume contraction at x.

ADI measures the magnitude of anisotropy deformation, i.e., it describes how much stretching occurs along one or two directions in a three-dimensional space. The formula for computing ADI is: ADI=(λ1−λ2λ2)2+(λ2−λ3λ3)2. If the volume changes uniformly along each direction (isotropically), then ADI equals zero. The larger the ADI value, the more anisotropic the deformation.

SRI is a measure of the nature of anisotropy, i.e., it measures whether the volume change is predominant along one direction (rod-like) or two directions (slab-like). SRI is computed using the formula: SRI=tan−1λ3(λ1−λ2)λ2(λ2−λ3)π/2. Note that the value of SRI ranges from zero to one where zero corresponds to a slab-like shape and one corresponds to a rod-like shape. SRI is undefined when ADI is zero.

Figure 2 illustrates the relationships between ADI and SRI. The α and β axes are defined as following: α=J−1|J−1|(λ2λ3−1), and β=J−1|J−1|(λ1λ2−1) [16]. Notice that the terms in the parentheses for both α