International Journal of Biomedical Imaging

Volume 2012 (2012), Article ID 285136, 17 pages

doi:10.1155/2012/285136

Research Article

Improving Intensity-Based Lung CT Registration Accuracy Utilizing Vascular Information

Kunlin Cao,¹ Kai Ding,² Joseph M. Reinhardt,³ and Gary E. Christensen4,5

¹Biomedical Image Analysis Lab, GE Global Research Center, Niskayuna, NY 12309, USA

²Department of Radiation Oncology, Virginia Commonwealth University, Richmond, VA 23298, USA

³Department of Biomedical Engineering, The University of Iowa, Iowa City, IA 52242, USA

⁴Department of Electrical and Computer Engineering, The University of Iowa, Iowa City, IA 52242, USA

⁵Department of Radiation Oncology, The University of Iowa, Iowa City, IA 52242, USA

Received 18 June 2012; Revised 14 September 2012; Accepted 28 September 2012

Academic Editor: Ayman El-Baz

Copyright © 2012 Kunlin Cao et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Accurate pulmonary image registration is a challenging problem when the lungs have a deformation with large distance. In this work, we present a nonrigid volumetric registration algorithm to track lung motion between a pair of intrasubject CT images acquired at different inflation levels and introduce a new vesselness similarity cost that improves intensity-only registration. Volumetric CT datasets from six human subjects were used in this study. The performance of four intensity-only registration algorithms was compared with and without adding the vesselness similarity cost function. Matching accuracy was evaluated using landmarks, vessel tree, and fissure planes. The Jacobian determinant of the transformation was used to reveal the deformation pattern of local parenchymal tissue. The average matching error for intensity-only registration methods was on the order of 1 mm at landmarks and 1.5 mm on fissure planes. After adding the vesselness preserving cost function, the landmark and fissure positioning errors decreased approximately by 25% and 30%, respectively. The vesselness cost function effectively helped improve the registration accuracy in regions near thoracic cage and near the diaphragm for all the intensity-only registration algorithms tested and also helped produce more consistent and more reliable patterns of regional tissue deformation.

1. Introduction

Image registration is used to find the spatial correspondence between two images and plays an important role in pulmonary image analysis. In a sequence of pulmonary scans, image registration provides the spatial locations of corresponding voxels. The computed correspondences describe the motion of the lung between a pair of images at the voxel level. Registration of lung volumes across time or across modalities has been utilized to establish lung atlases [1], estimate regional ventilation and local lung tissue expansion [2–5], assess lobar slippage during respiration [6, 7], and measure pulmonary function change following radiation therapy [8].

Lung registration methods are typically intensity-based [2, 5, 9–13] or feature-based [14–18]. Intensity-based methods consider intensity patterns of the whole lung regions to define similarity measures. They take advantage of the strong contrast between the lung parenchyma and the chest wall, and between the parenchyma, the blood vessels and larger airways. Commonly used intensity-based methods include minimizing intensity difference [3, 10], maximizing mutual information or normalized cross correlation [5, 9], and preserving tissue volume or lung mass [12, 13]. Since intensity-based methods do not use anatomical knowledge, these methods can get stuck in local minima resulting in mismatches of important anatomical structures such as bifurcations of smaller airway and vessel branches. On the other hand, feature-based methods define transformations utilizing corresponding features derived from original images. They usually utilize corresponding landmarks and local intensity patterns [14, 16, 17] and surfaces correspondences [14, 15, 18]. However, due to the sparsity of features, good alignment of features can not guarantee satisfactory matching accuracy for all lung regions.

Since registration methods using either intensity-only or feature-only registration have their limitations, it is desirable to design lung registration methods that utilize both intensity and feature information [19–23]. It has been shown that hybrid registration methods that combine intensity and feature information can help improve matching accuracy. Most of these methods incorporate distributed landmark pairs selected at airway or vessel branch points, that were identified manually or semiautomatically [24], centerline of the airway and vessel tree structures, and lung surface information. These feature extractions may be difficult and time-consuming tasks. Therefore, fast feature extractions and effective methods to utilize feature information are needed to improve registration accuracy.

In this paper, we couple intensity and feature information together to match 3D lung CT images acquired during breath hold at two different levels of inflation. We propose a feature-based similarity criterion utilizing the information of vessel locations and shapes in the registration process. This vesselness preserving cost function is added to four existing intensity-based similarity costs, and comparison experiments show that this criterion helps improve the registration accuracy. Higher matching accuracy makes the postanalysis of regional tissue mechanical properties more plausible and reliable.

Our preliminary work on the effectiveness of using a vesselness matching similarity term was described in [25–27]. The work presented in this paper provides a complete description of our vesselness matching approach. This paper extends our previous work by describing how to choose weighting factors for the different cost terms, describes a multiresolution optimization scheme, and provides more validation. This paper studies the effect of the vessel matching when used with various intensity similarity metrics such as the sum of squared intensity difference, sum of squared tissue volume difference, mutual information, and normalized cross correlation. This paper also examines the effect of using a linear-elastic constraint to regularize the displacement fields.

2. Material and Methods

2.1. Image Data Sets

In this study, six pairs of volumetric CT data sets from six human subjects in the supine orientation were collected on a Siemens Sensation 64 multidetector CT scanner. Each image pair was acquired during breath-holds near functional residual capacity (FRC) and total lung capacity (TLC) in the same scanning session. For subject 1, data were acquired at functional residual capacity (FRC) with 21.8% of the vital capacity (VC) and total lung capacity (TLC) with 95.6% of the VC. For subject 2, data were acquired at FRC with 30.5% of the VC and TLC with 89.6% of the VC. For subject 3, data were acquired at FRC with 26.3% of the VC and TLC with 95.7% of the VC. For subject 4, data were acquired at FRC with 11.0% of the VC and TLC with 68.9% of the VC. For subject 5, data were acquired at FRC with 25.8% of the VC and TLC with 92.9% of the VC. For subject 6, data were acquired at FRC with 26.5% of the VC and TLC with 102.0% of the VC. Each volumetric data set was acquired at a section spacing of 0.5~0.6 mm and a reconstruction matrix of . In-plane pixel spatial resolution was approximately 0.6 mm × 0.6 mm. The parenchyma regions in the FRC and TLC data sets were segmented using the method described in [28]. Figure 1 gives an illustration of pulmonary CT images with renderings of the lung segmentations.

(a)

(b)

Figure 1: Pulmonary CT images acquired at breath hold (a) maximum exhale and (b) maximum inhale with renderings of the lung segmentations (green objects).

2.2. Image Registration

The goal of image registration is to estimate a transformation that defines the pointwise correspondences between two images and . More formally, let define the domain of the images and . In this work, the transformation is assumed to be a diffeomorphism. Let denote the set of all diffeomorphisms from to . The optimal transformation is estimated by minimizing a cost function that consists of a similarity cost function and a regularizing term, that is,

The similarity cost function describes the characteristics of two images that should agree for corresponding image points. For example, these characteristics often correspond to matching intensity characteristics and features. The regularization cost function is used to enforce desired properties of the transformation such as minimum distortion. The constant is used to balance the influence of the regularization cost with respect to the similarity cost. In general, the similarity and the regularization costs can be decomposed into linear combinations of more specific cost functions.

2.2.1. Intensity-Based Similarity Cost Functions

For intensity-based image registration, it is usually assumed that intensities of corresponding voxels are related to each other in some way. Many criteria to construct the intensity relationship between corresponding points have been suggested as the cost function for aligning two images. Examples of intensity-based cost functions are the mean square difference (MSD), correlation coefficient, mutual information, pattern intensity, and gradient correlation [29, 30]. In this work, we investigated three commonly used intensity-based cost functions and one intensity-based cost function designed for matching lung CT images.

Sum of Squared Difference (SSD)

Minimizing the intensity difference at corresponding points between two images is an intuitive method to register grayscale images. A simple and common cost function is the sum of squared difference (SSD) defined by

where and are the template and target image intensity functions, respectively. denotes the union of lung regions in target image and deformed template image. The underlying assumption of SSD is that the image intensity at corresponding points between two images should be similar. This is true when registering images of the same modality. However, considering the change in CT intensity as air inspired and expired during the respiratory cycle, the grayscale range is different within the lung region in two CT images acquired at different inflation levels. To balance this grayscale range difference, intensity normalization is needed. For example, a histogram matching procedure [31] can be used before SSD registration to modify the histogram of template image so that it is similar to that of target image.

Mutual Information (MI)

Mutual information (MI) similarity cost function can accommodate intensity difference between two images and is therefore well-suited to accommodate the CT intensity change during inflation and deflation of the lung. Mutual information expresses the amount of information that one image contains about the other one. Analogous to the Kullback-Leibler measure, the negative mutual information cost of two images is defined as [9, 32]

where is the joint intensity distribution of transformed template image and target image ; and are their marginal distributions, respectively. The histogram bins of and are indexed by and . The experiments of MI-driven registration use histogram bins to estimate joint distribution.

Normalized Cross Correlation (NCC)

Normalized cross correlation can be used for multimodality registration problems since it is insensitive to a constant multiplicative factor between the images. This cost function measures the pixel-wise cross-correlation between image intensities normalized by the square root of the autocorrelation of each image. Mathematically, the negative normalized cross correlation measure is given by [33]

where the negative sign was added so that the optimal transformation is found by minimization. When the factor of the intensity patterns from two images is a constant, the measure equals −1. Misalignments between the images will result in decrease of the normalized cross correlation, and thus, increase of the similarity cost .

Sum of Squared Tissue Volume Difference (SSTVD)

The sum of squared tissue volume difference (SSTVD) cost function [13] accounts for the variation of intensity in the lung CT images during respiration. Assume that lung is a mixture of two materials: air and tissue/blood (nonair). Then the Hounsfield units (HU) in lung CT images is a function of tissue and air content. From the HU of CT lung images, the regional tissue volume and air volume can be estimated following the air-tissue mixture model by Hoffman and Ritman [34]. Let be the volume element at location . Then the tissue volume within the volume element can be estimated as , where we assume that and . The intensity similarity cost function SSTVD is defined as [13]

The Jacobian of a transformation estimates the regional volume changes resulted from mapping corresponding regions. Thus, the tissue volumes in image and are related by .

2.2.2. Feature-Based Similarity Cost Function

Feature information extracted from the grayscale image is important to help guide the registration process. During the respiration cycle, blood vessels keep their tubular shapes and tree structures. Therefore, the shape and spatial information of vessels can be utilized to help improve the registration accuracy. In CT images, blood vessels have larger intensity values than that of parenchyma tissues. This intensity difference between parenchyma and blood vessels can effectively help intensity-based registration. However, the diameter of vessel becomes smaller as the blood vessel branches. The small blood vessels are difficult to visualize because of their low intensity contrast. Therefore, grayscale information of the small vessels give almost no contribution to the intensity-based registration. In order to better utilize the information of blood vessel locations, we utilize the vesselness measure (VM) computed from intensity images rather than using their grayscales directly.

Sum of Squared Vesselness Measure Difference (SSVMD)

The vesselness measure is based on the analysis of eigenvalues of the Hessian matrix of image intensity. The eigenvalues, which are geometrically interpreted as principal curvatures, can be used to indicate the shape of underlying object. In 3D lung CT images, isotropic structures such as parenchyma tissues (dark) are associated with three similar nonzero positive eigenvalues. Tubular structures such as blood vessels (bright) are associated with one negligible eigenvalue and two similar nonzero negative eigenvalues [35]. Ordering the eigenvalues of a Hessian matrix by magnitude , the Frangi’s vesselness function [35] is defined as

with , where distinguishes between plate-like and tubular structures, accounts for the deviation from a blob-like structure, and differentiates between tubular structures and noise. , , and control the sensitivity of the vesselness measure. The experiments in this paper use , , and .

The Hessian matrix is computed by convolving the intensity image with second and cross derivatives of the Gaussian function. In a multiscale analysis, the response of the vesselness filter will achieve the maximum at a scale, which approximately matches the size of vessels to detect. Therefore, the vesselness measure is estimated by computing (6) for a range of scales and selecting the maximum response: , where is the standard deviation of the Gaussian function [36].

The vesselness image is rescaled to and can be considered as a probability-like estimate of vesselness features. Larger vesselness value indicates that the underlying object is more likely to be a vessel structure, as shown in Figure 2. Let and represent the vesselness measures of images and at location , respectively. In order to match similar vesselness patterns between two images, the sum of squared vesselness measure difference (SSVMD) is proposed as

Mismatch from vessel to tissue structures will result in a larger SSVMD cost.

(a)

(b)

Figure 2: The vesselness images calculated from lung CT grayscale images. (a) A coronal slice of TLC data and (b) its vesselness measure. Vesselness measure is computed in multiscale analysis with and rescaled to .

2.2.3. Elastic Regularization

Enforcing constraints on the transformation helps generate physiologically more meaningful registration results. Continuum mechanical models such as linear elasticity [11, 37, 38] and viscous fluid [37, 39] can be used to regularize the transformation. A common way to constraint the deformation is applying a differential operator on the transformation and formulating an additive cost term in the objective cost function [11, 22, 40–46]. In our registration algorithms, a linear-elastic constraint was used to regularize the displacement fields , where . This regularization term is formed as

where can be any nonsingular linear differential operator [47]. Here the linear elasticity operator is formed as , where and .

Using the linear elasticity differential operator can help smooth the transformation and help eliminate abrupt changes in the displacement fields. The linear elasticity operator is used in this work to help avoid the transformation from folding onto itself. However, it cannot prevent the Jacobian of the transformation from going negative, that is, destroying the image topology under the transformation [48]. Additional constraints on the displacement parameters are applied in the optimization method.

2.2.4. Total Cost Function

Finally, the total cost is defined as a linear combination of the intensity-based costs, vesselness measure preserving cost, and regularization constraint

can be one of the four intensity-based similarity cost functions: , , , or . Constants and are weights to adjust the significance of the three terms.

2.3. Transformation Parameterization

The transform defines how points from the template image are mapping to their corresponding points in the target image . In three dimensional space, let define a voxel coordinate in the image domain of the target image . The transformation is a vector-valued function defined on the voxel lattice of target image, and gives the corresponding location in template image to the point .

The lung is composed of nonhomogenous soft tissue. Lung tissue expansion varies in different lung regions. Since lung expansion is nonuniform, nonrigid transformations are required to model the lung motion across different inflation levels. To represent the locally varying geometric distortions, the transformation can be represented by various forms of basis function, such as Fourier transform, thin-plate splines, and B-splines. B-splines are well suited for image registration and are able to capture the local nonrigid deformation between two images [40, 45]. Considering the computational efficiency and accuracy requirement, the cubic B-spline based parameterization was chosen to represent the transformation.

Let be the coefficients of the th control point on the spline lattice along each direction. The transformation is represented as

where is a separable convolution kernel. is the uniform cubic B-spline basis function defined as

2.4. Optimization Method

Most registration cost functions can be minimized using standard optimization techniques. There are several existing methods in numerical analysis such as the partial differential equation (PDE) solvers used to solve for elastic and fluid transformations, steepest gradient descent method, conjugate gradient method, and so forth. Similar to [9, 13], our similarity cost functions were optimized using a limited memory quasi-Newton minimization method with bounds (L-BFGS-B) [49] algorithm. It is well suited for optimization with a high-dimensional parameter space. In addition, this algorithm allows bound constraints on the independent variables.

The bound constraints are applied on B-Spline coefficients to guarantee the local injectivity (one-to-one property) of the transformation [50], that is, the transformation maintains the topology of two images. According to their analysis, the displacement fields are locally injective over the domain if the B-Spline coefficients satisfy the conditions that , , and , where , , and are the B-Spline grid sizes along each direction, and is a constant approximately equal to 2.479772335.

2.5. Registration Accuracy Assessment

Validation and evaluation of image registration accuracy is an important task to quantify the performance of registration algorithms. Due to the absence of a “gold standard” to judge a registration algorithm, multiple evaluation methods are needed to evaluate the performance of image registration algorithms.

2.5.1. Landmark Matching Accuracy

Intrasubject CT images of the lung contain identifiable landmarks such as airway-tree and vascular-tree branch points. For each pair of FRC and TLC data, 100–150 distinctive landmark pairs were selected as branch points of the vascular tree using a semiautomatic method [24]. The landmarks were selected so that they were dispersed throughout the lungs. An example of the landmark distribution is shown as green points in Figure 3. The Euclidean distance between the registration-predicted landmark position and its true position is defined as landmark error. Let and be the location of landmark on template image and target image , respectively. The landmark error is calculated as .

(a)

(b)

Figure 3: Distribution of landmarks (green points) selected at vessel-tree branch points on (a) FRC and (b) TLC scans of one subject. Vessel trees are marked as red curves. Different lobes are marked using different colors.

2.5.2. Vessel Matching Accuracy

Vessels in the lung keep their tubular shape and tree structures during the respiratory process. The vascular tree provides rich spatial and shape information in parenchyma regions. Therefore, evaluating the alignment on vessel trees is an important approach to validate the matching accuracy at the lung feature level.

The registration accuracy on the vessel tree was evaluated by vessel matching distance, which is calculated as the distance between a point on the target vessel tree and its closet point on warped template vessel tree. Mathematically, this distance can be stated as the vessel positioning error (VPE)

for a given point in , where and , respectively, are the set of all points in the vessel trees extracted from image and , respectively, and defines the Euclidean distance. The vessel positions used for validation are segmented using vessel segmentation algorithm [51]. Examples of the vessel tree extractions are shown as red curves in Figure 3.

2.5.3. Fissure Alignment Distance

The human lungs are divided into five independent compartments called lobes. Lobar fissures are defined as the division between adjacent lung lobes and represent important physical boundaries within the lungs. Therefore, alignment of fissures is an important evaluation criterion. Fissure locations are extracted from the images by segmenting the lobes using [52] and then identifying voxels adjacent to two lobe segmentations. The fissure positioning error (FPE) is defined as the minimum distance between a point on the deformed fissure and the closest point on the corresponding target fissure

for a given point in , where and , respectively, are the set of all points in the fissure in image and , respectively. Examples of the lobe segmentations are shown in different colors in Figure 3.

2.5.4. Assessment of Lung Function by the Jacobian Determinane

The lung tissue deformation pattern can be used as an index to assess lung function. In this work, the Jacobian determinant of the transformation field derived by image registration is used to estimate the local tissue deformation [53].

The Jacobian determinant (often simply called the Jacobian) [11, 48, 54] is a measurement to estimate the pointwise expansion and contraction during the deformation. The Jacobian of the transformation is defined as

Using a Lagrangian reference frame, a Jacobian value of one corresponds to zero expansion or contraction. Local tissue expansion corresponds to a Jacobian greater than one, and local tissue contraction corresponds to a Jacobian less than one.

2.6. Preprocessing

Preprocessing starts by identifying the lung regions in all images using the method described in [28]. Images and masks are downsampled by a factor of 2 in each dimension to reduce computation time. For each pair of data sets, FRC images are used as the target image, and TLC images are used as the template image.

In order to evaluate how the vesselness cost function affects the registration algorithm results, we performed registration experiments using different similarity costs on parenchyma region in each pair of data sets for comparison. There were four registration methods driven by intensity-only similarity cost functions described in Section 2.2.1, and the same four registration methods that included the feature-based similarity cost SSVMD. After registration, the results from each method were evaluated and compared through matching accuracy on landmarks, vessels, fissures, and underlying transformation properties.

3. Results

3.1. Tuning Weighing Parameters

We designed experiments to discover good parameter settings on intensity-based cost, feature-based cost , and regularization term in (9). Here we discuss our approach for selecting registration parameters for . Parameter settings for registration algorithms using other intensity-based cost functions, for example, , , and can be tuned in the same way.

Table 1 and Figure 4 show the results for 20 CT-to-CT registration experiments, as the weighting values (ved) and (smooth) are varied. The cost functions were values averaged on the results from six subjects. The values of and ranged from 0 to 2 and 0 to 0.5, respectively.

Table 1: Registration experimental parameters and cost values.

Experiment
SSVMD
REG

Jacobian min
Jacobian 1/max

CT01
No
No
0.0
0
34118
57879
172605
0.25
0.11

CT02
No
Yes
0.0
0.05
32751
49291
111770
0.38
0.12

CT03
No
Yes
0.0
0.1
32059
46156
86780
0.44
0.13

CT04
No
Yes
0.0
0.2
33571
44504
63530
0.50
0.13

CT05
No
Yes
0.0
0.4
36736
44853
43827
0.58
0.15

CT06
Yes
No
0.5
0
31288
29492
153068
0.33
0.12

CT07
Yes
Yes
0.5
0.05
30275
28611
110499
0.40
0.12

CT08
Yes
Yes
0.5
0.1
30074
27326
87990
0.44
0.13

CT09
Yes
Yes
0.5
0.2
31194
26812
63612
0.50
0.14

CT10
Yes
Yes
0.5
0.4
33426
28111
46094
0.57
0.14

CT11
Yes
No
1
0
30873
23963
161940
0.36
0.11

CT12
Yes
Yes
1
0.05
29991
22560
115510
0.40
0.11

CT13
Yes
Yes
1
0.1
29923
22328
94784
0.44
0.13

CT14
Yes
Yes
1
0.2
30956
22261
70079
0.49
0.13

CT15
Yes
Yes
1
0.4
32830
23576
50549
0.57
0.15

CT16
Yes
No
1.5
0
31713
20082
181465
0.34
0.11

CT17
Yes
Yes
1.5
0.05
31046
20059
124563
0.39
0.12

CT18
Yes
Yes
1.5
0.1
30871
20270
103340
0.43
0.12

CT19
Yes
Yes
1.5
0.2
31407
20024
76930
0.49
0.13

CT20
Yes
Yes
1.5
0.4
32982
21413
55329
0.56
0.14

CT21
Yes
No
2
0
32532
18907
192035
0.34
0.11

CT22
Yes
Yes
2
0.05
31566
18396
134317
0.39
0.12

CT23
Yes
Yes
2
0.1
31709
18742
109223
0.41
0.12

CT24
Yes
Yes
2
0.2
32172
19218
82619
0.49
0.13

CT25
Yes
Yes
2
0.4
33252
19946
59525
0.56
0.14

(a) SSTVD cost

(b) SSVMD cost

(c) REG cost

Figure 4: Plots of (a) SSTVD, (b) SSVMD, and (c) REG costs using different parameter settings. Data are averaged through results from six subjects.

Experiment CT01 corresponds to unconstrained estimation, in which the transformation was estimated only according to the intensity similarity cost. This experiment produced relatively the worse registration results as evident by the large values of , , and in the respective tables.

Experiments CT02, CT03, CT04, and CT05 demonstrate the effect of estimating the transformations without minimizing the vesselness similarity cost while varying the weight of the linear elastic cost. Minimizing the linear elastic cost is good for optimizing the other two similarity costs and , as we can see from Figures 4(a) and 4(b). values larger than 0.2 are not recommended since they may cause the to increase dramatically. Increasing the constraint weights results in a worse intensity match between images.

Experiments CT06, CT11, CT16, and CT21 demonstrate the effect of using vesselness similarity cost function without enforcing the linear elasticity constraint. The values for these experiments are much lower than the previous cases since they are being minimized. The intensity similarity costs also decreased using registration with vesselness constraint, especially when is in the range of .

The remaining experiments show the effect of jointly estimating the transformations, while varying the weights on both the vesselness similarity cost and the linear elasticity constraint. These experiments show that it is possible to find a set of parameters that produce better results using both constraints than only using one or none.

The optimal set of parameters should be chosen to provide a good intensity match and vesselness match, while producing less spatial distortion as measured by an acceptable level of linear elastic cost. From the experiments, we observe that and are good for minimizing the three costs at the same time. In this work, weighting parameters were set to and when using SSTVD as intensity cost function. Parameters in registration using other three intensity cost functions are set in the same way.

3.2. Multiresolution Optimization Scheme

A spatial multiresolution optimization procedure from coarse to fine was used to improve speed, accuracy, and robustness of the registration. In the experiments, the spatial multiresolution strategy proceeds from low to high resolution, starts at one-eighth of the spatial resolution, and increases by a factor of two until the full resolution is reached. Meanwhile, a hierarchy of B-Spline grid spacings from large to small was also used. The finest B-spline grid space used in the experiments was 4 mm. An example multiresolution scheme design for minimizing the total cost function is listed in Table 2. The image spatial resolution and B-spline grid spacing were refined alternatively.

Table 2: An example multiresolution scheme.

Image resolution
B-spline grid size
Max. iteration

1/8
128 mm
200

64 mm
200

1/4
32 mm
200

16 mm
200

1/2
16 mm
50

8 mm
50

1
8 mm
20

4 mm
20

Figure 5 lists the cost values at each iteration of one registration. At 1/8 and 1/4 image resolutions, registration speed is fast and runs for more iterations. The optimization is stopped before reaching the maximum iterations if the total cost change is nominal between consecutive iterations. Global shapes are matched at these two levels. At 1/2 image resolution, the inner structures get clearer and are aligned roughly. Registration at full resolution further adjusts the local structures matching. During the registration procedure, the Jacobian values are checked to make sure that they remain positive.

(a)

(b)

(c)

(d)

Figure 5: Cost functions versus iteration associated with experiment CT13. Cost function values are recorded at (a) 1/8, (b) 1/4, (c) 1/2, and (d) full image resolution.

3.3. Landmark Matching Accuracy

The transformation determined from registration can be used to track the landmark movements. The original average landmark error is  mm with a maximum landmark error of 72.79 mm. Table 3 shows the mean and standard deviation of landmark errors through all six subjects after using different registration methods. Figure 6 shows the box plot of landmark errors.

Table 3: Landmark errors (mm) using six subjects.

Without SSVMD
With SSVMD

Avg.
Max
Avg.
Max

SSD

15.97

4.56

MI

16.72

5.11

NCC

15.18

6.99

SSTVD

15.38

4.1

Figure 6: Box plot of landmark errors for six subjects after using eight registration methods. Results from methods with SSD, MI, NCC, and SSTVD metric along contain outliers beyond the error range.

3.4. Vessel Matching Accuracy

The original average vessel positioning error was  mm. Table 4 shows the vessel positioning errors for six subjects after using different registration methods. The average errors and standard deviations were all decreased after adding the vesselness constraint.

Table 4: Vessel positioning errors (mm) averaged over six subjects.

SSD
MI
NCC
SSTVD

Without SSVMD

With SSVMD

Figure 7 shows the distance map on FRC vessel tree from one subject after using eight different registration methods. These results show that large errors between the deformed source and target vessel trees are reduced after adding the SSVMD constraint.

(a) SSD

(b) SSD + SSVMD

(c) MI

(d) MI + SSVMD

(e) NCC

(f) NCC + SSVMD

(g) SSTVD

(h) SSTVD + SSVMD

Figure 7: Vessel positioning errors (mm) on target vessel tree. Results are generated from eight registration methods. Arrows denote regions of large discrepancies between the deformed source and target vessel trees. Note that the errors in these regions were reduced after adding the SSVMD constraint to the registration algorithms.

3.5. Fissure Alignment Distance

For each pair of FRC and TLC images, the parenchyma regions were segmented into five lobes, and the three fissures were identified, where the lobe segmentations touched each other. The average fissure positioning error was  mm before registration. The mean and standard deviation of fissure positioning error over all three fissures after using eight registration methods are shown in Table 5. The average fissure positioning error across the six subjects was 9.20 mm and decreased to around 1.50 mm and 1.00 mm without and with SSVMD cost, respectively. The fissure positioning matching accuracy improved by approximately 30% after adding SSVMD cost.

Table 5: Fissure positioning errors (mm) using six subjects.

SSD
MI
NCC
SSTVD

Without SSVMD

With SSVMD

The distance map on fissure planes from one subject after using eight different registration methods is shown in Figure 8. Notice that adding SSVMD helped improve registration accuracy in the lung regions near the thoracic cage.

(a) SSD

(b) SSD + SSVMD

(c) MI

(d) MI + SSVMD

(e) NCC

(f) NCC + SSVMD

(g) SSTVD

(h) SSTVD + SSVMD

Figure 8: Fissure positioning error (mm) on fissure planes. Arrows denote regions of large discrepancies between the deformed source and target fissure planes.

3.6. Assessment of Lung Function by Jacobian

Registration methods producing similar matching accuracy on the feature locations may result in different underlying parenchymal tissue functions. In order to reveal the lung tissue deformation pattern, the Jacobian of the transformation field is used to estimate the local tissue deformation [4]. The Jacobian determinant at a given point gives important information about the behavior of transformation near that point. Figure 9 shows the Jacobian maps resulted from eight registration methods. Arrows denote regions that show different deformation patterns using intensity-only registration methods, but they are more similar after adding vesselness cost function.

(a) SSD

(b) SSD + SSVMD

(c) MI

(d) MI + SSVMD

(e) NCC

(f) NCC + SSVMD

(g) SSTVD

(h) SSTVD + SSVMD

Figure 9: The color-coded Jacobian maps of a sagittal slice resulted from eight registration methods. Blue and purple regions have larger lung deformation, while red and orange regions are deforming less. Arrows denote regions which show different deformation patterns using intensity-only registration methods.

4. Discussion

The experiments presented in this paper were designed to evaluate the performance of the vesselness constraint when it was added to intensity-based registration algorithms. The experimental results of tuning weighting parameters in the total cost function suggest that using both vesselness and smoothing constraint can help minimize the similarity cost, as shown in Figure 4. Figure 5 shows that the multiresolution scheme is important and useful for solving complex registration problem efficiently. This is because the registration is first performed at a coarse scale, where the images have much fewer voxels, which are fast and can help eliminate local minima.

Table 3 and Figure 6 demonstrate that adding the SSVMD cost function reduced the mean landmark errors of the four basic registration methods. Landmarks with large errors, shown as outliers in the box plot, are aligned much better when SSVMD is used. Figure 7 reflects the fact that the SSVMD constraint helps improve matching accuracy over all four basic methods on small vessels, around lung boundaries and in the region near diaphragm. The reason for this is that blood vessels in those regions are usually small and have low intensity contrast, and thus they contribute little to conventional intensity similarity criteria. The vesselness measurement enhances blood vessel information and strengthens contribution of small vessels to registration process when using the SSVMD similarity cost. Figure 8 indicates that the SSVMD not only helps match vessel structures, but also helps improve registration accuracy in other regions, such as positions on the fissure planes near the thoracic cage.

Good matching accuracy on the feature locations does not guarantee that the parenchymal tissue is correctly aligned. Rather than evaluates the alignment accuracy, the Jacobian evaluates the quality of the transformation properties. It reveals how well the transformation preserves topology and measures the differential lung volume change. In Figure 9, the left column shows that the Jacobian maps generated by the four registration methods without SSVMD have a similar ventral to dorsal gradient as expected since the subjects were imaged in the supine orientation. However, the local tissue deformation patterns derived from these methods are different even in the methods pair SSD and SSTVD, which have similar landmark errors as shown in Table 3. This is consistent with the findings that while the intermethod variability on the landmark error is small, there may be discriminating difference in the Jacobian maps [55]. The Jacobian map from SSTVD method shows more local structure in the dorsal region. The right column shows that adding the SSVMD constraint produces Jacobian images that are much more similar across different registration methods and reveal more detailed deformation patterns especially near vessel locations. Generally, vessels have smaller volume changes comparing with parenchymal tissues during breathing cycles. The four Jacobian maps produced using registration methods with SSVMD are similar, which may imply that the derived local deformation patterns are more reliable.

Comparing the four intensity-only registration methods, registration driven by SSTVD achieved lower landmark error and lower vessel and fissure positioning errors than other three methods driven by SSD, MI, and NCC. Figures 7 and 8 reflect that SSTVD-driven registration resulted in more accurate matching within the region near the thoracic cage. The reason may be that SSTVD cost function contains a local Jacobian factor, which can constrain incorrect large displacement and help prevent distortion near the thoracic cage. After adding SSVMD on the four intensity-based cost functions, these methods generated similar transformations. SSVMD helps improve matching accuracy in regions near the thoracic cage and near the diaphragm. In our experiments, registration method using SSTVD + SSVMD resulted in the smallest matching errors. Therefore, we may consider that SSTVD + SSVMD is the best similarity cost function to register lung CT images according to our evaluation.

The effectiveness of the vesselness preserving metric was tested on a variety of lung CT data sets as part of the grand challenge “Evaluation of Methods for Pulmonary Image Registration 2010” (EMPIRE10) [56, 57]. More than 20 individual registration algorithms from different groups were applied to 30 pairs of lung CT scans in the EMPIRE10 challenge. Besides our tissue volume and vesselness preserving method, one other method called the Robust TreeReg [58] also combined intensity and feature information. The Robust TreeReg algorithm performed a robust tree registration (RTR) and added correspondences between bifurcations of the vessel tree to the voxel-based mutual information driven registration. For this registration challenge, our tissue volume and vesselness preserving method had better performance than the Robust TreeReg method. This result may imply that the vesselness measure provides more feature information than the bifurcation landmarks of the vessel tree. In addition, our tissue volume (or mass) and vesselness preserving method was shown to improve matching results compared to methods that only incorporated mass preservation [59].

5. Conclusions

This paper presented nonrigid registration algorithms driven by commonly used intensity-based criteria for lung registration, a feature-based vesselness constraint, and a linear elastic smoothing constraint. Results were presented to show that adding the SSVMD constraint to existing similarity cost functions such as SSD, MI, NCC, and SSTVD reduces landmark error and improves overlap on vascular tree and fissure planes. The purpose of adding the vesselness cost in registration process is that it can help correct the mismatches of small vessels and their surrounding lung tissues. Using the SSVMD constraint was shown to produce a more detailed expansion pattern for local tissue, especially near vessel locations. Also, the expansion patterns were similar across different registration methods. This demonstrates that using the SSVMD constraint not only helps match on feature structures, but also helps align corresponding parenchymal tissues providing a more reliable pattern of local lung tissue deformation. In this paper, registration method preserving both tissue volume and vesselness measurement performed best on matching 3D lung CT data according to our evaluation.

Acknowledgments

The authors would like to thank Dr. K. Murphy and Dr. B. van Ginneken for providing the software iX for generating and annotating landmarks. This work was supported by Grant nos. HL079406, HL64368, HL080285, EB004126, and CA129022.

References

B. Li, G. E. Christensen, E. A. Hoffman, G. McLennan, and J. M. Reinhardt, “Establishing a normative atlas of the human lung: intersubject warping and registration of volumetric CT images,” Academic Radiology, vol. 10, no. 3, pp. 255–265, 2003.

T. Guerrero, K. Sanders, E. Castillo et al., “Dynamic ventilation imaging from four-dimensional computed tomography,” Physics in Medicine and Biology, vol. 51, no. 4, pp. 777–791, 2006.

G. E. Christensen, J. H. Song, W. Lu, I. E. Naqa, and D. A. Low, “Tracking lung tissue motion and expansion/compression with inverse consistent image registration and spirometry,” Medical Physics, vol. 34, no. 6, pp. 2155–2163, 2007.

J. M. Reinhardt, K. Ding, K. Cao, G. E. Christensen, E. A. Hoffman, and S. V. Bodas, “Registration-based estimates of local lung tissue expansion compared to xenon CT measures of specific ventilation,” Medical Image Analysis, vol. 12, no. 6, pp. 752–763, 2008.

T. A. Sundaram and J. C. Gee, “Towards a model of lung biomechanics: pulmonary kinematics via registration of serial lung images,” Medical Image Analysis, vol. 9, no. 6, pp. 524–537, 2005.

K. Ding, Y. Yin, K. Cao et al., “Evaluation of lobar biomechanics during respiration using image registration,” in Proceedings of the 12th International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI '09), vol. 5761, pp. 739–746, 2009.

Y. Yin, E. A. Hoffman, and C. L. Lin, “Lung lobar slippage assessed with the aid of image registration,” in Proceedings of the 13th International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI '10), vol. 13, pp. 578–585, Springer, 2010.

K. Ding, J. E. Bayouth, J. M. Buatti, G. E. Christensen, and J. M. Reinhardt, “4DCT-based measurement of changes in pulmonary function following a course of radiation therapy,” Medical Physics, vol. 37, no. 3, pp. 1261–1272, 2010.

D. Mattes, D. R. Haynor, H. Vesselle, T. K. Lewellen, and W. Eubank, “PET-CT image registration in the chest using free-form deformations,” IEEE Transactions on Medical Imaging, vol. 22, no. 1, pp. 120–128, 2003.

D. Sarrut, V. Boldea, S. Miguet, and C. Ginestet, “Simulation of four-dimensional CT images from deformable registration between inhale and exhale breath-hold CT scans,” Medical Physics, vol. 33, no. 3, pp. 605–617, 2006.

G. E. Christensen and H. J. Johnson, “Consistent image registration,” IEEE Transactions on Medical Imaging, vol. 20, no. 7, pp. 568–582, 2001.

V. Gorbunova, P. Lo, H. Ashraf, A. Dirksen, M. Nielsen, and M. de Bruijne, “Weight preserving image registration for monitoring disease progression in lung CT,” in Proceedings of the 11th International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI '08), vol. 5242, pp. 863–870, 2008.

Y. Yin, E. A. Hoffman, and C. L. Lin, “Mass preserving nonrigid registration of CT lung images using cubic B-spline,” Medical Physics, vol. 36, no. 9, pp. 4213–4222, 2009.

L. Fan, Chang Wen Chen, J. M. Reinhardt, and E. A. Hoffman, “Evaluation and application of 3D lung warping and registration model using HRCT images,” in Medical Imaging 2001: Physiology and Function from Multidimensional Images, Proceedings of SPIE, pp. 17–22, February 2001.

M. R. Kaus, K. K. Brock, V. Pekar, L. A. Dawson, A. M. Nichol, and D. A. Jaffray, “Assessment of a model-based deformable image registration approach for radiation therapy planning,” International Journal of Radiation Oncology Biology Physics, vol. 68, no. 2, pp. 572–580, 2007.

M. M. Coselmon, J. M. Balter, D. L. McShan, and M. L. Kessler, “Mutual information based CT registration of the lung at exhale and inhale breathing states using thin-plate splines,” Medical Physics, vol. 31, no. 11, pp. 2942–2948, 2004.

S. Krishnan, K. C. Beck, J. M. Reinhardt et al., “Regional lung ventilaton from volumetric CT scans using image warping functions,” in Proceedings of the 2nd IEEE International Symposium on Biomedical Imaging: Macro to Nano, pp. 792–795, Washington, DC, USA, April 2004.

V. Gorbunova, S. Durrleman, P. Lo, X. Pennec, and M. de Bruijne, “Curve- and surface-based registration of lung ct images via currents,” in Proceedings of the 2nd International Workshop on Pulmonary Image Analysis, pp. 15–25, 2009.

Y. Yin, E. A. Hoffman, K. Ding, J. M. Reinhardt, and C. L. Lin, “A cubic B-spline-based hybrid registration of lung CT images for a dynamic airway geometric model with large deformation,” Physics in Medicine and Biology, vol. 56, no. 1, pp. 203–218, 2011.

K. Cao, G. E. Christensen, K. Ding, and J. M. Reinhardt, “Intensity-and-landmarkdriven, inverse consistent, B-Spline registration and analysis for lung imagery,” in Proceedings of the 2nd International Workshop on Pulmonary Image Analysis, pp. 137–148, 2009.

V. Gorbunova, S. Durrleman, P. Lo, X. Pennec, and M. De Bruijne, “Lung CT registration combining intensity, curves and surfaces,” in Proceedings of the 7th IEEE International Symposium on Biomedical Imaging (ISBI '10), pp. 340–343, April 2010.

H. J. Johnson and G. E. Christensen, “Consistent landmark and intensity-based image registration,” IEEE Transactions on Medical Imaging, vol. 21, no. 5, pp. 450–461, 2002.

C. V. Stewart, Y. L. Lee, and C. L. Tsai, “An uncertainty-driven hybrid of intensity-based and feature-based registration with application to retinal and lung CT images,” in Proceedings of the 7th International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI '04), pp. 870–877, September 2004.

K. Murphy, B. van Ginneken, J. P. W. Pluim, S. Klein, and M. Staring, “Semi-automatic reference standard construction for quantitative evaluatoin of lung CT registration,” in Proceedings of the 11th International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI '08), vol. 5242, pp. 1006–1013, 2008.

K. Cao, K. Ding, G. E. Christensen, and J. M. Reinhardt, “Tissue volume and vesselness measure preserving nonrigid registration of lung CT images,” in Medical Imaging 2010: Image Processing, Proceedings of SPIE, usa, February 2010.

K. Cao, K. Du, K. Ding, J. M. Reinhardt, and G. E. Christensen, “Regularized nonrigid registration of lung CT images by preserving tissue volume and vesselness measure,” in Grand Challenges in Medical Image Analysis, 2010.

K. Cao, K. Ding, G. E. Christensen, M. L. Raghavan, R. E. Amelon, and J. M. Reinhardt, “Unifying vascular information in intensity-based nonrigid lung ct registration,” in Proceedings of the 4th International Workshop on Biomedical Image Registration, vol. 6204 of Lecture Notes in Computer Science, pp. 1–12, Springer, 2010.

S. Hu, E. A. Hoffman, and J. M. Reinhardt, “Automatic lung segmentation for accurate quantitation of volumetric X-ray CT images,” IEEE Transactions on Medical Imaging, vol. 20, no. 6, pp. 490–498, 2001.

G. P. Penney, J. Weese, J. A. Little, P. Desmedt, D. L. G. Hill, and D. J. Hawkes, “A comparison of similarity measures for use in 2-D-3-D medical image registration,” IEEE Transactions on Medical Imaging, vol. 17, no. 4, pp. 586–595, 1998.

M. Holden, D. L. G. Hill, E. R. E. Denton et al., “Voxel similarity measures for 3-D serial MR brain image registration,” IEEE Transactions on Medical Imaging, vol. 19, no. 2, pp. 94–102, 2000.

Wikipedia, “Histogram matching,” Hypertext Document, 2012

P. Thévenaz and M. Unser, “Spline pyramids for inter-modal image registration using mutual information,” in Wavelet Applications in Signal and Image Processing V, vol. 3169 of Proceedings of SPIE, pp. 236–247, July 1997.

T. S. Yoo, Insight Into Images: Principles and Practice for Segmentation, Registration, and Image Analysis, AK Peters, 2004.

E. A. Hoffman and E. L. Ritman, “Effect of body orientation on regional lung expansion in dog and sloth,” Journal of Applied Physiology, vol. 59, no. 2, pp. 481–491, 1985.

A. F. Frangi, W. J. Niessen, K. L. Vincken, and M. A. Viergever, “Multiscale vessel enhancement filtering,” in Proceedings of the Medical Image Computing and Computer-Assisted Intervention Conference (MICCAI '98), vol. 1496, pp. 130–137, 1998.

A. Enquobahrie, L. Ibanez, E. Bullitt, and S. Aylward, “Vessel enhancing diffusion filter,” The Insight Journal—MICCAI Open Science Workshop, 2007.

G. E. Christensen, S. C. Joshi, and M. I. Miller, “Volumetric transformation of brain anatomy,” IEEE Transactions on Medical Imaging, vol. 16, no. 6, pp. 864–877, 1997.

G. E. Christensen, R. D. Rabbitti, and M. I. Miller, “3D brain mapping using a deformable neuroanatomy,” Physics in Medicine and Biology, vol. 39, no. 3, pp. 609–618, 1994.

G. E. Christensen, R. D. Rabbitt, and M. I. Miller, “Deformable templates using large deformation kinematics,” IEEE Transactions on Image Processing, vol. 5, no. 10, pp. 1435–1447, 1996.

J. Kybic, Elastic image registration using parametric deformation models [Ph.D. thesis], Swiss Federal Institute of Technology, Lausanne, Switzerland, 2001.

D. Shen and C. Davatzikos, “HAMMER: hierarchical attribute matching mechanism for elastic registration,” IEEE Transactions on Medical Imaging, vol. 21, no. 11, pp. 1421–1439, 2002.

X. Geng, D. Kumar, and G. E. Christensen, “Transitive inverse-consistent manifold registration,” in Proceedings of the International Conference on Information Processing in Medical Imaging, vol. 3564 of Lecture Notes in Computer Science, pp. 468–479, Springer, 2005.

M. Vaillant and J. Glaunès, “Surface matching via currents,” in Information Processing in Medical Imaging, pp. 381–392, Springer, 2005.

J. Kybic and M. Unser, “Fast parametric elastic image registration,” IEEE Transactions on Image Processing, vol. 12, no. 11, pp. 1427–1442, 2003.

D. Rueckert, “Nonrigid registration using free-form deformations: application to breast mr images,” IEEE Transactions on Medical Imaging, vol. 18, no. 8, pp. 712–721, 1999.

G. E. Christensen, “Inverse consistent registration with object boundary constraints,” in Proceedings of the 2nd IEEE International Symposium on Biomedical Imaging: Macro to Nano, pp. 591–594, IEEE, Arlington, Va, USA, April 2004.

M. Miller, A. Banerjee, G. Christensen et al., “Statistical methods in computational anatomy,” Statistical Methods in Medical Research, vol. 6, no. 3, pp. 267–299, 1997.

G. E. Christensen, R. D. Rabbitt, M. I. Miller et al., “Topological properties of smooth anatomic maps,” in Information Processing in Medical Imaging, Y. Bizais, C. Braillot, and R. D. Paola, Eds., vol. 3, pp. 101–112, 1995.

R. H. Byrd, P. Lu, J. Nocedal, and C. Zhu, “A limited memory algorithm for bound constrained optimization,” SIAM Journal on Scientific Computing, vol. 16, no. 5, pp. 1190–1208, 1995.

Y. Choi and S. Lee, “Injectivity conditions of 2D and 3D uniform cubic B-spline functions,” Graphical Models, vol. 62, no. 6, pp. 411–427, 2000.

H. Shikata, E. A. Hoffman, and M. Sonka, “Automated segmentation of pulmonary vascular tree from 3D CT images,” in Medical Imaging 2004: Physiology, Function, and Structure from Medical Images, A. A. Amini and A. Manduca, Eds., vol. 5369 of Proceedings of SPIE, pp. 107–116, 2004.

S. Ukil and J. M. Reinhardt, “Anatomy-guided lung lobe segmentation in X-ray CT images,” IEEE Transactions on Medical Imaging, vol. 28, no. 2, pp. 202–214, 2009.

J. M. Reinhardt, K. Ding, K. Cao, G. E. Christensen, E. A. Hoffman, and S. V. Bodas, “Registration-based estimates of local lung tissue expansion compared to xenon CT measures of specific ventilation,” Medical Image Analysis, vol. 12, no. 6, pp. 752–763, 2008.

R. C. Buck, Advanced Calculus, McGraw-Hill Book Company, St. Louis, Miss, USA, 3rd edition, 1978.

S. Kabus, T. Klinder, K. Murphy, B. van Ginneken, C. Lorenz, and J. P. W. Pluim, “Evaluation of 4D-CT lung registration,” in Proceedings of the 12th International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI '09), pp. 747–754, 2009.

K. Murphy and B. van Ginneken, “Evaluation of methods for pulmonary image registration,” Hypertext Document, 2010

K. Murphy, B. van Ginneken, and Reinhardt, “Evaluation of registration methods on thoracic CT: The EMPIRE10 challenge,” IEEE Transactions on Medical Imaging, vol. 30, no. 11, pp. 1901–1920, 2011.

D. Loeckx, D. Smeets, J. Keustermans et al., “3D lung registration using splineMIRIT and robust tree registration (RTR),” in Grand Challenges in Medical Image Analysis, 2010.

V. Gorbunova, J. Sporring, P. Lo, A. Dirksen, and M. de Bruijne, “Mass preserving image registration: results of evaluation of methods for pulmonary image registration 2010 challenge,” in Grand Challenges in Medical Image Analysis, 2010.