LUNG CT REGISTRATION COMBINING INTENSITY, CURVES AND SURFACES

LUNG CT REGISTRATION COMBINING INTENSITY, CURVES AND SURFACES Vladlena Gorbunova 1 ∗, Stanley Durrleman2,3, Pechin Lo1, Xavier Pennec2, Marleen de Bruijne1,4 1Department of Computer Science, Univercity of Copenhagen, Denmark 2INRIA Sophia Antipolis, Asclepios, France 3Centre de Mathematique et Leurs Application, ENS Cachan, France 4Biomedical Imaging Group Rotterdam, Erasmus MC, Rotterdam, the Netherlands ABSTRACT In this paper we propose a registration method that com- bines intensity information with geometrical information in the form of curves and surfaces derived from lung CT im- ages. Vessel tree centerlines and lung surfaces were extracted from segmented structures. First, a current-based registration was applied to align the pulmonary vessel tree and the lung surfaces. Subsequently, the resulting deformation field was used to constrain an intensity-based registration method. We applied the combined registration on a set of image pairs, ex- tracted at the end exhale and the end inhale phases of 4D- CT scans. The proposed combined registration was com- pared to intensity-based registration, using a set of manually selected landmarks. The proposed registration decreases the mean and the standard deviation of the target registration er- rors for all 5 cases to on average 1.47±1.05 mm, compared to the intensity-based registration without constraint 1.74±1.31 mm. Index Terms— Image registration, BSplines, currents, lung CT. 1. INTRODUCTION The ultimate goal of any registration algorithm is to es- tablish dense point-to-point correspondence between two im- ages. Generally, registration of lung CT images is a difficult problem due to the possible large variation between the scans. Scans of the same patient taken at maximum inspiration, can have more than 0.5 liter difference in lung volume. The reg- istration of end exhale and end inhale phases of 4D-CT lung images is an even more difficult problem due to the large and non-uniform deformations during the breathing cycle [1]. Image registration methods can be divided into two groups of methods: intensity-based and feature-based. A feature-based method establishes deformations based on low- dimensional features, derived from the original images, while ∗ This work is financially supported by the Danish Council for Strategic Research under the Programme Commission for Nanoscience and Technol- ogy, Biotechnology and IT (NABIIT), the Netherlands Organization for Sci- entific Research (NWO), and AstraZeneca, Lund, Sweden. Authors would like to thank Jon Sporring, Copenhagen University, Department of Computer Science, eScience Center, for fruitful discussions. intensity-based method considers intensity information over complete image. The state-of-the art registration methods for lung CT images are mainly intensity-based approaches [2] because the feature-based methods generally produce less accurate results [3]. Recently, Li et al. [4] developed an image registration al- gorithm where the intensity-based registration was improved with the subsequent bio-mechanical simulation of lung in- flation. Results showed an improvement in both accuracy of registration and physical plausibility of the deformation field for the combined approach. We previously developed a feature-based algorithm for registering lung CT images and compared it to intensity-based registration [5]. The overall accuracy of the feature-based algorithm was slightly worse than that of the intensity-based algorithm, but in 35 % of landmarks the feature-based registration outperformed the intensity-based method. The results inspired us to investi- gated how the intensity-based registration can be improved with the results of the feature-based registration. The direct combination of two completely different reg- istration methods is usually not possible, particularly if the underlying deformation models are different. For exam- ple, in parametric non-rigid registration, deformation fields are commonly modeled with b-spline functions, while in non-parametric methods deformation fields are usually mod- eled using partial differential equations. Furthermore, in landmark-based registration deformations are modeled by thin-plate splines or radial basis functions. We propose a solution to this problem and instead of combining the models we constrain the final results of the registration - the defor- mation fields - in a least square sense. While feature-based methods can more accurately estimate deformation fields of the features, the intensity-based method can benefit from its results and improve the overall accuracy of alignment further away from the features. The study [6] presents a similar solution to this problem, the registration algorithm which in- tegrates intensity-based and feature-based methods. The cost function incorporates difference in intensities and difference in the distances to the annotated surfaces. In this paper, we combine the previously developed 340978-1-4244-4126-6/10/$25.00 ©2010 IEEE ISBI 2010 feature-based algorithm with the B-spline intensity regis- tration algorithm and evaluate performance on 5 image pairs with manually annotated landmarks. 2. METHOD The section briefly recalls the feature-based registration method and the non-rigid intensity-based registration method. We propose to combine both approaches as it is described in Section 2.3. 2.1. Current-based registration In our previous work, we developed a feature-based regis- tration, where the vessel centerlines and lung surfaces, were used to establish correspondence between lung CT scans [5]. Both vessel centerlines and lung surfaces were represented in a framework of currents and aligned using the metric on currents. The current μ for a vessel centerline C is repre- sented by tangential vectors attached to the centerline points and for a triangulated surface S it is represented by normal directions attached to the centers of each face. Norm of a current μ(C), μ(S) is defined via a path integral, in case of curves, or flux integral for surfaces [7]. The cost function between anatomical lung structures in a fixed image Cf , Sf and a moving image Cm, Sm is defined as a weighted sum of the similarity measures between currents for the vessel cen- terlines Cf , Cm, the similarity between currents for surfaces the Sf , Sm, and a regularization term: E(Cf , Sf ;Cm, Sm) = γC ||μ(Cf )− μ(φ(Cm))||2W + γS ||μ(Sf )− μ(φ(Sm))||2W + γφ ∫ 1 0 ||vt||2V dt. (1) Diffeomorphic transformation φ of curves and surfaces was modeled in the framework of large deformation diffeomor- phic matching [7], where deformation of each feature point is defined by a velocity vector field vt = φ′t. The smooth velocity field vt is described via a Gaussian kernel with stan- dard deviation σV , where σV determines the typical scale of the deformations [8]. The smoothness of the currents is deter- mined by the parameter σW [8]. 2.2. Intensity-based registration via BSplines In this paper we used a multi-resolution image registra- tion framework similar to the framework developed in [9]. First, lung regions were extracted from the CT images and the background value was set to 0 HU. Images were aligned with affine transform TA. Subsequently, a series of 3 B-Spline transforms T i=1..3B−Spline with decreasing grid size was applied to the affinely registered images. Thus, the final deforma- tion is a composition of the affine transform and 3 levels of B-Spline transforms: Tfinal(x) = T 3B−Spline ◦ ... ◦ T 1B−Spline ◦ TA(x), (2) where x is a point in the moving image. We use the sum of squared intensity differences as the similarity measure be- tween the images, Eint(If , Im;T ) = 1 |Ω| ||If (x)− Im(T (x))|| 2 L2 , (3) where If is the fixed image, Im is the moving image andΩ the region of intersection. Each level was optimized separately using a stochastic gradient descent optimizer. 2.3. Combined registration We propose to constrain the intensity-based registration of Section 2.2 with the deformation field obtained from the current-based registration of Section 2.1. We constrain b- spline deformation field �Dbsp to match the given final defor- mation field �Dcurr by minimizing the L2 distance between the deformations. Since the current-based registration uses anatomical lung features to establish the correspondence, the deformation field in locations close to the extracted features is expected to be more reliable than further away from the features. Thus, we propose to incorporate a spatially vary- ing weight w(x) ∈ [0; 1], x ∈ Ω into the constrain between the deformation fields, which defines the trade off between matching intensity and deformations for every voxel x. The combined cost function then consists of the sum of squared intensity differences similarity function and constraint on the deformation field: E(If , Im;T ) = Eint + γEdef = 1 |Ω| ∫ Ω (1− w(x)) ||If (x)− Im(T (x))||2dx + λ |Ω| ∫ Ω w(x)|| �Dbsp(x)− �Dcur(x)||2dx, (4) where the coefficient λ compensates for the difference in units of the two terms. The deformation field �Dbsp(x) is a vector field defined as �Dbsp(x) = T (x) − x. Using vector notation, the gradient of the cost function (4) can be computed explic- itly: DaE(If , Im;T ) = − 2|Ω| ∫ Ω (1− w(x)) [If (x)− Im(T (x))] × [DxIm DaT ] dx − 2λ|Ω| ∫ Ω w(x)( �Dbsp(x)− �Dcur(x))TDaTdx. (5) The above method is naturally extended to an iterative ap- proach. After a level of the combined registration, the current- based registration is restarted with the deformed currents. The next level of combined registration starts from the final trans- form coefficients of the previous level and the new deforma- tion field obtained from the current-based registration. Using the described scheme we can iterate the current-based and the combined registration gradually improving the result. 341 3. EXPERIMENTS We conducted experiments on the five publicly available image pairs extracted at the end exhale and end inhale phases of the 4D-CT scans [10]. The study also provides 300 man- ually placed landmarks for each image pair. The landmarks were uniformly distributed over the lungs. In-plane resolution of the images varied from 0.97×0.97 mm to 1.16×1.16 mm and slice thickness was 2.5 mm. For each pair, an image ex- tracted at end inhale phase of 4D CT image was registered to an image extracted at end exhale phase. Lung fields, main bronchi and vessel tree were segmented as described in [11]. First, we applied the current-based registration [5] to register vessel trees and lung surfaces and computed the final deformation field for the whole image region. Then we applied the proposed registration, where the intensity-term was combined with the constraint on the de- formation fields as in Eq. (4). We iterated the two registration methods for the total number of iterations N = 2. For the first iteration, the parameters of the current-based registration were set to σ1W = 5 mm, σ 1 V = 25 mm and γ 1 φ = 10 −4. For the second iteration we decreased the smoothness kernel σ2W = 2.5 mm, σ 2 V = 25 mm and increase the γ 2 φ = 10 parameter in order to preserve more details of the currents and establish a locally accurate correspondence. Finally, we compared the results of the proposed com- bined registration to the registration with only the inten- sity term Eq.(3) and to the iterative registration where only current-based method is used. The current-based registra- tion was applied with the same parameters as in combined approach but the next iteration started from the results of the previous current-based registration. Internal parameters for the intensity-based and the combined registration were identical. (a) (b) Fig. 1. An example of spatially varying weights w(x) for the first 1(a) and the second 1(b) iteration of the combined registration. The coefficient λ in the Eq.(4) was set to 102 and 5 × 103 for the first and the second iterations respectively. The weights w(x) for the combined registration was constructed as follows. The lung surfaces were extracted from the seg- mented lungs. Then we erased the lung surfaces and vessel centerlines near the hilum area by first dilating the left and right main bronchus with a disk element of radius 20 voxels in axial plane and then deleting the constructed dilation from the lung surfaces and vessel centerlines. For the first iteration, we used lung surface alone and for the second iteration both lung surfaces and vessel centerlines. We computed the distance map to the constructed geometrical structures and evaluated the Gaussian kernel with the size κ1w = 2.0, κ 2 w = 5.0 mm on the distance image. Fig. 1 shows an example of a coronal slice of a weight image for the first 1(a) and the second 1(b) iteration. 4. RESULTS Visual comparison of the intensity-based registrations and the combined registration is presented in Fig. 2. Deformed images were interpolated using linear interpolation. Fig. 2. Right column shows every fourth slice from 36- 48 of the difference images between the fixed image and the moving image deformed after the combined registration. Left column shows corresponding difference images from the intensity-based registration. Difference images of the case 5 are shown in intensity window [−250; 250] HU. The overall accuracy of the image registration method was 342 Table 1. The mean and standard deviation of target registration error at the landmark positions in [mm] before the registration (Original); the current-based registration at each iteration (Curr It#); registration with combined cost (Comb It#); registration with intensity-cost (Intensity); after applying current-based registration (Curr). N Original Curr (It1) Comb (It1) Curr (It2) Comb (It2) Intensity Curr 1 3.89 ± 2.78 1.49 ± 0.75 1.16 ± 0.57 1.26 ± 0.73 1.15 ± 0.60 1.18 ± 0.57 1.44 ± 0.72 2 4.34 ± 3.90 2.26 ± 2.03 1.21 ± 0.64 1.15 ± 0.57 1.12 ± 0.55 1.26 ± 0.68 1.72 ± 1.38 3 6.94 ± 4.05 3.39 ± 3.09 1.79 ± 1.09 1.52 ± 0.87 1.46 ± 0.83 1.91 ± 1.15 2.97 ± 2.96 4 9.83 ± 4.86 3.90 ± 3.42 2.09 ± 1.48 1.72 ± 1.18 1.72 ± 1.16 2.12 ± 1.52 3.30 ± 2.61 5 7.48 ± 5.51 4.24 ± 3.34 2.14 ± 1.70 1.89 ± 1.53 1.92 ± 1.54 2.23 ± 1.79 3.52 ± 2.91 Average 6.50 3.06 1.68 1.51 1.47 1.74 2.59 defined as the mean Euclidean distance between the land- marks, target registration error (TRE), in millimeters. The mean and the standard deviation of TRE before registration, after registration with only the intensity term, after iterative registration using only the current-based method and after the proposed combined registration are reported in the Table 1. 5. DISCUSSION In this paper we presented a general framework for com- bining two registration methods. We combined the previously developed feature- and intensity-based registration using the constraint on deformation fields. We assumed that feature-based registration results in more accurate alignment of small, unclear structures, like small vessels where the gradient of the image is weak. Thus an intensity-based registration may result in a less accurate regis- tration of those structures. While both feature- and intensity- based methods implicitly use the intensity for registration, for the feature-based registration original intensities are less important. Segmentation process uses intensity and various derivatives of the intensity and results in a binary vessel tree. Thus large and small vessels are assigned the same value in feature-based registration whereas original intensities of those differ significantly. We supplement intensity information with the deformation field near the anatomical structures. The spa- tially varying weight defines both accuracy and location of the constraint. The maximum weight of 1 is at the lung border and the vessel centerline and decays elsewhere, thus implies the perfect fit of the deformation fields at the location of the segmented structures. But the actual effect of the constraint propagates within the support of the closest b-spline basis functions. The final solution brings minimum both to the sum of squared intensity differences cost far from an anatomical structure and the differences in the deformation fields close to it. Results show that both the feature-based registration and the intensity-based registration perform less accurate that the combined approach. Restarting the feature-based registration from the results of the combined registration result in bet- ter feature-based registration. Moreover the next iteration of the combined registration also improves it. We can conclude that the intensity-based registration is flexible enough to es- tablish the accurate transformation but lacks information near the lung border and small vessels. 6. REFERENCES [1] G. E. Christensen et al., “Tracking lung tissue motion and expansion/compression with inverse consistent image registra- tion and spirometry,” Med. Phys., vol. 34(6), pp. 2155–2164, 2007. [2] I. Sluimer et al., “Computer analysis of computed tomography scans of the lung: A survey.,” IEEE Trans. on Med. Imag., vol. 25, pp. 385–405, 2006. [3] S. Kabus et al., “Evaluation of 4D-CT lung registration.,” in MICCAI (1), 2009, vol. 5761 of LNCS, pp. 747–754. [4] P. Li et al., “Combination of intensity-based image registration with 3D simulation in radiation therapy,” Phys. Med. Biol, vol. 53, pp. 4621–4637, 2008. [5] V. Gorbunova et al., “Curve- and surface-based registration of lung CT images via currents,” in Second International Work- shop on Pulmonary Image Analysis, 2008. [6] H. Jingfeng et al., “Feature constrained non-rigid image regis- tration,” in 18th Symposium on Simulationstechnique, 2005. [7] Glaune`s et al, “Large deformation diffeomorphic metric curve mapping,” Int. J. Comput. Vision, vol. 80, no. 3, pp. 317–336, 2008. [8] S. Durrleman et al, “Inferring brain variability from diffeomor- phic deformations of currents: an integrative approach,” Med. Image Anal., vol. 12, no. 5, pp. 626–637, 2008. [9] V. Gorbunova et al., “Weight preserving image registration for monitoring disease progression in lung CT.,” in MICCAI, 2008, vol. 5242, pp. 863–870. [10] R. Castillo et al, “A framework for evaluation of deformable image registration spatial accuracy using large landmark point sets.,” Phys. Med. Biol., vol. 54, no. 7, pp. 1849–1870, 2009. [11] P. Lo et al., “Vessel-guided airway segmentation based on voxel classification,” in First International Workshop on Pul- monary Image Analysis, 2008. 343