1 Introduction

Pulmonary disease is one of the major causes of morbidity and mortality around the world [1, 2]. For example, the recent global outbreak of COVID-19 has killed tens of thousands of people in just a few months. Early diagnosis of pulmonary disease with computed tomography (CT) technique is crucial for making treatment decisions. In non-invasive detection and diagnosis of pulmonary disease, accurate lung segmentation is often a prerequisite for assessing the disease severity, it ensures that disease detection is not confounded by regions outside lungs [2]. However, inner structures of thoracic CT images are usually various with different textures and pixel densities. Additionally, intensities of pathological images are inhomogeneous and it is difficult to provide a reliable generic solution for a wide spectrum of lung abnormalities, which cause difficulties in lung segmentation.

Fig. 1
figure 1

Pathological lung segmentation. A pathological thoracic CT image (a) with its segmented lung mask (b)

Lung segmentation refers to the computer-based process of identifying the boundaries of lungs from surrounding thoracic tissue on CT images [3], see Fig. 1 for a reference. In this paper, the key issue we want to tackle is to segment pathological lungs from thoracic CT images accurately overcoming external distractions of lung diseases and abnormalities. Thus, we proposed a novel pathological segmentation algorithm based on random forest, deep convolutional network, and multi-scale superpixels. Our contributions in the paper are:

  • We propose a novel pathological lung segmentation algorithm combining three principal processes: feature extraction, classification fusion, and contour correction, which can generate more complete segmentations.

  • We put forward an effective classification fusion method based on a fractional-order gray correlation, which can produce more accurate fusion results in multi-scale classifications.

  • We present a new lung segmentation refinement approach based on a divide-and-conquer strategy of contour correction of left lungs and region repair of right lungs, which contributes to generating more accurate lung segmentations.

The remainder of this paper is organized as follows: in Sect. 2, we introduce the related research in lung segmentation, and in Sect. 3, we make a detailed description of our lung segmentation algorithm. In Sect. 4, we provide a set of experimental results and in Sect. 5 and Sect. 6, we make a discussion and summarize our algorithm. The basic pipeline of our algorithm is illustrated in Fig. 2.

Fig. 2
figure 2

The pipeline of our algorithm. A DCNN model is trained with samples captured from CT images and their ground truths. An input image is segmented into three groups of superpixels: \(P_1\), \(P_2\), and \(P_3\), with three different scales, respectively. The deep features, texture, and intensity features are extracted from each group of superpixels, which are classified with random forests (RFs). We fused the classification results (\(R_1\), \(R_2\), and \(R_3\)) and refined the fused segmentation result with a divide-and-conquer strategy and information propagation mechanism

2 Related Work

Many methods have been proposed to segment lungs from thoracic CT images in the past decades. In [4], three thresholding based approaches, connected threshold, neighborhood connected, and threshold level set were performed on lung segmentation. Prabin et al. [5] segmented thoracic CT images using a region growing algorithm along with a combination of supervised contextual clustering technique. Mansoor et al. [6] segmented pathological lungs from CT scans by combining region-based segmentation with a local descriptor based classification. Chen et al. [7] segmented pathological lungs from 3D low-dose CT images using an eigenspace sparse shape composition by integrating a sparse shape composition with an eigenvector space shape prior model. Revathi et al. [8] introduced a pathological lung identification system, where, FC was used for segmenting lungs with a diverse range of lung abnormalities, RF was then applied to refine the segmentation by identifying pathology and non-pathology tissues according to the features extracted from the gray-level co-occurrence matrix (GLCM), gray level run length matrices, histograms and so on. Soliman et al. [9] segmented pathological lungs from CT images based on a 3D joint Markov-Gibbs random field model which integrated the first-order visual appearance model, the second-order spatial interaction model, and a shape prior model. Hua et al. [10] used a graph search driven by a cost function combining the intensity, gradient, boundary smoothness, and rib information for pathological thoracic CT image segmentation. Hosseini-Asl et al. [11] proposed a nonnegative matrix factorization (NMF)-based pathological lung segmentation approach, which included three stages: removing image background from CT images by a conventional 3D region growing method, modeling visual appearance of the remaining chest-lung image with an NMF technique, extracting 3D lung voxels by a two-step data clustering and a region map cleaning approach. An automated lung segmentation algorithm was developed combining unsupervised and supervised techniques in [12]. The method combined an unsupervised MRF technique to provide an initial estimate of lung borders, a supervised texture analysis scheme based on an SVM classifier was applied on searching border regions and distinguishing lung tissues from their surrounding tissues. Liu et al. [13] segmented lungs from CT images with a random forest (RF) classifier, where texture and intensity features extracted from superpixels were taken as the input of RF Classifier. Meng et al. [14] presented a lung segmentation algorithm based on anatomical knowledge and a snake model. By setting the snake model’s initial curve on human’s rib edge in the thoracic CT image, their model could capture the concavities locating on lung boundaries well. Abdollahi et al. [15] segmented multi-scale initial lungs from CT images using a linear combination of discrete Gaussian approach and a Markov-Gibbs random field (GMRF) model. The initial segmentations were fused together using a Bayesian approach to obtain the final segmentation of lung regions. The potential benefits of deep learning techniques in image analysis have also generated remarkable results [16, 17]. Harrison et al. [1] presented a bottom-up deep-learning-based approach unaffected by any variations in lung shape. In the method, a deep model, progressive holistically-nested network, was used to produce finer detailed lung masks. Park et al. [18] employed a two-dimensional U-Net for lung parenchyma segmentation. Anthimopoulos et al. [19] used a deep purely convolutional neural network for the semantic segmentation of interstitial lung disease (ILD) patterns. Lung CT images of arbitrary sizes were taken as inputs and produced corresponding label maps.

In lung segmentation, manually crafted features such as shape, color, and/or texture play an important part in some algorithms. Nevertheless, the features are complementary in medical images and cannot make representations of high-level problem domain concepts [20]. Deep learning is often used to improve the descriptive ability in feature representation [21,22,23]. For example, Hong et al. [21] used deep neural networks to compute features for face pose estimation. Zhang et al. [22] extracted deep features with contractive autoencoders for unsupervised dimension reduction. Convolutional neural network (CNN) has an ability of representation learning, where input information is extracted and filtered layer by layer. The convolutional layers of a CNN can be seen as feature extractors, which generate local features of image patches in each layer, and the features are combined to produce a global deep feature vector in the last. Recent work demonstrated that 2D CT slices are expressive enough for segmenting complex organs [24]. Consequently, we fuse the deep features and low-level traditional features to character different regions of thoracic CT images and extract lung regions from CT images with RF according to the fused features. Furthermore, the lung segmentation is refined by contour correcting and region repairing.

3 Methods

The superpixel algorithm [25] can group pixels into perceptually meaningful atomic regions with similar features, such as intensity, texture, and so on. Approximately equally-sized superpixels with boundaries aligning to local image edges are more suitable for being taken as classification units than isolated image patches. Consequently, we take superpixels as classification units of deep convolutional neural network (DCNN) and RF classifiers for preserving lung contours well. Our lung segmentation algorithm mainly includes four steps: superpixel segmenting, feature extracting, RF classification, and segmentation refining.

3.1 Superpixel Segmenting

After preprocessed with morphological filters [26], thoracic CT images are segmented into a group of superpixels with the simple linear iterative clustering (SLIC) approach [25]. SLIC employs k-means clustering approach to efficiently generate superpixels with approximately equal sizes. The scales of superpixels in an image depend on a user-specified parameter, P. Assume the size of a thoracic CT image is Z, the scale of a superpixel is roughly \(\sqrt{Z/P}\). Hence, smaller P produces larger-scaled superpixels, which is helpful for segmenting inhomogeneous intensities images involving pathologic tissues with strong anti-noise interference ability. Conversely, bigger P generates smaller-scaled superpixels, which can capture detailed regions. Here, a thoracic CT image is segmented by SLIC with three Ps (\(P_1\), \(P_2\), \(P_3\)), respectively. As a consequence, we obtain three groups of superpixels with three different scales.

3.2 Feature Extraction

We map the superpixels into the thoracic CT image and extract deep, texture, and intensity features from superpixels.

3.2.1 Deep Feature Extraction

CNN is one of deep neural networks with multilayer perceptrons and usually consists of an input and an output layers, as well as multiple hidden layers typically including a series of convolutional layers, pooling layers, fully connected layers, and normalization layers [27]. Here, we construct an 8-layer DCNN model for deep feature extraction:

Fig. 3
figure 3

Our DCNN consists of three convolutional layers, three pooling layers, one fully connected layer and one output layer

The outputs of the convolutional layers are all activated by a Sigmoid function,

Three kinds of regions, lungs, pleural tissues, and backgrounds, are captured from thoracic CT images according to their ground truths. They are resized and used for training and testing of our DCNN model. Finally, the output of the last mean-pooling layer, a feature vector of 40 elements, is extracted and taken as our deep features.

3.2.2 Texture Feature Extraction

GLCM is a widely used texture statistical analysis and measurement tool. Its statistical characteristics are superior to those of fractal dimension, Markov model, Gabor filter and so on [28, 29]. We extract 46 features from GLCM, such as energy, contrast, homogeneity, correlation and so on. Some of them are defined as \(Energy=\sum _{I,j=0}^{N-1}{P_{i,j}^2}\), \(Contrast=\sum _{I,j=0}^{N-1}{P_{i,j}(i-j)^2}\), \(Homogeneity=\sum _{I,j=0}^{N-1} \frac{ P_{i,j}}{1+(i-j)^2}\), \(Correlation=\sum _{I,j=0}^{N-1} \frac{ P_{i,j}(i-\mu )}{\sigma ^2}\), where \(\mu \) and \(\sigma \) are the mean and standard deviations of GLCM.

Besides GLCM, moment invariants, \(\eta \)s, have been extensively used to characterize the patterns of images in a variety of applications [30].

$$\begin{aligned} \eta _{pq}=\int _{-\infty }^{\infty } \int _{-\infty }^{\infty } x^p f(x,y) dxdy, \quad p,q=0,1,2,\ldots \end{aligned}$$
(1)

Hu [31] introduced seven moment invariants \(\phi _1\sim \phi _7\).

Moment invariants are useful properties, which are robust to image scaling, translation, and rotation. Collectively, we obtain a total of 54 texture features after adding image entropy [32] into texture features.

3.2.3 Intensity Features

Since color information is mainly distributed in low-order moments, color moments have been proved an effective tool in representing the color distribution in images [33]. Stricker et al. [34] introduced three color moments to represent image color distribution including first-order moments (\(\mu \)), second-order moments (\(\delta \)), and third-order moments (s) defined as \(\mu _i=\frac{1}{N}\sum _{j=1}^{N}p_{i,j}\), \(\delta _i=\sqrt{\frac{1}{N}\sum _{j=1}^{N}(p_{i,j}-\mu _i)^2}\), \(s_i=\root 3 \of {\frac{1}{N}\sum _{j=1}^{N}(p_{i,j}-\mu _i)^3}\), where \(p_{i,j}\) stands for the ith color value of pixel i, N and \(\mu \) respectively represent the total number of pixels in the image and the mean value of ith color channel.

Here, we convert a CT image into a grayscale image and extract intensity features from it, and in total, we concatenate 97 features listed in Table 1.

Table 1 Features characterizing thoracic CT images

3.3 Random Forest Classification

We capture the maximum inscribed square patch from each superpixel as a region of interest (ROI), resize, and take them as inputs of our DCNN model to extract deep features. Combining with the deep, texture, and intensity features extracted from superpixels, Random forest (RF) [35] classifies the superpixels into three classes: lungs, pleural tissues, and image backgrounds.

The lung segmentation results in the initial classification with an RF are usually incomplete due to the existence of pathologic lung tissues. Because superpixels with different scales can cluster regions with different scales and details, we adopt the RF to classify multi-scale superpixels according to their features and fuse different segmentation results to get a relatively complete lung segmentation result. However, some image patches located at the same positions are assigned different classes in multiple classifications. In order to deal with the issue, we adopt a fusion technology based on a gray correlation algorithm instead of a simple and crude addition fusion.

The similarity of two data series can be determined by their slopes in corresponding periods according to the theory of gray correlation. If the slopes of two series are equal or similar in each period, the two series have a large correlation coefficient and vice versa. Thus, image patch similarity in terms of grayscale and texture is calculated by the gray correlation of two vectors.

Assume two data series are \(\left\{ x_{0}^{(0)}(k), k=1,2, \ldots n\right\} \) and \(\left\{ x_{1}^{(0)}(k), k=1,2, \ldots n\right\} \), the similarity between them is calculated as follows:

  1. 1.

    Fractional accumulation [36]

    We reduce the noise affect by fractional accumulation which is expressed as

    $$\begin{aligned} x^{\left( \frac{p}{q}\right) }(k)=\sum _{i=1}^{k} C_{k-i+\frac{p}{q}-1}^{k-i} x^{(0)}(i) \end{aligned}$$
    (2)

    where \(\frac{p}{q}\left( 0<\frac{p}{q}<1\right) \) is fractional accumulation operator. Let \(C_{\frac{p}{q}-1}^{0}=1, C_{k}^{k+1}=0, k=0,1, \ldots , n-1\), and

    $$\begin{aligned} C_{k-i+\frac{p}{q}^{k-1}}^{k-i}=\frac{\left( k-i+\frac{p}{q}-1\right) \left( k-i+\frac{p}{q}-2\right) \cdots \left( \frac{p}{q}+1\right) \frac{p}{q}}{(k-i) !} \end{aligned}$$
    (3)

    Then \(X^{\left( \frac{p}{q}\right) }=\left( x^{\frac{k}{q}}(1), x^{\frac{p}{q}}(2), \ldots , x^{\left( \frac{k}{q}\right) }(n)\right) \).

  2. 2.

    Initialization

    The purpose of initialization is to make each sequence comparable.

    $$\begin{aligned}&Y_{0}: y_{0}^{\left( \frac{p}{q}\right) }(k)=\left\{ x_{0}^{\left( \frac{p}{q}\right) }(k)- {\bar{x}}_{0}^{\left( \frac{p}{q}\right) } \right\} .\end{aligned}$$
    (4)
    $$\begin{aligned}&Y_{1}: y_{1}^{\left( \frac{p}{q}\right) }(k)=\left\{ x_{1}^{\left( \frac{p}{q}\right) }(k) -{\bar{x}}_{1}^{\left( \frac{p}{q}\right) } \right\} . \end{aligned}$$
    (5)

    where \({\bar{x}}^{(\frac{p}{q})}\) is the mean value of \(x^{(\frac{p}{q})}(k)\).

  3. 3.

    Inverse accumulation

    Inverse accumulation is to find slops of the curve at each time point. Here, we use the absolute values of difference of adjacent data in data series to reflect the overall distribution trend of data series.

    $$\begin{aligned} \begin{array}{lll} \alpha ^{(1)}\left( y_{(0)}^{\left( \frac{p}{q}\right) }(k+1)\right) =\left| y_{(0)}^{\left( \frac{p}{q}\right) }(k+1)-y_{(0)}^{\left( \frac{p}{q}\right) }(k)\right| \\ {\alpha ^{(1)}\left( y_{(1)}^{\left( \frac{p}{q}\right) }(k+1)\right) =\left| y_{(1)}^{\left( \frac{p}{q}\right) }(k+1)-y_{(1)}^{\left( \frac{p}{q}\right) }(k)\right| \quad k=1,2, \ldots n-1}\end{array} \end{aligned}$$
    (6)
  4. 4.

    Correlation coefficient calculation

    The correlation coefficient is calculated by Eq. 7.

    $$\begin{aligned} \xi (k+1)=\frac{1}{1+\left| \alpha ^{(1)}\left( y_{0}^{\left( \frac{p}{q}\right) }(k+1)\right) -\alpha ^{(1)}\left( y_{1}^{\left( \frac{p}{q}\right) }(k+1)\right) \right| } \quad k=1,2, \ldots n-1\nonumber \\ \end{aligned}$$
    (7)
  5. 5.

    Fractional-order gray correlation degree calculation

    Fractional-order gray correlation degree of two data series is defined as Eq. 8.

    $$\begin{aligned} e=-\frac{1}{n-1} \sum _{k=2}^{n} \xi (k) \end{aligned}$$
    (8)

With the fractional-order gray correlation degree, we can distinguish any two image patches and assign them an appropriate class in classification. The segmentation fusion is described detailedly in Algorithm 1.

figure a

In Fig. 4, we set the numbers of superpixels, Ps, in SLIC method to 200, 300, and 500, respectively, which generates three scales of superpixels. Because lung tissues and image backgrounds have low intensities, which are very different from pleural tissues with high intensities, we group lung tissues and image backgrounds into one class and pleural tissues into another class in the output of our algorithm. By fusing the different classification results with Algorithm 1, we obtain a relatively complete segmentation.Additionally, because we use three small-sized Ps, the time cost is much less than that with a big-sized P. For example, in Fig. 4, the time cost of the fusion algorithm with \(P_1=200\), \(P_2=300\), and \(P_3=500\) is 0.68 times of that with \(P=2000\). At the same time, the fused segmentation is better than one-time segmentation.

We obtain the initial lung segmentation results by extracting lung regions from the fusion results of RFs with morphological operations.

Fig. 4
figure 4

Segmentation fusion. We segment lungs from a pathologic thoracic CT image (a), fuse the classification results of an RF with \(P_1=200\) (b), \(P_2=300\) (c), and \(P_3=500\) (d). Our fusion result (e) is more completed compared with one-time segmentation with \(P=2000\) (f)

3.4 Segmentation Refinement

Pathologic lungs usually have similar densities with surrounding tissues, which often leads to incomplete segmentation of lung regions. Additionally, the left and right lungs closed to each other usually connected in some segmentation results. In order to tackle these problems, we utilize a lung contour correction approach and lung region repair method to refine segmentations after lung separation, see Fig. 5 for a reference.

Fig. 5
figure 5

Segmentation refinement. It covers three main steps: lung separation, contour correction of left lung and region repair of the right lung

As shown in Fig. 6, in some pathological lung CT images, the right lung regions usually have larger defects and the left lung regions are relatively more complete. Accordingly, we adopt a divide-and-conquer strategy to refine segmentation results by combing contour correction of left lung contours and region repair of the right lungs. We first presented a contour correction approach based on SUSAN operator for left lung contour correction. Because the left lung and the right lung are approximately symmetrical and they usually have similar contours in a thoracic CT image, the right lung region is repaired supervised by left lung contours. It should be noted that the segmentation refinement methods for the left and right lungs can be switched in the other general images. Anyway, the first step of our segmentation refinement is to separate two connected lungs.

Fig. 6
figure 6

A set of lung segmentation results. The right lung regions usually have larger defects and the left lung regions are usually more complete in some initial lung segmentation results

3.4.1 Lung Separation

In a thoracic CT slice sequence \((I_1\cdots , I_k, \cdots I_L)\), left and right lungs are closed to each other and may connected from \(I_k\) in a segmentation sequence \((I_k\cdots , I_{k+1}, \cdots I_{k+n})\). Here, we utilize a lung separation line propagation approach to separate the connected lungs by using the separation line of its former one.

3.4.2 Contour Correction of Left Lung

The Smallest Univalue Segment Assimilating Nucleus (SUSAN) algorithm is a famous corner detection technique. In order to repair lung contour concaves caused by pathological abnormalities, we used a contour correction approach based on SUSAN operator [37] for lung contour correcting, which covers three stages: candidate corner detecting, validate corner filtering and corner connecting. We create a copy, \(M_1\), of lung mask M, remove the right lung region from \(M_1\) and retain left lung. The main steps of contour correction are described in Algorithm 2.

Fig. 7
figure 7

Segmentation refinements. For initial segmentations (a), we extract left lungs and correct their lung contours with a corner detection-based method (b). We capture a segment of left part contours from the left lungs, stretch them (c) to supervise right lung correction (e). With contour correction, we obtain the final corrected results (f)

figure b

3.4.3 Region Repair of Right Lung

The left lung and right lung in a thoracic CT image usually have similar contours. Accordingly, instead of using the whole mask of the left lung to repair the right lung, we take a part of corrected left lung contour as the contour mask of the right lung. We divide the repair process into three stages: extraction of reference contour, repair of the right part, and repair of the left part of the lung.

  1. a.

    Extraction of reference contour

    1. (1)

      Create a copy, \(M_2\), of lung mask M. Remove the left lung region from \(M_2\) and retain the right lung.

    2. (2)

      Extract single-pixel wide contour, E, of the corrected left lung in \(M_1\).

    3. (3)

      Calculate bounding box, \(B_l\), of the left lung and capture the left part of E, \(E_l\), according to the top and bottom points of \(B_l\), see Fig. 7 for a reference.

    4. (4)

      Extract the middle segment (reference contour) of \(E_l\).

  2. b.

    Repair of right part

    1. (1)

      Flip \(E_l\) horizontally, move it to the right until \(E_l\) reaches the right boundary of right lung \(M_2\).

    2. (2)

      Stretch \(E_l\) from its two endpoint according to its slope, respectively.

    3. (3)

      Merge \(M_2\) with E: \(M_2=M_2\cap E\).

    4. (4)

      Fill holes of \(M_2\) with morphological operations.

  3. c.

    Repair of left part

    1. (1)

      Calculate bounding box, \(B_r\), of the right lung \(M_2\).

    2. (2)

      Extract left part of right lung according to the top and bottom points of \(B_r\).

    3. (3)

      Correct the left part of right lung contours with Algorithm 2 in \(M_2\).

After the contour correction of left lung and region repair of the right lung, we merge their results and obtain the final refined segmentation: \(M=M_1\cap M_2\).

4 Experiments

Interstitial lung disease (ILD) is the leading cause of mortality and characterized by widespread fibrotic and inflammatory abnormalities of lungs [38]. The images used for training and testing our algorithm come from ILDs database [39]. It contains high-resolution computed tomography (HRCT) image series with pathologically proven diagnoses of ILDs. The images which are difficult to segment by conventional methods are used to test the performance of our segmentation approach.

In the following, we first introduce metrics in our algorithm evaluation and then evaluate the influence of preset superpixel numbers. Finally, we perform our algorithm on a set of pathological lung CT images and analyze the experiment results.

4.1 Evaluation Method

TNR, TPR, PPV, ACC, Error (Er) are used to evaluate our lung segmentation performance [40]. Four classical metrics [41], over-segmentation rate (OR), under-segmentation rate (UR), Dice similarity coefficient (DSC) and Jaccard’s similarity index (JSI) are also utilized for evaluating the performance of our algorithm. Larger values of TNR, TPR, PPV, ACC, DSC and JSI and smaller values of Er, UR and OR indicate more accurate segmentation of lung images.

4.2 The Influence of Preset Superpixel Number

The preset superpixel number, P, affects the efficiency and accuracy of our segmentation algorithm greatly. In order to detect a suitable group of Ps, we randomly select 80 CT images affected with sarcoidosis to test our algorithm efficiency and accuracy with different Ps. We test four groups of Ps: \(G_1\) (\(P_1=100\), \(P_2=200\), \(P_3=300\)), \(G_2\) (\(P_1=300\), \(P_2=400\), \(P_3=500\)), \(G_3\) (\(P_1=500\), \(P_2=600\), \(P_3=700\)), \(G_4\) (\(P_1=700\), \(P_2=800\), \(P_3=900\)) and the corresponding average time costs are depicted in Fig. 8. For example, the average time coats on \(P_1\)=100, \(P_2\)=200 and \(P_3\)=300 are 14.2s, 21.8s, and 27.8s, respectively. The total time cost of \(G_1\) is 64.93s. The segmentation is implemented in MATLAB environment on a computer with Intel i7-6500U CPU and 16GB of RAM. It can be seen, the time cost increases with the increase of P.

In Fig. 9, we show the performance of our algorithm with four groups of Ps. We can see that the segmentation accuracies in terms of DSC and JSI increase with an increase of Ps. Nevertheless, the DSCs are all above 97%. In order to reduce the time cost while keeping a high segmentation accuracy, \(P_i\) is ranged from 100 to 700.

Fig. 8
figure 8

Time costs (s) of our algorithm on four groups of Ps

Fig. 9
figure 9

Performance of our algorithm in terms of DSC and JSI on four groups of Ps (G1\(\sim \)G4)

4.3 Result Analysis

We chose a challenging data set of thoracic CT images containing density pathologies in varying degrees of severity from ILDs database to test the performance of our lung segmentation algorithm: images affected with ground glass (GG), fibrosis (F), nodules (N), reticulation (R) and PCP (P). The number of the tested images in each type thoracic CT images ranged from 40 to 80, and we list the average performances of our algorithm in terms of TNR, TPR, PPV, ACC and Er in Table 2, and DSC, JSI, OR and UR in Table 3. It can be seen that the segmentation accuracy of R images has the highest accuracy such as TPR because of relatively homogeneous intensities in the images. Conversely, P images have the lowest accuracy among our segmentations due to the severe uneven distribution of gray scales.

Table 2 Metric values of segmentation performance of our algorithm in terms of TNR, TPR, PPV, ACC and Er
Table 3 Metric values of segmentation performance of our algorithm in terms of DSC, JSI, OR and UR

The average metric values (AVE) in terms of TNR, TPR, PPV, ACC and Er is shown in Fig. 10, and average DSC, JSI, OR and UR in Fig. 11. It can be seen that ACC and DSC are 0.9880 and 0.9645, respectively. The average PR curve of the five types of pathologic lung images is shown in Fig. 12.

Fig. 10
figure 10

Average metric values of segmentation performance of our algorithm in terms of TNR, TPR, PPV, ACC

Fig. 11
figure 11

Average metric values of segmentation performance of our algorithm in terms of DSC, JSI, OR and UR

Fig. 12
figure 12

Average Precision-Recall curve

Fig. 13
figure 13

Error distance

Fig. 14
figure 14

Comparison of corrected and initial segmentation results of CT images of PCP

Fig. 15
figure 15

Lung segmentations with different methods. ae are original lung CT images, lung segmentation results with U-Net,random forest and our algorithm, respectively

Fig. 16
figure 16

Comparison of segmentation accuracy of different methods. Our algorithm achieves higher accuracy than other methods in terms of DSC, JSI, TPR, and ACC

We randomly select 270 images from the above mentioned five types of images and the error distance defined as \((FN+FP)/(FN+TP)\) is shown in Fig. 13. Error distances are mainly concentrated in a range of (0, 0.2). The P images usually be destructed by pneumocystis pneumonia (PCP) and the intensities are inhomogeneous. Fig. 14 depicts a comparison of the segmentation accuracy of initial segmentation and refined segmentation of lungs in P images. It can be seen that the accuracies of the refined segmentation results in terms of TPR, ACC, DSC, and JSI are all higher than those of the initial segmentation results.

In Fig. 15, we show a group of segmentation results with some state-of-art methods [13, 18]. It can be seen that the other methods cannot deal with some intensity inhomogeneous images well. Conversely, our algorithm can segment more complete lungs from pathological lung CT images. In Fig. 16, we compare our algorithm with the other methods on a set of CT images affected with fibrosis and PCP. The results show that our algorithm achieves higher accuracies in terms of DSC, JSI, TPR, and ACC, which implies the segmentation ability of deep learning networks [18] on pathological lung CT images is hindered by a limited size of a training set (1000). The segmentation accuracy of a single-scale RF method [13] is also smaller than ours on the whole. Our algorithm is relatively more robust on lung segmentation than other methods because of the multi-scale classification fusion and effective post-processing operations.

5 Discussion

The fusion algorithm can generate high segmentation accuracy. However, it may generate a little bit higher OR shown in Fig. 11 due to two main factors. One factor arises from superpixel segmentation. We take superpixels as the input of the RF classifier. In order to accelerate our algorithm, we use three Ps with smaller values, which generates a group of multi-scaled superpixels with the SLIC approach. However, a bigger superpixel is difficult to capture narrow regions and sometimes brings about high ORs, see Fig. 17b for a reference. Accordingly, OR of fusion result is high by simply adding the multi-scale segmentation results together. In our algorithm, the correlation approach can solve this tissue and reduce OR. Another factor lies in region repair of right lung supervised by the left lung contour, and the over-segmentation of the left lung may propagate to the right lung. In order to overcome the issue, we capture a middle segment of the left lung contour and stretch it according to its local slope to match the right lung contour, which reduces the probability of wrong contour propagation, as well as OR.

Fig. 17
figure 17

An example of over-segmentation by our algorithm. For a thoracic CT image (a) with narrow lung regions, bigger-scale superpixels used in our multi-scale segmentation may result in over-segmentation (b) compared with the corresponding ground truth (c)

6 Conclusion

Pathologic thoracic CT image segmentation is a challenging issue for the presence of inhomogeneous intensities and various abnormalities. In this paper, we introduce a novel pathologic lung segmentation algorithm based on the DCNN model and random forest. In our algorithm, two fusion operations enhance the performance of our segmentation algorithm. First, the fusion of features of deep, texture and intensity enriches the classification information, which contributes to the classification of random forest. Second, the fusion of classification results of random forest based on multi-scale superpixels can produce more accurate results than one-time segmentation with relatively low time cost due to small-sized Ps. In order to improve segmentation accuracy, we introduce a divide-and-conquer strategy for refining segmentation results by correcting and repairing the left and right lungs. In our future work, we will improve the fusion techniques, augment training samples, and further enhance the performance of our lung algorithm on pathologic lung segmentation.