Modified and Optimized Method for Segmenting Pulmonary Parenchyma in CT Lung Images, Based on Fractional Calculus and Natural Selection

3.1 Two-Dimensional Histogram

The 2D Otsu method is based on the input image gray level intensity function f(x, y) and local mean level g(x, y). The input image size m×n and the number of pixels N (N=m×n) and f(x, y) gray level vary from 0 to L−1 levels. If the f(x, y), which is the gray level, is i and its neighborhood mean value is j, then the pair is (i, j). The total number of (i, j) pairs in existence is O_ij and its probability [11] is

where i, j=0, 1, …, L−1.

The top view of the 2D histogram q_ij is shown in Figure 1, in which the x coordinates represent f(x, y) gray levels and the y coordinates represent g(x,y) local mean levels. This histogram is divided into four sections by vector (u, v), where 0<u, v≤L−1, because of homogeneity. The background or the object is distributed along the diagonal quadrants 1 and 3 because the background and object gray levels of a pixel and their local mean are similar. There is no similarity in the pixel gray levels of edge, which is in between the objects and the background. Therefore, the pixels near the edges and noise are distributed along the off-diagonal quadrants 2 and 4.

Figure 1:

Two-Dimensional Histogram.

(A) Projection in the plane x0y. (B) 2D histogram of the CT image 3 (CT-3).

3.2 Two-Dimensional Otsu

The input image pixel is divided into two classes: object (C₀) and background (C₁) by threshold vector (u, v). ω₀ and ω₁ are probabilities of C₀ and C₁, respectively:

and the respective mean levels are

and the total mean from the 2D histogram is

The 2D histograms of off-diagonal elements are neglected.

The variance among classes of discrete matrix is defined as

The trace S_B measured in between class [22] is

The maximum value of trS_B taken as optimal threshold vector (u′, v′) is

This iterative 2D Otsu method has more complexity, of about O (L⁴), as discussed in Refs. [8, 27, 25]. Therefore, this paper introduces DPSO and FODPSO to resolve this complexity. Figure 2 shows the 2D Otsu algorithm results. One corner of the segmented lung region circular structured dilation is applied as a smoothing filter to avoid losing nodules that are adhering to the walls of the lung.

Figure 2:

Segmentation of lung parenchyma in thoracic CT image by 2D Otsu.

(A) CT slice image 1 (CT-1). (B) 2D histogram. (C) Segmentation done by 2D Otsu. (D) Dilation image. (E) Segmented lung parenchyma. (F) Enhanced image.

3.3 PSO

The discussion on PSO is needed to understand the DPSO and FODPSO algorithms. The PSO was initially proposed by Eberhart and Kennedy [5]. This algorithm is based on the swarm intelligence concept. Imagine that the number of crying birds is proportional to the amount of food present in the current location. This makes the birds communicate with each other and leads the swarm to move to the location having the highest concentration of food. This swarm movement depends mainly on three rules: (i) flying in the same course, (ii) reaching the location where there is the highest concentration of food, and (iii) moving toward the bird that cries the loudest [6].

In the PSO method, the candidate solutions are called particles. Here, 50 particles are initialized as a random swarm, each having (u and v) unknowns. With the effect of its best solution and its neighbor, the best solution particle adjusts its velocity (vtn) and momentum (xtn). Using the Otsu function [Eq. (10)], the particle fitness is evaluated. The algorithm executes by storing and replacing the local best (x⌣tn) of the particles, the global best (g⌣tn), and the neighborhood best (n⌣tn). The particle velocity and position are updated through Eqs. (11) and (12):

The constants w, ρ₁, ρ₂, and ρ₃ assign weights to the inertial impact; the g⌣tn,x⌣tn, and n⌣tn decide the new velocity individually. The particles use the 2D histogram as a search space for g⌣tn. The g⌣tn value is used for segmenting the CT image, as shown in Table 1.

3.4 DPSO

The main drawback in PSO is that the algorithm may get trapped in a local optimum solution, such that it works only for some problems. The DPSO was proposed by Tillett et al. [23] as a better model of using PSO. The DPSO is the extension of PSO with the basic principle of “survival of the fittest”. This helps PSO to escape from the local optimum. In DPSO, multiple PSO algorithms run as multiple swarms. The results of each swarm are compared and the best one extends its life and the stagnated one is deleted.

The DPSO rules are as follows: the stagnated particles will be deleted, which leads to particle population decrease. If there are insufficient particles in swarms, then the swarm will be deleted and a new swarm will be initialized and follow the main loop steps as shown in Table 2. The 2D Otsu optimized by DPSO is executed on lung CT images, and the results are superior when compared to 2D Otsu with PSO.

3.5 FODPSO

The FODPSO was introduced by Couceiro et al. [3] as an extension model for DPSO, in which the convergence rate of the algorithm is controlled by fractional calculus.

The fractional calculus method attracts several researchers and is applied in several scientific fields. The Grunwald–Letnikov definition of fractional differential with coefficient α∈∁ of the signal x(t) is

where r is the truncation order, T is the sampling time, and γ is the gamma function. With influence of Eq. (11) as w=1, assuming T=1, similar to Pires et al. [15]:

where r≥4, as complexity increases with r. Then, Eqs. (11), (13), and (14) are rewritten as

where α=1 (without memory) is a case of DPSO. Hence, particles disregard their previous activity and get trapped in a local solution (i.e. exploitation). The larger the α value, the nearer to 1, which helps the particles to diversify and get new solutions (i.e. exploration); however, if the exploration is more, then the algorithm consumes more time for the global solution. As Ghamisi et al. [6] showed, experimentation on fractional coefficient leads to α=0.6; thus, balancing between exploration and exploitation is considered in the FODPSO algorithm. The detailed steps of lung parenchyma segmentation in a thoracic CT image by using fractional calculus are shown in Table 3.

4.1 Regional Non-Uniformity (NU)

The segmentation is assessed using the NU technique. The NU is defined as

where F_T signifies the foreground threshold and B_T signifies the background threshold of the image. σ² corresponds to the variance of the complete image and σf2 corresponds to the foreground variance. The estimation of regional NU measures the degree of segmentation carried out. If the NU result is around 0, then the segmentation is righteous, and if near 1 it represents an unrighteous segmentation [8].

4.2 Correlation

The correlation of two images signifies how firmly connected the images are. The quality of segmentation can be proved by correlation. In this study, the correlation is done between the ground truth contour image and the segmentation resultant image [18].

4.3 Peak Signal-to-Noise Ratio (PSNR) and Structural Similarity Index Matrix (SSIM)

Along with the above evaluation functions, the CPU execution time and image quality evaluation by PSNR and SSIM are also considered [16].

where f and g are the original and segmented images, and σf2 and σg2 are variances. σ_fg is a covariance, μ_f and μ_g are average values, and C₁=(K₁L)² and C₂=(K₂L)² with L=256, K₁=0.01, and K₂=0.03.

5.1 Experimental Setup

The proposed algorithm is verified with the simulations, and MATLAB R2015a is used in a Windows 10 OS, i5 processor system. The performance of the DPSO and FODPSO algorithms was compared with the traditional 2D Otsu. The initial parameters of the algorithms are shown in Table 4. The segmentation accuracy of the proposed method is compared with the widely used region growing algorithm based on seed points and the k-means clustering segmentation methods.

In experimentation, as many as 40 cases are taken, among which four cases are reported here. Lung CT images of size 512×512, from Lung Image Database Consortium datasets [1] are generated by GE Medical Systems models (Light-Speed 16 and VCT models). In this experimentation, with the help of two radiologists, manual contouring is done on all CT slice images and these contour images are considered as ground truth for reference. Lung region segmentation is done using the optimized Otsu algorithm by DPSO and FODPSO, and then the segmented binary mask is used for pulmonary parenchyma extraction. The results are depicted in Figure 3.

Figure 3:

Segmentation of lung parenchyma in thoracic CT image using optimization techniques.

(A) CT images (CT-1,2,3,4). (B) Lung segmentation done by 2D Otsu optimized by PSO. (C) Lung segmentation done by 2D Otsu optimized by DPSO. (D) Lung segmentation done by 2D Otsu optimized by FODPSO. (E) Enhanced resultant image of FODPSO.

5.2 Segmented Output

The proposed methods of lung parenchyma segmentation are as shown in Figure 3 by PSO, DPSO, and FODPSO. Figure 4 illustrates the segmentation results of region growing, k-means clustering, and the proposed method. As shown in Figure 4D, the segmentation result includes nodules adhering to the inner wall of the lung and accurate lung parenchyma segmentation.

Figure 4:

Comparison of segmentation results.

(A) CT images (CT-5,6). (B) Lung segmentation done by region growing. (C) Lung segmentation done by k-means clustering method. (D) Lung segmentation done by the proposed method (optimized by FODP).

5.3 Computation Time

The average computing time from 20 runs in the conventional and proposed methods are shown in Table 5.

5.4 Evaluation Results and Discussion

As Table 5 depicts, when there is large gray level in the image, 2D Otsu has poor practicability and is time consuming, because complexity depends on L (gray levels). The optimal threshold calculation in 1D histogram will be in the range [0, L], while a 2D histogram has the range [0, L×L], which makes time complexity about O (L⁴) for processing a complete 2D histogram plane. As discussed here, all optimization techniques are random population and stochastic based. The PSO optimization technique results in higher computing time compared to DPSO and FODPSO. Here, the PSO population is fixed, such as 50 particles mean 50 different solutions for evaluation within the same swarm. As DPSO and FODPSO are composed of multiple swarms, these techniques take less time than PSO. In the tabulated values, the average computing time for threshold calculation to segment lung parenchyma by the traditional 2D Otsu method has been decreased to a greater extent by using the proposed DPSO and FODPSO. When considering more number of CT slices for CAD, the proposed methods will diminish an enormous amount of image analysis time.

As discussed above in performance assessment by NU, correlation, PSNR, and SSIM were applied to the all CT slice images. The results are given in Tables 6 and 7 . Table 7 lists the performance evaluation results of the algorithms. The region growing and k-means clustering algorithms have limitations and result in inaccurate lung segmentation as shown in Figure 4B and C. In Table 6 values, the region NU is nearer to 1 because the traditional and optimized segmentation results are similar. These tabular values verify the efficiency of the proposed DPSO and FODPSO algorithms.

As discussed above in performance assessment by region NU, the results are given in Table 6 and the value is nearer to 1 because the traditional and optimized segmentation results are similar. These tabulated values confirm the efficiency of the proposed DPSO and FODPSO algorithms. The correlation, PSNR, and SSIM were applied to the all CT slice images. Table 7 lists the performance evaluation results of the algorithms for comparison. As compared to the proposed method for segmenting lung parenchyma by 2D Otsu optimized by FODPSO, region growing and k-means clustering algorithms have limitations and result in inaccurate lung segmentation, as shown in Figures 4 and 5.

Figure 5:

Comparison of Performance Evaluation Results Between Region Growing, k-Means Clustering, and Our Algorithm.

5.5 Inference from Evaluation Results

From the above resultant figures and tabulated values, we can infer that traditional 2D Otsu and the proposed method have the same results. However, compared to traditional methods, the proposed method removes an immense amount of computation time.