Quantification of regional deformation of the lungs by non-rigid registration of three-dimensional contrast-enhanced magnetic resonance imaging

Introduction

Assessment of lung function is vital for the diagnosis of a variety of pathological conditions, such as emphysema and pulmonary embolism. The pulmonary function test (PFT) has been a standard method for evaluation of lung function; however, PFT is a global measurement and does not provide an assessment of regional lung function, which is desirable for evaluating heterogeneous conditions such as chronic obstructive pulmonary disease (COPD).

Many methods have been proposed for evaluation of regional lung function using different imaging modalities, including single-photon emission computed tomography (SPECT) (1), positron emission tomography (PET) (2), computed tomography (CT) (3-7), and magnetic resonance imaging (MRI) (8-13). Both SPECT and PET are based on the application of radioactively labeled tracers, directly provide three-dimensional (3D) lung ventilation images, and are widely accepted as standard methods for measurement of regional pulmonary function, but are subject to relatively poor spatial resolution. For CT-based methods, aligning two volumes (14-16) or multivolume CT (17,18) using non-rigid image registration offers a way of obtaining 3D pulmonary function at high-resolution. After alignment of two or more volumes, the regional lung ventilation can be measured using information of regional intensity changes at different volumes or spatial Jacobian of the transformation obtained by image registration, which has been demonstrated to have significant correlations with PFT, PET measurements as well as SPECT measurement (7). The cost of CT-based methods is the exposure to ionizing radiation, which limits their values for frequent follow-up examinations in children and young adults.

Recent advances in MRI permit lung ventilation assessment using hyperpolarized noble gasses, such as ³He (19), as contrast agents. Research has been proposed to study pulmonary mechanics and kinematics using two-dimensional (2D) MRI (20). This allows estimation of regional lung tissue mechanics but is limited to 2D information. We propose a method for quantification of 3D regional deformation of the lung parenchyma by registering two volumetric contrast-enhanced MR angiography (CE-MRA) images. We hypothesize that the greatly enhanced conspicuity of pulmonary blood vessels due to gadolinium contrast injection will facilitate non-rigid image registration of two CE-MRA images acquired at end-inspiration and end-expiration, respectively. The goal is to develop a technique to evaluate lung regional deformation via registration of 3D CE-MRA images which avoids ionizing radiation while achieving acceptable accuracy and can provide accurate 3D lung regional deformation information.

Methods

Image acquisition

The study was approved by the Institutional Review Board and was compliant with the Health Insurance Portability and Accountability Act. Three pairs of 3D CE-MRA data sets, acquired from three healthy volunteers, were used in this study. For each volunteer, two breath-held CE-MRA scans were acquired 5 min after gadolinium contrast injection: one at end-expiration and one at end-inspiration. The images were acquired using a 3D spoiled gradient echo sequence with the same flip angle =30 degrees and parallel imaging (2X GRAPPA) on a 1.5 T scanner (Avanto, Siemens Healthcare, Erlangen, Germany). Detailed parameters of the three pairs of datasets were summarized in Table 1. They were modified to make the breath-holds time less than 25 seconds for all volunteers. For volunteer 1 and volunteer 3, the end-inspiration scans were acquired immediately after the end-expiration scans, while for volunteer 2, the end-expiration scan was acquired immediately after the end-inspiration scan.

Table 1 Parameters for CE-MRA scans of each volunteer
Full table

Image registration

Registration algorithm

Non-rigid image registration is widely used in medical image applications, including registration of lung CT images with application to measuring regional lung function. Various non-rigid registration algorithms have been proposed and evaluated (21). The Elastix tool (22) for the non-rigid image registration was used in this work, and the applied formulation of the image registration algorithm was adopted as the following parameterized optimization problem:

where I_r and I_s denote the reference and source image respectively, T_u is the spatial transformation relating the two and vector u is its parameters, S() is the similarity measure that defines the alignment quality, P is a penalty term that regularizes the transformation T_u, and k is a parameter that weighs the regularity against similarity.

For registration of each dataset, the end-expiration image was always set as the reference image, and the end-inspiration image was set as the source image. The main components of image registration applied in this work included lung mask selection, spatial transformations, similarity measure, penalty term, optimization algorithm and hierarchical strategy, which are described as follows.

Lung mask selection

A rectangle box was manually selected as a mask in the reference image (end-expiration image) for each dataset, and no mask was selected in the source image. The registration algorithm was only applied in the mask area of the reference image. The mask was required to cover the whole lungs in both right-to-left (RL) and inferior-to-superior (IS) dimensions. A limited area outside of the lungs but inside the mask in these two directions are acceptable (Figure 1A), where the yellow rectangle is drawn on the maximum intensity projections (MIPs) of the end-expiration scan represents the selected mask edges in RL and IS directions.

Figure 1 Maximum intensity projections (MIPs) of lung contrast-enhanced MR angiography (CE-MRA) (zoomed, same windowing level) acquired from volunteers 2 in anterior-to-posterior (AP) and inferior-to-superior (IS) directions. With contrast enhancement, numerous blood vessels were observed in the images. The selected landmarks were labeled using numbers at their projection positions, and the same number represents the same branching point of blood vessels.

In the posterior area of the lungs, there are denser and brighter blood vessel signals in end-expiration scans than that in end-inspiration scans, which may cause larger misalignment of blood vessels in other areas of the lungs if this area is included in the lung mask. Therefore, in the anterior-to-posterior (AP) direction, a few coronal slices that covered the last posterior area of the lungs were excluded when manually selecting the mask (Figure 1B), where the yellow rectangle is drawn on the MIPs of the end-expiration scan represents the selected mask edges in RL and AP directions.

Spatial transformations and hierarchical strategy

Cubic B-splines were used as non-rigid transformations for registration due to its computational efficiency and performance in capturing non-rigid deformation (23,24). An important parameter that influences registration accuracy in this transformation is the grid spacing of the B-splines control points. A multi-gird strategy (22) was applied with B-splines transformation, which requires the setting of the B-spine control point spacing to accommodate different strategies. The grid spacing was set to 120, 60, 36, 18, 12, and 6 mm in each direction for each resolution, respectively. The larger grid spacing is used to initially capture the coarse deformations, and then progressively smaller grid spacing was used to gradually capture smaller structures.

The registration process was initiated using smoothed and down-sampled images to increase the robustness of successful registration. To ensure compatibility with the multi-grid approach applied in transformation, six resolutions were used for the registration. The images were down-sampled (22) with a factor of 16, 8, 4, 2, 1 and 1 in RL and IS directions and a factor of 8, 4, 2, 1, 1 and 1 in AP direction, for each resolution. Table 2 lists the corresponding parameters.

Table 2 Grid spacing and down-sample factors for resolution levels R1–R6
Full table

Similarity measure and penalty term

Normalized cross correlation (NCC) (25) was selected as the similarity measure, which can compensate for global intensity differences between end-expiration and end-inspiration images for mono-modal image registration. The bending energy term was applied to restrict sharp deviations of the transformation to avoid undesired deformation (24). The weight for the bending energy term was set to 0.05 in this work.

Optimization algorithm

Iterative procedures were applied to optimize Eq. [1]. The adaptive stochastic gradient descent (ASGD) (26) algorithm was used, and the iteration number was set to 2,000 for each resolution.

Data analysis

To validate and quantitatively analyze the accuracy of image registration results, branching points of pulmonary blood vessels were selected manually as landmarks for validation and analysis. Pairs of landmarks were chosen in the lung area for each pair of the 3D images by an experienced observer. Two other observers were shown the chosen landmarks on the end-expiration scans and asked to locate the corresponding landmarks on the end-inspiration scans. Each landmark position in end-inspiration scans was calculated by averaging the three positions corresponding to the same landmark in the end-expiration scan selected by the three observers. These sets of averaged landmarks were treated as “ground truth” sets of landmarks for end-inspiration scans.

The landmark distance before registration, which is the displacement of a pulmonary blood vessel branching point from end-expiration to end-inspiration, was calculated using positions of the selected end-expiration landmarks and the corresponding averaged end-inspiration landmarks. The landmark distance after registration was calculated using the landmarks positions in the source image (end-inspiration scan) and the predicted landmarks positions by applying the calculated transformation to the landmarks in the reference image (end-expiration scan).

Results

Image acquisition

The MIPs of end-expiration scans and end-inspiration scans acquired from volunteer 2 are shown in Figure 1 as an example. Figure 1A and Figure 1C show the AP direction MIPs of the end-expiration and end-inspiration scans, respectively. Figure 1B and Figure 1D show the IS direction MIPs of the end-expiration and end-inspiration scans, respectively. With gadolinium contrast enhancement, numerous blood vessels signals were enhanced in the images, providing significant contrast information for non-rigid image registration in the lung area.

Validation

Landmarks validation

The selected landmarks for end-expiration scan of volunteer 2 were displayed using numbers as labels in Figure 1A,B, and the corresponding “ground truth” sets of landmarks for end-inspiration images were labeled using the same numbers shown in Figure 1C,D. The same number represents the same branching point of blood vessels.

The average of landmarks distances between end-expiration and end-inspiration scans before registration are listed in Table 3, which shows that there are primarily large distances in IS and AP directions. Average landmarks distances after registration are summarized in Table 4. Comparing Tables 3,4, it is clear that the original landmarks distances are greatly reduced by non-rigid registration, and the average landmarks distances in three pairs of datasets are reduced from 17.9, 20.3 and 16.3 mm, to 1.0, 1.6 and 1.2 mm, respectively. After registration, the average distances error of each pair of datasets were less than 0.6 mm in the RL direction, less than 0.9 mm in the IS direction, and less than 1.2 mm in the AP direction.

Table 3 Average landmarks distances before registration at different directions
Full table

Table 4 Average landmarks distances after registration at different directions
Full table

Visual inspection

Figure 2 shows the alignment of blood vessels in the lung area before and after registration using MIPs from datasets A (volunteer 1). The red channel in each color image (Figure 2) represents the reference image (the end-expiration scan); the green channel in (Figure 2A-C) represents the source image (the end-inspiration scan); the green channel in (Figure 2D-F) represents the warped image obtained by non-rigid registration; and the value of blue channels are all set to be zero. Figure 2A,D shows the composited MIP images in AP direction. Figure 2B,E shows the composited MIP images of the right lung in RL direction, and Figure 2C,F depict the composited MIP CE-MRA images of the left lung in RL direction before and after registration. All the images in Figure 2 are shown in the same window level. With these settings, the well-aligned structures in the lung area show yellow color, and the places with separate red and green colors represent misaligned structures. The alignment of blood vessel structures in each direction in the lung area was greatly improved after registration.

Figure 2 Composite maximum intensity projections (MIPs) of lung contrast-enhanced MR angiography (CE-MRA) images (zoomed, same windowing level) of the first volunteer before registration and after registration, where the red channel in each RGB image represents the reference image, the green channel in image (A-C) represent the source image, and the green channel in image (D-F) represent the warped image obtained by non-rigid image registration. (A,D) the composite MIP images in AP direction; (B,F) the composite MIP images of the right lung in RL direction; (C,F) the composite MIP images of the left lung in right-to-left (RL) direction. The area with yellow color represents good alignment in this area.

Figure 3 compares the reference image of datasets A and the corresponding warped image obtained from registration in two slices (coronal and sagittal). A coronal slice of the reference image is shown in Figure 3A, the corresponding slice of the warped image is shown in Figure 3B, and the slice of the source image at the same coronal position is shown in Figure 3C. A sagittal slice of the reference image is shown in Figure 3D, the corresponding slice of the warped image is shown in Figure 3E, and the slice of the source image at the same coronal position is also shown in Figure 3F. The excellent conspicuity of the pulmonary blood vessels, even at distal branches, helped the non-rigid registration algorithm achieve accurate registration results.

Figure 3 A comparison between the reference image and the warped image of the first volunteer. (A) A coronal slice from the reference image; (B) a sagittal slice of the right lung from the reference image; (C) a coronal slice from the warped image obtained by non-rigid registration; (D) a sagittal slice of the right lung from the warped image; (E) a coronal slice from the source image; and (F) a sagittal slice of the right lung from the source image.

Lung regional deformation estimation

Figure 4 shows the sampled deformation vector fields obtained by non-rigid registration of end-expiration and end-inspiration scan images. The sampled deformation vector fields are overlaid on one coronal slice and one axial slice of each end-expiration scan for illustration using green vectors. Green vectors indicate corresponding lung regional movement from end-expiration to end-inspiration breath-hold.

Figure 4 Sampled deformation vector fields overlaid on two slices (coronal and axial) of the end-expirations contrast-enhanced MR angiography (CE-MRA) scans of three volunteers, where the green vectors indicate the lung regional movement due to inspiration. (A-C) one coronal slice of each volunteer, respectively; (D-F) one axial slice of each volunteer, respectively.

Discussion

We investigated the potential quantification of lung regional deformation by non-rigid registration of two 3D CE-MRA images acquired at end-expiration and end-inspiration. The 3D CE-MRA lung imaging experiment demonstrated that, with gadolinium contrast, pulmonary blood vessel signal was greatly enhanced, providing sufficient contrast and structural information for non-rigid image registration to work in the lungs, which are well known to have lower MRI signal in the absence of contrast agents. The non-rigid registration methods were applied using the B-splines transformation, the NCC similarity measure, combined with the bending energy penalty term, and the multi-resolution approach with refined grid spacing. To quantitatively validate the accuracy of non-rigid image registration, pairs of anatomical landmarks at various places in the pulmonary vascular tree were identified in end-expiration and end-inspiration scans. Quantitative results demonstrated that, with carefully selected parameters, each pair of 3D CE-MRA images could be aligned well with average landmarks registration error less than a voxel size for each dataset, therefore providing reliable lung regional deformation estimation.

In this work, we observed that end-expiration scans had denser and brighter signals in the posterior area of the lungs. We speculate that this might be due to the gravity of the blood, the heart and adjacent structures (27), which could compress the posterior area of the lung when the patient is in the supine posture and result in higher density in those areas. As we used NCC as similarity measure in the registration, it can compensate for global intensity differences between end-expiration and end-inspiration images. Therefore, our methods can work well for global intensity change due to global volume change between end-inspiration and end-expiration scans. Furthermore, we expect our method to work well in the presence of signal intensity changes due to phased array coil profile alterations or coil signal intensity correction/normalization methods, as long as these changes are smooth in space. However, in the presence of large regional signal changes between end-inspiration and end-expiration, such as in the posterior area of the lungs, the registration accuracy may be reduced if these regions were included. This is a limitation of our technique. As a proof of concept study, only three subjects were studied in this work, which is a limitation; further studies, especially in patient populations, are warranted to evaluate the effectiveness of our method. In our study, the breath-holding time was approximately 25 seconds, which could be problematic for patients with lung disease. We used a modest 2× parallel imaging and further acceleration factor may be used for patients who cannot hold their breath for 25 seconds.

Our approach may potentially be useful for evaluating regional expansion/compression of the lung parenchyma, an important index for regional pulmonary ventilation assessment, and be useful for evaluating heterogeneous conditions such as COPD. Similar approaches have been studied based on CT images, and to our knowledge, our study was the first 3D contrast-enhanced MRI-based approach to evaluate 3D lung regional deformations. Aside from assessing lung ventilation using our approach, contrast-enhanced MRI also has the capability to measure the pulmonary perfusion. Hence, we speculate that once our technique is clinically validated, the lung ventilation measurements obtained from our technique may complement perfusion assessment from conventional first-pass contrast-enhanced MR (28) such that a perfusion-ventilation mismatch may be diagnosed in a variety of lung diseases.

Acknowledgements

None.

Footnote

Conflicts of Interest: The authors have no conflicts of interest to declare.

Ethical Statement: The study was approved by the Institutional Review Board (No. 11-002490) and written informed consent was obtained from all volunteers.

References

1. Petersson J, Sánchez-Crespo A, Larsson SA, Mure M. Physiological imaging of the lung: single-photon-emission computed tomography (SPECT). J Appl Physiol (1985) 2007;102:468-76. [Crossref] [PubMed]
2. Harris RS, Schuster DP. Visualizing lung function with positron emission tomography. J Appl Physiol (1985) 2007;102:448-58. [Crossref] [PubMed]
3. Chae EJ, Seo JB, Goo HW, Kim N, Song KS, Lee SD, Hong SJ, Krauss B. Xenon ventilation CT with a dual-energy technique of dual-source CT: initial experience. Radiology 2008;248:615-24. [Crossref] [PubMed]
4. Castillo R, Castillo E, Martinez J, Guerrero T. Ventilation from four-dimensional computed tomography: density versus Jacobian methods. Phys Med Biol 2010;55:4661-85. [Crossref] [PubMed]
5. Sase S, Nakano H, Suzuki H, Honda M. Subtraction lung image for evaluating pulmonary ventilation in xenon-enhanced CT. Med Phys 2010;37:4464-74. [Crossref] [PubMed]
6. Du K, Bayouth JE, Ding K, Christensen GE, Cao K, Reinhardt JM. Reproducibility of intensity-based estimates of lung ventilation. Med Phys 2013;40:063504. [Crossref] [PubMed]
7. Yamamoto T, Kabus S, Lorenz C, Mittra E, Hong JC, Chung M, Eclov N, To J, Diehn M, Loo BW Jr, Keall PJ. Pulmonary ventilation imaging based on 4-dimensional computed tomography: comparison with pulmonary function tests and SPECT ventilation images. Int J Radiat Oncol Biol Phys 2014;90:414-22. [Crossref] [PubMed]
8. Kauczor HU, Chen XJ, van Beek EJ, Schreiber WG. Pulmonary ventilation imaged by magnetic resonance: at the doorstep of clinical application. Eur Respir J 2001;17:1008-23. [Crossref] [PubMed]
9. Edelman RR, Hatabu H, Tadamura E, Li W, Prasad PV. Noninvasive assessment of regional ventilation in the human lung using oxygen-enhanced magnetic resonance imaging. Nat Med 1996;2:1236-9. [Crossref] [PubMed]
10. Wielpütz M, Kauczor HU. MRI of the lung: state of the art. Diagn Interv Radiol 2012;18:344-53. [PubMed]
11. Bauman G, Puderbach M, Deimling M, Jellus V, Chefd'hotel C, Dinkel J, Hintze C, Kauczor HU, Schad LR. Non-contrast-enhanced perfusion and ventilation assessment of the human lung by means of fourier decomposition in proton MRI. Magn Reson Med 2009;62:656-64. [Crossref] [PubMed]
12. Fain SB, Korosec FR, Holmes JH, O'Halloran R, Sorkness RL, Grist TM. Functional lung imaging using hyperpolarized gas MRI. J Magn Reson Imaging 2007;25:910-23. [Crossref] [PubMed]
13. Sá RC, Asadi AK, Theilmann RJ, Hopkins SR, Prisk GK, Darquenne C. Validating the distribution of specific ventilation in healthy humans measured using proton MR imaging. J Appl Physiol (1985) 2014;116:1048-56. [Crossref] [PubMed]
14. Coselmon MM, Balter JM, McShan DL, Kessler ML. Mutual information based CT registration of the lung at exhale and inhale breathing states using thin-plate splines. Med Phys 2004;31:2942-8. [Crossref] [PubMed]
15. Yim Y, Hong H, Shin YG. Deformable lung registration between exhale and inhale CT scans using active cells in a combined gradient force approach. Med Phys 2010;37:4307-17. [Crossref] [PubMed]
16. Gorbunova V, Sporring J, Lo P, Loeve M, Tiddens HA, Nielsen M, Dirksen A, de Bruijne M. Mass preserving image registration for lung CT. Med Image Anal 2012;16:786-95. [Crossref] [PubMed]
17. Yang D, Lu W, Low DA, Deasy JO, Hope AJ, El Naqa I. 4D-CT motion estimation using deformable image registration and 5D respiratory motion modeling. Med Phys 2008;35:4577-90. [Crossref] [PubMed]
18. Pennati F, Salito C, Baroni G, Woods J, Aliverti A. Comparison between multivolume CT-based surrogates of regional ventilation in healthy subjects. Acad Radiol 2014;21:1268-75. [Crossref] [PubMed]
19. Fain S, Schiebler ML, McCormack DG, Parraga G. Imaging of lung function using hyperpolarized helium-3 magnetic resonance imaging: Review of current and emerging translational methods and applications. J Magn Reson Imaging 2010;32:1398-408. [Crossref] [PubMed]
20. Gee J, Sundaram T, Hasegawa I, Uematsu H, Hatabu H. Characterization of Regional Pulmonary Mechanics from Serial MRI Data. In: Dohi T, Kikinis R. editors. Medical Image Computing and Computer-Assisted Intervention — MICCAI 2002. Berlin: Springer Berlin Heidelberg, 2002:795-9.
21. Murphy K, van Ginneken B, Reinhardt JM, Kabus S, Ding K, Deng X, Cao K, Du K, Christensen GE, Garcia V, Vercauteren T, Ayache N, Commowick O, Malandain G, Glocker B, Paragios N, Navab N, Gorbunova V, Sporring J, de Bruijne M, Han X, Heinrich MP, Schnabel JA, Jenkinson M, Lorenz C, Modat M, McClelland JR, Ourselin S, Muenzing SE, Viergever MA, De Nigris D, Collins DL, Arbel T, Peroni M, Li R, Sharp GC, Schmidt-Richberg A, Ehrhardt J, Werner R, Smeets D, Loeckx D, Song G, Tustison N, Avants B, Gee JC, Staring M, Klein S, Stoel BC, Urschler M, Werlberger M, Vandemeulebroucke J, Rit S, Sarrut D, Pluim JP. Evaluation of registration methods on thoracic CT: the EMPIRE10 challenge. IEEE Trans Med Imaging 2011;30:1901-20. [Crossref] [PubMed]
22. Klein S, Staring M, Murphy K, Viergever MA, Pluim JP. elastix: a toolbox for intensity-based medical image registration. IEEE Trans Med Imaging 2010;29:196-205. [Crossref] [PubMed]
23. Yin Y, Hoffman EA, Lin CL. Mass preserving nonrigid registration of CT lung images using cubic B-spline. Med Phys 2009;36:4213-22. [Crossref] [PubMed]
24. Kanai T, Kadoya N, Ito K, Onozato Y, Cho SY, Kishi K, Dobashi S, Umezawa R, Matsushita H, Takeda K, Jingu K. Evaluation of accuracy of B-spline transformation-based deformable image registration with different parameter settings for thoracic images. J Radiat Res 2014;55:1163-70. [Crossref] [PubMed]
25. Penney GP, Weese J, Little JA, Desmedt P, Hill DL, Hawkes DJ. A comparison of similarity measures for use in 2-D-3-D medical image registration. IEEE Trans Med Imaging 1998;17:586-95. [Crossref] [PubMed]
26. Klein S, Pluim JP, Staring M, Viergever MA. Adaptive stochastic gradient descent optimisation for image registration, Int J Compu. Vis 2009;81:227-39.
27. Hopkins SR, Henderson AC, Levin DL, Yamada K, Arai T, Buxton RB, Prisk GK. Vertical gradients in regional lung density and perfusion in the supine human lung: the Slinky effect. J Appl Physiol (1985) 2007;103:240-8. [Crossref] [PubMed]
28. Tsai SY, Wu MT, Lin YR, Hsieh KS, Lin CC, Huang TY, Chung HW, Pan JY, Huang YL, Pan HB, Yang CF. Temporal correlation-based dynamic contrast-enhanced MR imaging improves assessment of complex pulmonary circulation in congenital heart disease. Magn Reson Med 2006;56:517-26. [Crossref] [PubMed]