- Email: [email protected]
Automatic lung nodule matching on sequential CT images
Computers in Biology and Medicine 38 (2008) 623 – 634 www.intl.elsevierhealth.com/journals/cobm
Automatic lung nodule matching on sequential CT images Helen Hong a,1 , Jeongjin Lee b,∗ , Yeny Yim b,2 a Division of Multimedia Engineering, College of Information and Media, Seoul Women’s University, 126 Gongreung-Dong, Nowon-Gu, Seoul 139-774, Korea b School of Computer Science and Engineering, Seoul National University, San 56-1 Shinlim 9-dong, Kwanak-gu, Seoul 151-742, Korea
Received 16 December 2005; accepted 29 February 2008
Abstract We propose an automatic segmentation and registration method that provides more efficient and robust matching of lung nodules in sequential chest computed tomography (CT) images. Our method consists of four steps. First, the lungs are extracted from chest CT images by the automatic segmentation method. Second, gross translational mismatch is corrected by optimal cube registration. This initial alignment does not require extracting any anatomical landmarks. Third, the initial alignment is step-by-step refined by hierarchical surface registration. To evaluate the distance measures between lung boundary points, a three-dimensional distance map is generated by narrow-band distance propagation, which drives fast and robust convergence to the optimal value. Finally, correspondences of manually detected nodules are established from the pairs with the smallest Euclidean distances. Experimental results show that our segmentation method accurately extracts lung boundaries and the registration method effectively finds the nodule correspondences. 䉷 2008 Elsevier Ltd. All rights reserved. Keywords: Surface registration; Affine transformation; Distance map; Distance measure; Optimization; CT; Lung nodule
1. Introduction Chest computed tomography (CT) is a well-established means of diagnosing pulmonary metastasis in oncology patients and evaluating lung disease progression and regression during treatment [1,2]. With ever-improving resolution and availability of CT scanners, the sensitivity for detection of pulmonary nodules has improved [3–5]. Recently, CT techniques have been applied to screening for lung cancer in high-risk individuals and are very promising for identifying early lung cancer [6–9]. With sequential follow-up CT scans, early changes in nodule size and number can be assessed [10]. However, it is often very difficult for radiologists to identify subtle changes, particularly in lesions that involve overlap with anatomic structures such as vessels, heart and diaphragm. Furthermore, thin-section CT scans of the whole thorax generates a large data set—typically 300–500 images of 1 mm section ∗ Corresponding author. Tel.: +82 2 3010 5050; fax: +82 2 476 0090.
E-mail addresses: [email protected] (H. Hong), [email protected] (J. Lee), [email protected] (Y. Yim). 1 Tel.: +82 2 970 5756; fax: +82 2 970 5981. 2 Tel.: +82 2 880 1860. 0010-4825/$ - see front matter 䉷 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.compbiomed.2008.02.010
thickness—and requires radiologists to spend a considerable amount of time for interpreting the images. Thus, the automatic matching of corresponding pulmonary nodules in sequential follow-up CT scans would be very useful for radiologists to improve follow-up study. Several methods have been suggested for the matching of lung nodules in sequential follow-up chest CT scans. In the method of Betke et al. [11–13], anatomical landmarks such as the sternum, vertebrae and tracheal centroids are used for initial global registration. Such initial alignment is then refined by an iterative closest point (ICP) process. Most of the computation time for the ICP process is in finding the point correspondences of lung boundaries obtained from two time-interval CT scans, as the Euclidean distance between each point pairs must be computed. Hong et al. [14] proposed a multi-level method for surface registration to cope with the problems of Betke et al. [11–13]. The multi-level method first reduces the original number of lung boundary points and aligns them using the ICP algorithm. The optimal solution for aligning the coarse surface models is then applied to match finer surface models. In addition, they proposed a midpoint approach to define the point correspondence instead of using the point with
624
H. Hong et al. / Computers in Biology and Medicine 38 (2008) 623 – 634
the smallest Euclidean distance as in the original ICP algorithm. However, the midpoint approach has a trade-off between accuracy and efficiency, because additional processing time is needed to find the second closest point and to compute the midpoint. Mullaly et al. [15] developed a multi-criterion nodule segmentation and registration method. Their registration method matches nodules based on affine registration of lung boundaries, combined with nodule size, location, shape and exclusivity measures. The method requires additional measures for geometric criteria. Gurcan et al. [16] developed an automated global matching of temporal thoracic helical CT scans. This method uses anatomical information from coronal maximum intensity projection (MIP) images such as ribs without requiring any anatomical landmark identification or organ segmentation. However, it would be difficult to align correctly since the method uses only limited information obtained by MIP images. Kusanagi et al. [17] proposed a method to find pairs of corresponding nodules by representing the nodule position as a relative coordinate in the lung region. The circumscribed quadrangles are set to left and right lung regions, respectively, in axial and coronal sections. The distance from the quadrangle origin to the nodule center is measured along the direction of the x-, y- and z-axes. The nodule position is represented by three components as the ratio of the distance and the length of quadrangle edge in each direction. However, in their experiments, false correspondences often occurred at the peripheral areas of the apex and bottom of the lungs because of respiration. Kubo et al. [18] proposed a method to find corresponding slice positions from previous and current images for comparative readings of CT images taken for mass screening. This method is based on template matching of the lungs, heart, aorta and blood vessels that finds nodule pairs in accordance with the relative position of each nodule in the lung region and their distance. Current approaches still need more progress in computational efficiency and accuracy for investigating changes of lung nodules in sequential chest CT scans. In this paper, we propose an automatic segmentation and registration method that provides more efficiency and robustness for matching lung nodules in successive chest CT scans. Our method consists of four steps. First, the lungs are extracted from the chest CT scans by the automatic segmentation method. Second, gross translational mismatch is corrected by optimal cube registration. This initial registration does not require extracting any anatomical landmarks. Third, the initial alignment is step-by-step refined by hierarchical surface registration. To evaluate the distance measures between lung boundary points, a three-dimensional (3D) distance map is generated by narrow-band distance propagation, which drives fast and robust convergence to the optimal value. Finally, nodule correspondences are established from the pairs with the smallest Euclidean distances. In our study, nodules are detected manually by an expert radiologist beforehand. The organization of the paper is as follows. In Section 2, we discuss how to extract the lung parenchyma from chest CT images and how to correct the geometrical mismatches in sequential chest CT images. In Section 3, experimental results show how our segmentation method accurately extracts the lungs and how our registration method rapidly aligns lung
nodules in two time-interval chest CT images. This paper is concluded in Section 4 with a brief discussion of the results. 2. Automatic lung nodule matching For the registration of the current CT scans, called the target volume, with the previous CT study, called the template volume, we apply the pipeline shown in Fig. 1. As the examinee in the screening for lung cancer is scanned in the course of the end-inspiration pause, most of the lung shape of the examinee is similar in each scan. From this, we assume that each CT scan is acquired at the maximal inspiration and the data set includes the thorax from the trachea to below the diaphragm. Based on these assumptions and our experience, we found that affine transformation would be sufficient for the matching of lung nodules in sequential chest CT scans. In our method, the centroids of pulmonary nodules are detected manually by an expert radiologist beforehand. The template and target scans from the same patient were read in the same session. Our method determines nodule correspondences automatically. Therefore, our method aids the radiologist in the evaluation of lung nodules for interval growth, which includes changes in nodule size, and investigation of the patient’s reaction to treatment regimens.
Fig. 1. The pipeline of automatic lung nodule matching.
H. Hong et al. / Computers in Biology and Medicine 38 (2008) 623 – 634
625
Fig. 2. The pipeline of lung boundary segmentation.
2.1. Lung boundary segmentation Lung segmentation is a precursor to the procedure for finding the corresponding lung nodules in successive chest CT images. With the introduction of spiral CT and more recently multidetector row CT techniques, it is critical to develop an efficient segmentation method that requires minimal or no human interaction. In this section, we describe an automatic segmentation method for accurately identifying the lungs in chest CT images. As shown in Fig. 2, the method consists of the following three steps: (1) the lung extraction step to identify the lungs, (2) the airway extraction step to segment trachea and left and right mainstem bronchi and (3) the subtraction step to separate the airway from the lungs. The goal of the lung extraction step is to separate the voxels of the lungs from the surrounding structures. Generally, grayscale thresholding and 3D region growing are used to identify the lungs [19,22–24]. As these methods based on difference in attenuation values may produce holes in highdensity vessels within the lungs, such unwanted interior cavities should be filled by morphological operators. However, it is difficult to decide the size of the structural elements of the morphological operators to eliminate the holes while distorting the lung boundary as little as possible. We solve this limitation by using the inverse operation of the two-dimensional (2D) region growing (iRG) method used to automatically segment the thorax from the background air and the lungs from the thorax, as shown in Figs. 3(b) and (c). Generally, in the CT image, air appears with a mean density of approximately −1000 Hounsfield units (HU), most lung tissue is in the range of −910 to −500 HU, while the chest wall, blood and bone is much more dense [19–21]. However, the tissue density in CT images varies between subjects according to radiation dose of the CT scanner. Instead of using a fixed threshold in 2D iRG, we use optimal thresholding [25] to automatically select a threshold value, Toptimal . Optimal thresholding allows us to adapt to variations in the chest CT
Fig. 3. The results of lung extraction: (a) chest CT image, (b) thorax extraction from the surrounding air, (c) lung separation from the thorax and (d) binary image construction.
images. For this process, we assume that the image contains only two principal brightness regions: (1) high-density regions, Rhigh , within the chest wall structures and (2) low-density regions, Rlow , in the lungs. The threshold is selected through an iterative procedure. We first select an initial threshold T0 and apply the initial threshold to the image to separate the pixels into high- and low-density regions. The new threshold Ti for the i-th step is the average of the mean densities of the high two regions i−1 , low i−1 . This threshold update procedure is repeated until there is no change in the threshold. To segment the thorax from the background air, the initial seed point is selected at (0, 0) in the upper left-hand corner of each 2D slice. Background air, which surrounds the body, is extracted by 2D region growing [26] and then the thorax is segmented by the inverse operation. To separate the lungs
626
H. Hong et al. / Computers in Biology and Medicine 38 (2008) 623 – 634
Fig. 4. The removal of unwanted regions: (a) the result of 2D iRG and (b) connected component labeling at low resolution eliminates bowel gas (indicated by square and displayed by enlarged image).
Fig. 5. The effect of median filtering (a) the result of airway extraction without median filtering and (b) the result of airway extraction with median filtering.
from the thorax, the seed point is chosen at a pixel of the thorax with a higher density value than the pre-selected threshold value, Toptimal by optimal thresholding. As a result, thorax that has similar density value to the seed point is extracted. By inverting the result, we obtain the segmented lungs and airways without unwanted interior cavities and distortion of the lung boundaries. Binary images are then constructed as shown in Fig. 3(d). The 3D connected component labeling [26] is applied to ensure that non-pulmonary structures, such as bowel gas, are not erroneously identified as lung regions, as shown in Fig. 4. To reduce memory and processing time, we perform 3D connected component labeling at low resolution and restore the result in the original images. As the densities of airways are similar to those of the lungs, the lungs resulting from the lung extraction step still contain the airways. Thus, the airway extraction step segments the trachea and left and right mainstem bronchi using our branch-based 3D region growing and subtracts the results from the results of the lung extraction step. Due to junctions between the lungs and the airways, the general 3D region growing for airway extraction may create leakage into the lung parenchyma, as shown in Fig. 5(a). To prevent this occurrence, as preprocessing we apply a 2D median filter with a 3 × 3 mask to each slice, which can make these junctions with weak contrast thin and separate the lungs and the airways, as shown in Fig. 5(b). Then branchbased 3D region growing with 26-connectivity is applied to the filtered images. The seed point is automatically searched for in the large, circular, air-filled region near the center of the first few slices in the data set. When the region growing proceeds from the trachea to the narrow airways, it often stops because of noise or insufficient resolution of the image. For this reason,
the region growing conditions have to be changed adaptively according to the narrow airway densities. Thus, we partition the airway into two parts, trachea and left and right mainstem bronchi, based on the branching point of the airways, called the carina. For the upper part, we obtain a threshold value using optimal thresholding. For the lower part, the region growing condition is defined using Eq. (1): − hI (x, y, z)
(1)
where I (x, y, z) is a current density value, is the average density of the upper part and is the standard deviation of the upper part; h is a parameter that controls the ease of region growing and is experimentally determined in the range from −0.4 to 0.4. When the region growing is finished, the high-density airway wall is often not included and unwanted cavities are still remained. We apply 2D binary dilation and closing with a 3 × 3 mask to each slice as a postprocessing step. We perform 2D dilation instead of 3D dilation to prevent anisotropic dilation along the z direction due to large slice thickness of CT data set. After the trachea and left and right mainstem bronchi are extracted, the results are subtracted from the results of the lung extraction step, as shown in Fig. 6. 2.2. Initial registration using an optimal cube According to the imaging protocol and the patient’s posture, the position of the lungs between the template and target volumes can be quite different. For the efficient registration of such volumes, an initial gross correction is usually applied.
H. Hong et al. / Computers in Biology and Medicine 38 (2008) 623 – 634
627
where P (x, y, z) is the center point of the optimal cube, Wa and Da are the ranges of the optimal cube along x- and y-axes in the axial section, and Hc is the range of the optimal cube along z-axis in the coronal section. The processing time for initial registration using an optimal cube is dramatically reduced, as it does not require any anatomical landmark detection. In addition, this method leads subsequent surface registration to robust convergence to the optimal value, as the search space is limited to near the optimal value. Fig. 6. The result of lung boundary segmentation (a) the result of the lung extraction step, (b) the result of the airway extraction step and (c) lungs extracted by image subtraction.
Fig. 7. Initial registration using an optimal cube: (a) optimal cube in axial section, (b) optimal cube in coronal section, (c) initial position of the optimal cube and (d) initial registration of the optimal cube.
Several landmark-based registrations have been used for initial gross correction [12–14]. To achieve initial alignment of the lungs, these landmark-based registrations require the detection of landmarks and point-to-point registration of corresponding landmarks. These processes greatly degrade the performance of the whole system. To minimize the processing time and maximize the effectiveness of the initial registration, we propose a simple method of global alignment using the circumscribed boundary of the lungs, as shown in Fig. 7. A bounding volume, which includes left and right lungs, is enlarged by a predefined distance, d, to be an optimal cube. The initial registration of the two volumes is accomplished by aligning the centers of the optimal cubes as in Eq. (2): Wa Da Hc P (x, y, z) = + d, + d, +d (2) 2 2 2
2.3. Hierarchical surface registration using narrow-band distance propagation In surface registration using a distance map, the calculation of distance from a surface boundary to a certain point can be done using a preprocessed distance map based on chamfer matching [27–29], a technique that finds the best fit between two sets of edge points by minimizing a generalized distance between them. Edge points of one image are transformed by a parametric transformation that describes how the images can be geometrically distorted in order to perform the registration. Chamfer matching reduces the generation time for a distance map using an approximated distance transformation compared to a Euclidean distance transformation [30,31]. However, the computation time for the distance is still expensive by the twostep distance transformation of forward and backward masks. In particular, when the initial alignment almost corrects the gross translational mismatch, the generation of a 3D distance map of the whole volume is unnecessary. From this observation, we propose a hierarchical surface registration using narrow-band distance propagation for the efficient generation of a distance map. To generate a 3D distance map, we approximate the global distance computation with repeated propagation of local distances within a predefined narrow-band. To approximate Euclidean distances, we consider 26-neighbor relations for 1-distance propagation, as shown in Eq. (3). The distance value tells how far it is from surface boundary points. The narrow-band distance propagation shown in Fig. 8 is applied to surface boundary points only in the template volume. We can generate a 3D distance map very quickly as pixels are propagated only in the direction of increasing distance to the maximum neighborhood: D(i) = min min (D(j ) + 1), D(i) (3) j ∈26-neighbors(i) where i and j are current and the 26-neighbor positions of the 3D distance map, respectively. The distance measure in Eq. (4) is used to determine the degree of resemblance of surface boundaries of the template and target volumes. The surface boundary of the target volume is superimposed on the distance map of the template volume, as shown in Fig. 9. The average of the distance values that the surface boundary hits is the measure of correspondence between the template and target volumes. The average of absolute distance difference (AADD) reaches a minimum when the
628
H. Hong et al. / Computers in Biology and Medicine 38 (2008) 623 – 634
Fig. 8. Generation of the 3D distance map using narrow-band distance propagation: (a) lung boundary, (b) distance 1 propagation, (c) distance 2 propagation and (d) distance dmax propagation.
the images are coarse and thus not very accurate. At high resolution, the images are accurate, but perhaps noisy and computations may be slow. Consequently, only final adjustments are made at high resolution. 2.4. Nodule registration for refinement
Fig. 9. The computation of the distance measure.
surface boundary points of the template and target volumes are correctly aligned: AADD =
N 1 |D(i)| N
(4)
i=1
where D(i) is the distance value of the i-th distance map in the template volume and N is the number of surface boundary points in the target volume. As the search space of our distance measure is limited to the surrounding surface boundaries, Powell’s multi-dimensional method [32–34] is sufficient for evaluating AADD instead of using a more powerful optimization algorithm. We calculate nine parameters, three translations, three rotations and three shears during the iterative minimization process. For minimization with nine parameters, three scaling factors, Sx , Sy , Sz , are directly computed with the ratio of the acquired voxel sizes. We define orders for Powell’s multi-dimensional optimization as (Tx , Ty , Rz , Rx , Ry , Tz , SHx , SHy , SHz ) where Tx , Ty , Tz , Rx , Ry , Rz , SHx , SHy , SHz are the coefficients of translation, rotation and shear about the three x-, y- and z-axes, respectively [35,36]. The search starts at the lowest resolution and the results from the low-resolution matching guides the computation at finer levels. Optimal local positions with the smallest distances in coarse resolution are used as initial positions for matching in the next finer resolution level. When the search at the zero level is completed, global optimal transformations correspond to the smallest distance. Generally, at low resolution, the computations are fast. In addition, noise and irrelevant small features are averaged out. On the other hand,
After the initial registration and subsequent hierarchical surface registration, lung nodules between template and target volumes may be mismatched as both registrations are methods for aligning lung boundaries rather than lung nodules. Thus, additional processing is required for finding accurate nodule correspondences between template and target volumes. The nodules of the template volume are transformed into target volume with the same affine transformation that matches the lung boundaries in the hierarchical surface registration step. As shown in Eq. (5), we define a corresponding mapping function, ND, based on the Euclidean distances between the centerof-mass points of nodules in the two volumes. In particular, point xi in X is the point corresponding to yi in Y if the Euclidean distance between T (xi ) and yi is the shortest among all distances between yi and any transformed point in X. Correspondences are established from the pairs with the smallest Euclidean distances: ND(yi ) = xi = arg min yi − T (xk ) xk ∈X
(5)
where xi and yi are the i-th nodules in template and target volumes, respectively. 3. Experimental results Our method has been applied to 20 pairs of successive chest CT scans with pulmonary nodules from patients who were referred to undergo thin-section (2 mm) or thick-section (5, 7, and 8 mm) CT scanning for the evaluation of pulmonary nodules. The intervals of time between initial and follow-up studies ranged from one to nine months. Twenty chest CT scans were acquired using a Philips Mx8000 16-slice multi-detector row CT scanner. Each image had a matrix size of 512 × 512 pixels with in-plane resolutions ranging from 0.51 to 0.70 mm. The number of images per scan ranged from 41 to 454.
H. Hong et al. / Computers in Biology and Medicine 38 (2008) 623 – 634
629
Fig. 10. The information of lung nodules in the 20 patients: (a) numbers and (b) average of the closest neighbor Euclidean distances between nodules.
In particular, subjects 1–5 had thin-section CT scanning with 351 image slices on average. The number of images and in-plane resolutions were different between the initial and follow-up scans of 18 patients. Fig. 10(a) shows the distribution of numbers of lung nodules in 20 patients. Eighty-six nodules were identified. Eighteen of 20 patients had one to nine pulmonary nodules and the other two patients had 14 and 21 nodules on CT images. There were no changes in the number of pulmonary nodules between initial and follow-up CT scans. Fig. 10(b) shows the average of the closest neighbor Euclidean distances between nodules in each CT images. The average of these distances for 20 patients in initial and follow-up CT scans was 77.45 and 76.70 mm, respectively. The lung volumes for the initial and follow-up CT scans of 20 patients are shown in Fig. 11. During the scanning, patients were supposed to be at maximal inspiration. Thus, lung volumes in the initial and follow-up scans should be similar for patients who were able to maintain maximal inspiration throughout both scans. The change in volume was 6.5% on average. There are three patients with relatively large changes, 15.4% on average, which indicate considerable differences in inspiration between scans. For the other 17 patients, an average change of only 4.9% was measured. The performance of our method was evaluated from the aspects of visual inspection, accuracy and processing time.
Fig. 11. The lung volumes from the initial and follow-up scans of the 20 patients.
Fig. 12 shows the results of lung boundary segmentation for subjects 1, 2 and 3. The second and third column shows segmented lung regions in 2D while the fourth column shows the shaded-surface display of the corresponding lungs. These results show that the proposed method accurately segments lung boundaries with large curvature. Fig. 13 shows coronal views of the lungs for subject 1, 2 and 4 with pulmonary nodules.
630
H. Hong et al. / Computers in Biology and Medicine 38 (2008) 623 – 634
Fig. 12. The results of lung boundary segmentation in subjects 1, 2 and 3. (a) Transverse CT images of the chest, (b) segmented lung boundary, (c) the transverse CT image enclosed by the segmented lung boundary and (d) shaded-surface display of the segmented lungs.
Fig. 13. The results of lung nodule matching in subjects 1, 2, and 4: (a) initial misalignment, (b) after optimal cube registration as initial registration, (c) after hierarchical surface registration and (d) after nodule registration.
H. Hong et al. / Computers in Biology and Medicine 38 (2008) 623 – 634
631
Fig. 14. Evaluation of the accuracy of corresponding nodules using the AED error in each subject.
Fig. 15. The reduction of AADD error and std. of ADD error per iteration in the hierarchical surface registration for the 20 patients at fine resolution: (a) AADD error and (b) standard deviation of ADD error.
632
H. Hong et al. / Computers in Biology and Medicine 38 (2008) 623 – 634
Fig. 16. Evaluation of the accuracy of corresponding lung boundaries using AADD error in each subject.
Fig. 17. The comparison of total processing time (a) without multi-resolution techniques and (b) with multi-resolution techniques.
Fig. 13(a) shows the initial misalignment which has translational, rotational and distorted differences. Fig. 13(b) shows the impact of our initial alignment. The positional differences
between corresponding lung nodules in the initial and followup CT scans shown in Fig. 13(a) are much reduced by the optimal cube registration, as shown in Fig. 13(b). This initial
H. Hong et al. / Computers in Biology and Medicine 38 (2008) 623 – 634
alignment is further refined by hierarchical surface registration until the lung boundaries of the initial and follow-up CT scans are aligned exactly, as in Fig. 13(c). From the exact matching of lung boundaries, lung nodule correspondences for each subject with nodules in the follow-up CT scans (dark gray) and nodules transformed from initial CT scans (light gray) into the follow-up CT scans are shown. Nodule correspondences for 86 (100%) out of 86 nodules were correctly established in the 20 patients. In Fig. 14, the results of the lung nodule alignment in the 20 patients are reported on a per-center-of-mass point basis using the average Euclidean distance (AED) between corresponding nodules in the initial and follow-up CT scans before and after alignment. Our initial registration reduces the initial AED error by 57.5% on average. The most substantial AED reduction in error in the initial misalignment is obtained with the initial registration and subsequent surface registration; 81.7% on average. Since the corresponding lung boundaries are almost aligned, the corresponding nodules are at or near the optimal position. After nodule registration, the average AED error of the initial misalignment is reduced from 30.1 to 1.8 mm. Fig. 15(a) shows how the error, the average of the AADD, is reduced by the hierarchical surface registration. In addition, Fig. 15(b) shows how the standard deviation of the absolute distance difference is reduced by our method. Since the positional difference is almost aligned by the initial registration, subsequent hierarchical surface registration rapidly converges to the optimal position. Following this, the optimal position for all subjects in the hierarchical registration is determined within 32 iterations on average. The average AADD error of the 20 patients is reduced by 1.4 voxels. To compare the accuracy of our method with conventional methods using a distance map, we evaluated the AADD error of the corresponding lung boundaries of the 20 patients using the hierarchical surface registration process shown in Fig. 16. The average AED error of the method based on Euclidean distance (Method 1), the method based on chamfer matching (Method 2) and our method (Method 3) were 1.58, 1.45 and 1.43 voxels, respectively. In particular, in subjects 1–5 with large data sets, the average reduction rate of the AADD error of Method 3 was relatively greater at 15.7% (Method 1) and 2.6% (Method 2). The total processing time is summarized in Fig. 17, where the execution time is measured for the initial registration, surface registration and nodule registration in comparison with Method 1 and Method 2. All our evaluations were performed on an Intel Pentium IV PC with a 3.0 GHz CPU and 2.0 GB of main memory. Due to optimal cube registration as an initial alignment, our distance map generation time is much faster than that of Method 1 and Method 2. As shown in Fig. 17(a), the average processing times of Methods 1, 2 and 3 without multiresolution techniques were 145.1, 95.5 and 74.8 s, respectively. On the other hand, the average processing times with multiresolution technique were 97.8, 49.1 and 37.1 s, respectively. With multi-resolution techniques, the total processing time is reduced by 56% on average. In particular, our method is 2.6 times faster than Method 1 and 1.3 times faster than Method 2 when
633
applying multi-resolution techniques. In addition, the processing time is much reduced when our method is applied to a large data set as in subjects 1–5. 4. Summary and conclusions In this paper, we have described an accurate and fast method for matching lung nodules in sequential chest CT scans. Our lung boundary segmentation extracts accurate lung boundaries through three steps: the lung extraction, airway extraction and subtraction steps. The 2D iRG can segment the thorax from the background air and separate the lungs from the thorax without holes such as unwanted interior cavities. The connected component labeling at low resolution can eliminate non-pulmonary structures such as bowel gas and reduce memory and processing time. Branch-based 3D region growing can successfully segment the trachea and the left and right mainstem bronchi using adaptive growing conditions. Our three-stage matching algorithm can find nodule correspondences efficiently. Using optimal cube registration, initial gross correction can be done rapidly and effectively without detecting anatomical landmarks. In the hierarchical surface registration, distance measure, using a 3D distance map generated by narrow-band distance propagation, allows rapid and robust convergence to the optimal position. For the refinement of the nodule correspondences, paring with the smallest Euclidean distances is performed. Experimental results show that our segmentation method gives accurate lung boundaries. In particular, lungs with large curvature or complicated shapes are accurately extracted. Our registration method was able to correctly find nodule correspondences in 20 patients. After nodule registration, the average AED error of initial misalignment is reduced from 30.1 mm down to 1.8 mm. In particular, in subjects obtained from thinsection scanning, the reduction rate of the AADD error in our method is relatively greater at 15.7% and 2.6% compared to methods based on Euclidean distance and chamfer matching. In addition, our method is 2.6 times faster than Euclidean distancebased registration and 1.3 times faster than chamfer matchingbased registration when applying multiresolution techniques. Acknowledgments This work was supported in part by Grant no. 10888 from the Seoul Research and Business Development Program and in part by Grant no. R01-2006-000-11244-0 from the Basic Research Program of the Korea Science and Engineering Foundation. The authors are grateful to Dr. Jin Mo Goo at the Seoul National University Hospital of Seoul, Korea, for providing the chest CT data sets shown in this paper and giving advice valuable to our research. References [1] E.R. Heitzman, The role of computed tomography in the diagnosis and management of lung cancer: an overview, Chest 89 (4) (1986) 237S–241S. [2] M.A. Howe, B.H. Gross, CT evaluation of the equivocal pulmonary nodule, Comput. Radiol. 11 (1987) 61–67.
634
H. Hong et al. / Computers in Biology and Medicine 38 (2008) 623 – 634
[3] M. Remy-Jardin, J. Remy, F. Giraud, C.H. Marquette, Pulmonary nodules: detection with thick-section spiral CT versus conventional CT, Radiology 187 (1993) 513–520. [4] F. Fischbach, F. Knollmann, V. Griesshaber, T. Frund, E. Akkol, R. Felix, Detection of pulmonary nodules by multiscale computed tomography: improved detection rate with reduced slice thickness, Eur. Radiol. 13 (2003) 2378–2383. [5] J.P. Ko, H. Rusinek, E.L. Jacobs, M. Betke, G. McGuinness, Small pulmonary nodules: volume measurement at chest CT—phantom study, Radiology 228 (2003) 864–870. [6] C.I. Henschke, D.I. McCauley, D.F. Yankelevitz, et al., Early lung cancer action project: overall design and findings from baseline screening, Lancet 354 (1999) 99–105. [7] K.T. Bae, J.S. Kim, Y.H. Na, K.G. Kim, J.H. Kim, Pulmonary nodules: automated detection on CT images with morphologic matching algorithm—preliminary results, Radiology 236 (1) (2005) 286–294. [8] M.F. McNitt-Gray, E.M. Hart, N. Wyckoff, J.W. Sayre, J.G. Goldin, D.R. Aberle, A pattern classification approach to characterizing solitary pulmonary nodules imaged on high resolution CT: preliminary results, Med. Phys. 26 (1999) 880–888. [9] T.W. Way, L.M. Hadjiiski, B. Sahiner, H.P. Chan, P.N. Cascade, E.A. Kazerooni, N. Bogot, C. Zhou, Computer-aided diagnosis of pulmonary nodules on CT scans: segmentation and classification using 3D active contours, Med. Phys. 33 (2006) 2323–2337. [10] D.F. Yankelevitz, A.P. Reeves, W.J. Kostis, B. Zhao, C.I. Henschke, Small pulmonary nodules: volumetrically determined growth rates based on CT evaluation, Radiology 217 (2000) 251–256. [11] J.P. Ko, M. Betke, Chest CT: automated nodule detection and assessment of change over time—preliminary experience, Radiology 218 (2001) 267–273. [12] M. Betke, H. Hong, J.P. Ko, Automatic 3D registration of lung surfaces in computed tomography scans, in: Proceedings of the Medical Image Computing and Computer-Assisted Intervention (MICCAI), 2001, pp. 725–733. [13] M. Betke, H. Hong, J.P. Ko, Automatic 3D registration of lung surfaces in computed tomography scans, CS Technical Report 2001-004, Boston University, 2001. [14] H. Hong, M. Betke, S. Teng, Multilevel 3D registration of lung surfaces in computed tomography scans—preliminary experience, in: Proceedings of the International Conference on Diagnostic Imaging and Analysis (ICDIA), 2002, pp. 90–95. [15] W. Mullaly, M. Betke, H. Hong, J. Wang, K. Mann, J.P. Ko, Multicriterion 3D segmentation and registration of pulmonary nodules on CT: a preliminary investigation, in: Proceedings of the International Conference on Diagnostic Imaging and Analysis (ICDIA), 2002, pp. 176–181. [16] M.N. Gurcan, R.C. Hardie, S.K. Rogers, D.E. Dozer, B.H. Allen, J.W. Hoffmeister, Automated global matching of temporal thoracic helical CT studies: feasibility study, Int. Congr. Ser. 1256 (2003) 1031–1036. [17] T. Kusanagi, Y. Mekada, J. Hasegawa, J. Toriwaki, M. Mori, H. Natori, Correspondence of lung nodules in sequential chest CT images for quantification of the curative effect, Int. Congr. Ser. 1256 (2003) 983–989. [18] M. Kubo, T. Yamamoto, Y. Kawata, N. Niki, K. Eguchi, H. Ohmatsu, R. Kakinuma, M. Kaneko, M. Kusumoto, N. Moriyama, K. Mori, H. Nishiyama, CAD system for the assistance of comparative reading for lung cancer using retrospective helical CT images, Proc. SPIE Med. Imaging 4322 (2001) 1788–1795. [19] S. Hu, E.A. Hoffman, J.M. Reinhardt, Accurate lung segmentation for accurate quantitation of volumetric X-ray CT images, IEEE Trans. Med. Imaging 20 (6) (2001) 490–498. [20] M.S. Brown, M.F. McNitt-Gray, N.J. Mankovich, J.G. Goldin, J. Hiller, L.S. Wilson, D.R. Aberle, Method for segmenting chest CT image data using an anatomical model: preliminary results, IEEE Trans. Med. Imaging 16 (6) (1997) 828–839. [21] M. Wu, J. Chang, A.A. Chiang, J. Lu, H. Hsu, W. Hsu, C. Yang, Use of quantitative CT to predict postoperative lung function in patients with lung cancer, Radiology 191 (1994) 257–262.
[22] S.G. Armato III, H. MacMahon, Automated lung segmentation and computer-aided diagnosis for thoracic CT scans, Int. Congr. Ser. 1256 (2003) 977–982. [23] S.G. Armato III, W.F. Sensakovic, Automated lung segmentation for thoracic CT: impact on computer-aided diagnosis, Acad. Radiol. 11 (2004) 1011–1021. [24] J.K. Leader, B. Zheng, R.M. Rogers, F.C. Sciurba, A. Perez, B.E. Chapman, S. Patel, C.R. Fuhrman, D. Gur, Automated lung segmentation in X-ray computed tomography, Acad. Radiol. 20 (2003) 1224–1236. [25] M. Sonka, V. Hlavac, R. Boyle, Image Processing, Analysis and Machine Vision, Pacific Grove, PWS, CA, 1999. [26] R.C. Gonzalez, R.E. Woods, Digital Image Processing, second ed., Prentice-Hall, Englewood Cliffs, NJ. [27] G. Borgefors, Hierarchical chamfer matching: a parametric edge matching algorithm, IEEE Trans. Pattern Recognition Mach. Intell. 10 (6) (1988) 849–865. [28] O. Cuisenaire, Distance transformations: fast algorithms and applications to medical image processing, Ph.D. Thesis, 1999. [29] B. Fang, W. Hsu, M.L. Lee, Techniques for temporal registration of retinal images, in: Proceedings of the International Conference on Image Processing (ICIP), 2004, pp. 1089–1092. [30] G. Borgefors, Distance transformations in digital images, Comput. Vision Graphics Image Processing 34 (3) (1986) 344–371. [31] T. Saito, J.-I. Toriwaki, New algorithms for Euclidean distance transformation of an n-dimensional digitized picture with applications, Pattern Recognition 27 (11) (1994) 1551–1565. [32] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical recipes in C—the art of scientific computing, Cambridge University Press, Cambridge, 1992. [33] G.K. Matsopoulos, N.A. Mouravliansky, K.K. Delibasis, K.S. Nikita, Automatic retinal image registration scheme using global optimization techniques, IEEE Trans. Inf. Technol. Biomed. 3 (1) (1999) 47–60. [34] N. Mouravliansky, G.K. Matsopoulos, K. Delibasis, K.S. Nikita, Automatic retinal registration using global optimization techniques, in: Proceedings of the 20th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, vol. 20(2), 1998, pp. 567–570. [35] J.L. Bernon, V. Boudousq, J.F. Rohmer, M. Fourcade, M. Zanca, M. Rossi, D. Mariano-Goulart, A comparative study of Powell’s and Downhill Simplex algorithms for a fast multimodal surface matching in brain imaging, Comput. Med. Imaging Graphics 25 (2001) 287–297. [36] F. Maes, A. Collignon, D. Vandermeulen, G. Marchal, P. Suetens, Multimodality image registration by maximization of mutual information, IEEE Trans. Med. Imaging 16 (2) (1997) 187–198.
Helen Hong received BS degree (1994), MS degree (1996) and PhD degree (2001) in Computer Engineering at the Ewha Womans University, Seoul, Korea. From 2001 to 2003, she was at the Seoul National University as a post doctor. In 2003, she was at the INFINITT Co. Ltd. as a principal researcher. From 2003 to 2005, she was at the Seoul National University as a research professor. Since 2006, she is a professor in the Division of Multimedia Engineering at the Seoul Women’s University, Seoul, Korea. Her current research interests include image registration, medical image processing and analysis, volume visualization and computer graphics. Jeongjin Lee received BS degree (2000) in Department of Mechanical Engineering and MS degree (2002), PhD degree (2008) in School of Computer Science and Engineering from Seoul National University, Korea. Since 2007, he is a research professor in Department of Radiology University of Ulsan, College of Medicine, Korea. His current research interests include medical image registration, medical image segmentation, virtual colonoscopy, deformable modeling, and computer animation. Yeny Yim received BS degree (2001) in Computer Engineering at the Kyungpook National University, Daegu, Korea. Now, she is a PhD candidate of School of Computer Science and Engineering at the Seoul National University, Seoul, Korea. Her current research interests include medical image segmentation and registration.
Automatic lung nodule matching on sequential CT images Helen Hong a,1 , Jeongjin Lee b,∗ , Yeny Yim b,2 a Division of Multimedia Engineering, College of Information and Media, Seoul Women’s University, 126 Gongreung-Dong, Nowon-Gu, Seoul 139-774, Korea b School of Computer Science and Engineering, Seoul National University, San 56-1 Shinlim 9-dong, Kwanak-gu, Seoul 151-742, Korea
Received 16 December 2005; accepted 29 February 2008
Abstract We propose an automatic segmentation and registration method that provides more efficient and robust matching of lung nodules in sequential chest computed tomography (CT) images. Our method consists of four steps. First, the lungs are extracted from chest CT images by the automatic segmentation method. Second, gross translational mismatch is corrected by optimal cube registration. This initial alignment does not require extracting any anatomical landmarks. Third, the initial alignment is step-by-step refined by hierarchical surface registration. To evaluate the distance measures between lung boundary points, a three-dimensional distance map is generated by narrow-band distance propagation, which drives fast and robust convergence to the optimal value. Finally, correspondences of manually detected nodules are established from the pairs with the smallest Euclidean distances. Experimental results show that our segmentation method accurately extracts lung boundaries and the registration method effectively finds the nodule correspondences. 䉷 2008 Elsevier Ltd. All rights reserved. Keywords: Surface registration; Affine transformation; Distance map; Distance measure; Optimization; CT; Lung nodule
1. Introduction Chest computed tomography (CT) is a well-established means of diagnosing pulmonary metastasis in oncology patients and evaluating lung disease progression and regression during treatment [1,2]. With ever-improving resolution and availability of CT scanners, the sensitivity for detection of pulmonary nodules has improved [3–5]. Recently, CT techniques have been applied to screening for lung cancer in high-risk individuals and are very promising for identifying early lung cancer [6–9]. With sequential follow-up CT scans, early changes in nodule size and number can be assessed [10]. However, it is often very difficult for radiologists to identify subtle changes, particularly in lesions that involve overlap with anatomic structures such as vessels, heart and diaphragm. Furthermore, thin-section CT scans of the whole thorax generates a large data set—typically 300–500 images of 1 mm section ∗ Corresponding author. Tel.: +82 2 3010 5050; fax: +82 2 476 0090.
E-mail addresses: [email protected] (H. Hong), [email protected] (J. Lee), [email protected] (Y. Yim). 1 Tel.: +82 2 970 5756; fax: +82 2 970 5981. 2 Tel.: +82 2 880 1860. 0010-4825/$ - see front matter 䉷 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.compbiomed.2008.02.010
thickness—and requires radiologists to spend a considerable amount of time for interpreting the images. Thus, the automatic matching of corresponding pulmonary nodules in sequential follow-up CT scans would be very useful for radiologists to improve follow-up study. Several methods have been suggested for the matching of lung nodules in sequential follow-up chest CT scans. In the method of Betke et al. [11–13], anatomical landmarks such as the sternum, vertebrae and tracheal centroids are used for initial global registration. Such initial alignment is then refined by an iterative closest point (ICP) process. Most of the computation time for the ICP process is in finding the point correspondences of lung boundaries obtained from two time-interval CT scans, as the Euclidean distance between each point pairs must be computed. Hong et al. [14] proposed a multi-level method for surface registration to cope with the problems of Betke et al. [11–13]. The multi-level method first reduces the original number of lung boundary points and aligns them using the ICP algorithm. The optimal solution for aligning the coarse surface models is then applied to match finer surface models. In addition, they proposed a midpoint approach to define the point correspondence instead of using the point with
624
H. Hong et al. / Computers in Biology and Medicine 38 (2008) 623 – 634
the smallest Euclidean distance as in the original ICP algorithm. However, the midpoint approach has a trade-off between accuracy and efficiency, because additional processing time is needed to find the second closest point and to compute the midpoint. Mullaly et al. [15] developed a multi-criterion nodule segmentation and registration method. Their registration method matches nodules based on affine registration of lung boundaries, combined with nodule size, location, shape and exclusivity measures. The method requires additional measures for geometric criteria. Gurcan et al. [16] developed an automated global matching of temporal thoracic helical CT scans. This method uses anatomical information from coronal maximum intensity projection (MIP) images such as ribs without requiring any anatomical landmark identification or organ segmentation. However, it would be difficult to align correctly since the method uses only limited information obtained by MIP images. Kusanagi et al. [17] proposed a method to find pairs of corresponding nodules by representing the nodule position as a relative coordinate in the lung region. The circumscribed quadrangles are set to left and right lung regions, respectively, in axial and coronal sections. The distance from the quadrangle origin to the nodule center is measured along the direction of the x-, y- and z-axes. The nodule position is represented by three components as the ratio of the distance and the length of quadrangle edge in each direction. However, in their experiments, false correspondences often occurred at the peripheral areas of the apex and bottom of the lungs because of respiration. Kubo et al. [18] proposed a method to find corresponding slice positions from previous and current images for comparative readings of CT images taken for mass screening. This method is based on template matching of the lungs, heart, aorta and blood vessels that finds nodule pairs in accordance with the relative position of each nodule in the lung region and their distance. Current approaches still need more progress in computational efficiency and accuracy for investigating changes of lung nodules in sequential chest CT scans. In this paper, we propose an automatic segmentation and registration method that provides more efficiency and robustness for matching lung nodules in successive chest CT scans. Our method consists of four steps. First, the lungs are extracted from the chest CT scans by the automatic segmentation method. Second, gross translational mismatch is corrected by optimal cube registration. This initial registration does not require extracting any anatomical landmarks. Third, the initial alignment is step-by-step refined by hierarchical surface registration. To evaluate the distance measures between lung boundary points, a three-dimensional (3D) distance map is generated by narrow-band distance propagation, which drives fast and robust convergence to the optimal value. Finally, nodule correspondences are established from the pairs with the smallest Euclidean distances. In our study, nodules are detected manually by an expert radiologist beforehand. The organization of the paper is as follows. In Section 2, we discuss how to extract the lung parenchyma from chest CT images and how to correct the geometrical mismatches in sequential chest CT images. In Section 3, experimental results show how our segmentation method accurately extracts the lungs and how our registration method rapidly aligns lung
nodules in two time-interval chest CT images. This paper is concluded in Section 4 with a brief discussion of the results. 2. Automatic lung nodule matching For the registration of the current CT scans, called the target volume, with the previous CT study, called the template volume, we apply the pipeline shown in Fig. 1. As the examinee in the screening for lung cancer is scanned in the course of the end-inspiration pause, most of the lung shape of the examinee is similar in each scan. From this, we assume that each CT scan is acquired at the maximal inspiration and the data set includes the thorax from the trachea to below the diaphragm. Based on these assumptions and our experience, we found that affine transformation would be sufficient for the matching of lung nodules in sequential chest CT scans. In our method, the centroids of pulmonary nodules are detected manually by an expert radiologist beforehand. The template and target scans from the same patient were read in the same session. Our method determines nodule correspondences automatically. Therefore, our method aids the radiologist in the evaluation of lung nodules for interval growth, which includes changes in nodule size, and investigation of the patient’s reaction to treatment regimens.
Fig. 1. The pipeline of automatic lung nodule matching.
H. Hong et al. / Computers in Biology and Medicine 38 (2008) 623 – 634
625
Fig. 2. The pipeline of lung boundary segmentation.
2.1. Lung boundary segmentation Lung segmentation is a precursor to the procedure for finding the corresponding lung nodules in successive chest CT images. With the introduction of spiral CT and more recently multidetector row CT techniques, it is critical to develop an efficient segmentation method that requires minimal or no human interaction. In this section, we describe an automatic segmentation method for accurately identifying the lungs in chest CT images. As shown in Fig. 2, the method consists of the following three steps: (1) the lung extraction step to identify the lungs, (2) the airway extraction step to segment trachea and left and right mainstem bronchi and (3) the subtraction step to separate the airway from the lungs. The goal of the lung extraction step is to separate the voxels of the lungs from the surrounding structures. Generally, grayscale thresholding and 3D region growing are used to identify the lungs [19,22–24]. As these methods based on difference in attenuation values may produce holes in highdensity vessels within the lungs, such unwanted interior cavities should be filled by morphological operators. However, it is difficult to decide the size of the structural elements of the morphological operators to eliminate the holes while distorting the lung boundary as little as possible. We solve this limitation by using the inverse operation of the two-dimensional (2D) region growing (iRG) method used to automatically segment the thorax from the background air and the lungs from the thorax, as shown in Figs. 3(b) and (c). Generally, in the CT image, air appears with a mean density of approximately −1000 Hounsfield units (HU), most lung tissue is in the range of −910 to −500 HU, while the chest wall, blood and bone is much more dense [19–21]. However, the tissue density in CT images varies between subjects according to radiation dose of the CT scanner. Instead of using a fixed threshold in 2D iRG, we use optimal thresholding [25] to automatically select a threshold value, Toptimal . Optimal thresholding allows us to adapt to variations in the chest CT
Fig. 3. The results of lung extraction: (a) chest CT image, (b) thorax extraction from the surrounding air, (c) lung separation from the thorax and (d) binary image construction.
images. For this process, we assume that the image contains only two principal brightness regions: (1) high-density regions, Rhigh , within the chest wall structures and (2) low-density regions, Rlow , in the lungs. The threshold is selected through an iterative procedure. We first select an initial threshold T0 and apply the initial threshold to the image to separate the pixels into high- and low-density regions. The new threshold Ti for the i-th step is the average of the mean densities of the high two regions i−1 , low i−1 . This threshold update procedure is repeated until there is no change in the threshold. To segment the thorax from the background air, the initial seed point is selected at (0, 0) in the upper left-hand corner of each 2D slice. Background air, which surrounds the body, is extracted by 2D region growing [26] and then the thorax is segmented by the inverse operation. To separate the lungs
626
H. Hong et al. / Computers in Biology and Medicine 38 (2008) 623 – 634
Fig. 4. The removal of unwanted regions: (a) the result of 2D iRG and (b) connected component labeling at low resolution eliminates bowel gas (indicated by square and displayed by enlarged image).
Fig. 5. The effect of median filtering (a) the result of airway extraction without median filtering and (b) the result of airway extraction with median filtering.
from the thorax, the seed point is chosen at a pixel of the thorax with a higher density value than the pre-selected threshold value, Toptimal by optimal thresholding. As a result, thorax that has similar density value to the seed point is extracted. By inverting the result, we obtain the segmented lungs and airways without unwanted interior cavities and distortion of the lung boundaries. Binary images are then constructed as shown in Fig. 3(d). The 3D connected component labeling [26] is applied to ensure that non-pulmonary structures, such as bowel gas, are not erroneously identified as lung regions, as shown in Fig. 4. To reduce memory and processing time, we perform 3D connected component labeling at low resolution and restore the result in the original images. As the densities of airways are similar to those of the lungs, the lungs resulting from the lung extraction step still contain the airways. Thus, the airway extraction step segments the trachea and left and right mainstem bronchi using our branch-based 3D region growing and subtracts the results from the results of the lung extraction step. Due to junctions between the lungs and the airways, the general 3D region growing for airway extraction may create leakage into the lung parenchyma, as shown in Fig. 5(a). To prevent this occurrence, as preprocessing we apply a 2D median filter with a 3 × 3 mask to each slice, which can make these junctions with weak contrast thin and separate the lungs and the airways, as shown in Fig. 5(b). Then branchbased 3D region growing with 26-connectivity is applied to the filtered images. The seed point is automatically searched for in the large, circular, air-filled region near the center of the first few slices in the data set. When the region growing proceeds from the trachea to the narrow airways, it often stops because of noise or insufficient resolution of the image. For this reason,
the region growing conditions have to be changed adaptively according to the narrow airway densities. Thus, we partition the airway into two parts, trachea and left and right mainstem bronchi, based on the branching point of the airways, called the carina. For the upper part, we obtain a threshold value using optimal thresholding. For the lower part, the region growing condition is defined using Eq. (1): − hI (x, y, z)
(1)
where I (x, y, z) is a current density value, is the average density of the upper part and is the standard deviation of the upper part; h is a parameter that controls the ease of region growing and is experimentally determined in the range from −0.4 to 0.4. When the region growing is finished, the high-density airway wall is often not included and unwanted cavities are still remained. We apply 2D binary dilation and closing with a 3 × 3 mask to each slice as a postprocessing step. We perform 2D dilation instead of 3D dilation to prevent anisotropic dilation along the z direction due to large slice thickness of CT data set. After the trachea and left and right mainstem bronchi are extracted, the results are subtracted from the results of the lung extraction step, as shown in Fig. 6. 2.2. Initial registration using an optimal cube According to the imaging protocol and the patient’s posture, the position of the lungs between the template and target volumes can be quite different. For the efficient registration of such volumes, an initial gross correction is usually applied.
H. Hong et al. / Computers in Biology and Medicine 38 (2008) 623 – 634
627
where P (x, y, z) is the center point of the optimal cube, Wa and Da are the ranges of the optimal cube along x- and y-axes in the axial section, and Hc is the range of the optimal cube along z-axis in the coronal section. The processing time for initial registration using an optimal cube is dramatically reduced, as it does not require any anatomical landmark detection. In addition, this method leads subsequent surface registration to robust convergence to the optimal value, as the search space is limited to near the optimal value. Fig. 6. The result of lung boundary segmentation (a) the result of the lung extraction step, (b) the result of the airway extraction step and (c) lungs extracted by image subtraction.
Fig. 7. Initial registration using an optimal cube: (a) optimal cube in axial section, (b) optimal cube in coronal section, (c) initial position of the optimal cube and (d) initial registration of the optimal cube.
Several landmark-based registrations have been used for initial gross correction [12–14]. To achieve initial alignment of the lungs, these landmark-based registrations require the detection of landmarks and point-to-point registration of corresponding landmarks. These processes greatly degrade the performance of the whole system. To minimize the processing time and maximize the effectiveness of the initial registration, we propose a simple method of global alignment using the circumscribed boundary of the lungs, as shown in Fig. 7. A bounding volume, which includes left and right lungs, is enlarged by a predefined distance, d, to be an optimal cube. The initial registration of the two volumes is accomplished by aligning the centers of the optimal cubes as in Eq. (2): Wa Da Hc P (x, y, z) = + d, + d, +d (2) 2 2 2
2.3. Hierarchical surface registration using narrow-band distance propagation In surface registration using a distance map, the calculation of distance from a surface boundary to a certain point can be done using a preprocessed distance map based on chamfer matching [27–29], a technique that finds the best fit between two sets of edge points by minimizing a generalized distance between them. Edge points of one image are transformed by a parametric transformation that describes how the images can be geometrically distorted in order to perform the registration. Chamfer matching reduces the generation time for a distance map using an approximated distance transformation compared to a Euclidean distance transformation [30,31]. However, the computation time for the distance is still expensive by the twostep distance transformation of forward and backward masks. In particular, when the initial alignment almost corrects the gross translational mismatch, the generation of a 3D distance map of the whole volume is unnecessary. From this observation, we propose a hierarchical surface registration using narrow-band distance propagation for the efficient generation of a distance map. To generate a 3D distance map, we approximate the global distance computation with repeated propagation of local distances within a predefined narrow-band. To approximate Euclidean distances, we consider 26-neighbor relations for 1-distance propagation, as shown in Eq. (3). The distance value tells how far it is from surface boundary points. The narrow-band distance propagation shown in Fig. 8 is applied to surface boundary points only in the template volume. We can generate a 3D distance map very quickly as pixels are propagated only in the direction of increasing distance to the maximum neighborhood: D(i) = min min (D(j ) + 1), D(i) (3) j ∈26-neighbors(i) where i and j are current and the 26-neighbor positions of the 3D distance map, respectively. The distance measure in Eq. (4) is used to determine the degree of resemblance of surface boundaries of the template and target volumes. The surface boundary of the target volume is superimposed on the distance map of the template volume, as shown in Fig. 9. The average of the distance values that the surface boundary hits is the measure of correspondence between the template and target volumes. The average of absolute distance difference (AADD) reaches a minimum when the
628
H. Hong et al. / Computers in Biology and Medicine 38 (2008) 623 – 634
Fig. 8. Generation of the 3D distance map using narrow-band distance propagation: (a) lung boundary, (b) distance 1 propagation, (c) distance 2 propagation and (d) distance dmax propagation.
the images are coarse and thus not very accurate. At high resolution, the images are accurate, but perhaps noisy and computations may be slow. Consequently, only final adjustments are made at high resolution. 2.4. Nodule registration for refinement
Fig. 9. The computation of the distance measure.
surface boundary points of the template and target volumes are correctly aligned: AADD =
N 1 |D(i)| N
(4)
i=1
where D(i) is the distance value of the i-th distance map in the template volume and N is the number of surface boundary points in the target volume. As the search space of our distance measure is limited to the surrounding surface boundaries, Powell’s multi-dimensional method [32–34] is sufficient for evaluating AADD instead of using a more powerful optimization algorithm. We calculate nine parameters, three translations, three rotations and three shears during the iterative minimization process. For minimization with nine parameters, three scaling factors, Sx , Sy , Sz , are directly computed with the ratio of the acquired voxel sizes. We define orders for Powell’s multi-dimensional optimization as (Tx , Ty , Rz , Rx , Ry , Tz , SHx , SHy , SHz ) where Tx , Ty , Tz , Rx , Ry , Rz , SHx , SHy , SHz are the coefficients of translation, rotation and shear about the three x-, y- and z-axes, respectively [35,36]. The search starts at the lowest resolution and the results from the low-resolution matching guides the computation at finer levels. Optimal local positions with the smallest distances in coarse resolution are used as initial positions for matching in the next finer resolution level. When the search at the zero level is completed, global optimal transformations correspond to the smallest distance. Generally, at low resolution, the computations are fast. In addition, noise and irrelevant small features are averaged out. On the other hand,
After the initial registration and subsequent hierarchical surface registration, lung nodules between template and target volumes may be mismatched as both registrations are methods for aligning lung boundaries rather than lung nodules. Thus, additional processing is required for finding accurate nodule correspondences between template and target volumes. The nodules of the template volume are transformed into target volume with the same affine transformation that matches the lung boundaries in the hierarchical surface registration step. As shown in Eq. (5), we define a corresponding mapping function, ND, based on the Euclidean distances between the centerof-mass points of nodules in the two volumes. In particular, point xi in X is the point corresponding to yi in Y if the Euclidean distance between T (xi ) and yi is the shortest among all distances between yi and any transformed point in X. Correspondences are established from the pairs with the smallest Euclidean distances: ND(yi ) = xi = arg min yi − T (xk ) xk ∈X
(5)
where xi and yi are the i-th nodules in template and target volumes, respectively. 3. Experimental results Our method has been applied to 20 pairs of successive chest CT scans with pulmonary nodules from patients who were referred to undergo thin-section (2 mm) or thick-section (5, 7, and 8 mm) CT scanning for the evaluation of pulmonary nodules. The intervals of time between initial and follow-up studies ranged from one to nine months. Twenty chest CT scans were acquired using a Philips Mx8000 16-slice multi-detector row CT scanner. Each image had a matrix size of 512 × 512 pixels with in-plane resolutions ranging from 0.51 to 0.70 mm. The number of images per scan ranged from 41 to 454.
H. Hong et al. / Computers in Biology and Medicine 38 (2008) 623 – 634
629
Fig. 10. The information of lung nodules in the 20 patients: (a) numbers and (b) average of the closest neighbor Euclidean distances between nodules.
In particular, subjects 1–5 had thin-section CT scanning with 351 image slices on average. The number of images and in-plane resolutions were different between the initial and follow-up scans of 18 patients. Fig. 10(a) shows the distribution of numbers of lung nodules in 20 patients. Eighty-six nodules were identified. Eighteen of 20 patients had one to nine pulmonary nodules and the other two patients had 14 and 21 nodules on CT images. There were no changes in the number of pulmonary nodules between initial and follow-up CT scans. Fig. 10(b) shows the average of the closest neighbor Euclidean distances between nodules in each CT images. The average of these distances for 20 patients in initial and follow-up CT scans was 77.45 and 76.70 mm, respectively. The lung volumes for the initial and follow-up CT scans of 20 patients are shown in Fig. 11. During the scanning, patients were supposed to be at maximal inspiration. Thus, lung volumes in the initial and follow-up scans should be similar for patients who were able to maintain maximal inspiration throughout both scans. The change in volume was 6.5% on average. There are three patients with relatively large changes, 15.4% on average, which indicate considerable differences in inspiration between scans. For the other 17 patients, an average change of only 4.9% was measured. The performance of our method was evaluated from the aspects of visual inspection, accuracy and processing time.
Fig. 11. The lung volumes from the initial and follow-up scans of the 20 patients.
Fig. 12 shows the results of lung boundary segmentation for subjects 1, 2 and 3. The second and third column shows segmented lung regions in 2D while the fourth column shows the shaded-surface display of the corresponding lungs. These results show that the proposed method accurately segments lung boundaries with large curvature. Fig. 13 shows coronal views of the lungs for subject 1, 2 and 4 with pulmonary nodules.
630
H. Hong et al. / Computers in Biology and Medicine 38 (2008) 623 – 634
Fig. 12. The results of lung boundary segmentation in subjects 1, 2 and 3. (a) Transverse CT images of the chest, (b) segmented lung boundary, (c) the transverse CT image enclosed by the segmented lung boundary and (d) shaded-surface display of the segmented lungs.
Fig. 13. The results of lung nodule matching in subjects 1, 2, and 4: (a) initial misalignment, (b) after optimal cube registration as initial registration, (c) after hierarchical surface registration and (d) after nodule registration.
H. Hong et al. / Computers in Biology and Medicine 38 (2008) 623 – 634
631
Fig. 14. Evaluation of the accuracy of corresponding nodules using the AED error in each subject.
Fig. 15. The reduction of AADD error and std. of ADD error per iteration in the hierarchical surface registration for the 20 patients at fine resolution: (a) AADD error and (b) standard deviation of ADD error.
632
H. Hong et al. / Computers in Biology and Medicine 38 (2008) 623 – 634
Fig. 16. Evaluation of the accuracy of corresponding lung boundaries using AADD error in each subject.
Fig. 17. The comparison of total processing time (a) without multi-resolution techniques and (b) with multi-resolution techniques.
Fig. 13(a) shows the initial misalignment which has translational, rotational and distorted differences. Fig. 13(b) shows the impact of our initial alignment. The positional differences
between corresponding lung nodules in the initial and followup CT scans shown in Fig. 13(a) are much reduced by the optimal cube registration, as shown in Fig. 13(b). This initial
H. Hong et al. / Computers in Biology and Medicine 38 (2008) 623 – 634
alignment is further refined by hierarchical surface registration until the lung boundaries of the initial and follow-up CT scans are aligned exactly, as in Fig. 13(c). From the exact matching of lung boundaries, lung nodule correspondences for each subject with nodules in the follow-up CT scans (dark gray) and nodules transformed from initial CT scans (light gray) into the follow-up CT scans are shown. Nodule correspondences for 86 (100%) out of 86 nodules were correctly established in the 20 patients. In Fig. 14, the results of the lung nodule alignment in the 20 patients are reported on a per-center-of-mass point basis using the average Euclidean distance (AED) between corresponding nodules in the initial and follow-up CT scans before and after alignment. Our initial registration reduces the initial AED error by 57.5% on average. The most substantial AED reduction in error in the initial misalignment is obtained with the initial registration and subsequent surface registration; 81.7% on average. Since the corresponding lung boundaries are almost aligned, the corresponding nodules are at or near the optimal position. After nodule registration, the average AED error of the initial misalignment is reduced from 30.1 to 1.8 mm. Fig. 15(a) shows how the error, the average of the AADD, is reduced by the hierarchical surface registration. In addition, Fig. 15(b) shows how the standard deviation of the absolute distance difference is reduced by our method. Since the positional difference is almost aligned by the initial registration, subsequent hierarchical surface registration rapidly converges to the optimal position. Following this, the optimal position for all subjects in the hierarchical registration is determined within 32 iterations on average. The average AADD error of the 20 patients is reduced by 1.4 voxels. To compare the accuracy of our method with conventional methods using a distance map, we evaluated the AADD error of the corresponding lung boundaries of the 20 patients using the hierarchical surface registration process shown in Fig. 16. The average AED error of the method based on Euclidean distance (Method 1), the method based on chamfer matching (Method 2) and our method (Method 3) were 1.58, 1.45 and 1.43 voxels, respectively. In particular, in subjects 1–5 with large data sets, the average reduction rate of the AADD error of Method 3 was relatively greater at 15.7% (Method 1) and 2.6% (Method 2). The total processing time is summarized in Fig. 17, where the execution time is measured for the initial registration, surface registration and nodule registration in comparison with Method 1 and Method 2. All our evaluations were performed on an Intel Pentium IV PC with a 3.0 GHz CPU and 2.0 GB of main memory. Due to optimal cube registration as an initial alignment, our distance map generation time is much faster than that of Method 1 and Method 2. As shown in Fig. 17(a), the average processing times of Methods 1, 2 and 3 without multiresolution techniques were 145.1, 95.5 and 74.8 s, respectively. On the other hand, the average processing times with multiresolution technique were 97.8, 49.1 and 37.1 s, respectively. With multi-resolution techniques, the total processing time is reduced by 56% on average. In particular, our method is 2.6 times faster than Method 1 and 1.3 times faster than Method 2 when
633
applying multi-resolution techniques. In addition, the processing time is much reduced when our method is applied to a large data set as in subjects 1–5. 4. Summary and conclusions In this paper, we have described an accurate and fast method for matching lung nodules in sequential chest CT scans. Our lung boundary segmentation extracts accurate lung boundaries through three steps: the lung extraction, airway extraction and subtraction steps. The 2D iRG can segment the thorax from the background air and separate the lungs from the thorax without holes such as unwanted interior cavities. The connected component labeling at low resolution can eliminate non-pulmonary structures such as bowel gas and reduce memory and processing time. Branch-based 3D region growing can successfully segment the trachea and the left and right mainstem bronchi using adaptive growing conditions. Our three-stage matching algorithm can find nodule correspondences efficiently. Using optimal cube registration, initial gross correction can be done rapidly and effectively without detecting anatomical landmarks. In the hierarchical surface registration, distance measure, using a 3D distance map generated by narrow-band distance propagation, allows rapid and robust convergence to the optimal position. For the refinement of the nodule correspondences, paring with the smallest Euclidean distances is performed. Experimental results show that our segmentation method gives accurate lung boundaries. In particular, lungs with large curvature or complicated shapes are accurately extracted. Our registration method was able to correctly find nodule correspondences in 20 patients. After nodule registration, the average AED error of initial misalignment is reduced from 30.1 mm down to 1.8 mm. In particular, in subjects obtained from thinsection scanning, the reduction rate of the AADD error in our method is relatively greater at 15.7% and 2.6% compared to methods based on Euclidean distance and chamfer matching. In addition, our method is 2.6 times faster than Euclidean distancebased registration and 1.3 times faster than chamfer matchingbased registration when applying multiresolution techniques. Acknowledgments This work was supported in part by Grant no. 10888 from the Seoul Research and Business Development Program and in part by Grant no. R01-2006-000-11244-0 from the Basic Research Program of the Korea Science and Engineering Foundation. The authors are grateful to Dr. Jin Mo Goo at the Seoul National University Hospital of Seoul, Korea, for providing the chest CT data sets shown in this paper and giving advice valuable to our research. References [1] E.R. Heitzman, The role of computed tomography in the diagnosis and management of lung cancer: an overview, Chest 89 (4) (1986) 237S–241S. [2] M.A. Howe, B.H. Gross, CT evaluation of the equivocal pulmonary nodule, Comput. Radiol. 11 (1987) 61–67.
634
H. Hong et al. / Computers in Biology and Medicine 38 (2008) 623 – 634
[3] M. Remy-Jardin, J. Remy, F. Giraud, C.H. Marquette, Pulmonary nodules: detection with thick-section spiral CT versus conventional CT, Radiology 187 (1993) 513–520. [4] F. Fischbach, F. Knollmann, V. Griesshaber, T. Frund, E. Akkol, R. Felix, Detection of pulmonary nodules by multiscale computed tomography: improved detection rate with reduced slice thickness, Eur. Radiol. 13 (2003) 2378–2383. [5] J.P. Ko, H. Rusinek, E.L. Jacobs, M. Betke, G. McGuinness, Small pulmonary nodules: volume measurement at chest CT—phantom study, Radiology 228 (2003) 864–870. [6] C.I. Henschke, D.I. McCauley, D.F. Yankelevitz, et al., Early lung cancer action project: overall design and findings from baseline screening, Lancet 354 (1999) 99–105. [7] K.T. Bae, J.S. Kim, Y.H. Na, K.G. Kim, J.H. Kim, Pulmonary nodules: automated detection on CT images with morphologic matching algorithm—preliminary results, Radiology 236 (1) (2005) 286–294. [8] M.F. McNitt-Gray, E.M. Hart, N. Wyckoff, J.W. Sayre, J.G. Goldin, D.R. Aberle, A pattern classification approach to characterizing solitary pulmonary nodules imaged on high resolution CT: preliminary results, Med. Phys. 26 (1999) 880–888. [9] T.W. Way, L.M. Hadjiiski, B. Sahiner, H.P. Chan, P.N. Cascade, E.A. Kazerooni, N. Bogot, C. Zhou, Computer-aided diagnosis of pulmonary nodules on CT scans: segmentation and classification using 3D active contours, Med. Phys. 33 (2006) 2323–2337. [10] D.F. Yankelevitz, A.P. Reeves, W.J. Kostis, B. Zhao, C.I. Henschke, Small pulmonary nodules: volumetrically determined growth rates based on CT evaluation, Radiology 217 (2000) 251–256. [11] J.P. Ko, M. Betke, Chest CT: automated nodule detection and assessment of change over time—preliminary experience, Radiology 218 (2001) 267–273. [12] M. Betke, H. Hong, J.P. Ko, Automatic 3D registration of lung surfaces in computed tomography scans, in: Proceedings of the Medical Image Computing and Computer-Assisted Intervention (MICCAI), 2001, pp. 725–733. [13] M. Betke, H. Hong, J.P. Ko, Automatic 3D registration of lung surfaces in computed tomography scans, CS Technical Report 2001-004, Boston University, 2001. [14] H. Hong, M. Betke, S. Teng, Multilevel 3D registration of lung surfaces in computed tomography scans—preliminary experience, in: Proceedings of the International Conference on Diagnostic Imaging and Analysis (ICDIA), 2002, pp. 90–95. [15] W. Mullaly, M. Betke, H. Hong, J. Wang, K. Mann, J.P. Ko, Multicriterion 3D segmentation and registration of pulmonary nodules on CT: a preliminary investigation, in: Proceedings of the International Conference on Diagnostic Imaging and Analysis (ICDIA), 2002, pp. 176–181. [16] M.N. Gurcan, R.C. Hardie, S.K. Rogers, D.E. Dozer, B.H. Allen, J.W. Hoffmeister, Automated global matching of temporal thoracic helical CT studies: feasibility study, Int. Congr. Ser. 1256 (2003) 1031–1036. [17] T. Kusanagi, Y. Mekada, J. Hasegawa, J. Toriwaki, M. Mori, H. Natori, Correspondence of lung nodules in sequential chest CT images for quantification of the curative effect, Int. Congr. Ser. 1256 (2003) 983–989. [18] M. Kubo, T. Yamamoto, Y. Kawata, N. Niki, K. Eguchi, H. Ohmatsu, R. Kakinuma, M. Kaneko, M. Kusumoto, N. Moriyama, K. Mori, H. Nishiyama, CAD system for the assistance of comparative reading for lung cancer using retrospective helical CT images, Proc. SPIE Med. Imaging 4322 (2001) 1788–1795. [19] S. Hu, E.A. Hoffman, J.M. Reinhardt, Accurate lung segmentation for accurate quantitation of volumetric X-ray CT images, IEEE Trans. Med. Imaging 20 (6) (2001) 490–498. [20] M.S. Brown, M.F. McNitt-Gray, N.J. Mankovich, J.G. Goldin, J. Hiller, L.S. Wilson, D.R. Aberle, Method for segmenting chest CT image data using an anatomical model: preliminary results, IEEE Trans. Med. Imaging 16 (6) (1997) 828–839. [21] M. Wu, J. Chang, A.A. Chiang, J. Lu, H. Hsu, W. Hsu, C. Yang, Use of quantitative CT to predict postoperative lung function in patients with lung cancer, Radiology 191 (1994) 257–262.
[22] S.G. Armato III, H. MacMahon, Automated lung segmentation and computer-aided diagnosis for thoracic CT scans, Int. Congr. Ser. 1256 (2003) 977–982. [23] S.G. Armato III, W.F. Sensakovic, Automated lung segmentation for thoracic CT: impact on computer-aided diagnosis, Acad. Radiol. 11 (2004) 1011–1021. [24] J.K. Leader, B. Zheng, R.M. Rogers, F.C. Sciurba, A. Perez, B.E. Chapman, S. Patel, C.R. Fuhrman, D. Gur, Automated lung segmentation in X-ray computed tomography, Acad. Radiol. 20 (2003) 1224–1236. [25] M. Sonka, V. Hlavac, R. Boyle, Image Processing, Analysis and Machine Vision, Pacific Grove, PWS, CA, 1999. [26] R.C. Gonzalez, R.E. Woods, Digital Image Processing, second ed., Prentice-Hall, Englewood Cliffs, NJ. [27] G. Borgefors, Hierarchical chamfer matching: a parametric edge matching algorithm, IEEE Trans. Pattern Recognition Mach. Intell. 10 (6) (1988) 849–865. [28] O. Cuisenaire, Distance transformations: fast algorithms and applications to medical image processing, Ph.D. Thesis, 1999. [29] B. Fang, W. Hsu, M.L. Lee, Techniques for temporal registration of retinal images, in: Proceedings of the International Conference on Image Processing (ICIP), 2004, pp. 1089–1092. [30] G. Borgefors, Distance transformations in digital images, Comput. Vision Graphics Image Processing 34 (3) (1986) 344–371. [31] T. Saito, J.-I. Toriwaki, New algorithms for Euclidean distance transformation of an n-dimensional digitized picture with applications, Pattern Recognition 27 (11) (1994) 1551–1565. [32] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical recipes in C—the art of scientific computing, Cambridge University Press, Cambridge, 1992. [33] G.K. Matsopoulos, N.A. Mouravliansky, K.K. Delibasis, K.S. Nikita, Automatic retinal image registration scheme using global optimization techniques, IEEE Trans. Inf. Technol. Biomed. 3 (1) (1999) 47–60. [34] N. Mouravliansky, G.K. Matsopoulos, K. Delibasis, K.S. Nikita, Automatic retinal registration using global optimization techniques, in: Proceedings of the 20th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, vol. 20(2), 1998, pp. 567–570. [35] J.L. Bernon, V. Boudousq, J.F. Rohmer, M. Fourcade, M. Zanca, M. Rossi, D. Mariano-Goulart, A comparative study of Powell’s and Downhill Simplex algorithms for a fast multimodal surface matching in brain imaging, Comput. Med. Imaging Graphics 25 (2001) 287–297. [36] F. Maes, A. Collignon, D. Vandermeulen, G. Marchal, P. Suetens, Multimodality image registration by maximization of mutual information, IEEE Trans. Med. Imaging 16 (2) (1997) 187–198.
Helen Hong received BS degree (1994), MS degree (1996) and PhD degree (2001) in Computer Engineering at the Ewha Womans University, Seoul, Korea. From 2001 to 2003, she was at the Seoul National University as a post doctor. In 2003, she was at the INFINITT Co. Ltd. as a principal researcher. From 2003 to 2005, she was at the Seoul National University as a research professor. Since 2006, she is a professor in the Division of Multimedia Engineering at the Seoul Women’s University, Seoul, Korea. Her current research interests include image registration, medical image processing and analysis, volume visualization and computer graphics. Jeongjin Lee received BS degree (2000) in Department of Mechanical Engineering and MS degree (2002), PhD degree (2008) in School of Computer Science and Engineering from Seoul National University, Korea. Since 2007, he is a research professor in Department of Radiology University of Ulsan, College of Medicine, Korea. His current research interests include medical image registration, medical image segmentation, virtual colonoscopy, deformable modeling, and computer animation. Yeny Yim received BS degree (2001) in Computer Engineering at the Kyungpook National University, Daegu, Korea. Now, she is a PhD candidate of School of Computer Science and Engineering at the Seoul National University, Seoul, Korea. Her current research interests include medical image segmentation and registration.