3D segmentation of lungs with juxta-pleural tumor using the improved active shape model approach

2.1Marker points extraction

First, the traditional method is used to segment the pulmonary region of the samples in training set and obtained lung boundary contour. The segmentation steps include extra-pulmonary region exclusion by iteration threshold and region growing segmentation methods, trachea elimination, internal cavity filling, liver exclusion and left and right lungs separation [16].

Second, extract marker points on left or right lung boundary contour, respectively according to Shenshen et al. [17]. In this way, there are 5 main marker points and 90 marker points in the left or right lung boundary (as shown in Fig. 3).

Figure 3.

Extracted marker points. Main marker points are shown in red; other marker points are shown in green. (a) normal lung (b) pathological lungs with juxta-pleural tumor.

2.2Shape model established

Aligned lung boundaries in training sample dataset by procrustes analysis method. In order to remove the influence of the noise contours such as Fig. 3b, PCA is replaced by RPCA.

The boundary contours of images in the training data set is d1,d2,…,dn, represented in vector form as D=(d1,d2,…,dn)T. di=((xi⁢0,…,xi⁢j⁢…,xi⁢m)T,(yi⁢0,…,yi⁢j,…,yi⁢m)T)T is marker points’ coordinate, in which, (xi⁢j,yi⁢j) is the coordinates of the j marker point of the i⁢t⁢h image in the training data set.

Firstly, D is processed by low rank.

where X is the decomposed low rank component and E is the sparse component. Equation (2.2) is a non-convex function with no unique solution. So in order to let Eq. (2.2) be a convex function,  rank(X) is replaced by the kernel norm and the L0 norm is replaced by the L1 norm.

where λ is a regularized parameter which is greater than 0, and δ is a very small positive number. ||||F is Frobenius norm. Equation (2.2) can be written in the form of Eq. (3) by Lagrange multiplier method.

where μ is a small positive number. The accelerated approximation gradient [18] is used to solve Eq. (3). Then, the mean vector and covariance matrix of X are obtained. The mean vector is mx=E⁢{X}. Covariance matrix is Σ=E⁢{(X-mx)⁢(X-mx)T}. Secondly, the covariance matrix Σ is decomposed by low rank.

where Φ is the decomposed low rank component and γ is the sparse component. Equation (2.2) is a non-convex function with no unique solution. So in order to let Eq. (2.2) be a convex function, rank(X) is replaced by the kernel norm, and the L0 norm is replaced by the L1 norm.

Equation (2.2) can be written in the form of Eq. (6). The accelerated approximation gradient [18] is used to solve Eq. (6).

Calculating the eigenvalues of Φ, and arranged them in descending order and the eigenvectors p=[p1;p2;…;pt]T corresponding to top-rank t, t is set 6 in this paper, eigenvalues b=[b1,b2,…,bt] are selected to construct the shape model X=mx+p⁢b, where mx is the mean shape vector. The shape model is shown in Fig. 4. Figure 4a is the shape expressed X by p=[p1]T, Fig. 4b is the shape expressed X by p=[p1;p2]T, Fig. 4c is the shape expressed X by p=[p1;p2;p3]T, Fig. 4d is the shape expressed X by p=[p1;p2;p3;p4]T.

Figure 4.

The mean shape mx of training sample set including noise samples is shown in blue, the deformed shape generated by the shape model based on RPCA of ASM is shown in red.

2.3Model solution

The accelerated approximation gradient [18] is used to solve Eq. (2.2).

The accelerated approximation gradient [18] is used to solve Eq. (2.5.2).

2.4Local appearance feature model construction

In this paper, a local appearance feature model is established based on PCA using gray features. The gray-level appearance model that describes the typical image structure around each marker point is obtained from pixel profiles, sampled (using linear interpolation) around each marker point, perpendicular to the boundary contour. On either side k pixels are sampled (step size = 1), which gives profiles of length 2⁢k+1, in this paper, k=8.

For each marker point j (j∈[1,95]), gray values of the sampled points i(i∈[1,17]) which is on the profile perpendicular to the contours of every images in the training data set constitutes the local appearance gray feature matrix gj. The PCA model of local appearance gray feature matrix gj corresponding to each landmark is established. The eigenvectors gj.𝐸𝑣𝑒𝑐𝑡𝑜𝑟, eigenvalues gj.𝐸𝑣𝑒𝑐𝑡𝑜𝑟𝑠 and mean vectors gj.𝐸𝑚𝑒𝑎𝑛 of the PCA model are calculated.

2.5Target contour search

2.5.1Recognition of abnormal marker point

Because of the adhesion between tumor and chest wall, concave lung contour is obtained based on the traditional ASM method, as shown in Fig. 5 in blue. The marker points on the contour notch are the abnormal marker points. These abnormal marker points are updated to the boundary between lung parenchyma and tumor. Therefore, the different search function of abnormal marker points should be given for abnormal marker points.

Firstly, lung contour with juxta-pleural tumor 𝑆ℎ𝑎𝑝𝑒1⁢(x1,y1,…⁢xk,yk) is aligned based on marker points and the mean contour of training samples 𝑆ℎ𝑎𝑝𝑒2⁢(w1,z1,…⁢wk,zk). The alignment steps are: translation components deletion, scaling and rotation transformation. The formula for deleting translation components is following: (X1-∑i=1kXi/K,y1-∑i=1kyi/k). The scaling scale formula is: ((x-x¯)/s,(y-y¯)/s), where, s=𝑠𝑞𝑟𝑡⁢(((x1-x¯)2+(y1-y¯)2+…)/k). The formula of rotation angle in rotation transformation is as follows:

θ=arctan⁡∑i=1k(wi⁢yi-xi⁢zi)∑i=1k(wi⁢xi+yi⁢zi),

the alignment results are shown in Fig. 5.

Figure 5.

Align the mean contour of the training samples which is shown in red and lung contour with juxta-pleural tumor which is shown in blue.

Figure 6.

The corresponding points distance between the lung contour with juxta-pleural tumor and the mean contour of the training samples.

If the distance of the corresponding marker point between the lung contour with juxta-pleural tumor and the mean contour of the training samples is bigger than the threshold value ε, then the marker point is an abnormal marker point, which is shown in Fig. 6.

3.1Experiment data set

First, the training dataset is established. The training dataset consists of 135 multi-slice spiral CT lung scan images collected from the First Affiliated Hospital of Guangzhou Medical University. And among them, sixteen scans are included juxta-pleural tumor on the left side and fourteen scans are included juxta-pleural tumor on the right side. The size of the lung scan images is 512*512*100 ∼ 200 layers, and their thickness is 0.6 mm*0.6 mm*2 mm (or 3 mm). The test dataset consists of 122 scans from EMPIRE 10 and LIDC-IDRI [19]. It had the gold standard of lung contour and among them, 52 scans are included juxta-pleural tumor on the left side, 30 scans are included juxta-pleural tumor on the right side. Figure 7 is three scans from training dataset. Figure 7a is CT scan without any tumor on the lung wall. Figure 7b and c are scans which include juxta-pleural tumors. The juxta-pleural tumor in Fig. 7b is located in the upper lobe and the juxta-pleural tumor in Fig. 7c is located in the lower lobe. Figure 7b and c can be regarded as noise samples in the training data set, because their lung contours described in green lines are different from those in Fig. 7a.

Figure 7.

Training samples including noise samples, the five main marking points are red, the marking points are green. (a) marker points in lung contour with no juxta-pleural tumor. (b) and (c) marker points in lung contour with juxta-pleural tumor.

3.2Parameter setting

In low rank algorithm, λ= 0.129, δ= 0.9. Maximum number of iterations is 10000. Tolerance for stopping criterion is 10-7. In PCA algorithm, number of searching iterations is 500. Number of contour points interpolated between the major landmarks is 20.

Figure 8.

Mean shape of training sample set including noise samples is shown in blue, the deformed shape based on RPCA of ASM is shown in red, when t=1∼4 respectively. The shape model is shown in Fig. 9 using ASM of PCA for the training data set with noise samples. The code of traditional ASM was used [7].

3.3Comparisons of generated shape models

The shape model is shown in Fig. 8 using ASM of RPCA for the training data set with noise samples. It can be seen that the larger t is, the more similar the shape obtained is to the average shape. And the code of low rank was also used [20]. The shape model is shown in Fig. 10 using ASM of PCA for the training data set without noise samples. When t= 1, the shape is represented only by the eigenvector corresponding to the maximum eigenvalue. It can be seen that Fig. 8a is most similar to the mean shape and Fig. 10a is more similar to the mean shape, while Fig. 9a is least similar to the mean shape. When t= 2, the shape is represented by the eigenvectors corresponding to the first two eigenvalues. It can be seen that Fig. 7b and Fig. 9b are more similar to the mean shape, while Fig. 8b is the least similar to the mean shape. When t= 4, the shape is represented by the eigenvectors corresponding to the first two eigenvalues. It can be seen that Fig. 7b and 9b are similar to the mean shape, while Fig. 8b is the least similar to the mean shape. The red shapes of (a), (b), (c) in Fig. 8 which is the trained shape model are no longer within the normal lung shape range, while the red shapes of Figs 7 and 10 have been within the normal lung shape range. It is shown that the shape model established in Fig. 8 is affected by noise samples. Because of the application of RPCA, the shape model established in Fig. 7 is not affected by the noise samples, which is basically the same as the shape model established in training data set without noise samples in Fig. 10. It can be seen that, to a certain degree, RPCA can eliminate the influence of noise samples on the shape model.

Figure 9.

Mean shape of training sample set including noise samples is shown in blue, the deformed shape based on PCA of ASM is shown in red, when t=1∼4 respectively.

Figure 10.

Mean shape of training sample set without noise samples is shown in blue, the deformed shape based on PCA of ASM is shown in red, when t=1∼4 respectively.

3.4Evaluation of segmentation results

The contour containment region obtained by the segmentation algorithm is compared with the contour containment region manually drawn by the doctor (as the golden standard), and the coverage degree overlap is calculated.

While R𝑔𝑜𝑙𝑑 and R𝑛𝑒𝑤 refer to the gold standard pixel set and the pixel set obtained by segmentation algorithm respectively. When the training sample set contains noise samples, 122 lung images of the test set are segmented. The results are as follows: the average overlap degree of lung region segmentation based on RPCA ASM is 94.5%, while that of traditional ASM [4] is 69.5%. Compared with the method of Shao et al. [15], the average rate of segmentation overlap is 85% because the influence of training set noise samples is not considered, as shown in the Fig. 11. The traditional ASM [4] is compared with the RPCA based ASM segmentation effect for juxta-pleural tumor lung, as shown in Fig. 12. That is because traditional ASM give the same search function to all marker points, while RPCA give different search functions to the marker points which locate in the place of tumor and pleura connection.

Figure 11.

The result of lung segmentation. (a) [15]’s method (b) ASM based on RPCA.

Figure 12.

Lung segmentation processing with juxta-pleural tumor by the ASM based on RPCA.

Figure 13.

Compares the lung segmentation results with juxta-pleural tumor based on ASM based on PCA and ASM based on RPCA; (a) ∼ (c) the segmentation results of ASM based on RPCA; (d) ∼ (f) to the segmentation results of traditional methods.

Figure 14.

Lung segmentation result with tumor by ASM based on RPCA.

Figure 15.

3D visualization results. (a) ASM based on PCA; (b) ASM based on RPCA.

The segmentation results were compared with the ground truth according to TPF (true positive fraction) under three situations, shown as Table 1. The first situation is the train dataset doesn’t contain noise samples and using traditional ASM; the second situation is the train dataset includes noise samples and using traditional ASM; the third situation is the train dataset includes noise samples and using RPCA-based ASM. From Table 1, it can be seen that RPCA-based ASM is largely unaffected by the training set noise images. The relationship between noise samples’ number in the training dataset and segmentation accuracy for RPCA-based ASM is shown in Table 2. It can be seen that with the increase of the number of noise samples, the accuracy of segmentation is decreasing. Figure 12 shows the initial position, middle position and final position of the lung contour searching process by the proposed method. Figure 13 compares the ASM based on RPCA and traditional ASM method in lung segmentation with juxta-pleural tumor. It can be seen that the proposed method is also useful for large juxta-pleural tumors. Figure 14 shows the effect of the lung segmentation with tumor by the ASM based on RPCA. Figure 15 shows 3D visualization results of ASM based on PCA and ASM based on RPCA. In 3D visualization results, ASM based on PCA cannot include juxta-pleural tumors into the lung parenchyma, so a hole is shown.

3.5Computer time

The proposed algorithm was implemented in C++ and Matlab Lib and test on a PC with Intel i7-8550U CPU and 8GB of RAM. The average computing time of the algorithm was 5 ± 29 s.