Maximal Contrast Adaptive Region Growing for CT Airway Tree Segmentation

Maximal Contrast Adaptive Region Growing for CT Airway Tree Segmentation Carlos S. Mendoza , Bego˜na Acha, and Carmen Serrano Universidad de Sevilla Av. de los Descubrimientos s/n 41092 Sevilla, Spain {csanchez1,bacha,cserrano}@us.es Abstract. In this paper we propose a fully self-assessed adaptive region growing airway segmentation algorithm. We rely on a standardized and self-assessed region-based approach to deal with varying imaging conditions. Initialization of the algorithm requires prior knowledge of trachea location. This can be provided either by manual seeding or by automatic trachea detection in upper airway tree image slices. The detection of the optimal parameters is managed internally using a measure of the varying contrast of the growing region. Extensive validation is provided for a set of 20 chest CT scans. Our method exhibits very low leakage into the lung parenchyma, so even though the smaller airways are not obtained from the region growing, our fully automatic technique can provide robust and accurate initialization for other methods. Key words: CT, airways, segmentation, multi-tolerance, adaptive region growing, median filter 1 Introduction CT scans are frequently used for pulmonary disorder assessment [1]. Pathologies that could affect sufficient lung function include tumors, pulmonary embolism, atelectasis, pneumonia, emphysema, asthma, bronchiectasis, and many others. Certain lung diseases can be diagnosed based on airway wall thickness measurements, diameter, branching geometry and rate of tapering. CT is currently the only readily accessible, relatively noninvasive technique that is capable of providing airway tree quantitative structural data in vivo [2]. Traditionally, analysis of CT chest scans was performed manually by skilled radiologist who recognized areas of abnormal airway properties in consecutive slices of the examined scan. However, analyzing about 400 slices covering the This work was supported by ”Fundaci´on Reina Mercedes” from ”Hospital Universitario Virgen del Roc´ıo” (Sevilla), and ”Consejer´ıa de Salud de la Junta de Andaluc´ıa”. Carlos S. Mendoza was supported by a doctoral scholarship financed by Universidad de Sevilla. EXACT'09 -285- chest area is very tedious and too cumbersome for everyday clinical use. Moreover, manual analysis performed by the radiologist was only qualitative estimation of airway abnormalities without accurate quantification of pathological changes. Even though airway tree abnormalities can be detected based on 2D slices, the ability to extract a full 3D model of the airway tree from a 3D image has several key advantages. For example, slice based measurements can be inaccurate if the airway is not perpendicular to the slice. Also, information available from the slices is deprived of useful context, making it harder for a radiologist to keep track of the generation number of an airway, the structure of nearby airways or the overall shape of a segment [3]. There has been a number of efforts to try to delineate the airway tree in chest CT scans. Airway tree segmentation is a complex task, mainly due to inhomogeneous grey level of the voxels located inside the bronchial lumen, artefacts caused by blood vessels adjacent to airway walls and changes of intensity levels along airway walls. Many airway segmentation methods use region-growing algorithms, which attempt to separate air and soft tissue voxels using an HU (Hounsfield Unit) threshold [4–8]. Region growing is fast and assumes no prior knowledge of the shape or size of the airways. Choosing an appropriate global HU threshold is difficult, however, as the lungs are filled with air and misclassifying a single wall voxel can allow the segmentation to leak into the lung parenchyma. Other methods make use of grey-level morphological operators [9–12], or front wave propagation schemes [13–15], to impose structural properties derived from a priori anatomical knowledge. Region growing is the preferred method for initializing several of the aformentioned algorithms, being a fast and intuitive technique [3,10,11]. Some recent segmentation techniques perform region growing stages on images derived from the intensity of the scan, like the morphological gradient [16], or the posterior of a classification stage [17]. Thus, fully automatic region growing approaches become mandatory, specially those that can provide as much information as possible for the following algorithmic stage, with very low leakeage. Few fully automatic region growing methods have proven successfull in these tasks over a significant variety of imaging devices. One of the most common techniques for threshold determination was proposed by Mori et al. [6]. Their approach, widely known as explosion-controlled region growing, tries to determine optimal thresholds by detection of sudden volume increase in the segmented region. The difficulty lies in defining how much volume increase is considered explosion, as compared to normal region growth. In our work, we also propose a self-assessed region growing algorithm. Our assessment strategy is founded on a previous normalization stage, and makes use of an iteratively computed contrast measure. In our method there is no need to define any confidence margins, as the segmented region is determined precisely by the evolution of the proposed contrast measure. -286- EXACT'09 Apart from manual trial-and-error adaptive threshold selection [18], some self-assessed adaptive region growing strategies have been proposed in the past, in a variety of application contexts. In their work [19], Hojjatoleslami and Kittler proposed a method based on finding the global maxima for two different contrast measures which they computed iteratively, as intensity-decreasing pixels were added to the segmented region. The success of the assessment was founded on the assumption that maximal contrast occurred on region boundaries, which is a reformulation of approaches assuming that the variation of the gray values within regions is smaller than across regions, an inherent assumption in all region growing techniques [20]. Unfortunately, the exhaustivity of their per-pixel approach entailed very low computational efficiency, aggravated by their multiple complex peripheral measures. Revol-Muller et al. [21] used morphological measures to assess the multiplier of the adaptive range in region growing. Instead of computing their assessment function for every pixel addition to the region, they sampled the function for an evenly-spaced set of values. In our method we propose an assessment function based on a simple measure of the evolving contrast for the region growing sequence. To make this approach computationally feasible in 3D, we produce only evenly-spaced samples of this function, along the values of the assessed parameter defined on the normalized dynamic range of the image. This sampling strategy dramatically reduces computational complexity while preserving most critical values. Globally, our goal is to provide fast automatic segmentation based only on the location of the upper part of the trachea. Our automatic segmentation algorithm extracts a large fraction of the visible airways with few false positive branches. For some applications, however, even this accuracy rate can be insufficient. Our method can then be considered as an initialization stage for other refining strategies, in a wide variety of imaging devices, and with no significant leakage, as suggested by our experimental results. 2 Method 2.1 Airway Lumen Intensity Model Since our goal is providing a mechanism for airway segmentation with minimal user intervention, we have established a model that takes into account their intensity distribution in CT images. We model then our object of interest as a connected region whose pixel intensities are sampled from a Gaussian distribution with unknown mean and standard deviation. We know that our region of interest is surrounded by other tissues derived from other, sometimes adjacent, intensity distributions. Although common in the literature, this assumption for the intensities is rarely met in practice, in the sense that the intensity distributions of tissues are only approximately Gaussian, as can be inferred from direct observation of histograms. Besides, partial overlap between adjacent distributions often occurs. For dealing with these inconveniences we propose the use of an assessment function, that EXACT'09 -287- is to be evaluated along a sequence of region growing stages (region growing sequence). 2.2 Segmentation Algorithm Normalization and Denoising Since our method was conceived for images from a wide range of scanners and acquisition protocols, we have developed a normalizing stage that accounts for such variability. As we will introduce later on, for the self-assessed region growing stage of the algorithm, we require the input intensity dynamic range to be normalized with respect to some parameter estimates of the objective intensity distribution. In the following equations in which we describe the normalization process, Nis a cubic neighborhood of radius Raround the seed, xis a voxel position, f(x) is the intensity for voxel at x,¯ fNis the mean intensity estimate in Nand |N|is the cardinality of N. Moreover, σfNis the estimated standard deviation for intensities in N,Kis a constant parameter, and f(x), f(x) are the input and output intensities for the non-linear mapping described below. ¯ fN=1 |N| xk∈N f(xk),(1) σfN=1 |N| xk∈Nf(xk)−¯ fN2,(2) f(x)=⎛ ⎝1 + exp ⎛ ⎝−f(x)−¯ fN KσfN 3⎞ ⎠⎞ ⎠ −1 .(3) In a first step we proceed by Gaussian distribution maximum-likelihood (ML) estimation of the mean and standard deviation as in (1-2), and then perform a non-linear normalization using a sigmoidal transfer function centered on the estimated mean as in (3). The width of the sigmoidal window extends Kσ around the center ¯ fNof the mapping. For K= 3 the width of the window would be enough to map 99.7% of the samples, of a Gaussian distribution with similar mean and standard deviation. Greater values of Kensures robust mapping for the estimated distribution (that of the tissues of interest). The sigmoidal mapping has been chosen because of its smoothness, and its ability to focus the output dynamic range on a given input intensity range of interest. Finally, we perform non-linear denoising using an in-slice bidimensional median filter with kernel radius Γ. Other denoising schemes would be valid, always keeping in mind that stronger smoothing usually involves loss of smaller airways. Self-Assessed Region Growing Departing from a normalized and filtered version of the image under study, whose intensities lie in the range [0,1], we apply our self-assessed contrast-maximizing algorithm. The initial region that -288- EXACT'09 needs to be provided (the upper trachea) can be obtained by means of manual seeding, or by automatic detection using one of the many methods available in the literature [22, 23]. For this initial implementation we decided to simply provide the seeds manually, as described below. Considering an initial region R0defined by several seeds along the upper trachea, we provide the analytical description of the i-th iteration of the algorithm: 1. Update multiplier ki=k0+iΔk 2. Compute, in last iteration grown region Ri−1, ML estimates for the mean (available from last iteration) and standard deviation ( ¯ fRi−1,σ f Ri−1) 3. For every candidate voxel xci−1being 26-connected to Ri−1,xci−1∈Riif fxci−1∈¯ fRi−1±kiσf Ri−1(4) 4. Compute the assessment function Oi¯ fRi,¯ fPiusing the intensity average ¯ fRiin Riand the intensity average ¯ fPiin the external perimeter Piof Ri according to (1) and the following eqs.: Pi={xci}∩RiC,(5) Oi¯ fRi,¯ fPi= ¯ fPi−¯ fRi ¯ fPi+¯ fRi (6) 5. If Oi−1was a local maximum, when compared to Oi−2and Oi(only when i≥2), then the algorithm stops and the output is Ri−1. Otherwise another iteration takes place Of all aforementioned parameters only k0and Δk are critical for the performance of the algorithm. k0affects computational efficiency requiring a greater number of iterations before a local maximum of O¯ fRi,¯ fPiis actually found. Therefore, its fine tuning for a specific scanner, could save some computational time. From observation of the region growing sequence, we conclude that these first iterations are typically very fast, so the improvement is frequently negligible. In what concerns Δk, the choice must guarantee that the assessment function is being sampled adequately in order to detect its local variations. Since the estimates for the mean and standard deviation are continually updated as the region grows, the estimates become increasingly close to the theoretical values. We argue that setting Δk below one tenth of 3 (which is the theoretical value multiplying the standard deviation of a Gaussian distribution for 99.7% of its samples to be included in a range of that width around the mean) is enough for the segmentation process to be able not to miss the available local maxima of the assessment function. This claim is supported by our experimental results. 3 Results We have implemented our algorithm using open source medical image processing libraries, more precisely the Insight Toolkit [24] for algorithm development, and EXACT'09 -289- the command line executable module infrastructure provided by 3DSlicer for fast prototyping, calibration, evaluation, and manual segmentation on real images for further validation [25]. The algorithm that we will validate, uses the following parameter values: R=2,K =12,Γ =1,k 0= 1 and Δk =0.1. The values for these parameters were determined from other non-thoracic CT images, and are intended to suit any imaging conditions. No thoracic CT scans were use for the tuning of these parameters, as the algorithm was initially conceived for general purpose segmentation. The algorithm was implemented and executed in a 2 GHz Intel Core 2 Duo Windows PC with 2 GB RAM, and the average running time was 129 ±27 s. For generation of the presented airway segmentations, the initial region was provided using three manually selected seeds along the upper trachea. The segmentation process took between 1 and 2 minutes for each dataset, including the reading/writing of the images. The evaluation of this algorithm has taken place in comparison with other 14 algorithms using a set of 40 CT scans, with different acquisition parameters. First 20 datasets were used for training (unnecessary in our approach), and last 20 for testing. The produced segmentations were centrally evaluated by a team of trained observers. The objective was to compare performance. For this purpose, a ground truth was constructed from all available segmentations from the different algorithms, and all results were subsequently evaluated with respect to this ground truth. Evaluation of each individual segmentation is performed in the following steps: 1. Airway segments or branches were extracted from submitted airway tree segmentations using a fast marching based algorithm [14]. 2. Airway segments were evaluated visually on a set of extracted slices from both a reoriented view and a reformatted view with straightened airway centerlines. 3. Each segment was scored as ”correct” or ”wrong”, by at least two observers. The criterion used is whether the extracted airway segment indeed belongs to the airway tree; the exact airway shape and dimensions are not taken into account. The ground truth is then defined as the union of all valid airway segments from all submitted segmentations. The following measurements are computed and used for comparing the submitted results: 1. Branch count: The number of branches that were detected correctly. A branch is considered detected as long as the length of the centerlines is more than 1 mm. 2. Branch detected: The fraction of branches that were detected, with respect to the branches present in the ground truth. 3. Tree length: The sum of the length of the centerlines of all correctly detected branches. -290- EXACT'09 Table 1. Evaluation measures for the twenty cases in the test set. Branch Branch Tree Tree length Leakage Leakage False count detected length detected count volume positive (%) (cm) (%) (mm3) rate (%) CASE21 89 44.7 46.0 41.6 0 0.0 0.00 CASE22 54 14.0 37.0 11.2 8 2085.7 30.67 CASE23 33 11.6 27.3 10.5 0 0.0 0.00 CASE24 49 26.3 43.2 26.6 0 0.0 0.00 CASE25 83 35.5 63.5 25.2 0 0.0 0.00 CASE26 22 27.5 15.8 24.0 0 0.0 0.00 CASE27 35 34.7 26.0 32.1 0 0.0 0.00 CASE28 56 45.5 40.5 37.0 0 0.0 0.00 CASE29 74 40.2 44.4 32.2 0 0.0 0.00 CASE30 44 22.6 30.2 19.8 0 0.0 0.00 CASE31 77 36.0 53.7 30.6 3 31.1 0.35 CASE32 80 34.3 62.4 28.6 2 314.3 2.81 CASE33 83 49.4 56.9 38.7 0 0.0 0.00 CASE34 266 58.1 189.8 53.1 2 39.3 0.18 CASE35 112 32.6 78.4 25.3 0 0.0 0.00 CASE36 59 16.2 54.1 13.1 0 0.0 0.00 CASE37 46 24.9 39.2 22.0 0 0.0 0.00 CASE38 35 35.7 26.9 40.5 0 0.0 0.00 CASE39 93 17.9 73.8 18.0 4 65.8 0.95 CASE40 40 10.3 30.9 8.0 0 0.0 0.00 Mean 71.5 30.9 52.0 26.9 0.9 126.8 1.75 Std. dev. 51.7 13.1 36.5 11.8 2.0 466.5 6.84 Min 22 10.3 15.8 8.0 0 0.0 0.00 1st quartile 40 17.9 30.2 18.0 0 0.0 0.00 Median 58 33.4 43.8 25.9 0 0.0 0.00 3rd quartile 89 44.7 63.5 38.7 2 39.3 0.35 Max 266 58.1 189.8 53.1 8 2085.7 30.67 4. Tree length detected: The fraction of tree length in the ground truth that was detected correctly. 5. Leakage count: The number of unconnected groups of ”correct” regions that were neighboring with a ”wrong” region. Indicates how easy/difficult it is to manually separate leakages from the correctly detected branches. 6. Leakage volume: The volume of regions that were wrongly detected. 7. False positive rate: The fraction of the volume of regions that were detected wrongly over the volume of all detected regions. Trachea was excluded from the branch length and branch count related measurements. For the leakage based measures, both trachea and main bronchi were excluded. In Table 1 we provide all the computed experimental measures. We also present 3D renderings for the 20 datasets used for validation, in Fig. 1. These results prove that our algorithm is able to produce reasonably complete segmen- EXACT'09 -291- (a)Case 21 (b)Case 22 (c)Case 23 (d)Case 24 (e)Case 25 (f)Case 26 (g)Case 27 (h)Case 28 (i)Case 29 (j)Case 30 (k)Case 31 (l)Case 32 (m)Case 33 (n)Case 34 (o)Case 35 (p)Case 36 (q)Case 37 (r)Case 38 (s)Case 39 (t)Case 40 Fig. 1. (a)-(t). Segmented airways for test cases 21 through 40 -292- EXACT'09 tations, with very limited leakage. In most of the available test sets, the maximal contrast condition has produced fully automatic segmentations which compare fairly with those obtained using more sophisticated methods. The results were obtained with the same parameters for all different acquisition varieties, in a reduced time frame. 4 Concluding Remarks and Future Work Our experimental results suggest that our modified region growing strategy, when used as an initialization stage, might benefit many airway segmentation algorithms currently available. Also, in those techniques which exploit region growing approaches applied over images obtained through processing of the image intensities, the use of our maximal-constrast stopping condition could be useful for further automatization. Our immediate future line of work will include further study of our novel region-growing stopping criterion. We intend to develop an exhaustive comparison between previously available explosion detection and our newer approach. Further along the way, we would like to develop a hybrid method, using our contrast-assessed region growing for initialization, and adding some refining stages using morphological operations and surface interpolation. This will provide for further comparison between the different initialization strategies. References 1. Sluimer, I., Schilham, A., Prokop, M., Ginneken, B.V.: Computer analysis of computed tomography scans of the lung: A survey. IEEE Transactions on Medical Imaging 25(4) (2006) 385–405 2. Coxson, H.O., Rogers, R.M.: Quantitative computed tomography of chronic obstructive pulmonary disease. Academic Radiology 12(11) (2005) 1457–1463 3. Szymczak, A., Vanderhyde, J.: Airway segmentation by topology-driven local thresholding. In: Progress in Biomedical Optics and Imaging - Proceedings of SPIE. Volume 6914. (2008) 4. Mori, K., Hasegawa, J.I., Suenaga, Y., Toriwaki, J.I.: Automated anatomical labeling of the bronchial branch and its application to the virtual bronchoscopy system. IEEE Transactions on Medical Imaging 19(2) (2000) 103–114 5. Summers, R.M., Feng, D.H., Holland, S.M., Sneller, M.C., Shelhamer, J.H.: Virtual bronchoscopy: Segmentation method for real-time display. Radiology 200(3) (1996) 857–862 6. Mori, K., Hasegawa, J., Toriwaki, J., Anno, H., Katada, K.: Recognition of bronchus in three-dimensional X-ray CT images with application to virtualized bronchoscopy system. Proceedings of the 13th International Conference on Pattern Recognition 3(1996) 528–532 7. Kiraly, A.P., Higgins, W.E., Hoffman, E.A., McLennan, G., Reinhardt, J.M.: 3D human airway segmentation for virtual bronchoscopy. In: Proceedings of SPIE - The International Society for Optical Engineering. Volume 4683. (2002) 16–29 EXACT'09 -293-