Fully Automated Segmentation of Lung Parenchyma Using Break and Repair Strategy

1 Introduction

Lung-related cancer is one of the most predominant diseases. In comparison to other cancers like breast cancer, lymphoma, melanoma, and prostate cancer, the survival rate of patients can be improved by early detection and diagnosis [20, 25]. Computed tomography (CT) is widely used to analyze lung nodule size variations and structure. Nowadays, thinner CT slices make more number of frames, which results in more analysis data, and this burdens radiologists. Computer-aided diagnosis/detection (CAD) is more helpful in assisting radiologists in lung parenchyma CT rendering [27].

In the lung CAD system, extracting lung parenchyma from a thoracic CT image is done to reduce analysis time. Reduction in misdiagnosis and improving the efficiency of the CAD system can be done through effective lung parenchyma segmentation. If lung parenchyma segmentation will be a procedural part in the CAD system for the analysis and assessment of lung pulmonary nodule structure, this will influence the precision of the entire lung CAD system.

As of late, researchers have made advances in lung segmentation techniques [2, 7, 19, 23]. In general, these techniques can be divided into three classifications: (i) pattern classification technique [12, 16], (ii) region growing technique [1, 11], and (iii) thresholding technique [10, 13, 14]. The pattern classification technique requires an extensive number of test information preparations and necessities to separate the components; thus, its handling time is longer [21, 22, 26]. It cannot meet the continuous requirements of the CAD framework in clinical applications. The region growing technique is not fully automated. In this technique, high-density regions of the lung attached to the lung parenchyma edges are not included/excluded accurately. It cannot separate the two lung projections. In this technique, multiple seed points need to be chosen physically. It is difficult to choose parameters from texture of area developing and fusing, which may influence the stability. The thresholding technique is straightforward and quick. It cannot preclude the branchial/tracheal zones adequately, and it can exclude nodules that are connected to the inner surface of the lung parenchyma.

In lung CT images, a sort of tumor/nodule clings to the inner wall of the lung, like juxtapleural nodules. It is a vital requirement in the judgment of early-stage lung malignancy/infection because the gray level is fundamentally the same as fat present outside the lung. Owing to juxtapleural nodule gray level variations, frequent errors are made such as identifying the tumor as fat in lung parenchyma segmentation. This leads to misanalysis and impacts the finding of the ailments.

The nodules/tumors on the lung inner wall with the same intensity gray level as fat surrounding the lung turn to the sunken boundary of the lung. Because of the significance of this nodule/tumor in diagnosis, researchers are attempting to develop measures to enhance the detection technique. The region growing active contour and level set the strategy among the most commonly utilized ones [5, 6]. Nevertheless, active counter based on the Chan-Vese segmentation technique cannot precisely judge the extent of the deformity; thus, it is difficult to identify the tumor region while growing.

Considering the above-mentioned issues in lung segmentation, this paper introduces a fully automated procedure for lung parenchyma extraction. The implementation process includes two-dimensional (2D) optimal threshold technique for binarization; then, the fat surrounding the lung and trachea is removed by 2D morphological reconstruction. In the end, the extracted lung parenchyma boundary is modified in enhanced chain code and connectivity can be done using the Bresenham method. The entire procedure of the proposed fully automatic segmentation is detailed in Figure 1. The lung parenchyma segmentation proposed in this paper can successfully extract missing juxtapleural nodules. In comparison with the active contour method, the proposed segmentation incredibly enhances the accuracy and computational speed. Thoracic CT slices generated by different machines and settings are considered to check the robustness and accuracy of the proposed method of lung parenchyma segmentation.

Figure 1:

Flowchart of Automated Segmentation of Lung Parenchyma Using Break and Repair Strategy.

2 Materials and Methods

The proposed method is separated into two phases: (i) lung parenchyma segmentation based on the thoracic CT slice image and (ii) modification of the segmented lung parenchyma boundary using enhanced chain code and Bresenham methods. Each phase involves distinctive subphases. As shown in the Figure 2A, the thoracic CT image incorporates the lung parenchyma, trachea, fat, examination table, and outside area. To reduce the computation and processing time in the first phase, except for the lung parenchyma all areas are subtracted to avoid interference in diagnosis and to accelerate the accurate lung parenchyma extraction. The initial phase of segmentation is to remove regions outside of the lung parenchyma.

Figure 2:

(A) Thoracic CT Image. (B) 2D Histogram of (A).

2.1 Initial Segmentation

The density of the lung parenchyma is low because it is filled with air. The density of other tissues like fat, soft tissues, and bone is higher and different. Therefore, the Hounsfield units (HUs) in the thoracic CT image are proportional to density. The HU estimations of the lung in CT images are practically the same for various cases. In this paper, initially, 2D optimal thresholding is applied to segment the lung parenchyma. The 2D thresholding process is as follows:

1. Generate an L×L size 2D histogram as shown in Figure 2B by pairing input pixel gray level F_i and its local mean F_m , where L is gray level.

2. Divide the histogram by two classes based on object (C_o ) and background (C₁) through

   (1)ω0=∑i=0s∑j=0tqij; ω1=∑i=s+1L−1∑j=t+1L−1qij,

   where q_ij is the probability of the F_i and F_m pair, and ω₀ and ω₁ are probabilities of C₀ and C₁, respectively.

3. Calculate the variance in C_o and C₁ by μ₀ and μ₁, respectively, with total variance μ_T .

4. Use those parameters to obtain variance among class:

   (2)Trace T0(s, t)=ω0[(μ0i−μTi)2+(μ0j−μTj)2]+ω1[(μ1i−μTi)2+(μ1j−μTj)2].

5. Run step (iv) in iterations until T₀(s, t) reaches maximum and this value is taken as optimal threshold:

   (3)Trace T0(s′, t′)=max0≤s,t≤L−1{trace T0(s, t)}.

The above-mentioned 2D optimal thresholding gives better results in noisy images [18, 29]. The segmentation results are shown in Figure 3A (lung lobes and sections are in white).

Figure 3:

(A) Segmentation Result. (B) 2D Reconstruction. (C) Initial Lung Parenchyma Segmentation.

The optimal threshold is used to obtain the binary lung CT image. In the segmented image, tissues outside the lung parenchyma are in black and the lung parenchyma is in white. In the segmented image, the outside of the body and lung will be suppressed by 2D morphological reconstruction that touches the boundary of the image. The 2D morphological reconstruction image is shown in Figure 3B. The lung parenchyma is in white; however, it includes the trachea and its nodular branches. Using an area opening operation, the trachea and tissues present in between the lungs are removed. In this phase, except larger lobes, all objects are removed in the binary image, and resultant image is shown in Figure 3C. This lung-only binary image is used to label lung lobes by 2D connecting components; then, holes are filled to make a proper mask and morphological elements are applied to close uneven edges of the lung. The lung lobes are separated by extracting only the two largest blobs in ascending/descending form. Here, the segmentation and separation of the lung lobes are referred to as the “break” strategy.

2.2 Lung Boundary Modification

Lung parenchyma modification is necessary because some nodules or tumors adhere to the inner lung pleural wall. This is not included properly, as described in Section 2.1. Also, this is the main drawback of pulmonary parenchyma segmentation schemes. Tissues attached to the lung pleural wall are often considered as fat and are wiped out, because the lung border is not even and includes some curved structures. In this paper, an improved chain code algorithm is used to analyze the lung border and the Bresenham method is used to include missed nodules/tumor in the lung border.

2.2.1 Improved Chain Code

In the initial segmented lung parenchyma, the boundary includes concave and convex points. A chain code is used to symbolize convexity and concavity through determined length and directions. This concavity and convexity is analyzed using the Freeman chain code [9, 15]. Assume there are N Freeman chain code boundary points C(x), where (x=1, 2, …, N−1). C(x) is the chain code point and C(x+1) is the next chain code point, and so on. As shown in Figure 4, using this chain code, a relative chain code can be generated. The relative chain code R(x) is defined as the relation between C(x) and C(x−1):

Figure 4:

Eight-Direction Freeman Chain Code.

(4)R(x)=([C(x)−C(x−1)+8] mod 8)>4? R(x)−8:R(x).

The absolute chain code can be generated using relative chain code with forehead absolute chain code values:

A(0)=0,

(5)A(x)=A(x−1)+R(x),

where R(x) is the relative chain code and A(x) is the absolute chain code.

Then, the three-point sum is generated using absolute chain code current point A(x) using two preceding points:

(6)Sum(x)=A(x)+A(x−1)+A(x−2).

Finally, convexity and concavity are obtained by the three-point sum difference:

(7)Diff(x)=sum(x+3)−sum(x).

This three-point sum difference is proportional to the curvature of the boundary line because it is a two-direction sum difference. The concavity difference is positive when the boundary line chain code calculation is in the clockwise course; else, it is identified as convexity.

As indicated in the chain code difference of the boundary, we can reconstruct the boundary line with the relative difference limit [30]. We can recognize convexity and concavity as summarized below:

1. Obtain lung boundary from the binary image.

2. Obtain the Freeman chain code with the origin point in the lung boundary line; the upper left corner point is used as the initial point for subsequent processing.

3. Decide the chain code heading of pivot: clockwise or anticlockwise.

4. Obtain the absolute chain code, relative chain code, and three-point sum Diff(x) of the boundary line.

5. Set a limit as indicated by the attributes of the chain code contrast Diff(x) and then recognize the concavity and convexity points: concave limit Diff(x)>2 and convex limit Diff(x)<−1. These difference limits are empirical values. Figure 5 demonstrates the improved chain code method. The image contains an assortment of convex and concave points computations.

Figure 5:

(A, B) Segregation of Concave and Convex Points.

Figure 5A shows an input test image, and marked concave-convex points are shown in Figure 5B. Symbol “o” represents the concave point in the boundary line and symbol “×” represents the convex point. Four points represent the concave because the chain code difference is >2 and 13 convex points <−1 threshold, as from our literature survey this different threshold value can be judged [24]. According to the above discussions and from the results, the improved chain code is very accurate.

2.2.2 Modification of the Sunken Boundary

With the help of concave-convex points, the sunken boundary can be filled. The following steps are adopted to modify the sunken boundary:

1. Initiate the loop with the chain code starting location and search for the convex point and mark it as P, its location as (x₁, y₁), and get its position from the chain code.

2. Search for concave point Q and its location as (x₂, y₂) and its position from the chain code.

3. Spot successive convex point R and its location as (x₃, y₃) and the respective position from the chain code.

4. In step (iii), on successive search, if the next convex R is also a convex point, then update the R successive convex point location and its position.

5. Finally, connecting convex P and R for filling can be done by connecting the boundary line from the point next to P, that is p+1 to r+1, and boundary line points from p+2 to r+2, and so on, until the boundary line point reaches (p+r)/2, as shown in Figure 6.

Figure 6:

Delineation of Sunken Boundary Filling.

2.2.3 Filling the Sunken Boundary

The sunken boundary is filled in between P(x, y) and R(x, y). The connection of two points in 2D linear space is through the Bresenham method [4]. The Bresenham method is popular in connecting two points because of its accuracy and simplicity. Here, the (x₁, y₁) and (x₂, y₂) coordinates of point P and R, respectively, need to be connected. The procedure of interconnection is as follows.

The slope in between connecting two points can be 0–1 as

dx=x2−x1dy=y2−y1.

Then, the slope will be

(8)k=dydx, 0≤k≤1.

As shown in Figure 7, a point f₁(x_s , y_s ) is the nearest point to P(x, y). The one unit in x direction increase and in y direction about k, and the intersection point is f₂(x_s+1, y_k ). The successive pixels are connected using the relation

Figure 7:

Delineation of Bresenham Method Pixel Interconnection.

(9)xi+1=xi+1yi+1=(d≥0.5)? yi+1:yi,

where d is the distance between the slope line and pixel position. At the pinpoint, d is zero and it will be increased to d=d+k for every X direction pixel position subtracted to one d=d−1 when the Y direction pixel position is forward. The following modifications are applied for accurate filling of sunken boundaries.

1. Exchange dx and dy, x and Y when the absolute value of the slope is >1.0.

2. To control x and y, the sign of dx and dy can be used to increase or decrease by 1.

Here, modifying the sunken boundary of the lung lobes is referred to as “repair” strategy. To check the accuracy of this filling method, we used a previously discussed image and an image similar to the sunken structure of the lung parenchyma.

Figure 8 shows the results of the above steps to fill the sunken boundary of the structure as shown in Figure 5. Figure 9 shows the structure similar to lung lobes. Figure 9B shows the concave and convex points obtained by the improved chain code. In this example, three concave points can be seen in the sunken area. The interconnection of convex points is done to fill the sunken boundary. The results of this type of filling in all tests are accurate and satisfactory.

Figure 8:

Sunken Boundary Filling Result for Figure 5A.

Figure 9:

Structure Similar to Lung Lobes for Filling.

Figure 10 shows the right lung extracted from Figure 3. Figure 10B shows the marking of concave and convex points and the presence of sunken boundaries. These sunken boundaries are filled with the Bresenham algorithm and the resultant mask is shown in Figure 10D. Figure 10E shows the final extracted left lung with complete inclusion.

Figure 10:

Illustration of Extracted Lung Parenchyma Repairing.

3 Performance Evaluation Method

The performance of the proposed fully automated lung parenchyma segmentation method is evaluated and compared with manual segmentation. In this work, two radiologists helped in manually contouring of all slices of the thoracic CT scan image in every case, from the top portion of the lung to the bottom. The manual segmentation results were evaluated by experts and the results are considered as ground truth for original reference. The ground truth images are used to evaluate the proposed fully automated segmentation method. The accuracy of the fully automated segmentation assessed by the following five parameters as well as the accuracy of juxtapleural nodule inclusion are evaluated.

3.1 Volumetric Overlap Error (VOE%)

The rate of volume overlap is estimated with the help of the Jaccard coefficient. The Jaccard coefficient is the ratio of the total value of intersection and union of segmented region A and ground truth region B, as shown in Figure 11A. The VOE is defined as

Figure 11:

Delineation of Segmentation Accuracy.

(10)VOE%=(1−A∩BA∪B)× 100,

where A is the segmented region and B is the ground truth region. The overlap error rate varies from 0% to 100%, i.e. an evaluated value near 0% is perfect segmentation. In perfect segmentation, the segmented and ground truth regions completely overlap [17, 28]. VOE value near 100% represents poor segmentation because of less or no overlap between the two regions, leading to poor segmentation results.

4 Experiments and Results

To evaluate the performance and accuracy of the proposed method, here we consider thoracic CT slices from three different datasets. Each dataset contains CT slices generated by different machines or settings. Table 1 lists the significant parameters of the three different datasets. From the Lung Image Database Consortium and Image Database Resource Initiative (LIDC-IDRI) [3] dataset, there were 82 thoracic CT scans from 82 patients; 29 scans from the Lung TIME1 dataset; and 4 scans from the Lung TIME2 dataset [8]. In total, there were 115 scans and each scan contains 80–140 slices. From the same datasets, 21 thoracic CT scans containing juxtapleural nodules were considered.

The results of the proposed fully automated lung parenchyma segmentation were accurate for all images, which included 31 juxtapleural nodules from 21 thoracic CT scans. Figure 12 illustrates the segmentation and boundary modification results of five slices in column (A) input image, (B) initial segmentation result, (C) markings of concave-convex points on the boundary line, (D) sunken boundary filling using the Bresenham method, (E) segmentation mask, and (F) segmented lung lobe. The proposed method was compared with the most cited active contour segmentation method. In this active counter based on the Chan-Vese segmentation technique, the region of interest is marked by seed points, to segment the lung parenchyma. Figure 13 shows the segmentation results of two slices in comparison to the proposed method. In Figure 13, column (A) shows the input image, (B) segmentation results from active contour technique, and (C) segmentation results from the proposed method. The results showed that the segmentation of juxtapleural nodules by the proposed method was very accurate, compared the active contour technique with poor inclusion of nodules/tumors on the lung pleural wall.

Figure 12:

Experimental Results of Segmentation and Repairing.

Figure 13:

Segmentation Results Comparison between the (A) Input Thoracic CT Images, (B) Active-Contour Method and (C) Proposed Method.

The segmentation accuracy for all 115 scans was evaluated by comparing to ground truth contour using DSC in Eq. (11). The resultant similarity coefficient mean value was about 96±0.35%, compared to the active contour segmentation overlap mean of 82±3.55%, as depicted in Figure 14A. All test results from the proposed method were well in the satisfactory level except small non-overlapping on lobe edges. The juxtapleural nodule inclusion accuracy was evaluated using Eq. (15), and the resultant values are listed in Table 2. The proposed segmentation method was more accurate in inclusion of juxtapleural nodules for all cases. By analyzing all the experimental results, the proposed method is more satisfactory and stable for all cases in lung parenchyma segmentation. The proposed method is 100% accurate in sensing nodules adhering to the lung pleural wall. This resulted in 98% juxtapleural nodule inclusion; however, the complete lung parenchyma segmentation was 96% because of minor over-/undersegmentation on the hilar region. This lower evaluation value is not affected by tracing nodules/tumors adhering on the boundary, so it will not lead to a misdiagnosis.

Figure 14:

Comparison of Performance Evaluation Results.

The computational complexity of the initial segmentation relies on L (gray level). The optimal threshold calculation in a 2D histogram having range [0, L×L] that makes time complexity is about O (L⁴). However, the computation complexity of the proposed method is mainly dependent on the number of points N on the boundary line to generating chain code and interconnection by the Bresenham method. The filling of sunken boundaries is mainly dependent on the k number of filling operations and distance D between convex points; that is, the computational complexity is O(kD²) in the most pessimistic scenario O(N²). For the active contour method computation, complexity is mainly dependent on the number of iterations n to grow segmentation O(n²logn) [6]. The proposed algorithm was verified with the simulations, and MATLAB R2014a was used for programing under Windows 10 OS, Intel i5-3210M processor system. A comparison of execution times is listed in Table 2, and comparatively the execution time is less in the proposed method.

The quality of the proposed segmentation method was determined by three evaluation techniques, and resultant values are shown in Table 3 and Figure 14B. VOE, ASSD, and RMSD measurements ought to have low incentive for good quality segmentation. The values in Table 3 demonstrate that the proposed method of segmenting lung parenchyma is more accurate. The findings related to the visual investigation were obtained by the radiologist.

5 Conclusion

In this paper, the authors propose a new fully automated lung parenchyma segmentation method. The segmentation phase includes 2D optimal threshold selection and 2D morphological reconstruction for lung parenchyma segmentation. Then, lung parenchyma boundaries are modified using improved chain code and Bresenham pixel interconnection. A total of 21 thoracic CT slices having juxtapleural nodules and 115 lung pulmonary parenchyma scans were used to verify the robustness and accuracy of the proposed method. The proposed method was 100% sensitive to inclusion of nodules/tumors present on the lung pleural wall. The juxtapleural nodule segmentation accuracy was >98% and the lung parenchyma segmentation accuracy was >96%. The results of the proposed method were compared with the most cited active contour methods. The proposed method showed good performance in inclusion of juxtapleural nodules, low computational complexity, and no user interference. The proposed method meets the prerequisites of the lung CAD framework, and provides good performance in the processing of a sequence of images.