EXACT'09 -285-

Maximal Contrast Adaptive Region Growing

for CT Airway Tree Segmentation

Carlos S. Mendoza  , Begoña 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 condi-
tions. 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 grow-

ing, 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 mea-
surements, 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ón Reina Mercedes” from ”Hospital Univer-
sitario Virgen del Rocı́o” (Sevilla), and ”Consejerı́a de Salud de la Junta de An-
dalucı́a”. Carlos S. Mendoza was supported by a doctoral scholarship financed by
Universidad de Sevilla.
-286- EXACT'09

chest area is very tedious and too cumbersome for everyday clinical use. More-
over, manual analysis performed by the radiologist was only qualitative esti-
mation 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 inhomo-
geneous 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 afor-
mentioned 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 deter-
mine optimal thresholds by detection of sudden volume increase in the segmented
region. The difficulty lies in defining how much volume increase is considered ex-
plosion, 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.
EXACT'09 -287-

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 com-
puted iteratively, as intensity-decreasing pixels were added to the segmented re-
gion. 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]. Un-
fortunately, the exhaustivity of their per-pixel approach entailed very low com-
putational 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 com-
putational 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 al-
gorithm extracts a large fraction of the visible airways with few false positive
branches. For some applications, however, even this accuracy rate can be insuf-
ficient. 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 in-
tensities 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 ap-
proximately 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
-288- EXACT'09

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,
N is a cubic neighborhood of radius R around the seed, x is a voxel position,
f (x) is the intensity for voxel at x, f¯N is the mean intensity estimate in N and
|N | is the cardinality of N . Moreover, σfN is the estimated standard deviation
for intensities in N , K is a constant parameter, and f (x), f  (x) are the input
and output intensities for the non-linear mapping described below.
1 
f¯N = f (xk ) , (1)
|N |
xk ∈N


1   2
σfN = f (xk ) − f¯N , (2)
|N |
xk ∈N
⎛ ⎛ ⎞⎞−1
f (x) − ¯N
f
f  (x) = ⎝1 + exp ⎝−  Kσ ⎠⎠ . (3)
fN
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 es-
timated mean as in (3). The width of the sigmoidal window extends Kσ around
the center f¯N of 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 K ensures robust mapping for
the estimated distribution (that of the tissues of interest). The sigmoidal map-
ping 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 me-
dian 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
EXACT'09 -289-

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 R0 defined by several seeds along the upper tra-
chea, 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 , σfR )
i−1
3. For every candidate voxel xci−1 being 26-connected to Ri−1 , xci−1 ∈ Ri if
 
f  xci−1 ∈ f¯ Ri−1 ± ki σfR (4)
i−1

 
4. Compute the assessment function Oi f¯ Ri , f¯ Pi using the intensity average
f¯ Ri in Ri and the intensity average f¯ Pi in the external perimeter Pi of Ri
according to (1) and the following eqs.:

Pi = {xci } ∩ Ri C , (5)
 
   f¯ − f¯ Ri 
Oi f¯ Ri , f¯ Pi =  ¯ Pi  (6)
f Pi + f¯ Ri 
5. If Oi−1 was a local maximum, when compared to Oi−2 and 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 k0 and Δk are critical for the performance
of the algorithm. k0 affects computational efficiency
 requiring a greater number
of iterations before a local maximum of O f¯Ri , f¯Pi is actually found. There-
fore, 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
-290- EXACT'09

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, k0 = 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 seg-
mentation 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 sub-
mitted 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.
EXACT'09 -291-

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 mea-
surements. 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-
-292- EXACT'09

(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

EXACT'09 -293-

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 im-
age 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 compar-
ison 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 pro-
vide 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 ob-
structive 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 label-
ing 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 Pat-
tern 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
-294- EXACT'09

8. Singh, H., Crawford, M., Curtin, J.P., Zwiggelaar, R.: Automated 3D segmentation
of the lung airway tree using gain-based region growing approach. MICCAI 2
(2004) 975–982
9. Fetita, C.I., Prteux, F., Beigelman-Aubry, C., Grenier, P.: Pulmonary airways: 3-D
reconstruction from multislice CT and clinical investigation. IEEE Transactions
on Medical Imaging 23(11) (2004) 1353–1364
10. Graham, M.W., Gibbs, J.D., Higgins, W.E.: Robust system for human airway-tree
segmentation. In: Progress in Biomedical Optics and Imaging - Proceedings of
SPIE. Volume 6914. (2008)
11. Sonka, M., Park, W., Huffman, E.A.: Rule-based detection of intrathoracic airway
trees. IEEE Transactions on Medical Imaging 15(3) (1996) 314–326
12. Aykac, D., Huffman, E.A., McLennan, G., Reinhardt, J.M.: Segmentation and
analysis of the human airway tree from three-dimensional X-ray CT images. IEEE
Transactions on Medical Imaging 22(8) (2003) 940–950
13. Tschirren, J., Huffman, E.A., McLennan, G., Sonka, M.: Intrathoracic airway trees:
Segmentation and airway morphology analysis from low-dose CT scans. IEEE
Transactions on Medical Imaging 24(12) (2005) 1529–1539
14. Schlathölter, T., Lorenz, C., Carlsen, I.C., Renisch, S., Deschamps, T.: Simultane-
ous segmentation and tree reconstruction of the airways for virtual bronchoscopy.
In: Proceedings of SPIE - The International Society for Optical Engineering. Vol-
ume 4684 I. (2002) 103–113
15. Ginneken, B.V., Baggerman, W., Rikxoort, E.M.V.: Robust segmentation and
anatomical labeling of the airway tree from thoracic CT scans. Lecture Notes in
Computer Science (including subseries Lecture Notes in Artificial Intelligence and
Lecture Notes in Bioinformatics) 5241 LNCS(PART 1) (2008) 219–226
16. Fabijanska, A.: Two-pass region growing algorithm for segmenting airway tree from
MDCT chest scans. Computerized Medical Imaging and Graphics (2009) Article
in Press.
17. Lo, P., Bruijne, M.D.: Voxel classification based airway tree segmentation. In:
Progress in Biomedical Optics and Imaging - Proceedings of SPIE. Volume 6914.
(2008)
18. Adams, R., Bischof, L.: Seeded region growing. IEEE Transactions on Pattern
Analysis and Machine Intelligence 16(6) (1994) 641–647
19. Hojjatoleslami, S.A., Kittler, J.: Region growing: A new approach. IEEE Trans-
actions on Image Processing 7(7) (1998) 1079–1084
20. Haralick, R.M., Shapiro, L.G.: Image segmentation techniques. Computer Vision,
Graphics, & Image Processing 29(1) (1985) 100–132
21. Revol-Muller, C., Peyrin, F., Carrillon, Y., Odet, C.: Automated 3D region growing
algorithm based on an assessment function. Pattern Recognition Letters 23(1-3)
(2002) 137–150
22. Hu, S., Hoffman, E.A., Reinhardt, J.M.: Automatic lung segmentation for accu-
rate quantitation of volumetric X-ray CT images. IEEE Transactions on Medical
Imaging 20(6) (2001) 490–498
23. Sluimer, I., Prokop, M., Ginneken, B.V.: Toward automated segmentation of the
pathological lung in CT. IEEE Transactions on Medical Imaging 24(8) (2005)
1025–1038
24. Yoo, T.S., Ackerman, M.J., Lorensen, W.E., Schroeder, W., Chalana, V., Aylward,
S., Metaxas, D., Whitaker, R.: Engineering and algorithm design for an image
processing API: a technical report on ITK–the Insight Toolkit. Studies in health
technology and informatics 85 (2002) 586–592
EXACT'09 -295-

25. Pieper, S., Lorensen, B., Schroeder, W., Kikinis, R.: The NA-MIC Kit: ITK, VTK,
pipelines, grids and 3D Slicer as an open platform for the medical image computing
community. In: 2006 3rd IEEE International Symposium on Biomedical Imaging:
From Nano to Macro - Proceedings. Volume 2006. (2006) 698–701