Automated lung segmentation and smoothing techniques for inclusion of juxtapleural nodules and pulmonary vessels on chest CT images

Biomedical Signal Processing and Control 13 (2014) 62–70

Contents lists available at ScienceDirect

Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc

Automated lung segmentation and smoothing techniques for inclusion of juxtapleural nodules and pulmonary vessels on chest CT images Shengjun Zhou a , Yuanzhi Cheng a,∗ , Shinichi Tamura b a b

School of Computer Science and Technology, Harbin Instituteof Technology, Harbin 150001, PR China Center for Advanced Medical Engineering and Informatics, Osaka University, Suita 565-0871, Japan

a r t i c l e

i n f o

Article history: Received 29 May 2013 Received in revised form 7 February 2014 Accepted 24 March 2014 Keywords: Lung segmentation Juxtapleural nodule Pulmonary vessel Fuzzy c-means clustering Curvature threshold

a b s t r a c t Segmentation of the lung is often performed as an important preprocessing step for quantitative analysis of chest computed tomography (CT) imaging. However, the presence of juxtapleural nodules and pulmonary vessels, image noise or artifacts, and individual anatomical variety make lung segmentation very complex. To address these problems, a fast and fully automatic scheme based on iterative weighted averaging and adaptive curvature threshold is proposed in this study to facilitate accurate lung segmentation for inclusion of juxtapleural nodules and pulmonary vessels and ensure the smoothness of the lung boundary. Our segmentation scheme consists of four main stages: image preprocessing, thorax extraction, lung identification and lung contour correction. The aim of preprocessing stage is to encourage intra-region smoothing and preserve the inter-region edge of CT images. In the thorax extraction stage, the artifacts external to the patient’s body are discarded. Then, a fuzzy-c-means clustering method is used in the thorax region and all lung parenchyma is identified according to fuzzy membership value and connectivity. Finally, the segmented lung contour is smoothed and corrected with iterative weighted averaging and adaptive curvature threshold on each axis slice, respectively. Our method was validated on 20 datasets of chest CT scans containing 65 juxtapleural nodules. Experimental results show that our method can re-include all juxtapleural nodules and achieve an average volume overlap ratio of 95.81 ± 0.89% and an average mean absolute border distance of 0.63 ± 0.09 mm compared with the manually segmented results. The average processing time for segmenting one slice was 2.56 s, which is over 20 times faster than manual segmentation. © 2014 Elsevier Ltd. All rights reserved.

1. Introduction According to the American Cancer Society, lung cancer is the leading cause of cancer death in both women and men, killing more than 150,000 people per year—more than colon, breast, ovarian and prostate cancers combined in the U.S [1]. The existence of lung cancer is usually indicated by the identification of pulmonary nodules. Early detection of potentially cancerous pulmonary nodules may be a way to improve a patient’s chances for survival [2]. Computed tomography (CT) based computer aided diagnoise/detection (CAD) is the most commonly used diagnosis technique because of its high sensitivity of small pulmonary nodules and flexible representation of the human thorax [3]. Accurate segmentation and quantitative

∗ Corresponding author. Tel.: +86 451 86413316; fax: +86 451 86413309. E-mail address: [email protected] (Y. Cheng). http://dx.doi.org/10.1016/j.bspc.2014.03.010 1746-8094/© 2014 Elsevier Ltd. All rights reserved.

assessment of the lung from CT images are of great importance in routine clinical practice for assessing the existence and severity of specific diseases and for monitoring responses to therapies. However, accurate segmentation of the lung CT images with nodules is not an easy task because of the highly varied properties of pulmonary nodules such as shape, size, intensity, and location. In [4], Armato et al. showed that 5–17% of the lung nodules in their test data was missed due to the inaccurate preprocessing segmentation, depending on whether or not the segmentation algorithm was adapted specifically to the nodule detection task. Other problems leading to the inaccurate lung segmentation [5] can be summarized as (a) false positive inclusion of other organs outside the lung, (b) nonsmooth exclusion of the trachea and pulmonary vessels in the hilar regions, and c) posterior and anterior attachments between left and right lungs. Although radiologists are able to identify the lung boundary slice by slice, it is extremely time-consuming and prone to intra- and

S. Zhou et al. / Biomedical Signal Processing and Control 13 (2014) 62–70

inter-observer variability [6]. Therefore, the availability of an automated, reliable, and accurate computerized segmentation tool for this purpose could be of great value. Thresholding method combined with region growing [7,8] is commonly used to extract lung parenchyma, by which those tissues having higher gray level than the selected threshold are excluded from the thoracic region. Nevertheless, since the juxtapleural nodules have similar densities with their surrounding tissues, these juxtapleural nodules, together with their surrounding tissues, will be removed from the lung area and will never be recognized in the later stages. Furthermore, high density pulmonary vessels are also ruled out from the lung area, which bring about indentations and salience in the lung boundary near the mediastinum. To compensate for the lost juxtapleural nodules and ensure the smoothness of the lung boundary, several methods have been proposed to correct the lung contours. In [9], a rolling ball algorithm was applied to the lung contours segmented by thresholding to avoid the loss of the juxtapleural nodules. However, as is pointed out in [10], due to the extremely varied sizes of juxtapleural nodules, it is difficult to select a proper rolling ball radius that is suitable for all juxtapleural nodules. In [11], Gurcan et al. designed an indentation detection method for ameliorating the initial lung contour. Bellotti et al. [12] put forward a new active contour, called glued elastic band, to correct the contour of lung parenchyma. Both two above methods highly rely on the predefined threshold, unfortunately, it cannot be chosen automatically but empirically. Pu et al. [13] proposed an adaptive border marching algorithm (ABM) to correct the lung boundary. The problem with ABM is that the border marching with fixed threshold may failed in some slices. For example, the lung lobe containing juxtapleural nodule needs a smaller threshold while the lung lobe without juxtapleural nodule and containing convex region requires higher threshold. Varshini et al. [14] proposed an improved ABM algorithm with two thresholds to solve this problem. In [15,16], a segmentation-by-registration scheme is used to segment and correct the lung boundary. This type of methods works well on the lungs even containing (substantial) pathologic abnormalities. But the creation of a model used in registration is not a trivial issue due to the high density pathologies in varying degrees of severity. In addition to the abovementioned studies, a few algorithms [17,18] based on curvature information are used to correct the lung border with juxtapleural nodules and pulmonary vessels. The motivation for the curvature-based method comes from the observation that a rapid change in curvature usually indicates a nodule, large vessel or bronchus. But for initial rough contour it is difficult to compute an accurate curvature, since the smoothness of the contour has a great impact on curvature calculation due to the second derivatives involved. Generally, in order to smooth the lung contour, morphological erosion and dilation operations were commonly used techniques [16,19,20]. However, designing a proper structuring element (SE) is still a problem since a large size SE may cause over-segmentation, and vice versa. In this study, we propose an accurate lung segmentation method in thorax CT images and an efficient scheme for smoothing and

63

correcting the segmented lung boundary for including juxtapleural nodules and pulmonary vessels. Our approach begins with image filtering via nonlinear anisotropic diffusion. Then, a fuzzy-c-means clustering method is used to extract initial lung boundary. Finally, a correction scheme based on adaptive curvature threshold is proposed to re-include juxtapleural nodules and pulmonary vessels. To compute an accurate curvature, a scheme termed “iterative weighted averaging” is used to smooth the rough lung contour in advance. Comparing to the works discussed in the literature, the work presented in this paper is different in the following ways: (1) We develop a fully automatic approach for segmenting lung boundary in chest CT images which does not need any user interaction. (2) The proposed method can re-include all juxtapleural nodules and properly include the pulmonary vessels near the mediastinum. (3) Additionally, our method can produce a smooth lung boundary instead of using morphological smoothing. The remainder of this paper is organized as follows: Section 2 provides a detailed description of our method. The experimental results are presented with a brief discussion in Section 3. And the paper concludes in Section 4.

2. Description of the method Our segmentation method consists of four main stages: image preprocessing, thorax extraction, lung identification and lung contour correction, as displayed in Fig. 1. We explain the proposed method in detail.

2.1. Image preprocessing The purpose of preprocessing is to reduce image noise, which benefits the following clustering step, since image clustering algorithms can be sensitive to noise. Linear low-pass filtering gives poor results as it incurs even more edge blurring and detail loss. However, nonlinear anisotropic diffusion filtering can overcome this drawback by introducing an implicit edge detection step into the filtering process so as to encourage intra-region smoothing and preserve the inter-region edge. The basic nonlinear anisotropic diffusion is characterized by the partial differential equation model given by [21]

∂I = div(c(∇ I)∇ I) ∂t

(1)

where  . and div denote the L2 norm and the divergence of the associated quantities, respectively; ∇ represents the gradient of the diffusing image I; and c(.) is a diffusivity function, also referred

Fig. 1. Overview of our lung segmentation method.

64

S. Zhou et al. / Biomedical Signal Processing and Control 13 (2014) 62–70

to as the diffusion coefficients, which is chosen locally as a function of the modulus of the image intensity gradient 2

c(∇ I) = e−(∇ I/ω)

(2)

ω is referred as the diffusion constant and determines the filtering behavior. A commonly used discrete version of (1) is given by the forward Euler numerical approximation with the initial condition of the original image Ij(n1 +1) = Ij(n1 ) + [div(e

−(∇ I (n1 ) /ω)2 j

∇ Ij(n1 ) )]

(3)

where  is the step size of the independent variable t used to approximate ∂I/∂t; n1 is the discrete-time index (or the iteration number); and Ij(n1 ) stands for the image intensity at the position j and the iteration number n1 . For a 2D image I, j ranges from 1 to M × N, where M and N are the size of I. 2.2. Thorax extraction In this step, rough extraction of the thorax based on thresholding is performed. The goal of extracting thorax is to discard the artifacts external to the patient’s body to a certain extent. Accurate extraction of the thorax using thresholding is not an easy task due to the following reasons. First, fat and the bed where patient lies are close to each other and has similar Hounsfield unit [HU] number, which limits the use of thresholding methods. Second, partial volume effects caused by the limited resolution of the imaging system raise more difficulties. Fortunately, precise extraction is not necessarily required in our method since primary purpose of the thorax extraction is to reduce computational cost in the following clustering step. As shown in Fig. 2, fat, bone, muscles and vascular tree in the thorax have high CT numbers. On the contrary, lung parenchyma, airway and background have low CT numbers, so a mean value of two peaks (lung parenchyma, airway and fat) is used to binarize the images. After image binarization, 2D hole filling [22] algorithm is used and the thorax is extracted by searching the largest connected components.

right lungs are separated by analyzing the remaining connected components slice by slice, as shown in Fig. 3. 2.3.1. Fuzzy c-means clustering Fuzzy segmentation has been favored over hard segmentation in some medical image applications since partial volume effects and image noise reduce the accuracy of hard segmentation. To cope with the limits of hard segmentation, a fuzzy c-means (FCM) clustering method is employed to classify pixels into several tissue categories. FCM is a well-known clustering technique used in nonsupervised image segmentation for pixel classification [23]. In FCM method, a set of tissue classes is first determined. Each pixel is then classified by its membership values of tissue classes according to its intensity. Membership value of a certain class indicates the likelihood of the pixel belonging to that class. Each tissue class has a centroid. The objective of FCM clustering is to compute the membership values for each pixel so that they are clustered around the centroid of each class. FCM can be written as the minimization of the weighted inter-class sum of the squared error objective function JFCM . JFCM =

J C  

where vc is the centroid of class c and C is the total number of tissue classes; ujc is the membership value of pixel j for class c and requires

C

ujc ∈ [0, 1] subject to u = 1; J is the total number pixels in the c=1 jc image and Ij is the image intensity at the position j. The objective function is minimized when a large membership is assigned to a pixel close to the centroid of a class and a small membership value is assigned to a pixel far away from the centroid. This is a nonlinear problem and can be solved iteratively. During each iteration, a new set of membership functions and class centroids are computed. The following steps describe the FCM method: (1) Set the initial values for class centroids, vc , c = 1, ..., C. (2) Compute the membership values Ij − vc −2

C

l=1

The lung identification stage is implemented by three steps: airfilled regions including lung parenchyma and airways are firstly detected by fuzzy c-means method; secondly, the large airways (primarily the trachea and main bronchi) are eliminated from airfilled regions by region-growing algorithm; thirdly, the left and

(4)

c=1 j=1

ujc =

2.3. Lung identification

u2jc Ij − vc 2

Ij − vl −2

(5)

(3) Compute the new centroid for each class

J

u2 I j=1 jc j

vc = J

u2 j=1 jc

Fig. 2. Histogram of pixels in CT images.

(6)

S. Zhou et al. / Biomedical Signal Processing and Control 13 (2014) 62–70

65

Fig. 3. The steps of lung identification visualized on a representative slice. (a) Original image. (b) Image of air-filled regions. (c) Image of large airways eliminated. (d) Image of two lungs separated. The range between two vertical line indicates the search area.

(4) Repeat steps (2) and (3) until the algorithm converges. Convergence is reached when the maximum changes over all centroids between two iterations is less than a predefined small threshold value  (0.001 in current algorithm). 2.3.2. Air-filled regions detection The presented technique for detecting the air-filled regions is based on one 2D transverse slice using FCM. We design a scheme to propagate the segmentation to 3D space. First, the segmentation is computed on the top slice of air-filled regions localization. Then, the segmentation is propagated to subsequent slice using the centroid of current segmentation as initialization. The propagation stops at the bottom slice where lung area disappears. We have defined four classes in the chest CT images: air-filled regions (including lung parenchyma and airway), fat, bone, other tissues (including muscles and vascular tree), as depicted in Fig. 2. Four membership values are computed for each pixel. As is known that FCM method can converge faster and more efficiently if they are initialized close to the desired centroid. The initial value of each class centroid is estimated from histogram. Specifically, we calculate the first and second derivatives of the histogram, and the initial estimations are defined as the CT numbers in the histogram with zero-crossings in the first derivative and negative second derivatives. For example, in Fig. 2, four CT values (125 HU, 925 HU, 1059 HU and 1297 HU) are used as the initial centroids for four tissue classes (air-filled regions, fat, other tissues and bone), respectively. Once the fuzzy membership values are obtained, we determine a hard segmentation by assigning each pixel solely to the class that has the highest membership value for that pixel. After performing the hard segmentation, each of the separated regions is a closed, then, we use 2D hole filling algorithm [22] to detect the air-filled regions. 2.3.3. Large airways elimination After detecting the air-filled regions, the trachea and main bronchi still remain in the result. To ensure these structures do not contribute to the segmented lung regions, the trachea and main bronchi are extracted and eliminated from the segmented air-filled region. The location of the trachea is automatically identified by searching the large, circular and air-filled region near the center on the top axial slices of the scan. If no suitable region is detected in the top slices, the bottom slices of the scan are inspected to be able to handle cases that were scanned in a reverse direction. Then,

a seed point is identified based on the center-of-mass location of the segmented trachea region in this slice. From this seed point, a 2D region-growing technique [24] is used to expand the identified trachea region. Regions in the current slice provide potential seed point positions for the next slice. Region growing ceases when the size of the region on a new slice increases dramatically, which indicates the main bronchi has merged into the low-density lung parenchyma. 2.3.4. Left and right lung separation After the large airways elimination, we analyze the remaining connected components slice by slice. This results either in one (both left and right lungs merge together) or two (both two lungs are separated). If the lungs represent one connected component, the lungs are separated using a dynamic programming approach similar to the technique described by [16]. Dynamic programming is applied to find the maximum cost path through a graph with weights proportional to pixel gray-level of the original image. The maximum cost path corresponds to the junction line position. The search area for the dynamic programming is determined in three dimensional, starting from the center of gravity of the lung region. Let xcg be the x-coordinate of the center of gravity. The search area on axial slices is now defined as all pixels (x, y) for which xcg − 10 ≤ x ≤ xcg + 10 instead of xcg − 5 ≤ x ≤ xcg + 5 used in [16], as depicted in Fig. 3d. From our experiment, we found this region to be more accurate to detect the anterior (or posterior) junction line. Pixels outside this range are not considered since it is unlikely that the anterior (or posterior) junction line locates in this area. 2.4. Lung contour correction The juxtapleural nodules tend to be excluded from the lung segmentation result in the previous step due to its density that is similar with thorax; moreover, high density pulmonary vessels are also ruled out from segmentation result which give rise to indentations and salience in the lung boundary near the mediastinum. To solve these problems, a correction algorithm based on adaptive curvature threshold is proposed. Curvature is an important feature for the lung contour where juxtapleural nodules and pulmonary vessels locate since these positions often have larger curvature values. However, it is difficult to compute an accurate curvature for each point on the initial rough contour since the behavior of the curvature is quite sensitive to local

66

S. Zhou et al. / Biomedical Signal Processing and Control 13 (2014) 62–70

perturbations due to the second derivatives involved. In such case, one may wish to smooth the contour in advance. To do so, we first use a fast and stable scheme termed “iterative weighted averaging” to smooth the lung contour, then, a correction algorithm based on adaptive curvature threshold is described. 2.4.1. Iterative weighted averaging Iterative weighted averaging [25] is a commonly used method to smooth the signal by means of averaging the neighborhood. (0) Let Sj = (xj , yj ) be a two-dimensional signal (i.e., lung contour in each slice) at the position j before smoothing. The smoothed signal Sj(n2 +1) on (2s + 1) neighborhoods at the (n2 + 1)th iteration is simply Sj(n2 +1) =

i=s 

(n2 ) Sj+i di+s

(7)

i=−s

with i=s 

di+s = 1

(8)

i=−s

The coefficients di are set symmetrically relative to the center, that is ds+i = ds−i (i = − s, . . . , s .). The expression (8) is equivalent to

2s

d = 1. The choice of the coefficients di (i = 0, . . . , 2s .) is flexible, i=0 i as long as it satisfies the above two conditions. In particular, if s is set as 1, the iterative weighted averaging can be extended to include the nearest-neighbors

(n2 ) (n2 ) Sj(n2 +1) = Sj−1 d0 + Sj(n2 ) d1 + Sj+1 d2

(9)

with d0 + d1 + d2 = 1 Taking d0 = d2 = q and d1 = 1 −2q, where q is a constant, this reduces to (n2 ) (n2 ) Sj(n2 +1) = q(Sj−1 + Sj+1 ) + (1 − 2q)Sj(n2 )

(10)

rearrangement of terms leads to (n2 ) (n2 ) Sj(n2 +1) − Sj(n2 ) = q(Sj−1 − 2Sj(n2 ) + Sj+1 )

(11)

which is a discrete approximation of the linear diffusion equation

∂S = q∇ 2 S ∂t

(12)

In this paper, di is chosen as Newton–Cotes coefficient in the mathematical subfield of numerical analysis since it satisfies the aforementioned two conditions. The expression of Newton–Cotes coefficient di with the equal sample size h can be formulated as (−1)m−i di = mi!(m − i)!



0

m

 m

k=0

k= / i

(t − k)dt

i = 0, ..., 2s.

(13)

Fig. 5. An arc composed of three successive points on the contour. i − 1, i, i + 1 are three successive points on the arc, r is the radius of the arc, l is the Euclidean distance between i and i − 1, w is the height of the arc.

where m = 2s. Noted that di is irrelevant to h if Sj is sampled with the equal step size, which makes the iterative weighted averaging easy to perform. After lung identification in previous step, an edge-tracing algorithm [24] is used to track the lung region automatically for each slice and form a closed initial contour. Starting from the jth edge point on the initial contour at Sj ∈ R2 , all other edge points are stored in counter-clockwise direction. The edge points with the equal sample length h can be derived by linear interpolation. Given Sj(n2 ) in the n2 th iteration for point j, the (n2 + 1)th iteration based on three neighbors can be given below: Sj(n2 +1) =

1 (n2 ) 2 (n2 ) 1 (n2 ) + Sj + Sj+1 S 6 j−1 3 6

(14)

Fig. 4 shows the smoothing results using iterative weighted averaging on a representative lung contour after 5 iterations. Compared with Fig. 4a and b, we can see that the lung contour is well smoothed, which is useful for correction based on curvature threshold in the next step. 2.4.2. Curvature threshold based correction As is described before, the juxtapleural nodules and pulmonary vessels locate at the concave regions and the corresponding bending radii are much smaller compared with those of the whole lung parenchyma. Thus, these concave points with larger curvature values should be corrected, which is associated with three problems such as (1) how to represent the curvature of a point, (2) how to judge the convexity–concavity of a point, and (3) how to determine the curvature threshold for correction. Next, we will solve the above three problems one by one. For an arc composed of three successive points i − 1, i and i + 1 on the contour, let ri be the radius of the arc; i , the curvature of the arc; li , the Euclidean distance between i and i − 1/i + 1; and wi , the height of the arc, i.e., the Euclidean distance from i to the line connecting i − 1 and i + 1 (see Fig. 5). The relationship between ri , li and wi is formulated as below: wi =

li2 2ri

=

li2 2

i

(15)

In the correction algorithm, li corresponds to the Euclidean distance between edge points, so it has a constant value for the equidistant sample, denoted as l. The expression (15) implies that

Fig. 4. (a) Initial contour. (b) Smoothed contour (5 iterations).

S. Zhou et al. / Biomedical Signal Processing and Control 13 (2014) 62–70

67

Table 1 Parameter settings in our algorithm

Fig. 6. Judgment of the convexity–concavity for points on a segment of a closed contour. 1, 2, 3 are three points on the contour. F is the reference vector always pointing to the inside of the contour, W is the height vector of the current point. 1 and 3 are concave points as their reference and height vectors have opposite directions, 2 is a convex point as its reference and height vectors have the same direction.

the height wi of an arc is inverse proportional to the radius ri and proportional to the curvature i . Therefore, wi can be treated as a measure of ri and i . For simplicity, the height wi and radius ri of an arc are also viewed as the height and radius of point i, respectively. In order to judge the convexity–concavity of a point i on a closed contour, we define a reference vector Fi , which originates from point i, plumbs the line connecting i − 1 and i + 1, and always directs to the inside of the contour. If the height wi is extended to a vector Wi that originates from point i, then, the point where the height vector Wi has opposite direction with reference vector Fi is concave, otherwise is convex. For example, in Fig. 6, point 1 and 3 are concave and point 2 is convex. The convexity–concavity of a point i can be indicated by adding a sign to the value of wi and denoted as w˜ i .



w˜ i =

wi

if iisconvex,

−wi

if iisconcave.

Index

Parameters

Value

1 2

 and n1 in Section 2.1 h, n2 and l in Section 2.4

0.15, 3 0.6, 5, 0.8

correction, we need consider the shape feature of a lung contour. Generally, a normal lung contour on each slice is not only smooth but has a approximated crescent-shaped outline, so it can be partitioned into two parts according to convexity–concavity. For a normal contour Q, let Q1 and Q2 be the sequence of all concave and convex pixels, respectively. As Q2 has a smaller bending radius than that of Q1 on average, we have w(Q1 ) < w(Q) < w(Q2 ). Thus, for a normal contour, our correction algorithm will not execute when choosing w(Q) as the curvature correction threshold because Q1 has a larger bending radius than average and Q2 is convex. But for a diseased contour Q containing juxtapleural nodules or pulmonary vessels, the pixels where juxtapleural nodules and pulmonary vessels locate are concave and have much larger curvature than average, when choosing w(Q) as the curvature correction threshold, these pixels will be amended while other pixels with larger curvature but convex or concave but smaller curvature are well protected. 3. Experimental results and discussion In this section, we first discuss the parameter settings used in our method. Then, we compare the result of our method with that of manual segmentation by an expert. All our implementation are programmed in Matlab R2010a environment and operated on a personal computer, equipped with an Intel Core 2 E6500 processor at 2.93GHz and 4GB RAM. 3.1. Image data set

(16)

The curvature threshold is used as the criteria for selecting the points to be corrected. As is noted that juxtapleural nodules and pulmonary vessels commonly locate at the concave parts of the contour and have very small bending radius, so those concave points with very large curvature values should be refined. Let  be the curvature correction threshold; Q be the sequence of all pixels on a contour; and w(Q) be the average of height magnitudes for the sequence Q. Our correction algorithm can be summarized below: (1) For an initial contour Q, apply the iterative weighted averaging method to get a smoothed contour, then, compute the average of height magnitudes w(Q) and set  = w(Q). (2) Resample the contour Q with the equal step length l. (3) For each point i ∈ Q, compute the height magnitude wi , the height vector Wi , the reference vector Fi , and the concavity w˜ i . If w˜ i < 0 and wi > , then it will be amended: connecting its two neighbors i − 1 and i + 1 with a line and removing it from the sequence Q. (4) Repeat steps (2) and (3) until the length of the contour Q does not change. In this correction algorithm, the first step is used to smooth the initial contour and determine the proper curvature threshold for correction. The second step is used to keep the sample length l a constant and assure the correction algorithm stability. The third step is used for judgement and correction. Next, we will explain the selection of the curvature correction threshold in detail. To determine a proper curvature threshold for

For our study, 20 multidetector computed tomography (MDCT) thorax scans of patients with pulmonary nodules were utilized. The images were generated by 2 kinds of MDCT scanners (in 10 cases by Somatom Sensation 64 of the Siemens Medical System; in 10 cases by Brilliance 64 of the Philips Medical Systems). Each patient was imaged by a common protocol (120 kV/Auto mA, helical pitch: 1.35/1) without any contrast enhancement, and the images were created using a normal reconstruction kernel. The image size varied from 512 × 512 × 210 to 512 × 512 × 540 voxels. The in-plane resolution of the images ranged from 0.58 × 0.58 to 0.82 × 0.82 mm and the slice thickness from 0.6 to 1.0 mm. All the manual segmentation results by an expert are available at the time of experiment. 3.2. Parameter settings Several parameters have been used in our algorithm. Table 1 lists these parameters and the values used in current algorithm. The iteration number n1 and the time step  of nonlinear anisotropic diffusion filtering influenced the segmentations of the following FCM method. From numerous experiments, the optimal value of the number of iterations was determined as 3 and the time step was set as 0.15, at which the computational efficiency was maximized without sacrificing the segmentation accuracy. The iteration number n2 and the sample size h of iterative weighted averaging influenced the degree of contour smoothness and thereby influenced the computation of curvature threshold for correction. The constant distance l assured the stability of lung correction process. We performed a number of experiments on the effects of smoothness and correction results with respect to the parameter variations. In this paper, the optimal parameters n2 , h

68

S. Zhou et al. / Biomedical Signal Processing and Control 13 (2014) 62–70

Fig. 7. Examples of lung segmentation. (a) and (b): Examples where there are juxtapleural nodules. (c) and (d): Examples where there are pulmonary vessels. The top row shows the original images, the middle row shows the results of manual segmentation and the bottom row shows the segmentation results of our proposed method.

and l were respectively set as 5, 0.6 and 0.8 based on our numerous experiments with consideration of maximizing the performance of our correction algorithm. We have not found that slight changes of three parameters greatly affect the correction results.

3.3. Experiments on clinical data

MABD(Bauto , Bmanu ) = 0.5 ∗ [dmin (Bauto , Bmanu ) + dmin (Bmanu , Bauto )]

Based on the above parameters, our method was applied to 20 datasets. Among 20 datasets, in total of 307 lung nodules (maximum diameter range: 1.2–17 mm; mean: 3.9 mm) were determined by a consensus panel of two radiologists, of which 65 were juxtapleural nodules. The segmentation results were evaluated by re-inclusion ratio of juxtapleural nodules, accuracy evaluation and processing time. Fig. 7 shows the segmentation results on the representative 2D slices between manually edited ground truth and our proposed segmentation. From Fig. 7, we can draw that our method was successfully applied to re-include the juxtapleural nodules (Fig. 7a and b). By testing against the 20 datasets, our experiments show that all of the juxtapleural nodules were re-included in the segmented lung regions after lung contour correction irrespective of its size and shape. Also, the pulmonary vessels near the mediastinum can be properly included and the entire lung contour can be well smoothed (Fig. 7c and d). As is observed in Fig. 7, there are some differences in the hilar area between manual segmentation and our proposed method. However, as illustrated in [5], it remains unclear at what points the segmentations are supposed to cut through the major bronchi and vessels in the hilar area where they enter the lung. Manual delineations also show much variance around the hilum, depending on definitions (or a lack thereof) and personal preferences with respect to smoothness. Thus, we believe our method based on iterative weighted averaging and adaptive curvature threshold provides a good basis for lung segmentation and smoothing. Fig. 8 shows a 3D surface-rendered view of the segmented lungs before and after contour correction, respectively. The corresponding visualization result is obtained by Marching cubes algorithm [26]. The accuracy evaluation was quantitatively performed using the following four measures [16,17]: volumetric overlap fraction, mean absolute border distance, false positive error and false negative error. • For two binary volumes Vauto and Vmanu , the volumetric overlap fraction is defined as the volume of their intersection divided by the volume of their union, VOF(Vauto , Vmanu ) =

|Vauto |Vauto

 V |  manu Vmanu |

• Given two borders Bauto and Bmanu , let dmin (p, Bmanu ) be the minimum distance of a point p on the border Bauto to the border Bmanu , and dmin (Bauto , Bmanu ) be the average minimum distance from all points on the border Bauto to the border Bmanu . The mean absolute border distance between Bauto and Bmanu is defined as

(18) • False positive (negative) error is the ratio of the total number of automatically (manually) segmented voxels, which are not included in the manual (automatic) segmentation result to the total number of manually segmented voxels. Table 2 shows the comparison of our proposed segmentation to the manually defined ground truth for 20 different datasets using the above metrics. The average of mean absolute border distance was 0.63 ± 0.09 mm. The average of false positive error and false negative error was 1.89 ± 0.89% and 2.39 ± 0.73%, respectively. The volume of our proposed method and manually segmented volume was also shown in Table 2. The average of volumetric overlap fraction was 95.81± 0.89 %. Less than 5% in the volume overlap error is generally acceptable in clinical practice. Therefore, it shows that the proposed algorithm is accurate segmentation scheme of the lung in CT images for the volume measurement in clinical practice. The processing time for each step of our proposed method is summarized in Fig. 9. Except for the lung identification step, other three steps were run in parallel computing. The average total processing time per set of CT scans was 778.28 s and the average processing time to segment each slice was 2.56 s. About 1 min was required for the expert to manually segment two lungs in a slice.

(17) Fig. 8. Lung segmentation results before and after contour correction.

S. Zhou et al. / Biomedical Signal Processing and Control 13 (2014) 62–70

69

Table 2 Assessment of accuracy Subject

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Mean Std. Dev.

VOF

MABD

FP

FN

Vmanu

Vauto

(%)

(mm)

(%)

(%)

(cm3 )

(cm3 )

95.86 94.87 95.09 97.15 95.85 95.96 95.19 95.28 95.47 96.90 97.12 96.58 95.45 94.07 95.71 95.64 97.16 94.87 96.91 95.16 95.81 0.89

0.54 0.61 0.74 0.51 0.63 0.61 0.68 0.63 0.81 0.64 0.51 0.53 0.71 0.74 0.69 0.66 0.56 0.74 0.53 0.62 0.63 0.09

2.56 2.01 1.87 1.64 1.07 1.12 1.86 1.23 2.82 1.16 1.72 1.23 1.67 4.19 1.54 2.12 1.24 3.83 1.07 1.79 1.89 0.89

1.67 3.23 3.14 1.25 3.13 2.97 3.05 3.56 1.81 1.98 1.43 2.24 2.96 1.94 2.83 2.34 1.64 1.46 2.06 3.14 2.39 0.73

3915.3 3546.9 3113.3 3423.1 3857.6 3987.7 2959.2 3491.5 2629.8 2893.4 3788.4 4900.2 3348.7 3461.4 2895.6 3676.4 3491.7 3243.2 3929.2 4785.2

3937.2 3531.4 3099.4 3438.6 3825.9 3966.3 2945.6 3459.8 2657.4 2880.7 3796.2 4881.4 3332.3 3504.6 2872.5 3671.3 3484.8 3275.4 3918.3 4769.9

Fig. 9. The processing time for each step of our proposed method.

Thus, our method was high-efficiency. The mean calculation time in each step of our method was 96.53 s in image preprocessing, 20.72 s in the thorax extraction, 590.38 s in lung identification, and 70.66 s in lung contour correction for 304 slices on average. The processing time of the first three steps was proportional to the number of images for each dataset. The processing time exhausted on the lung contour correction step depended on the number of juxtapleural nodules and pulmonary vessels to be corrected per scan. In our method, adaptive curvature threshold is used for correcting the segmented lung contour to re-include juxtapleural nodules and pulmonary vessels. Compared with the works presented in [9,11,12], our threshold is automatically determined and varies on different slices. This help to avoid the erroneous inclusion caused by a single threshold. Moreover, since designing a proper SE is difficult, we introduce “iterative weighted averaging” scheme instead

of using morphological smoothing presented in [17,20], to get a smooth lung contour. By this way, the curvature threshold can be calculated more accurate. We recognize the proposed technique still needs further improvement. As the examples in Fig. 10 demonstrate, lung contour smoothing based on iterative weighted averaging will miss local sharp features of lung regions. If a small juxtapleural nodule is coincidently located at such “sharp” regions, it may be missed in the segmented lung volume. This limitation contributes to the majority of the undersegmentation in the evaluation of lung volume segmentation (Table 2). It is difficult to overcome this limitation or dilemma because the smoothness property of iterative weighted averaging will definitely lead to the missing of very details with high curvatures. Furthermore, it is necessary to smooth the initial rough contour for computing an accurate curvature since the behavior of the curvature is quite sensitive to local perturbations

Fig. 10. The local enlargement of the regions. (a)–(d) Indicated by the boxes in Fig. 7a–d.

70

S. Zhou et al. / Biomedical Signal Processing and Control 13 (2014) 62–70

Fig. 11. Example of the region merging result. (a) Result of lung identification. (b) Result of lung contour correction. (c) Result of (a) and (b) merged.

due to the second derivatives involved. However, as the examples in Fig. 11 show, if the final segmentation results are given by the region merging between lung identification and lung contour correction, the undersegmentation errors caused by this dilemma will be alleviate largely. 4. Conclusion In this paper, we have developed an accurate and fully automatic method for segmenting lung boundary in chest CT images and an efficient scheme for smoothing and correcting the segmented lung boundary for inclusion of juxtapleural nodules and pulmonary vessels. The effectiveness of our approach was demonstrated on 20 diseased lung scans and the results were compared with the manual segmentations by an expert. Experimental results show that our method can re-include all juxtapleural nodules and properly include the pulmonary vessels near the mediastinum. The proposed method can achieve an average mean absolute border distance of 0.63 ± 0.09 mm and an average volume overlap ratio of 95.81 ± 0.89%. Less than 5% in the volume overlap error is generally acceptable in clinical application. So our method is an accurate lung segmentation scheme from CT images for the volume measurement in clinical practice. Accurate segmentation of lung in presence of severe pathologies such as lung cancer is still a challenging task. Our future work will mainly focus on the segmentation algorithm of lung with severe abnormalities. Acknowledgment This work is supported by the National Natural Science Foundation of China under Grant No. 61073124; and the Fundamental Research Funds for the Central Universities under Grant No. HIT. IBRSEM. 201328. References [1] Cancer Facts & Figures, American Cancer Society: Cancer Statistics, 2012 http://www.cancer.org [2] C.I. Henschke, D.I. McCauley, D.F. Yankelevitz, D.P. Naidich, G. McGuinness, O.S. Miettinen, D.M. Libby, M.W. Pasmantier, J. Koizumi, N.K. Altorki, J.P. Smith, Early lung cancer action project: overall design and findings from baseline screening, Lancet 354 (9173) (1999) 99–105. [3] E.A. Hoffman, G. McLennan, Assessment of the pulmonary structure–function relationship and clinical outcomes measures: quantitative volumetric CT of the lung, Acad. Radiol. 4 (11) (1997) 758–776. [4] S.G. Armato, W.F. Sensakovic, Automated lung segmentation for thoracic CT: impact on computer-aided diagnosis, Acad. Radiol. 11 (9) (2004) 1011–1021. [5] I. Sluimer, A. Schilham, M. Prokop, B. van Ginneken, Computer analysis of computed tomography scans of the pulmonary: a survey, IEEE Trans. Med. Imaging 25 (4) (2006) 385–405.

[6] Q. Li, F. Li, K. Suzuki, J. Shiraishi, H. Abe, R. Engelmann, Y. Nie, H. MacMahon, K. Doi, Computer-aided diagnosis in thoracic CT, Seminars in Ultrasound, CT, MRI 26 (5) (2005) 357–363. [7] J.P. Ko, M. Betke, Chest CT: automated nodule detection and assessment of change over time – preliminary experience, Radiology 218 (1) (2001) 267–273. [8] Y. Lee, T. Hara, H. Fujita, S. Itoh, T. Ishigaki, Automated detection of pulmonary nodules in helical CT images based on an improved template – matching technique, IEEE Trans. Med. Imaging 20 (7) (2001) 595–604. [9] S.G. Armato, M.L. Giger, C.J. Moran, J.T. Blackburn, K. Doi, H. MacMahon, Computerized detection of pulmonary nodules on CT scans, Radiographics 19 (5) (1999) 1303–1311. [10] D.Y. Kim, J.H. Kim, S.M. Noh, J.W. Park, Pulmonary nodule detection using chest CT images, Acta Radiol. 44 (3) (2003) 252–257. [11] M.N. Gurcan, B. Sahiner, N. Petrick, H.P. Chan, E.A. Kazerooni, P.N. Cascade, L. Hadjiiski, Lung nodule detection on thoracic computed tomography images: preliminary evaluation of a computer-aided diagnosis system, Med. Phys. 29 (11) (2002) 2552–2558. [12] R. Bellotti, F. De Carlo, G. Gargano, S. Tangaro, D. Cascio, E. Catanzariti, P. Cerello, S.C. Cheran, P. Delogu, I. De Mitri, C. Fulcheri, D. Grosso, A. Retico, S. Squarcia, E. Tommasi, B. Golosio, A CAD system for nodule detection in low-dose lung CTs based on region growing and a new active contour model, Med. Phys. 34 (12) (2007) 4901–4910. [13] J. Pu, J. Roos, C.A. Yi, S. Napel, G.D. Rubin, D.S. Paik, Adaptive border marching algorithm: automatic lung segmentation on chest CT images, Comput. Med. Imaging Graph. 32 (6) (2008) 452–462. [14] P.R. Varshini, S. Baskar, S. Alagappan, An improved adaptive border marching algorithm for inclusion of juxtapleural nodule in lung segmentation of CT-images, IEEE Int. Conf. Image Process. (2012) 230–235. [15] I. Sluimer, M. Prokop, B. van Ginneken, Toward automated segmentation of the pathological lung in CT, IEEE Trans. Med. Imaging 24 (8) (2005) 1025– 1038. [16] E.M. van Rikxoort, B. de Hoop, M.A. Viergever, M. Prokop, B. van Ginneken, Automatic lung segmentation from thoracic computed tomography scans using a hybrid approach with error detection, Med. Phys. 36 (7) (2009) 2934– 2947. [17] Y. Yim, H. Hong, Correction of segmented lung boundary for inclusion of pleural nodules and pulmonary vessels in chest CT images, Comput. Biol. Med. 38 (8) (2008) 845–857. [18] Y. Yim, H. Hong, J. Beom Seo, N. Kim, E. Jin Chae, Y. Gil Shin, Correction of lung boundary using the gradient and intensity distribution, Comput. Biol. Med. 39 (3) (2009) 239–250. [19] S. Hu, E.A. Hoffman, J.M. Reinhardt, Automatic lung segmentation for accurate quantitation of volumetric X-ray CT images, IEEE Trans. Med. Imaging 20 (6) (2001) 490–498. [20] S. Ukil, J.M. Reinhardt, Smoothing lung segmentation surfaces in 3D X-ray CT images using anatomic guidance, in: Proceedings of SPIE Conference on Medical Imaging: Image Processing, 2004, pp. 1066–1075. [21] P. Perona, J. Malik, Scale-space and edge detection using anisotropic diffusion, IEEE Trans. Pattern Anal. Mach. Intell. 12 (7) (1990) 629–639. [22] R. Gonzalez, R. Woods, Digital Image Processing, 2nd ed., Prentice-Hall, Upper Saddle River, NJ, 2002. [23] S. Sivakumar, C. Chandrasekar, Lungs image segmentation through weighted FCM, in: IEEE International Conference on Recent Advances in Computing and Software Systems, 2012, pp. 109–113. [24] M. Sonka, V. Hlavac, R. Boyle, Image Processing, Analysis, and Machine Vision, Brooks/Cole Publishing, Pacific Grove, CA, 1999. [25] P. Saint-Marc, Adaptive smoothing: a general tool for early vision, IEEE Trans. Pattern Anal. Mach. Intell. 13 (6) (1991) 514–529. [26] K.S. Delibasis, G.K. Matsopoulos, N.A. Mouravliansky, K.S. Nikita, A novel and efficient implementation of the marching cubes algorithm, Comput. Med. Imaging Graph. 25 (4) (2001) 343–352.