- Email: [email protected]
A fast and accurate automatic lung segmentation and volumetry method for MR data used in epidemiological studies
Computerized Medical Imaging and Graphics 36 (2012) 281–293
Contents lists available at SciVerse ScienceDirect
Computerized Medical Imaging and Graphics journal homepage: www.elsevier.com/locate/compmedimag
A fast and accurate automatic lung segmentation and volumetry method for MR data used in epidemiological studies Tetyana Ivanovska a,∗ , Katrin Hegenscheid b , Rene Laqua b , Jens-P. Kühn b , Sven Gläser c , Ralf Ewert c , Norbert Hosten b , Ralf Puls b , Henry Völzke a a
Institute of Community Medicine, Ernst-Moritz-Arndt University, Greifswald, Germany Institut of Diagnostic Radiology und Neuroradiology, Ernst-Moritz-Arndt University, Greifswald, Germany c Klinik und Poliklinik für Innere Medizin, Bereich Pneumologie und Infektiologie, Ernst-Moritz-Arndt University, Greifswald, Germany b
a r t i c l e
i n f o
Article history: Received 6 March 2011 Received in revised form 30 September 2011 Accepted 6 October 2011 Keywords: Segmentation MRI Pulmonary imaging Volumetry
a b s t r a c t In modern epidemiological population-based studies a huge amount of magnetic resonance imaging (MRI) data is analysed. This requires reliable automatic methods for organ extraction. In the current paper, we propose a fast and accurate automatic method for lung segmentation and volumetry. Our approach follows a “coarse-to-fine” segmentation strategy. First, we extract the lungs and trachea excluding the main pulmonary vessels. This step is executed very fast and allows for measuring the volume of both structures. Thereafter, we start a refinement procedure that consists of three main stages: trachea extraction, lung separation, and filling the cavities on the final lung masks. After the trachea extraction step the volumes of both lungs without the main vessels can be measured. The final segmentation step results in the volumes of the left and right lungs including the vessels. The method has been tested by processing MR datasets from ten healthy participants. We compare our results with manually produced masks and obtain high agreement between the expert reading and our method: the True Positive Volume Fraction is more than 95%. The proposed automatic approach is fast and accurate enough to be applied in clinical routine for processing of thousands of participants. © 2011 Elsevier Ltd. All rights reserved.
1. Introduction Whereas high-resolution X-ray computed tomography (CT) is considered to be the gold standard for pulmonary imaging, its employment in an epidemiological cohort or repeated CT scans for research purposes in single patients are ethically not justified due to the ionizing nature of CT. On the contrary, magnetic resonance imaging (MRI) is a non-radiation based examination method, which gains an increasing popularity in research settings in healthy subjects, which is especially relevant for epidemiological studies. For example, in the frames of the Study of Health in Pomerania (SHIP) [8,30] whole-body MR imaging is taken. One of the research tasks of this project is to establish populationbased MRI reference values for lung volumes and to correlate MR
∗ Corresponding author. Tel.: +49 3834867554. E-mail addresses: [email protected], [email protected] (T. Ivanovska), [email protected] (K. Hegenscheid), [email protected] (R. Laqua), [email protected] (J.-P. Kühn), [email protected] (S. Gläser), [email protected] (R. Ewert), [email protected] (N. Hosten), [email protected] (R. Puls), [email protected] (H. Völzke). 0895-6111/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compmedimag.2011.10.001
findings with non-invasive pulmonary function tests, which assess static or dynamic lung volumes. The obtained lung volumes will be associated with such factors as the influence of genetic variations, endocrine metabolism, cardiac or pulmonary diseases, systemic inflammation, and mortality [21,7,23,24,6]. The results of dynamic pulmonary tests may be significantly influenced by several factors, whereas manual segmentation of MR images is a very laborious, observer-dependent, and time-consuming process, which is prone to errors. Therefore, fast, reliable, and reproducible automatic segmentation methods are required. However, lung imaging is a challenge in MRI because of the predominance of air within the lungs and associated susceptibility issues as well as the lower spatial and depth resolution when compared to the CT imaging and low signal to noise ratio of the inflated lung parenchyma. Moreover, in CT images pulmonary regions have the highest contrast when compared to other tissues, whereas in MR images all tissues, such as bones, muscles, fat are clearly distinguished. In this paper, a fast automatic method for segmenting the lung in MR images is presented. The method includes four main steps: First, an initial segmentation is performed, where the both lungs and trachea are extracted. The vessels are excluded from the parenchyma. Then, the trachea is eliminated and lungs are separated. Finally,
282
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293
the holes and cavities of the lung tissues are closed. Each step of the pipeline is thoroughly described and an analogy between the current approach and other methods is drawn. The paper is organized as follows. We give an overview of the existing methods in Section 2 and describe the pulmonary MR data created in SHIP in Section 3. Then we present our pipeline in Section 4 and compare the results against the manually established ground truth in Section 5. The findings are discussed in Section 6. Section 7 concludes the paper.
2. Related work In the literature, numerous segmentation strategies for CT data exist, whereas only few methods for MR data are addressed. Lelieveldt et al. developed a model-based method that allowed for extraction of lungs as well as other nearby organs such as the heart and diaphragm from MR data with different projection planes [14]. The authors reported that their method was designed to find the approximate borders to initialize subsequent more accurate local segmentation process. However, for our purposes it is enough to extract the lungs directly without the model-based presegmentation. Although numerous methods for CT data can be applied to MR images, there appear a number of difficulties which are specific for this imaging modality. In the following, we describe the main directions for lung segmentation from CT images and explain what problems arise on the correspondent processing steps with the MR data. In pulmonary CT images, the natural contrast between the low-density lungs and the surrounding high-density chest wall is usually used to guide the segmentation process. A standard approach applied by many authors [2,10,15,31,5] is binary thresholding. Here, either optimal Otsu thresholding [20] or iterative thresholding is applied [10]. In some previous works [11], manual thresholding is utilized. Shojaii et al. proposed to apply markerbased Watershed segmentation [29] with the defined internal and external markers [25]. In the case of MRI, the binary thresholding is not sufficient due to the different contrast for bones, air, and soft tissues. Thereafter, the lungs are separated from the large airways. The task of the trachea extraction is considered in many approaches [18,25] for CT data and it is solved mainly by either region growing [27] or morphological operations [26]. The start of the trachea is automatically identified by searching for a large, circular, air-filled region near the center of the first few slices in the data set. Regions in the current slice provide potential seed points positions for the next slice. The slice-by-slice growing procedure terminates when the size of the region on a new slice increases dramatically, indicating that the airways have merged into the low-density lung tissue [10]. Due to the lower spatial resolution and signal to noise ratio in MR images, such a criterion is not fully sufficient and needs to be extended. We introduce a procedure that extends this approach and overcomes the problem. Then the lungs need to be separated from each other. One of the most popular approaches is to apply morphological operations to detect the region of interest and then a 2D maximum cost path algorithm [27] is used to follow the junction line. For example, Hu et al. [10] and Li et al. [15] applied erosion and dilation to separate the lungs and then to restore the boundary as much as possible. A search region is found on a 2D slice, and it is propagated to successive slices due to the fact that the junction line position varies slowly through the data set. Li et al. [15] proposed to use the geodesic dilation, which is faster than the conditional dilation. The geodesic dilation terminates when no pixels can be added without changing the connectivity of the lung regions [15]. The smoothness
of the pulmonary anatomy in MRI is not guaranteed due to the lower spatial resolution. Hence, the 2D based method can result in a rather coarse boundary, and the application of a procedure that produces a 3D boundary is more beneficial in terms of accuracy. For example, Kuhnigk et al. [13] applied a 3D marker-based watershed transform with one automatically generated seed point per lung to obtain separate lung masks [13]. However, such a method will be more computationally expensive. We use an optimization step that allows for finding the boundary in 3D and at the same time keeping the computational costs relatively low. For smoothing and filling the cavities, practically in all previous works [15,10,25,1] a morphological closing method is applied to the lung boundary. This method is also referred to as the rolling ball filter. We developed a processing pipeline for MR data, which allows us to overcome the above-mentioned difficulties and leads to a fast and reliable segmentation, which is extremely important when thousands of datasets are processed.
3. Materials The design of the SHIP project is given elsewhere [8]. In brief, SHIP is a population-based project consisting of the two independent studies SHIP and SHIP-Trend. The examination programme of the 11 year examination follow-up of the original SHIP (SHIP-2) and SHIP-Trend was supplemented with whole-body MR imaging [30,8]. All MRI studies are performed on a 1.5-T MR imager (Magnetom Avanto; Siemens Medical Systems, Erlangen, Germany). Subjects are placed in the supine position, and five phased array surface coils are placed on the head, neck, abdomen, pelvis, and lower extremities. A spinal coil is embedded in the patient table. Each participant undergoes standardized WB-MRI including an axial T1-weighted VIBE (volume-interpolated breath hold examination) sequence (typical parameters: repetition time (TR) 31 ms, echo time (TE) 11 ms, flip angle 8◦ , voxel size 1.8 × 1.8 × 3 mm) and an axial T2-weighted HASTE (half-Fourier single-shot turbo spinecho) sequence (TR 5500 ms, TE 220 ms, flip angle 150◦ , voxel size 2.3 × 1.8 × 5 mm) for chest imaging. To avoid motion artifacts 21 s of breath hold are needed for the VIBE sequence. The HASTE sequence is taken in two steps and then the parts are composed together with the navigator control. In total, 40 s are needed for the HASTE sequence. Because of the small slice thickness (3 mm) and bigger spatial resolution when compared to the HASTE sequences the automatic lung segmentation and volumetry is performed using the VIBE sequences. For our experiments we used 10 randomly selected datasets of healthy participants from the SHIP-Trend programme. Each dataset consists of 88 slices with a resolution of 512 × 416.
4. Methods Our method consists of four steps: the main extraction and three refinement steps. First, the lung and trachea excluding the main vessels are extracted. Second, the trachea is separated from the lung tissue, and the cut is made on the main bronchi entering the lung parenchyma. Due to the fact that a participant holds the breath during the scan, which results in a very thin and faint boundary between the lungs, the main extraction step may fail to detect this boundary on certain slices. Third, we apply the lung separation procedure. Fourth, the lung boundaries are smoothed and the cavities in the lung tissue are filled, which allows us to obtain the completed lung parenchyma. In each processing step, we compute the volume of the segmented tissue.
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293
283
Fig. 1. Lungs and trachea segmented with K-Means.
4.1. Lung and trachea segmentation MR technique allows for better contrast within the soft tissues, which are rather distinctively separated in terms of intensity, so we assumed that any automatic clustering technique could be applied here. One of the classical machine learning methods, K-Means [17], for example, is a good candidate due to its speed and simplicity. KMeans aims to separate n observations (x1 , x2 , . . ., xn ) into k classes (k ≤ n) S = S1 , S2 , . . ., Sk so as to minimize the within-cluster sum of squares: E = argmin S
k
xj − i 2 ,
i=1 xj ∈Si
where i is the mean of the points in Si . We experiment with the number of clusters K and determine that setting K = 4 or 5 is enough to separate the lungs from the vessels, which are brighter than the lung parenchyma, but much darker than the other tissues. Such a value of K was set automatically for all datasets. It is known that K-Means method results are strongly dependent on the initial cluster selection [4]. To avoid the inadequate clustering, e.g., the case when the cluster centers fall too close together, we set the initial cluster centers as far as possible from each other and in each step of the algorithm iterative procedure, if a cluster contains less than 2% of all data voxels, we place its center to the mean of the two other centers accounting for the most voxels. In Fig. 1, the result of the segmentation is shown.
components, say, Ci = {c1i , . . . , cKi }, have been detected there. We stop at the slice i, if a connected component cˆ ∈ Ci , which satisfies two conditions, is found: (1) {∀c ∈ Ci , A(ˆc ) ≥ ATr }, (2) {∀c ∈ Ci , |xcˆ − 0.5W | = min(|xc − 0.5W |)},
(1)
where A(c) defines the area of the connected component c in pixels, W is the slice width, xc is the x coordinate of the center of mass of the connected component c, and ATr is the area threshold. The largest connected component, whose center of mass lies as close as possible to the slice center in the x direction, is selected as the initial trachea region. Thereafter the trachea is located. In Fig. 3, three consecutive slices are shown. The segmented tracheal regions are depicted in green. We traverse the trachea along the z axis analysing each time a pair of consecutive slices. Let us assume the trachea is already detected in slice i and is to be determined in slice j = i − 1. In slice i the trachea consists of K connected components, say, Ci = {c1i , . . . , cKi }, and their cumulative area is A(Ci ). The set of connected components Ci is orthogonally projected to slice j and M j j connected components Cj = {c1 , . . . , cM } are selected if the intersection with the projected connected component set Ci is non-zero: j j ∃cki , ∃cm : cki ∩ cm = / , where k ∈ [1, . . ., K] and m ∈ [1, . . ., M]. If a cumulative area A(Ci ) or A(Cj ) is smaller than a threshold ATr , then the processing stops. Otherwise, we compute for all pairs of interj j secting components (cki , cm ), where cki ∈ Ci and cm ∈ Cj , the area ratios: j
4.2. Trachea and main bronchi extraction
R
The extraction task can be divided into two parts: first, the start component of the trachea is located. Second, the trachea and main bronchi are found using a 3D component analysis procedure. Under “connectedness” we understand 8- and 26-neighborhood for 2D and 3D, respectively. The cut is made on the spots, where the main bronchi enter the lung tissue. Here, we would like to explain both steps in detail. In Fig. 2, the 3D and 2D views of the trachea initialization are shown. Starting from the slices, which anatomically lie in the direction of the larynx, we process each slice and find its connected components. Let us denote the processed slice i, and K connected
If R
=
j (c i ,cm ) k
j
(c i ,cm )
A(cm ) A(cki )
.
(2)
is greater than a preselected threshold Rconst , we assume
k
j
that the trachea component cm has merged with the lung tissue and a separation procedure is called for it. To separate the airways from the lung tissue, we compute a distance map [27] D. The distance map is the result of the distance transform. The distance transform (DT) is the transformation that generates a map D whose value in each pixel p is the smallest distance from this pixel to background Bkg: D(p) = min{d(p, q)|q ∈ Bkg},
(3)
284
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293
Fig. 2. The start trachea component is marked red and shown in 3D (left) and 2D (right). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
where d is the Eucledian distance, p is a foreground pixel, and q is a background pixel. Then the distance map is inverted: Dinv = − D, and the 2D Watershed [29,22] is applied to Dinv . The Watershed result is a set of the regions separated by black boundaries. Each region has a unique ID. The regions that have a non-zero intersection with the projected trachea component cki are taken and the ratios are calculated again. If the Watershed transform failed to detect the sought boundary, i.e., the correspondent ratio is still greater than the threshold, then the boundary is created artificially by projecting the trachea region from the current slice i. We permit at most two consecutive artificial boundaries, so that a “fake” trachea inside the lung tissue is not created. Thereafter, the same procedure is repeated for the subsequent slice pairs till i ≥ 1. This procedure allows for separation of the trachea and main bronchi from the lung tissue, and the cut is approximately made on the spot, where the main bronchi enter the lung tissue. When the extraction is finished, we exclude the trachea and main bronchi mask from the further processing. Any small regions between two lungs, which can appear after the extraction and represent noise, are also excluded. In Fig. 4, a 3D view of the extracted trachea is shown.
4.3. Lung separation After the trachea extraction, we analyse the remaining connected components in 3D and exclude all components with area, which is smaller than a selected threshold A3D . This results either in one (the lungs have merged together) or two (the lungs are separate) connected components. If the lungs represent one connected component, they need to be separated. Our approach allows for detection of a boundary in 3D on the slices, where the lungs have merged together. It consists of the following steps. First, the slices, where the lungs merged together are detected. Second, a subvolume including the Volume of Interest (VOI) is selected. Third, the distance map [27] is built and corrected for the computational costs reduction. Finally, 3D Watershed [29] is applied to the corrected and inverted distance map, and two separated lung masks are built. Here, we describe the steps of the procedure in more detail. There are two possible lung junctions: anterior and posterior. Since the processing of both cases is identical, we describe it only once using the anterior junction example. In Fig. 5, an example of the anterior junction on a 2D slice is shown. First, we analyse a 3D dataset and select the slices that contain only one connected component. Then, a subvolume is constructed from this list including all intermediate slices and two additional slices for the beginning and the end, e.g., if the slice list consists of slice ids 38, 41, 42, the selected subvolume will include the slices that have numbers 37 → 43. This subvolume can be separated into two subvolumes: 37 → 39 and 40 → 43, which allows to reduce the computational costs in the Watershed step. For simplicity, we refer to one (total) subvolume in the text. We reindex the subvolume slice numbers and set them 1, . . ., N. Let V = (S1 , S2 , . . ., SN ) denote the voxel sets (slices) of our subvolume. Thereafter, a new subvolume V = , . . . , S is formed, so that the first and the last slices S 2 N−1 mid are projected on the whole subvolume and the coinciding voxels are masked: Si = Si \ (S1 ∪ SN ) ,
Fig. 3. Three consecutive slice with the segmented large airway regions.
where i ∈ [2, . . . , N − 1]. Such a masking allows one to exclude the coinciding voxels. However, it contains the voxels, which correspond to the junction, since there is no junction in the first and last slices. Then, we detect 3D connected components. The connected component cj that represents the anterior junc tion is visible in the whole subvolume V mid . The junction lies in
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293
285
Fig. 4. The extracted Trachea from Fig. 1 is shown in an overlaid manner.
the upper third of each slice and is close to the slice center in x direction: (1) {∀c ∈ V mid , ycj ≤ 0.3H, (2) {∀c ∈ V mid , |xcj − 0.5W | = min(|xc − 0.5W |)}, where xc and yc are the coordinates of the center of mass of the connected component c, W is the slice width, and H is the slice height. We build a bounding box V bb around the selected connected component to assure that the junction is included. The bounding box V bb include the voxels that satisfy: {∀v ∈ V : lxcj − f1 ≤ v · x ≤ rxcj + f1 , tycj − f2 ≤ v · y ≤ bycj + f2 }, where lx, rx, ty, by are the left and right x coordinates and top and bottom y coordinates of the connected component cj , correspondingly. Moreover, f1 and f2 are preselected factor values. Thereafter, a 3D distance map [27] is computed and inverted. The minimal distance map values are the “deepest” points, lying inside the lung tissue. The junction does not lie deep and has high values on the inverted distance map. Hence, the distance values can be thresholded. Let Dinv denote the inverted 3D distance map and the thresholded distance map in each voxel p is defined as D inv (p) = max(Dinv (p), const),
(4)
where const is a constant. It reduces the computational costs of the Watershed transformation, which starts from the minima of the distance map. The constant value const is computed from the distance values of the connected component cj and the distance map minima. Thereafter, a 3D Watershed [29] is applied to the cor rected distance transformation D inv . In Fig. 6, a close up view of the watershed result where the connected lung parts are successfully separated is shown.
Fig. 5. A 2D slice from the dataset shown in Fig. 1 with the anterior junction.
in an overlaid manner. One can observe a smooth lung boundary with the filled cavities. In Fig. 8, the final lung volumes are overlaid with the initial dataset.
4.4. Filling lung cavities and holes The final processing step is to fill the cavities in the lung boundaries, obtain the complete lung masks, and compute the final lung volumes. First, the holes in the lung tissue are closed. Thereafter, we apply a morphological closing procedure to the lung boundary to cover the boundary cavities. A spherical 2D kernel of radius 9 is used. Each lung is processed separately to avoid artificial lung junctions. In Fig. 7, the resulting mask for the slice from Fig. 5 is shown
Fig. 6. A close up view of the separation result for the 2D slice shown in Fig. 5.
286
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293
Fig. 7. The smoothed lung tissue for the 2D slice shown in Fig. 5 is presented. The left and right lungs are shown in red and green colors, correspondingly. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
5. Results We tested our processing pipeline and extracted the volumes on two stages: the total volume of both lungs and trachea excluding vessels (after the initial segmentation) and the separate volumes of both lungs with the filled holes and closed cavities (after the final step) are computed. The initial step takes around 1 min and the whole pipeline takes around 3 min on Intel(R) Core i7 − 950 (Quadcore) @ 3.06 GHz computer with 6 Gb DDR3 for one dataset with a resolution 512 × 416 × 88. All sequences have been processed automatically without an additional parameter tuning for each participant. The following parameter values have been set: K = 4 is the number of classes for segmentation, ATr = 50 is the minimum area of a 2D trachea connected component, Rconst = 3 is the maximum area factor between the connected component areas on consecutive slices. Moreover, A3D = 5000 is the minimum size of a 3D lung connected component, lungMinThresh2D = 500 is the minimum size of a 2D lung connected component, f1 = f2 = 50 are the padding values for the bounding box, Kr = 9 is the kernel radius for the rolling ball filter. The ground truth is determined by manual identification of the lung boundaries. In total, we obtained manual readings from two experts who traced the lung boundaries. One image analyst is a
Fig. 8. The final lung volumes are overlaid with the initial dataset. The left and right lungs are shown in red and green colors, correspondingly. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
Fig. 9. A 2D slice from dataset 6 is shown. The region close to mediastinum defines the difference between the expert readings for the left lung mask.
radiologist (R.L.), and he determined all types of masks on 880 slices. The other expert is a very experienced pulmonologist (R.E.), and he determined two types of masks (separate left and right lungs) on 880 slices. The image analysts worked independently of each other. We have evaluated the difference between the manually delineated volumes and observed that the expert agreed on around 95% of the volumes, i.e., the average proportions of the overlapping voxels in both volumes are 93.87% with the standard deviation 1.52 for the left lung masks and 94.76% with the standard deviation 1.01 for the right lung masks. The correspondent mean deviations for the detected volumes are 3.89% for the left lung masks and 3.46% for the right lung masks. In Fig. 9, an example slice from dataset 6 is shown. The region of the left lung, which is close to the mediastinum is depicted in a close-up view. The region of expert variations is colored in green and red. Fig. 10 shows the expert readings in an overlaid fashion. As it can be observed, the experts marked the region near the mediastinum differently, which results in the variation interval (0, 6.5]% for the complete masks. The result from one the experts is compared with the lung masks detected by our pipeline. To evaluate the accuracy of the method, we applied a methodology described by Udupa et al. [28]. We computed the following measures: False Negative Volume Fraction (FNVF), which indicates the fraction of tissue, which was manually marked but missed by the automatic method; False Positive Volume Fraction (FPVF), which indicates the fraction of the tissue that was marked by the method but not marked but the expert; True Positive Volume Fraction (TPVF) is the number of overlapping voxels marked both by the method and the expert. Additionally we calculated the True Positive Volume Error (TPVE), which shows the deviations between the automatically detected volume and the manually delineated one. However, we have to mention that this measure does not regard the exact boundaries. We computed it, as our ultimate goal is to assess the volume of the lungs. Let Cauto and Cexp denote the binary masks produced by our pipeline and the expert delineation, correspondingly. Moreover, let Volauto and Volexp to be the volumes computed by the pipeline and
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293
287
Fig. 10. The left lung masks from two experts for the slice from Fig. 9 are shown in green and red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
the expert, correspondingly, then the FNVF, FPVF, TPVF, and TPVE are defined as follows:
Cexp − Cauto Cauto − Cexp , , FPVF = FNVF = Cexp Cexp Cauto ∩ Cexp Volauto − Volexp , TPVE = TPVF = . Cexp Volexp
Our evaluations are presented in Tables 1, 2, and 3 for the total volume of both lung with trachea, the volume of the right lung, and the volume of the left lung, correspondingly. The following i notations for the headings are used: “Volexp ” is the expert-detected volume detected by expert i; “Volauto ” means the automatically detected volume; Cexp , “FNVF”, “FPVF”, “TPVF”, “TPVE” are defined as described above.
Table 1 Comparison of the automatic results to the expert-defined ground truth for ten different data sets. Lungs and trachea without vessels. Here, the manual results are produced by one expert (radiologist). Lungs and trachea without vessels 1
ID
Cexp (voxel)
FNVF (%)
FPVF (%)
TPVF (%)
Volexp (ml)
Volauto (ml)
TPVE (%)
1 2 3 4 5 6 7 8 9 10 Mean StDev
1,068,926 1,423,881 1,304,202 1,690,847 2,035,239 1,275,871 1,503,466 1,311,256 1,609,496 1,550,048
3.41 2.04 5.56 1.95 3.77 8.53 1.74 0.36 2.83 2.32 3.25 2.32
0.08 0.42 0.48 4.27 1.18 0.47 1.1 1.81 1.51 0.8 1.21 1.2
96.59 97.96 94.44 98.05 96.23 91.47 98.26 99.64 97.17 97.68 96.75 2.32
2477.16 3299.74 3022.39 3918.42 4716.52 2956.74 3484.17 3038.74 3729.38 3592.13
2394.74 3246.2 2868.71 4009.72 4594.69 2718.44 3461.78 3082.79 3680.55 3537.42
3.33 1.62 5.08 2.33 2.58 8.06 0.64 1.45 1.31 1.52 2.79 2.24
288
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293
Table 2 Comparison of the automatic results to the expert-defined ground truth for ten different data sets. Right lung with vessels and filled holes. Here, the manual results are produced by one expert (radiologist). Right lung ID
Cexp (voxel)
FNVF (%)
FPVF (%)
TPVF (%)
1 2 3 4 5 6 7 8 9 10 Mean StDev
608,632 851,018 833,722 963,975 1,141,684 704,715 837,662 706,295 935,049 885,919
2.32 1.55 4.23 0.89 2.7 5.8 2.03 1.09 4.08 1.87 2.66 1.57
7.74 2.31 3.07 4.98 2.58 2.87 4.39 5.37 2.84 3.42 3.96 1.69
97.68 98.45 95.77 99.11 97.3 94.2 97.83 98.91 95.92 98.13 97.33 1.57
1
Volexp (ml)
Volauto (ml)
1410.46 1972.17 1932.09 2233.94 2645.77 1633.13 1941.22 1636.79 2166.91 2053.05
1486.82 1986.99 1909.84 2325.4 2642.72 1585.18 1984.36 1706.77 2140.16 2084.89
TPVE (%) 5.41 0.75 1.15 4.09 0.12 2.94 2.22 4.28 1.23 1.55 2.37 1.74
Table 3 Comparison of the automatic results to the expert-defined ground truth for ten different data sets. Left lung with vessels and filled holes. Here, the manual results are produced by one expert (radiologist). Left lung 1
ID
Cexp (voxel)
FNVF (%)
FPVF (%)
TPVF (%)
Volexp (ml)
Volauto (ml)
TPVE (%)
1 2 3 4 5 6 7 8 9 10 Mean StDev
507,028 636,388 552,772 766,651 973,428 635,371 704,294 668,718 774,810 709,396
3.35 2.58 5.32 2.56 2.15 5.12 1.34 1.35 3.02 2.55 2.94 1.36
2.16 2.94 2.69 5.61 2.91 1.72 2.73 3.55 3.3 2.06 2.97 1.08
96.65 97.42 94.68 97.44 97.85 94.88 98.66 98.65 96.98 97.45 97.06 1.36
1175 1474.78 1281.01 1776.66 2255.85 1472.43 1632.15 1549.71 1795.57 1643.97
1160.96 1480.12 1247.3 1830.8 2272.94 1422.35 1654.75 1583.81 1800.67 1635.9
1.19 0.36 2.63 3.05 0.76 3.4 1.38 2.2 0.28 0.49 1.58 1.16
6. Analysis and discussion 6.1. Comparison to other techniques
we will compare our results to the well-known CT segmentation techniques in a step-by-step manner and explain the benefits of our approach.
One of the most popular directions for segmentation of lung CT data are pixel-based methods [10,15,13,9]. Here, main steps are thresholding, extraction of large airways, lung separation, and smoothing. Our method belongs to the same group and the pipeline steps are adapted for MRI data. In the following section,
6.1.1. Thresholding A thoracic CT scan contains two main groups of pixels: highintensity pixels, which are located in the body, and low-intensity pixels that are in the lung and surrounding air [9]. There is a huge difference between the pixel values in these groups, so the
Fig. 11. A slice from the initial dataset (left) and result for Otsu thresholding (right) are shown side-by-side. Otsu thresholding fails to separate the lungs from the other tissues.
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293
image is thresholded. Here, the simplest thresholding at a fixed value often suffices [12,9]. Another possibility is to apply a binary automatic thresholding, such as Otsu thresholding [20] or the optimal thresholding [10] proposed by Hu et al. Since MR technique allows for a better contrast within the soft tissues, the binary thresholding, which is usually applied on high resolution CT data [10,15,13], would fail in separating lungs and airways from other tissues, which is documented in Fig. 11. As it can be observed in Fig. 12, the MR thorax data consists of three main intensity classes: muscles (gray), bones and fat (light gray), lungs and trachea (dark), which are well separated from each other. We applied two classical approaches, namely Multiclass Otsu thresholding [16] and K-Means [17] clustering. Then we extracted the largest connected component, which represents the lungs, trachea, and main bronchi and compared the results. In Fig. 13, the segmentation of both lungs and large airways excluding the vessels is shown in a 2D slice. The results were comparable in quality: TPKMeans = 97.95%, TPOtsu = 97.4%, but K-Means had lower computational costs.
289
Fig. 12. An slice from a lung dataset.
Fig. 13. A close up view of the segmentation results for a slice shown in Fig. 11 obtained with K-Means (left) and Otsu (right). The results are similar in quality, but K-Means is superior in terms of speed.
Fig. 14. ROIs from three consecutive slices. The trachea component, which has not reached the carina yet, is connected on one slice (middle) and then disconnected again (right).
290
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293
6.1.2. Extraction of large airways To analyse the lung volumes, the trachea and large airways must be detected and separated from the left and right lungs. Since the airways are empty cavities, the intensity of pixels in the area is low. In CT images, many techniques can be applied to extract the trachea and large airways. Heuberger et al. reported that a predefined threshold often suffices [9]. To remove the trachea and large airways, the authors searched for an area with an average intensity lower than 0.5Tc , where Tc is a threshold value, computed by the optimal thresholding method [10]. Such an approach is definitely not applicable for the MR data, since the trachea and large airways have practically the same intensity and noise level as the lungs. Li et al. treated the trachea and bronchi as small round regions, which were removed using a simple morphology technique [15]. First, all connected regions (black and white one) on each slice were labeled. Then, the regions with the area smaller than a predetermined threshold were inverted, i.e. operator NOT was applied to each pixel inside the regions. Hu et al. proposed to use a slice-byslice region growing procedure, which stops, when the size of the region on a new slice increases dramatically, indicating that the airways have merged into the low-density lung tissue [10]. These two approaches have some disadvantages. According to the procedure, described by Hu et al. [10], only the separated trachea parts will be erased before the first connection to the lung tissue is reached. The approach proposed by Li et al. [15] will delete the separated trachea regions in the whole volume. However, on the slices, where the trachea is connected to the lung, no separation can be done. For the MR data, which have lower spatial resolution and signal-to-noise ratio, the result of these approaches will add a significant number of voxels, that belong to the large airways, to the lung volumes. For example, in Figs. 14 and 15 the extracted lung mask is overlaid with the original data and two groups of regions of interest (ROIs) that contain the trachea and main bronchi are shown. In Fig. 14, the trachea has not reached the carina yet, but is merged with the lung tissue and separated again on the consecutive slice. Here, Hu’s procedure would stop on the “merged” spot and the rest of the large airways (including the trachea, carina and main bronchi) would be referred to the lung volume. Li’s approach will exclude the trachea in the left and right slices, but not in the middle one. In Fig. 15, the trachea is separated into the main bronchi that can be at least partially excluded from the lung tissue, which can be accomplished by neither of three above mentioned methods. Our approach is also based on a region growing procedure, but it processes the whole trachea and main bronchi and separates the merged parts from the lung tissue. This allows us to extract the significant parts of the large airways and exclude them from the lung volumes.
6.1.3. Lung separation Finally, the two lungs are identified and, if needed, separated. Several methods have been proposed for CT images. Brown et al.
Fig. 15. ROIs from three consecutive slices. The trachea component reaches the carina and splits into two bronchi components (top). The right main bronchus is connected to the lung tissue (middle and bottom). The connection spot is marked with red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
proposed to apply a 2D maximum cost path algorithm [27], where the cost of a pixel is its graylevel, can be used to follow the anterior junction line [3]. Hu et al. and Li et al. utilized the same approach, but applied different procedures for detection of the connection ROI [10,15]. The dynamic programming approach is applicable in the case of CT data, since the soft tissue of the junction has a slightly higher gray level than the lungs, and forms a continuous line in each axial slice of the connecting region [3]. These two criteria are usually not fulfilled in the case of the MR data due to the lower spatial resolution and partial volume effects. It is documented by Fig. 16, where ROIs from three consecutive slices are shown. The junction is not always continuous here and junction voxels have the same intensity values as the lung tissue. Moreover, a 2D based approach may lead to the result, when the lungs will overlap in 3D, being separated in 2D.
Fig. 16. ROIs from three consecutive slices with the anterior junction marked red. The junction voxels have practically the same intensity as the lung tissue voxels. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293
291
Fig. 17. Automatically detected result (left) consists of less parenchyma (more vessels are excluded) than the expert-defined result (right) for a slice shown in Fig. 12.
We utilize an 3D marker based Watershed transform, similarly to Kuhnigk et al. [13]. The appropriate selection of the ROI and markers allows us to reduce the computational costs for the Watershed transform (cf. Section 4). 6.2. Quality analysis As it can be taken in Table 1, our automatic procedure produces the results with high TPVF (96.75%) and rather low FNVF and FPVF (3.25% and 1.21%, correspondingly). However, it is obvious that FNVF is higher. The expert has to manually trace not only the lung and trachea boundaries, but also cross out the vessels, which are brighter than the lung tissue. This procedure is very laborious and time-consuming. The automatic method definitely detects more bright vessels than the expert does. Fig. 17 shows a 2D slice from dataset 6, which is an outlier among our test datasets with the FNVF equal to 8.53%. The TPVE is equal to 2.79%, which is exact enough for our research purposes. The findings presented in Tables 2 and 3 also show a high agreement with the expert reading. The TPVF for the separate right and left lungs (97.06% and 97.33%, correspondingly) is higher than in the case of the detection of the lungs and trachea, where the vessels have to be manually excluded. However, closing the vessels also increases FPVF (3.96% and 2.97% for the right and left lungs, correspondingly). In Fig. 18, a 2D example of a right lung mask from dataset 1 demonstrating this finding is shown. Since our automatic procedure does not distinguish between the different types of holes, i.e., whether it is a hole where the vessel is located or a hole, which appears during the scanning, some “fake” lung voxels can be added to the mask. Nevertheless, the achieved accuracy in True Positive Volume is enough and lies within the variation interval between the experts. Moreover, TPVE is equal 2.37% for the right lung and 1.58% for the left lung. These errors are not significant and can be neglected for our research tasks. 6.3. Performance analysis We analyse the performance of our approach using the test datasets. As it can be seen in Fig. 19, the total execution times lies around 150 CPU seconds. However, there are several outliers, namely, datasets 4, 7, 10. The maximum execution time reaches 452.27 for dataset 4. We selected eight main steps in our algorithmic pipeline (see Table 4) and analysed the performance for each of them. The results are documented in Fig. 20. The most expensive step with the biggest dispersion is the Lung separation, which includes a marker based 3D Watershed algorithm on a selected ROI.
Fig. 18. A slice from dataset 1 (top). Automatically detected result (bottom) covers both vessel locations and a “fake hole”, which appears due to the scanning, whereas the expert marks only the vessel holes (middle).
292
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293
300
we implemented an application, which includes a batch processing mode. The application is developed under MeVisLab framework [19]. In total, the application can be used in 2 modes: single-file step-by-step processing and automatic batch processing modes. The former mode is primarily designed for parameter selection and evaluation. With the latter mode the application can run in the background over thousands of participants using the same parameter settings for all of them.
250
7. Conclusions
200
We presented an approach for automatic lung extraction from MR data, which consists of one main step and three refinement procedures. For our research task, the volumes of both lungs with trachea and excluded vessels and separate full lung masks are of a particular interest. The accuracy and speed of the proposed method allow for using it in the clinical routine to evaluate thousands of participants and provide the further research background in such epidemiological studies as the SHIP Project.
500 Execution time 450
CPU sec
400 350
150 100
0
1
2
3
4
5
6
7
8
9
10
11
Dataset # Fig. 19. The total execution time for 10 test dataset in CPU seconds. The average time is lower than 3 min, if no lung separation is done.
Acknowledgement
Table 4 Algorithm steps description. Step
Description
Min time (CPU sec)
Max time (CPU sec)
1 2 3 4 5 6 7 8
Smoothing Clustering Save result Exclude trachea Separate lungs Save result Fill holes Compute volumes, save result
51.92 9.95 0.99 2.5 12.33 3.4 41.29 8.01
52.8 13.91 1.09 5.8 311.92 3.66 59.62 8.27
References
Here, the average value is 65.173 CPU seconds and the standard deviation equals 104.7. The performance of the algorithm strongly depends on the ROI size. We assume that the performance can be improved if a more efficient algorithm implementation is used. We leave this improvement for future work. 6.4. Application for medical research The algorithm in its current state can be applied in an epidemiological setting. To be conveniently embedded for medical research, 200
Execution time for each step
CPU sec
150
100
50
0
-50
1
2
3
4
5
6
7
The authors would like to thank Prof. Dr. Volkmar Liebscher (Institute of Mathematics and Informatics, Ernst-Moritz-Arndt University, Greifswald, Germany) for support and valuable discussions.
8
Pipeline steps Fig. 20. Mean and standard deviation for each pipeline step are shown using error bars. The most expensive step with the biggest dispersion is the lung separation. This steps depends on how many slices it must be done.
[1] Armato SG, Giger ML, Moran CJ, Blackburn JT, Doi K, MacMahon H. Computerized detection of pulmonary nodules on CT scans. Imaging and Therapeutic Technology 1999. [2] Armato SG, MacMahon H. Automated lung segmentation and computeraided diagnosis for thoracic ct scans. International Congress Series 2003;1256:977–82. [3] Brown MS, McNitt-Gray MF, Mankovich NJ, Goldin JG, Hiller J, Wilson LS, et al. Method for segmenting chest ct data using an anatomical model: preliminary results. IEEE Transactions on Medical Imaging 1997;16(6):828–39. [4] Duda RO, Hart PE, Stork DG. Pattern classification. Wiley; 2000. [5] El-Baz A, Farag AA, Falk R, Rocca RL. A unified approach for detection, visualization, and identification of lung abnormalities in chest spiral ct scans. International Congress Series 2003;1256:998–1004, cARS 2003. Computer Assisted Radiology and Surgery. Proceedings of the 17th International Congress and Exhibition. [6] Gläser S, Ittermann T, Koch B, Völzke H, Wallaschofski H, Nauck M, et al. Airflow limitation, lung volumes and systemic inflammation in a general population. European Respiratory Journal 2011, Jun 30 [Epub ahead of print]. [7] Gläser S, Friedrich N, Ewert R, Schäper C, Nauck M, Dörr M, et al. Association between serum insulin-like growth factor (igf) i and IGF binding protein-3 and lung function. Journal of Clinical Endocrinology and Metabolism 2009;94(7):2452–8. [8] Hegenscheid K, Kühn J, Völzke H, Biffar R, Hosten N, Puls R. Whole-body magnetic resonance imaging of healthy volunteers: pilot study results from the population-based ship study. Rofo 2009. [9] Heuberger J, Geissbuehler A, Mueller H. Lung ct segmentation for image retrieval using the insight toolkit. In: Proceedings of the 1st International Conference on Medical Imaging and Telemedicine (MIT ‘05). 2005. [10] Hu S, Hoffman EA, Reinhardt JM. Automatic lung segmentation for accurate quantitation of volumetric X-ray CT images. IEEE Transactions on Medical Imaging 2001;20:490–8. [11] Kalender WA, Fichte H, Bautz W, Skalej M. Semiautomatic evaluation procedures for quantitative ct of the lung. Journal of Computer Assisted Tomography; 1991. [12] Kemerink GJ, Lamers RJ, Pellis BJ, Kruize HH, van Engelshoven JM. On segmentation of lung parenchyma in quantitative computed tomography of the lung. Medical Physics 1998;25(12):2432–9. [13] Kuhnigk J, Hahn HK, Hindennach M, Dicken V, Krass S, Peitgen H. Lung lobe segmentation by anatomy-guided 3d watershed transform. In: Proceedings of SPIE, vol. 25032. 2005. p. 1270–3. [14] Lelieveldt BPF, Sonka M, Bolinger L, Scholz TD, Kayser HWM, et al. Anatomical modeling with fuzzy implicit surface templates: application to automated localization of the heart and lungs in thoracic mr volumes. Computer Vision and Image Understanding 2000:1–20. [15] Li W, Nie SD, Cheng JJ. A fast automatic method of lung segmentation in CT images using mathematical morphology. Berlin, Heidelberg: Springer; 2007, 2419–2422. [16] Liao P, Chen T, Chung P. A fast algorithm for multilevel thresholding. Journal of Information Science and Engineering 2001;17(5):713–27.
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293 [17] MacQueen JB. Some methods for classification and analysis of multivariate observations. In: Proceedings of 5-th Berkeley Symposium on Mathematical Statistics and Probability, vol. 1. 1967. p. 281–97. [18] Memon NA, Mirza AM, Gilani S. Limitations of lung segmentation techniques, vol. 27. US: Springer; 2009. p. 753–66, chapter 76. [19] MeVisLab. Software for medical image processing and visualization; 2011. http://www.mevislab.de. [20] Otsu N. A threshold selection method from gray-level histograms. IEEE Transactions on Systems, Man, and Cybernetics 1979;9(1):62–6. [21] Repapi E, Sayers I, Wain LV, Burton PR, Johnson T, Obeidat M, et al. Genomewide association study identifies five loci associated with lung function. Nature Genetics 2010;42:36–44. [22] Roerdink JBTM, Meijster A. The watershed transform: definitions, algorithms and parallelization strategies. Fundamenta Informaticae 2001;41:187–228. [23] Schroeder EB, Welch VL, Couper D, Nieto FJ, Liao D, Rosamond WD, et al. Lung function and incident coronary heart disease: the atherosclerosis risk in communities study. American Journal of Epidemiology 2003;158:1171–81. [24] Schuenemann HJ, Dorn J, Grant BJ, Winkelstein W, Trevisan M. Pulmonary function is a long-term predictor of mortality in the general population. Chest 2000;118:656–64.
293
[25] Shojaii R, Alirezaie J, Babyn P. Automatic lung segmentation in CT images using watershed transform. In: IEEE International Conference on Image Processing, 2005. ICIP 2005, vol. 2. 2005. p. 1270–3. [26] Soille P. Morphological image processing: principles and applications. Cambridge University Press; 1999. [27] Sonka M, Hlavac V, Boyle R. Image processing, analysis, and machine vision. Thomson; 2008. [28] Udupa JK, LaBlanc VR, Schmidt H, Imielinska C, Saha PK, Grevera GJ, et al. Methodology for evaluating image-segmentation algorithms, vol. 4684, SPIE, 2002, p. 266–77. http://link.aip.org/link/?PSI/4684/266/1. [29] Vincent L, Soille P. Watersheds in digital spaces: an efficient algorithm based on immersion simulations. IEEE Transactions on Pattern Analysis and Machine Intelligence 1991;13(6):583–98. [30] Völzke H, Alte D, Schmidt CO, Radke D, Lorbeer R, Friedrich N, et al. C. O. S. Cohort profile: the study of health in pomerania. International Journal of Epidemiology; 2010. [31] Zhao B, Gamsu G, Ginsberg MS, Jiang L, Schwartz LH. Automatic detection of small lung nodules on ct utilizing a local density maximum algorithm. Journal of Applied Clinical Medical Physics; 2003.
Contents lists available at SciVerse ScienceDirect
Computerized Medical Imaging and Graphics journal homepage: www.elsevier.com/locate/compmedimag
A fast and accurate automatic lung segmentation and volumetry method for MR data used in epidemiological studies Tetyana Ivanovska a,∗ , Katrin Hegenscheid b , Rene Laqua b , Jens-P. Kühn b , Sven Gläser c , Ralf Ewert c , Norbert Hosten b , Ralf Puls b , Henry Völzke a a
Institute of Community Medicine, Ernst-Moritz-Arndt University, Greifswald, Germany Institut of Diagnostic Radiology und Neuroradiology, Ernst-Moritz-Arndt University, Greifswald, Germany c Klinik und Poliklinik für Innere Medizin, Bereich Pneumologie und Infektiologie, Ernst-Moritz-Arndt University, Greifswald, Germany b
a r t i c l e
i n f o
Article history: Received 6 March 2011 Received in revised form 30 September 2011 Accepted 6 October 2011 Keywords: Segmentation MRI Pulmonary imaging Volumetry
a b s t r a c t In modern epidemiological population-based studies a huge amount of magnetic resonance imaging (MRI) data is analysed. This requires reliable automatic methods for organ extraction. In the current paper, we propose a fast and accurate automatic method for lung segmentation and volumetry. Our approach follows a “coarse-to-fine” segmentation strategy. First, we extract the lungs and trachea excluding the main pulmonary vessels. This step is executed very fast and allows for measuring the volume of both structures. Thereafter, we start a refinement procedure that consists of three main stages: trachea extraction, lung separation, and filling the cavities on the final lung masks. After the trachea extraction step the volumes of both lungs without the main vessels can be measured. The final segmentation step results in the volumes of the left and right lungs including the vessels. The method has been tested by processing MR datasets from ten healthy participants. We compare our results with manually produced masks and obtain high agreement between the expert reading and our method: the True Positive Volume Fraction is more than 95%. The proposed automatic approach is fast and accurate enough to be applied in clinical routine for processing of thousands of participants. © 2011 Elsevier Ltd. All rights reserved.
1. Introduction Whereas high-resolution X-ray computed tomography (CT) is considered to be the gold standard for pulmonary imaging, its employment in an epidemiological cohort or repeated CT scans for research purposes in single patients are ethically not justified due to the ionizing nature of CT. On the contrary, magnetic resonance imaging (MRI) is a non-radiation based examination method, which gains an increasing popularity in research settings in healthy subjects, which is especially relevant for epidemiological studies. For example, in the frames of the Study of Health in Pomerania (SHIP) [8,30] whole-body MR imaging is taken. One of the research tasks of this project is to establish populationbased MRI reference values for lung volumes and to correlate MR
∗ Corresponding author. Tel.: +49 3834867554. E-mail addresses: [email protected], [email protected] (T. Ivanovska), [email protected] (K. Hegenscheid), [email protected] (R. Laqua), [email protected] (J.-P. Kühn), [email protected] (S. Gläser), [email protected] (R. Ewert), [email protected] (N. Hosten), [email protected] (R. Puls), [email protected] (H. Völzke). 0895-6111/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compmedimag.2011.10.001
findings with non-invasive pulmonary function tests, which assess static or dynamic lung volumes. The obtained lung volumes will be associated with such factors as the influence of genetic variations, endocrine metabolism, cardiac or pulmonary diseases, systemic inflammation, and mortality [21,7,23,24,6]. The results of dynamic pulmonary tests may be significantly influenced by several factors, whereas manual segmentation of MR images is a very laborious, observer-dependent, and time-consuming process, which is prone to errors. Therefore, fast, reliable, and reproducible automatic segmentation methods are required. However, lung imaging is a challenge in MRI because of the predominance of air within the lungs and associated susceptibility issues as well as the lower spatial and depth resolution when compared to the CT imaging and low signal to noise ratio of the inflated lung parenchyma. Moreover, in CT images pulmonary regions have the highest contrast when compared to other tissues, whereas in MR images all tissues, such as bones, muscles, fat are clearly distinguished. In this paper, a fast automatic method for segmenting the lung in MR images is presented. The method includes four main steps: First, an initial segmentation is performed, where the both lungs and trachea are extracted. The vessels are excluded from the parenchyma. Then, the trachea is eliminated and lungs are separated. Finally,
282
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293
the holes and cavities of the lung tissues are closed. Each step of the pipeline is thoroughly described and an analogy between the current approach and other methods is drawn. The paper is organized as follows. We give an overview of the existing methods in Section 2 and describe the pulmonary MR data created in SHIP in Section 3. Then we present our pipeline in Section 4 and compare the results against the manually established ground truth in Section 5. The findings are discussed in Section 6. Section 7 concludes the paper.
2. Related work In the literature, numerous segmentation strategies for CT data exist, whereas only few methods for MR data are addressed. Lelieveldt et al. developed a model-based method that allowed for extraction of lungs as well as other nearby organs such as the heart and diaphragm from MR data with different projection planes [14]. The authors reported that their method was designed to find the approximate borders to initialize subsequent more accurate local segmentation process. However, for our purposes it is enough to extract the lungs directly without the model-based presegmentation. Although numerous methods for CT data can be applied to MR images, there appear a number of difficulties which are specific for this imaging modality. In the following, we describe the main directions for lung segmentation from CT images and explain what problems arise on the correspondent processing steps with the MR data. In pulmonary CT images, the natural contrast between the low-density lungs and the surrounding high-density chest wall is usually used to guide the segmentation process. A standard approach applied by many authors [2,10,15,31,5] is binary thresholding. Here, either optimal Otsu thresholding [20] or iterative thresholding is applied [10]. In some previous works [11], manual thresholding is utilized. Shojaii et al. proposed to apply markerbased Watershed segmentation [29] with the defined internal and external markers [25]. In the case of MRI, the binary thresholding is not sufficient due to the different contrast for bones, air, and soft tissues. Thereafter, the lungs are separated from the large airways. The task of the trachea extraction is considered in many approaches [18,25] for CT data and it is solved mainly by either region growing [27] or morphological operations [26]. The start of the trachea is automatically identified by searching for a large, circular, air-filled region near the center of the first few slices in the data set. Regions in the current slice provide potential seed points positions for the next slice. The slice-by-slice growing procedure terminates when the size of the region on a new slice increases dramatically, indicating that the airways have merged into the low-density lung tissue [10]. Due to the lower spatial resolution and signal to noise ratio in MR images, such a criterion is not fully sufficient and needs to be extended. We introduce a procedure that extends this approach and overcomes the problem. Then the lungs need to be separated from each other. One of the most popular approaches is to apply morphological operations to detect the region of interest and then a 2D maximum cost path algorithm [27] is used to follow the junction line. For example, Hu et al. [10] and Li et al. [15] applied erosion and dilation to separate the lungs and then to restore the boundary as much as possible. A search region is found on a 2D slice, and it is propagated to successive slices due to the fact that the junction line position varies slowly through the data set. Li et al. [15] proposed to use the geodesic dilation, which is faster than the conditional dilation. The geodesic dilation terminates when no pixels can be added without changing the connectivity of the lung regions [15]. The smoothness
of the pulmonary anatomy in MRI is not guaranteed due to the lower spatial resolution. Hence, the 2D based method can result in a rather coarse boundary, and the application of a procedure that produces a 3D boundary is more beneficial in terms of accuracy. For example, Kuhnigk et al. [13] applied a 3D marker-based watershed transform with one automatically generated seed point per lung to obtain separate lung masks [13]. However, such a method will be more computationally expensive. We use an optimization step that allows for finding the boundary in 3D and at the same time keeping the computational costs relatively low. For smoothing and filling the cavities, practically in all previous works [15,10,25,1] a morphological closing method is applied to the lung boundary. This method is also referred to as the rolling ball filter. We developed a processing pipeline for MR data, which allows us to overcome the above-mentioned difficulties and leads to a fast and reliable segmentation, which is extremely important when thousands of datasets are processed.
3. Materials The design of the SHIP project is given elsewhere [8]. In brief, SHIP is a population-based project consisting of the two independent studies SHIP and SHIP-Trend. The examination programme of the 11 year examination follow-up of the original SHIP (SHIP-2) and SHIP-Trend was supplemented with whole-body MR imaging [30,8]. All MRI studies are performed on a 1.5-T MR imager (Magnetom Avanto; Siemens Medical Systems, Erlangen, Germany). Subjects are placed in the supine position, and five phased array surface coils are placed on the head, neck, abdomen, pelvis, and lower extremities. A spinal coil is embedded in the patient table. Each participant undergoes standardized WB-MRI including an axial T1-weighted VIBE (volume-interpolated breath hold examination) sequence (typical parameters: repetition time (TR) 31 ms, echo time (TE) 11 ms, flip angle 8◦ , voxel size 1.8 × 1.8 × 3 mm) and an axial T2-weighted HASTE (half-Fourier single-shot turbo spinecho) sequence (TR 5500 ms, TE 220 ms, flip angle 150◦ , voxel size 2.3 × 1.8 × 5 mm) for chest imaging. To avoid motion artifacts 21 s of breath hold are needed for the VIBE sequence. The HASTE sequence is taken in two steps and then the parts are composed together with the navigator control. In total, 40 s are needed for the HASTE sequence. Because of the small slice thickness (3 mm) and bigger spatial resolution when compared to the HASTE sequences the automatic lung segmentation and volumetry is performed using the VIBE sequences. For our experiments we used 10 randomly selected datasets of healthy participants from the SHIP-Trend programme. Each dataset consists of 88 slices with a resolution of 512 × 416.
4. Methods Our method consists of four steps: the main extraction and three refinement steps. First, the lung and trachea excluding the main vessels are extracted. Second, the trachea is separated from the lung tissue, and the cut is made on the main bronchi entering the lung parenchyma. Due to the fact that a participant holds the breath during the scan, which results in a very thin and faint boundary between the lungs, the main extraction step may fail to detect this boundary on certain slices. Third, we apply the lung separation procedure. Fourth, the lung boundaries are smoothed and the cavities in the lung tissue are filled, which allows us to obtain the completed lung parenchyma. In each processing step, we compute the volume of the segmented tissue.
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293
283
Fig. 1. Lungs and trachea segmented with K-Means.
4.1. Lung and trachea segmentation MR technique allows for better contrast within the soft tissues, which are rather distinctively separated in terms of intensity, so we assumed that any automatic clustering technique could be applied here. One of the classical machine learning methods, K-Means [17], for example, is a good candidate due to its speed and simplicity. KMeans aims to separate n observations (x1 , x2 , . . ., xn ) into k classes (k ≤ n) S = S1 , S2 , . . ., Sk so as to minimize the within-cluster sum of squares: E = argmin S
k
xj − i 2 ,
i=1 xj ∈Si
where i is the mean of the points in Si . We experiment with the number of clusters K and determine that setting K = 4 or 5 is enough to separate the lungs from the vessels, which are brighter than the lung parenchyma, but much darker than the other tissues. Such a value of K was set automatically for all datasets. It is known that K-Means method results are strongly dependent on the initial cluster selection [4]. To avoid the inadequate clustering, e.g., the case when the cluster centers fall too close together, we set the initial cluster centers as far as possible from each other and in each step of the algorithm iterative procedure, if a cluster contains less than 2% of all data voxels, we place its center to the mean of the two other centers accounting for the most voxels. In Fig. 1, the result of the segmentation is shown.
components, say, Ci = {c1i , . . . , cKi }, have been detected there. We stop at the slice i, if a connected component cˆ ∈ Ci , which satisfies two conditions, is found: (1) {∀c ∈ Ci , A(ˆc ) ≥ ATr }, (2) {∀c ∈ Ci , |xcˆ − 0.5W | = min(|xc − 0.5W |)},
(1)
where A(c) defines the area of the connected component c in pixels, W is the slice width, xc is the x coordinate of the center of mass of the connected component c, and ATr is the area threshold. The largest connected component, whose center of mass lies as close as possible to the slice center in the x direction, is selected as the initial trachea region. Thereafter the trachea is located. In Fig. 3, three consecutive slices are shown. The segmented tracheal regions are depicted in green. We traverse the trachea along the z axis analysing each time a pair of consecutive slices. Let us assume the trachea is already detected in slice i and is to be determined in slice j = i − 1. In slice i the trachea consists of K connected components, say, Ci = {c1i , . . . , cKi }, and their cumulative area is A(Ci ). The set of connected components Ci is orthogonally projected to slice j and M j j connected components Cj = {c1 , . . . , cM } are selected if the intersection with the projected connected component set Ci is non-zero: j j ∃cki , ∃cm : cki ∩ cm = / , where k ∈ [1, . . ., K] and m ∈ [1, . . ., M]. If a cumulative area A(Ci ) or A(Cj ) is smaller than a threshold ATr , then the processing stops. Otherwise, we compute for all pairs of interj j secting components (cki , cm ), where cki ∈ Ci and cm ∈ Cj , the area ratios: j
4.2. Trachea and main bronchi extraction
R
The extraction task can be divided into two parts: first, the start component of the trachea is located. Second, the trachea and main bronchi are found using a 3D component analysis procedure. Under “connectedness” we understand 8- and 26-neighborhood for 2D and 3D, respectively. The cut is made on the spots, where the main bronchi enter the lung tissue. Here, we would like to explain both steps in detail. In Fig. 2, the 3D and 2D views of the trachea initialization are shown. Starting from the slices, which anatomically lie in the direction of the larynx, we process each slice and find its connected components. Let us denote the processed slice i, and K connected
If R
=
j (c i ,cm ) k
j
(c i ,cm )
A(cm ) A(cki )
.
(2)
is greater than a preselected threshold Rconst , we assume
k
j
that the trachea component cm has merged with the lung tissue and a separation procedure is called for it. To separate the airways from the lung tissue, we compute a distance map [27] D. The distance map is the result of the distance transform. The distance transform (DT) is the transformation that generates a map D whose value in each pixel p is the smallest distance from this pixel to background Bkg: D(p) = min{d(p, q)|q ∈ Bkg},
(3)
284
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293
Fig. 2. The start trachea component is marked red and shown in 3D (left) and 2D (right). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
where d is the Eucledian distance, p is a foreground pixel, and q is a background pixel. Then the distance map is inverted: Dinv = − D, and the 2D Watershed [29,22] is applied to Dinv . The Watershed result is a set of the regions separated by black boundaries. Each region has a unique ID. The regions that have a non-zero intersection with the projected trachea component cki are taken and the ratios are calculated again. If the Watershed transform failed to detect the sought boundary, i.e., the correspondent ratio is still greater than the threshold, then the boundary is created artificially by projecting the trachea region from the current slice i. We permit at most two consecutive artificial boundaries, so that a “fake” trachea inside the lung tissue is not created. Thereafter, the same procedure is repeated for the subsequent slice pairs till i ≥ 1. This procedure allows for separation of the trachea and main bronchi from the lung tissue, and the cut is approximately made on the spot, where the main bronchi enter the lung tissue. When the extraction is finished, we exclude the trachea and main bronchi mask from the further processing. Any small regions between two lungs, which can appear after the extraction and represent noise, are also excluded. In Fig. 4, a 3D view of the extracted trachea is shown.
4.3. Lung separation After the trachea extraction, we analyse the remaining connected components in 3D and exclude all components with area, which is smaller than a selected threshold A3D . This results either in one (the lungs have merged together) or two (the lungs are separate) connected components. If the lungs represent one connected component, they need to be separated. Our approach allows for detection of a boundary in 3D on the slices, where the lungs have merged together. It consists of the following steps. First, the slices, where the lungs merged together are detected. Second, a subvolume including the Volume of Interest (VOI) is selected. Third, the distance map [27] is built and corrected for the computational costs reduction. Finally, 3D Watershed [29] is applied to the corrected and inverted distance map, and two separated lung masks are built. Here, we describe the steps of the procedure in more detail. There are two possible lung junctions: anterior and posterior. Since the processing of both cases is identical, we describe it only once using the anterior junction example. In Fig. 5, an example of the anterior junction on a 2D slice is shown. First, we analyse a 3D dataset and select the slices that contain only one connected component. Then, a subvolume is constructed from this list including all intermediate slices and two additional slices for the beginning and the end, e.g., if the slice list consists of slice ids 38, 41, 42, the selected subvolume will include the slices that have numbers 37 → 43. This subvolume can be separated into two subvolumes: 37 → 39 and 40 → 43, which allows to reduce the computational costs in the Watershed step. For simplicity, we refer to one (total) subvolume in the text. We reindex the subvolume slice numbers and set them 1, . . ., N. Let V = (S1 , S2 , . . ., SN ) denote the voxel sets (slices) of our subvolume. Thereafter, a new subvolume V = , . . . , S is formed, so that the first and the last slices S 2 N−1 mid are projected on the whole subvolume and the coinciding voxels are masked: Si = Si \ (S1 ∪ SN ) ,
Fig. 3. Three consecutive slice with the segmented large airway regions.
where i ∈ [2, . . . , N − 1]. Such a masking allows one to exclude the coinciding voxels. However, it contains the voxels, which correspond to the junction, since there is no junction in the first and last slices. Then, we detect 3D connected components. The connected component cj that represents the anterior junc tion is visible in the whole subvolume V mid . The junction lies in
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293
285
Fig. 4. The extracted Trachea from Fig. 1 is shown in an overlaid manner.
the upper third of each slice and is close to the slice center in x direction: (1) {∀c ∈ V mid , ycj ≤ 0.3H, (2) {∀c ∈ V mid , |xcj − 0.5W | = min(|xc − 0.5W |)}, where xc and yc are the coordinates of the center of mass of the connected component c, W is the slice width, and H is the slice height. We build a bounding box V bb around the selected connected component to assure that the junction is included. The bounding box V bb include the voxels that satisfy: {∀v ∈ V : lxcj − f1 ≤ v · x ≤ rxcj + f1 , tycj − f2 ≤ v · y ≤ bycj + f2 }, where lx, rx, ty, by are the left and right x coordinates and top and bottom y coordinates of the connected component cj , correspondingly. Moreover, f1 and f2 are preselected factor values. Thereafter, a 3D distance map [27] is computed and inverted. The minimal distance map values are the “deepest” points, lying inside the lung tissue. The junction does not lie deep and has high values on the inverted distance map. Hence, the distance values can be thresholded. Let Dinv denote the inverted 3D distance map and the thresholded distance map in each voxel p is defined as D inv (p) = max(Dinv (p), const),
(4)
where const is a constant. It reduces the computational costs of the Watershed transformation, which starts from the minima of the distance map. The constant value const is computed from the distance values of the connected component cj and the distance map minima. Thereafter, a 3D Watershed [29] is applied to the cor rected distance transformation D inv . In Fig. 6, a close up view of the watershed result where the connected lung parts are successfully separated is shown.
Fig. 5. A 2D slice from the dataset shown in Fig. 1 with the anterior junction.
in an overlaid manner. One can observe a smooth lung boundary with the filled cavities. In Fig. 8, the final lung volumes are overlaid with the initial dataset.
4.4. Filling lung cavities and holes The final processing step is to fill the cavities in the lung boundaries, obtain the complete lung masks, and compute the final lung volumes. First, the holes in the lung tissue are closed. Thereafter, we apply a morphological closing procedure to the lung boundary to cover the boundary cavities. A spherical 2D kernel of radius 9 is used. Each lung is processed separately to avoid artificial lung junctions. In Fig. 7, the resulting mask for the slice from Fig. 5 is shown
Fig. 6. A close up view of the separation result for the 2D slice shown in Fig. 5.
286
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293
Fig. 7. The smoothed lung tissue for the 2D slice shown in Fig. 5 is presented. The left and right lungs are shown in red and green colors, correspondingly. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
5. Results We tested our processing pipeline and extracted the volumes on two stages: the total volume of both lungs and trachea excluding vessels (after the initial segmentation) and the separate volumes of both lungs with the filled holes and closed cavities (after the final step) are computed. The initial step takes around 1 min and the whole pipeline takes around 3 min on Intel(R) Core i7 − 950 (Quadcore) @ 3.06 GHz computer with 6 Gb DDR3 for one dataset with a resolution 512 × 416 × 88. All sequences have been processed automatically without an additional parameter tuning for each participant. The following parameter values have been set: K = 4 is the number of classes for segmentation, ATr = 50 is the minimum area of a 2D trachea connected component, Rconst = 3 is the maximum area factor between the connected component areas on consecutive slices. Moreover, A3D = 5000 is the minimum size of a 3D lung connected component, lungMinThresh2D = 500 is the minimum size of a 2D lung connected component, f1 = f2 = 50 are the padding values for the bounding box, Kr = 9 is the kernel radius for the rolling ball filter. The ground truth is determined by manual identification of the lung boundaries. In total, we obtained manual readings from two experts who traced the lung boundaries. One image analyst is a
Fig. 8. The final lung volumes are overlaid with the initial dataset. The left and right lungs are shown in red and green colors, correspondingly. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
Fig. 9. A 2D slice from dataset 6 is shown. The region close to mediastinum defines the difference between the expert readings for the left lung mask.
radiologist (R.L.), and he determined all types of masks on 880 slices. The other expert is a very experienced pulmonologist (R.E.), and he determined two types of masks (separate left and right lungs) on 880 slices. The image analysts worked independently of each other. We have evaluated the difference between the manually delineated volumes and observed that the expert agreed on around 95% of the volumes, i.e., the average proportions of the overlapping voxels in both volumes are 93.87% with the standard deviation 1.52 for the left lung masks and 94.76% with the standard deviation 1.01 for the right lung masks. The correspondent mean deviations for the detected volumes are 3.89% for the left lung masks and 3.46% for the right lung masks. In Fig. 9, an example slice from dataset 6 is shown. The region of the left lung, which is close to the mediastinum is depicted in a close-up view. The region of expert variations is colored in green and red. Fig. 10 shows the expert readings in an overlaid fashion. As it can be observed, the experts marked the region near the mediastinum differently, which results in the variation interval (0, 6.5]% for the complete masks. The result from one the experts is compared with the lung masks detected by our pipeline. To evaluate the accuracy of the method, we applied a methodology described by Udupa et al. [28]. We computed the following measures: False Negative Volume Fraction (FNVF), which indicates the fraction of tissue, which was manually marked but missed by the automatic method; False Positive Volume Fraction (FPVF), which indicates the fraction of the tissue that was marked by the method but not marked but the expert; True Positive Volume Fraction (TPVF) is the number of overlapping voxels marked both by the method and the expert. Additionally we calculated the True Positive Volume Error (TPVE), which shows the deviations between the automatically detected volume and the manually delineated one. However, we have to mention that this measure does not regard the exact boundaries. We computed it, as our ultimate goal is to assess the volume of the lungs. Let Cauto and Cexp denote the binary masks produced by our pipeline and the expert delineation, correspondingly. Moreover, let Volauto and Volexp to be the volumes computed by the pipeline and
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293
287
Fig. 10. The left lung masks from two experts for the slice from Fig. 9 are shown in green and red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
the expert, correspondingly, then the FNVF, FPVF, TPVF, and TPVE are defined as follows:
Cexp − Cauto Cauto − Cexp , , FPVF = FNVF = Cexp Cexp Cauto ∩ Cexp Volauto − Volexp , TPVE = TPVF = . Cexp Volexp
Our evaluations are presented in Tables 1, 2, and 3 for the total volume of both lung with trachea, the volume of the right lung, and the volume of the left lung, correspondingly. The following i notations for the headings are used: “Volexp ” is the expert-detected volume detected by expert i; “Volauto ” means the automatically detected volume; Cexp , “FNVF”, “FPVF”, “TPVF”, “TPVE” are defined as described above.
Table 1 Comparison of the automatic results to the expert-defined ground truth for ten different data sets. Lungs and trachea without vessels. Here, the manual results are produced by one expert (radiologist). Lungs and trachea without vessels 1
ID
Cexp (voxel)
FNVF (%)
FPVF (%)
TPVF (%)
Volexp (ml)
Volauto (ml)
TPVE (%)
1 2 3 4 5 6 7 8 9 10 Mean StDev
1,068,926 1,423,881 1,304,202 1,690,847 2,035,239 1,275,871 1,503,466 1,311,256 1,609,496 1,550,048
3.41 2.04 5.56 1.95 3.77 8.53 1.74 0.36 2.83 2.32 3.25 2.32
0.08 0.42 0.48 4.27 1.18 0.47 1.1 1.81 1.51 0.8 1.21 1.2
96.59 97.96 94.44 98.05 96.23 91.47 98.26 99.64 97.17 97.68 96.75 2.32
2477.16 3299.74 3022.39 3918.42 4716.52 2956.74 3484.17 3038.74 3729.38 3592.13
2394.74 3246.2 2868.71 4009.72 4594.69 2718.44 3461.78 3082.79 3680.55 3537.42
3.33 1.62 5.08 2.33 2.58 8.06 0.64 1.45 1.31 1.52 2.79 2.24
288
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293
Table 2 Comparison of the automatic results to the expert-defined ground truth for ten different data sets. Right lung with vessels and filled holes. Here, the manual results are produced by one expert (radiologist). Right lung ID
Cexp (voxel)
FNVF (%)
FPVF (%)
TPVF (%)
1 2 3 4 5 6 7 8 9 10 Mean StDev
608,632 851,018 833,722 963,975 1,141,684 704,715 837,662 706,295 935,049 885,919
2.32 1.55 4.23 0.89 2.7 5.8 2.03 1.09 4.08 1.87 2.66 1.57
7.74 2.31 3.07 4.98 2.58 2.87 4.39 5.37 2.84 3.42 3.96 1.69
97.68 98.45 95.77 99.11 97.3 94.2 97.83 98.91 95.92 98.13 97.33 1.57
1
Volexp (ml)
Volauto (ml)
1410.46 1972.17 1932.09 2233.94 2645.77 1633.13 1941.22 1636.79 2166.91 2053.05
1486.82 1986.99 1909.84 2325.4 2642.72 1585.18 1984.36 1706.77 2140.16 2084.89
TPVE (%) 5.41 0.75 1.15 4.09 0.12 2.94 2.22 4.28 1.23 1.55 2.37 1.74
Table 3 Comparison of the automatic results to the expert-defined ground truth for ten different data sets. Left lung with vessels and filled holes. Here, the manual results are produced by one expert (radiologist). Left lung 1
ID
Cexp (voxel)
FNVF (%)
FPVF (%)
TPVF (%)
Volexp (ml)
Volauto (ml)
TPVE (%)
1 2 3 4 5 6 7 8 9 10 Mean StDev
507,028 636,388 552,772 766,651 973,428 635,371 704,294 668,718 774,810 709,396
3.35 2.58 5.32 2.56 2.15 5.12 1.34 1.35 3.02 2.55 2.94 1.36
2.16 2.94 2.69 5.61 2.91 1.72 2.73 3.55 3.3 2.06 2.97 1.08
96.65 97.42 94.68 97.44 97.85 94.88 98.66 98.65 96.98 97.45 97.06 1.36
1175 1474.78 1281.01 1776.66 2255.85 1472.43 1632.15 1549.71 1795.57 1643.97
1160.96 1480.12 1247.3 1830.8 2272.94 1422.35 1654.75 1583.81 1800.67 1635.9
1.19 0.36 2.63 3.05 0.76 3.4 1.38 2.2 0.28 0.49 1.58 1.16
6. Analysis and discussion 6.1. Comparison to other techniques
we will compare our results to the well-known CT segmentation techniques in a step-by-step manner and explain the benefits of our approach.
One of the most popular directions for segmentation of lung CT data are pixel-based methods [10,15,13,9]. Here, main steps are thresholding, extraction of large airways, lung separation, and smoothing. Our method belongs to the same group and the pipeline steps are adapted for MRI data. In the following section,
6.1.1. Thresholding A thoracic CT scan contains two main groups of pixels: highintensity pixels, which are located in the body, and low-intensity pixels that are in the lung and surrounding air [9]. There is a huge difference between the pixel values in these groups, so the
Fig. 11. A slice from the initial dataset (left) and result for Otsu thresholding (right) are shown side-by-side. Otsu thresholding fails to separate the lungs from the other tissues.
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293
image is thresholded. Here, the simplest thresholding at a fixed value often suffices [12,9]. Another possibility is to apply a binary automatic thresholding, such as Otsu thresholding [20] or the optimal thresholding [10] proposed by Hu et al. Since MR technique allows for a better contrast within the soft tissues, the binary thresholding, which is usually applied on high resolution CT data [10,15,13], would fail in separating lungs and airways from other tissues, which is documented in Fig. 11. As it can be observed in Fig. 12, the MR thorax data consists of three main intensity classes: muscles (gray), bones and fat (light gray), lungs and trachea (dark), which are well separated from each other. We applied two classical approaches, namely Multiclass Otsu thresholding [16] and K-Means [17] clustering. Then we extracted the largest connected component, which represents the lungs, trachea, and main bronchi and compared the results. In Fig. 13, the segmentation of both lungs and large airways excluding the vessels is shown in a 2D slice. The results were comparable in quality: TPKMeans = 97.95%, TPOtsu = 97.4%, but K-Means had lower computational costs.
289
Fig. 12. An slice from a lung dataset.
Fig. 13. A close up view of the segmentation results for a slice shown in Fig. 11 obtained with K-Means (left) and Otsu (right). The results are similar in quality, but K-Means is superior in terms of speed.
Fig. 14. ROIs from three consecutive slices. The trachea component, which has not reached the carina yet, is connected on one slice (middle) and then disconnected again (right).
290
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293
6.1.2. Extraction of large airways To analyse the lung volumes, the trachea and large airways must be detected and separated from the left and right lungs. Since the airways are empty cavities, the intensity of pixels in the area is low. In CT images, many techniques can be applied to extract the trachea and large airways. Heuberger et al. reported that a predefined threshold often suffices [9]. To remove the trachea and large airways, the authors searched for an area with an average intensity lower than 0.5Tc , where Tc is a threshold value, computed by the optimal thresholding method [10]. Such an approach is definitely not applicable for the MR data, since the trachea and large airways have practically the same intensity and noise level as the lungs. Li et al. treated the trachea and bronchi as small round regions, which were removed using a simple morphology technique [15]. First, all connected regions (black and white one) on each slice were labeled. Then, the regions with the area smaller than a predetermined threshold were inverted, i.e. operator NOT was applied to each pixel inside the regions. Hu et al. proposed to use a slice-byslice region growing procedure, which stops, when the size of the region on a new slice increases dramatically, indicating that the airways have merged into the low-density lung tissue [10]. These two approaches have some disadvantages. According to the procedure, described by Hu et al. [10], only the separated trachea parts will be erased before the first connection to the lung tissue is reached. The approach proposed by Li et al. [15] will delete the separated trachea regions in the whole volume. However, on the slices, where the trachea is connected to the lung, no separation can be done. For the MR data, which have lower spatial resolution and signal-to-noise ratio, the result of these approaches will add a significant number of voxels, that belong to the large airways, to the lung volumes. For example, in Figs. 14 and 15 the extracted lung mask is overlaid with the original data and two groups of regions of interest (ROIs) that contain the trachea and main bronchi are shown. In Fig. 14, the trachea has not reached the carina yet, but is merged with the lung tissue and separated again on the consecutive slice. Here, Hu’s procedure would stop on the “merged” spot and the rest of the large airways (including the trachea, carina and main bronchi) would be referred to the lung volume. Li’s approach will exclude the trachea in the left and right slices, but not in the middle one. In Fig. 15, the trachea is separated into the main bronchi that can be at least partially excluded from the lung tissue, which can be accomplished by neither of three above mentioned methods. Our approach is also based on a region growing procedure, but it processes the whole trachea and main bronchi and separates the merged parts from the lung tissue. This allows us to extract the significant parts of the large airways and exclude them from the lung volumes.
6.1.3. Lung separation Finally, the two lungs are identified and, if needed, separated. Several methods have been proposed for CT images. Brown et al.
Fig. 15. ROIs from three consecutive slices. The trachea component reaches the carina and splits into two bronchi components (top). The right main bronchus is connected to the lung tissue (middle and bottom). The connection spot is marked with red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
proposed to apply a 2D maximum cost path algorithm [27], where the cost of a pixel is its graylevel, can be used to follow the anterior junction line [3]. Hu et al. and Li et al. utilized the same approach, but applied different procedures for detection of the connection ROI [10,15]. The dynamic programming approach is applicable in the case of CT data, since the soft tissue of the junction has a slightly higher gray level than the lungs, and forms a continuous line in each axial slice of the connecting region [3]. These two criteria are usually not fulfilled in the case of the MR data due to the lower spatial resolution and partial volume effects. It is documented by Fig. 16, where ROIs from three consecutive slices are shown. The junction is not always continuous here and junction voxels have the same intensity values as the lung tissue. Moreover, a 2D based approach may lead to the result, when the lungs will overlap in 3D, being separated in 2D.
Fig. 16. ROIs from three consecutive slices with the anterior junction marked red. The junction voxels have practically the same intensity as the lung tissue voxels. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293
291
Fig. 17. Automatically detected result (left) consists of less parenchyma (more vessels are excluded) than the expert-defined result (right) for a slice shown in Fig. 12.
We utilize an 3D marker based Watershed transform, similarly to Kuhnigk et al. [13]. The appropriate selection of the ROI and markers allows us to reduce the computational costs for the Watershed transform (cf. Section 4). 6.2. Quality analysis As it can be taken in Table 1, our automatic procedure produces the results with high TPVF (96.75%) and rather low FNVF and FPVF (3.25% and 1.21%, correspondingly). However, it is obvious that FNVF is higher. The expert has to manually trace not only the lung and trachea boundaries, but also cross out the vessels, which are brighter than the lung tissue. This procedure is very laborious and time-consuming. The automatic method definitely detects more bright vessels than the expert does. Fig. 17 shows a 2D slice from dataset 6, which is an outlier among our test datasets with the FNVF equal to 8.53%. The TPVE is equal to 2.79%, which is exact enough for our research purposes. The findings presented in Tables 2 and 3 also show a high agreement with the expert reading. The TPVF for the separate right and left lungs (97.06% and 97.33%, correspondingly) is higher than in the case of the detection of the lungs and trachea, where the vessels have to be manually excluded. However, closing the vessels also increases FPVF (3.96% and 2.97% for the right and left lungs, correspondingly). In Fig. 18, a 2D example of a right lung mask from dataset 1 demonstrating this finding is shown. Since our automatic procedure does not distinguish between the different types of holes, i.e., whether it is a hole where the vessel is located or a hole, which appears during the scanning, some “fake” lung voxels can be added to the mask. Nevertheless, the achieved accuracy in True Positive Volume is enough and lies within the variation interval between the experts. Moreover, TPVE is equal 2.37% for the right lung and 1.58% for the left lung. These errors are not significant and can be neglected for our research tasks. 6.3. Performance analysis We analyse the performance of our approach using the test datasets. As it can be seen in Fig. 19, the total execution times lies around 150 CPU seconds. However, there are several outliers, namely, datasets 4, 7, 10. The maximum execution time reaches 452.27 for dataset 4. We selected eight main steps in our algorithmic pipeline (see Table 4) and analysed the performance for each of them. The results are documented in Fig. 20. The most expensive step with the biggest dispersion is the Lung separation, which includes a marker based 3D Watershed algorithm on a selected ROI.
Fig. 18. A slice from dataset 1 (top). Automatically detected result (bottom) covers both vessel locations and a “fake hole”, which appears due to the scanning, whereas the expert marks only the vessel holes (middle).
292
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293
300
we implemented an application, which includes a batch processing mode. The application is developed under MeVisLab framework [19]. In total, the application can be used in 2 modes: single-file step-by-step processing and automatic batch processing modes. The former mode is primarily designed for parameter selection and evaluation. With the latter mode the application can run in the background over thousands of participants using the same parameter settings for all of them.
250
7. Conclusions
200
We presented an approach for automatic lung extraction from MR data, which consists of one main step and three refinement procedures. For our research task, the volumes of both lungs with trachea and excluded vessels and separate full lung masks are of a particular interest. The accuracy and speed of the proposed method allow for using it in the clinical routine to evaluate thousands of participants and provide the further research background in such epidemiological studies as the SHIP Project.
500 Execution time 450
CPU sec
400 350
150 100
0
1
2
3
4
5
6
7
8
9
10
11
Dataset # Fig. 19. The total execution time for 10 test dataset in CPU seconds. The average time is lower than 3 min, if no lung separation is done.
Acknowledgement
Table 4 Algorithm steps description. Step
Description
Min time (CPU sec)
Max time (CPU sec)
1 2 3 4 5 6 7 8
Smoothing Clustering Save result Exclude trachea Separate lungs Save result Fill holes Compute volumes, save result
51.92 9.95 0.99 2.5 12.33 3.4 41.29 8.01
52.8 13.91 1.09 5.8 311.92 3.66 59.62 8.27
References
Here, the average value is 65.173 CPU seconds and the standard deviation equals 104.7. The performance of the algorithm strongly depends on the ROI size. We assume that the performance can be improved if a more efficient algorithm implementation is used. We leave this improvement for future work. 6.4. Application for medical research The algorithm in its current state can be applied in an epidemiological setting. To be conveniently embedded for medical research, 200
Execution time for each step
CPU sec
150
100
50
0
-50
1
2
3
4
5
6
7
The authors would like to thank Prof. Dr. Volkmar Liebscher (Institute of Mathematics and Informatics, Ernst-Moritz-Arndt University, Greifswald, Germany) for support and valuable discussions.
8
Pipeline steps Fig. 20. Mean and standard deviation for each pipeline step are shown using error bars. The most expensive step with the biggest dispersion is the lung separation. This steps depends on how many slices it must be done.
[1] Armato SG, Giger ML, Moran CJ, Blackburn JT, Doi K, MacMahon H. Computerized detection of pulmonary nodules on CT scans. Imaging and Therapeutic Technology 1999. [2] Armato SG, MacMahon H. Automated lung segmentation and computeraided diagnosis for thoracic ct scans. International Congress Series 2003;1256:977–82. [3] Brown MS, McNitt-Gray MF, Mankovich NJ, Goldin JG, Hiller J, Wilson LS, et al. Method for segmenting chest ct data using an anatomical model: preliminary results. IEEE Transactions on Medical Imaging 1997;16(6):828–39. [4] Duda RO, Hart PE, Stork DG. Pattern classification. Wiley; 2000. [5] El-Baz A, Farag AA, Falk R, Rocca RL. A unified approach for detection, visualization, and identification of lung abnormalities in chest spiral ct scans. International Congress Series 2003;1256:998–1004, cARS 2003. Computer Assisted Radiology and Surgery. Proceedings of the 17th International Congress and Exhibition. [6] Gläser S, Ittermann T, Koch B, Völzke H, Wallaschofski H, Nauck M, et al. Airflow limitation, lung volumes and systemic inflammation in a general population. European Respiratory Journal 2011, Jun 30 [Epub ahead of print]. [7] Gläser S, Friedrich N, Ewert R, Schäper C, Nauck M, Dörr M, et al. Association between serum insulin-like growth factor (igf) i and IGF binding protein-3 and lung function. Journal of Clinical Endocrinology and Metabolism 2009;94(7):2452–8. [8] Hegenscheid K, Kühn J, Völzke H, Biffar R, Hosten N, Puls R. Whole-body magnetic resonance imaging of healthy volunteers: pilot study results from the population-based ship study. Rofo 2009. [9] Heuberger J, Geissbuehler A, Mueller H. Lung ct segmentation for image retrieval using the insight toolkit. In: Proceedings of the 1st International Conference on Medical Imaging and Telemedicine (MIT ‘05). 2005. [10] Hu S, Hoffman EA, Reinhardt JM. Automatic lung segmentation for accurate quantitation of volumetric X-ray CT images. IEEE Transactions on Medical Imaging 2001;20:490–8. [11] Kalender WA, Fichte H, Bautz W, Skalej M. Semiautomatic evaluation procedures for quantitative ct of the lung. Journal of Computer Assisted Tomography; 1991. [12] Kemerink GJ, Lamers RJ, Pellis BJ, Kruize HH, van Engelshoven JM. On segmentation of lung parenchyma in quantitative computed tomography of the lung. Medical Physics 1998;25(12):2432–9. [13] Kuhnigk J, Hahn HK, Hindennach M, Dicken V, Krass S, Peitgen H. Lung lobe segmentation by anatomy-guided 3d watershed transform. In: Proceedings of SPIE, vol. 25032. 2005. p. 1270–3. [14] Lelieveldt BPF, Sonka M, Bolinger L, Scholz TD, Kayser HWM, et al. Anatomical modeling with fuzzy implicit surface templates: application to automated localization of the heart and lungs in thoracic mr volumes. Computer Vision and Image Understanding 2000:1–20. [15] Li W, Nie SD, Cheng JJ. A fast automatic method of lung segmentation in CT images using mathematical morphology. Berlin, Heidelberg: Springer; 2007, 2419–2422. [16] Liao P, Chen T, Chung P. A fast algorithm for multilevel thresholding. Journal of Information Science and Engineering 2001;17(5):713–27.
T. Ivanovska et al. / Computerized Medical Imaging and Graphics 36 (2012) 281–293 [17] MacQueen JB. Some methods for classification and analysis of multivariate observations. In: Proceedings of 5-th Berkeley Symposium on Mathematical Statistics and Probability, vol. 1. 1967. p. 281–97. [18] Memon NA, Mirza AM, Gilani S. Limitations of lung segmentation techniques, vol. 27. US: Springer; 2009. p. 753–66, chapter 76. [19] MeVisLab. Software for medical image processing and visualization; 2011. http://www.mevislab.de. [20] Otsu N. A threshold selection method from gray-level histograms. IEEE Transactions on Systems, Man, and Cybernetics 1979;9(1):62–6. [21] Repapi E, Sayers I, Wain LV, Burton PR, Johnson T, Obeidat M, et al. Genomewide association study identifies five loci associated with lung function. Nature Genetics 2010;42:36–44. [22] Roerdink JBTM, Meijster A. The watershed transform: definitions, algorithms and parallelization strategies. Fundamenta Informaticae 2001;41:187–228. [23] Schroeder EB, Welch VL, Couper D, Nieto FJ, Liao D, Rosamond WD, et al. Lung function and incident coronary heart disease: the atherosclerosis risk in communities study. American Journal of Epidemiology 2003;158:1171–81. [24] Schuenemann HJ, Dorn J, Grant BJ, Winkelstein W, Trevisan M. Pulmonary function is a long-term predictor of mortality in the general population. Chest 2000;118:656–64.
293
[25] Shojaii R, Alirezaie J, Babyn P. Automatic lung segmentation in CT images using watershed transform. In: IEEE International Conference on Image Processing, 2005. ICIP 2005, vol. 2. 2005. p. 1270–3. [26] Soille P. Morphological image processing: principles and applications. Cambridge University Press; 1999. [27] Sonka M, Hlavac V, Boyle R. Image processing, analysis, and machine vision. Thomson; 2008. [28] Udupa JK, LaBlanc VR, Schmidt H, Imielinska C, Saha PK, Grevera GJ, et al. Methodology for evaluating image-segmentation algorithms, vol. 4684, SPIE, 2002, p. 266–77. http://link.aip.org/link/?PSI/4684/266/1. [29] Vincent L, Soille P. Watersheds in digital spaces: an efficient algorithm based on immersion simulations. IEEE Transactions on Pattern Analysis and Machine Intelligence 1991;13(6):583–98. [30] Völzke H, Alte D, Schmidt CO, Radke D, Lorbeer R, Friedrich N, et al. C. O. S. Cohort profile: the study of health in pomerania. International Journal of Epidemiology; 2010. [31] Zhao B, Gamsu G, Ginsberg MS, Jiang L, Schwartz LH. Automatic detection of small lung nodules on ct utilizing a local density maximum algorithm. Journal of Applied Clinical Medical Physics; 2003.