An automatic detection method for lung nodules based on multi-scale enhancement filters and 3D shape features

Abstract

In the computer-aided detection (CAD) system of pulmonary nodules, a high false positive rate is common because the density and the computed tomography (CT) values of the vessel and the nodule in the CT images are similar, which affects the detection accuracy of pulmonary nodules. In this paper, a method of automatic detection of pulmonary nodules based on multi-scale enhancement filters and 3D shape features is proposed. The method uses an iterative threshold and a region growing algorithm to segment lung parenchyma. Two types of multi-scale enhancement filters are constructed to enhance the images of nodules and blood vessels in 3D lung images, and most of the blood vessel images in the nodular images are removed to obtain a suspected nodule image. An 18 neighborhood region growing algorithm is then used to extract the lung nodules. A new pulmonary nodules feature descriptor is proposed, and the features of the suspected nodules are extracted. A support vector machine (SVM) classifier is used to classify the pulmonary nodules. The experimental results show that our method can effectively detect pulmonary nodules and reduce false positive rates, and the feature descriptor proposed in this paper is valid which can be used to distinguish between nodules and blood vessels.

1. Introduction

Lung cancer has one of the highest mortality rates of cancer throughout the world [1]. Most lung cancers manifest as solitary pulmonary nodules (SPNs) in the early stages, which are solitary, round, nodular lesions in the lung parenchyma with a maximum diameter of 30 mm. These lesions are not associated with lymph node enlargement, atelectasis, and pneumonia. Computed tomography (CT) has excellent density resolution for human lung lesions, and is the most effective and direct imaging method for the early diagnosis of lung cancer. With improvements in the imaging accuracy of the focus area, the thickness of the CT scan is decreasing, whereas CT image data is undergoing explosive growth. The massive amount of CT image data is increasing the workloads of doctors. Computer aided detection (CAD) systems can effectively reduce workloads and ensure the manual review of CT scans. However, a critical factor of the performance of lung cancer CAD systems is high false positive and low sensitivity, which will reduce the detection efficiency, and aggravate the burden on doctors. Therefore, a method for accurately detecting pulmonary nodules is significant for reducing the missing rates and false positive rates ensuring the early diagnosis and early treatment of lung cancer patients and improving their survival rates.

A lung cancer CAD system usually contains the following processes: lung segmentation, nodules detection, feature extraction and classification. In recent years, numerous scholars have investigated the detection methods of a lung cancer CAD system. Han et al. proposed a CAD system [2] based on hierarchical vector quantization (VQ) to detect pulmonary nodules in the early stage of cancer. High-level VQ was used to segment the lungs from CT images and low-level VQ was employed to detect and segment the nodules. Rule-based filtering was performed to select features for training SVM classifiers. Tabakhi et al. proposed a feature selection approach using ant colony optimization (ACO) [3], in which a cosine similarity measure (CSM) served as the filter evaluation function. Features with high pheromone values and low similarity values in the partially constructed subsets were added to the subsets. In 2012, Choi’s group applied 3D connected component labeling method to segment a lung body image and utilized an optimized multi-threshold and rule-based pruning technique to detect candidate nodules and extract their 2D and 3D features. The pulmonary nodules were classified by a genetic programming classifier (GPC), which can achieve a false positive rate of 5.45 and an accuracy of 94.1% [4]. In 2013, Choi’s group proposed a lung nodule detection method based on the theory of information entropy analysis [5]. They segmented a CT image into a three-dimensional image block and extracted the information entropy to segment suspected pulmonary nodules and employed a SVM to classify and identify nodule images. Compared with the previous experimental results, an accuracy of 95.28% and a false positive rate of 2.27 can be attained with this method. Stember [6] proposed a pseudo molecular motion model to measure the shape of an object. The model employs geometric characteristics to guide and visualize the internal motion of object and simulate and identify lung nodules. The experimental results revealed that the classification of round, oval and irregular graphics can achieve an accuracy rate of 97%. Ye X apply a fuzzy threshold to segment the lung and extract local shape features and partial dispersion information to describe pulmonary nodule images [7]. This method can detect solid nodules and ground glass nodules with an accuracy of 9 0.2% accuracy and a false positive rate of 8.2. Santos AM proposed a small nodule (diameter range of 2-10 mm) detection method that employs the Gauss mixture model [8]. The texture features measured by information entropy are extracted and the method can attain a classification accuracy rate of 88.4% and a false positive rate of 1.17.

Regarding the aspect of the feature description for pulmonary nodules, Lee et al. extracted global shape features, such as sphericity, compactness and narrow length of candidate nodules, to distinguish true and false positive nodules [9]. Based on a gray histogram, Messay et al. extracted the gray mean, variance, kurtosis and skewness as feature vectors for the classification of true and false positive nodules [10]. Skibbe H proposed a HOG operator based on the probability density function to detect 3D spherical objects [11]. Mikolajczyk K comparatively analyzed the advantages of a scale invariant feature transform (SIFT) and its variants to extract image features [12]. Zechao Li proposed a novel unsupervised feature selection algorithm by integrating cluster analysis and sparse structural analysis into a joint framework [13-15].

Although many scholars have researched pulmonary nodules detection, the accuracy and effectiveness of a detection system can not satisfy the clinical needs because the vascular density and CT values are similar to those of nodules and because the nodules and blood vessels usually cross within a lesion area. These conditions can cause false positives and low accuracy. Therefore, we consider the differences in the shapes and structures of nodules and vessels in 3D images and propose a method for automatic detection of pulmonary nodules based on multi-scale . Compared with existing methods, this method can effectively reduce false positives in the detection of solitary pulmonary nodules and achieve high sensitivity, specificity and accuracy.

2. Materials and methods

2.1 Materials

The sequence of the CT image data used in this paper is derived from the LIDC database. All data can be accessed at https://figshare.com/s/39ff2ec6469fdc354a6c. In the experiment, we selected data from 90 case, each case contains an average of 280 CT images, the thickness range for the CT layers is 0.5 mm-3.75 mm, and the image size is 512 * 512. The total number of nodules labeled by statistical experts is 165, and the diameter range is 3 mm-30 mm. The total number of suspected nodules detected by this method is 568.

2.2 Proposed Method

Our method primarily involves segmenting pulmonary parenchyma image sequences using a method of iterative threshold and a region growing algorithm. Two types of multi-scale enhancement filters are constructed for enhancing the 3D nodule images and vascular images and most of the blood vessel images in the nodular images are removed to obtain suspected nodule image. An 18-neighborhood region growing algorithm is used to extract lung nodules, and novel pulmonary nodule feature descriptors are proposed for classification by a SVM. A flow chart of the method is shown in Fig. 1.

Fig. 1. Flow chart of lung nodules detection  

2.2.1 Lung segmentation

In the CADs of pulmonary nodules, accurate segmentation of lung parenchyma is critical to the detection and classification of pulmonary nodules. First, binarize the lung image sequences with the best threshold iteration method. Second, extract the region of interest (ROI) of pleural and segment lung parenchyma with the region growing algorithm. Last, modify the lung contour using morphological methods to obtain the final segmentation mask and lung parenchyma image sequences. The specific steps are given as follows:

Step 1: Binarize for the lung image sequences.

Step 2: Extract ROI and the minimum bounding rectangle of the lung.

Step 3: Select seed points using the LRS and FCRS algorithm.

Step 4: Refine the lung contours with erosion and dilation.

Step 5: Acquire the final lung mask sequences.

1. Binarization and ROI extraction

In traditional medical image threshold segmentation, threshold selection is usually based on empirical values. Because the CT values of lung parenchyma and other tissues are similar, the selection of a fixed global threshold for image binarization is difficult. In this paper, we adopt an optimal iterative threshold algorithm to dynamically access the threshold. The steps are as follows:

Step 1: Initialize the threshold value T₀.

Step 2: Classify the global image into two pixels set using T: the set B of body voxels and the set N of nonbody voxels.

Step 3: Calculate the mean gray value μb of B and the mean gray value μn of N.

Step 4: Calculate the new threshold T according to Eq. (1)

T = (μ_b + μ_n)/2         (1)

Step 5: Repeat Step 2~4 until |T_n − T_n+1| is less than the preset parameter, where T_n represents the threshold at iteration n, and T_n+1 denotes the threshold at iteration n + 1. The threshold for the last iteration is T_opt.

In the CT image sequence, the CT value for the air is approximately -1000 Hounsfield units (HU). The CT value of most lung tissues ranges from -500HU to -910HU [16]. In this paper, -500 HU is the initial threshold T₀, and the optimal threshold T_opt is obtained by the threshold iterative algorithm, for image binarization. Because the lung parenchyma region in the CT image of human lungs is relatively fixed in a CT image, that is, it has position relative invariance, this paper also extracts the lung chest ROI of the CT image, which can effectively reduce the time of image processing. The results of the image binarization and ROI extraction of the original lung CT image are shown in Fig. 2 (d)-(f).

Fig. 2. Segmentation of the lung parenchyma image 

(a)-(c) are original lung images; (d)-(f) are the results of image binarization and ROI extraction; (g)-(i) are coarse segmented lung images; (j)-(l) are final lung parenchyma mask.

2. Lung image segmentation

In lung CT image sequences, the morphological changes vary significantly from the top to the bottom of the lung. The image of the top of the lung usually contains the trachea, whereas the image of the bottom of the lung usually exhibits a four-part diagonal distribution. To address the particularity of the lung image, this paper adopts a region growing algorithm, which can realize the adaptive segmentation of the CT image of the top of the lung, the middle of the lung and the bottom of the lung. The segmentation steps are listed as follows:

Step 1: Find the smallest circumscribed rectangle of the lung region in the chest ROI image.

Step 2: Select seed points of the left and right lung parenchyma images. LRS algorithm [17] is used for images of the top of the lung and the middle of the lung; FCRS algorithm [17] is used for images of the bottom of the lung.

Step 3: Use the improved regional growth method to segment the lung parenchyma images and remove the trachea and bronchi.

Step 4: Optimize the lung parenchyma mask image using corrosion and dilation operations in the morphology.

Step 5: Obtain image sequences of lung parenchyma. Using the LRS algorithm to find seed points of the top of the lung and the middle of the lung, the minimum circumscribed rectangular image is simultaneously scanned along the left and right sides of each pixel. When n continuous white pixels are observed along the Y direction, stop scanning, and establish the middle pixels are as seed points.

The process of the FCRS algorithm is to seek the seed points of the bottom of the image as shown as Fig. 3.

Fig. 3. Process of FCRS algorithm to seek the seed points of the bottom of the lung

Set the four vertexes of the circumscribed rectangular as the base points, and set x = x₁, y =y₁, mid_x = (x₁+x₂)/2 and mid_y = (y₁+y₂)/2 as the boundary lines to radiate the rotating radiation. When the upper left corner region (x₁, y₁) is set as the rotation point, the equation of the rotating radiation in the upper left corner is formula (2):

\(y-y_{1}=\tan (\theta)\left(x-x_{1}\right), \theta \in[0, \pi / 2) \text { and } x_{I}<=x<=m i d_{-} x \text { and } y_{1}<=y<=m i d\text_ y\)       (2)

where θ gradually increases from 0. For each pixel in the area of the upper left corner, when n continuous white pixels are observed along the rotating radiation direction, stop scanning, and record the distance between the first pixel of n consecutive white pixels and the rotation point distance d(θ). The pixel points that correspond to the minimum \(d(\theta), \theta \in[0, \pi / 2)\) are obtained as the lung parenchyma seed points of the upper left corner. If no points are eligible, no lung parenchyma seed point exists in this area. The method of scanning the lung parenchyma seed points in the upper right corner, lower left corner, and lower right corner is the same as the method of scanning in the upper left corner. 

After the seed point of lung parenchyma is obtained, the regional growth method is used to segment the lung parenchyma image. The lung parenchyma mask is repaired by morphological operation, and the final lung parenchyma image is obtained using the mask to segment the original image. The results of pulmonary parenchyma segmentation using the regional growth method are shown in Fig. 2 (g)-(i). Remove the tracheal bronchus and optimize the mask to obtain the final lung parenchymal mask as shown in Fig. 2 (j)-(l).

In addition, Fig. 4 gives the results of the 3D reconstruction of the lung parenchyma sequence images using the Visualization Toolkit (VTK).

Fig. 4. Front view, rear view and top view of the 3D lung image

2.2.2 Suspected nodules detection

After obtaining lung parenchyma 3D image sequences, we need to detect suspected pulmonary nodules. Shape-based selective filtering has been extensively employed in pulmonary nodule detection system. In this paper, a method for the detection of suspected pulmonary nodules based on multi-scale enhancement of a Hessian matrix is proposed. The 3D models of pulmonary nodules and blood vessels are constructed. Two types of multi-scale enhancement filters based on 3D shapes using a Hessian matrix are constructed to enhance the spherical objects and cylindrical objects, which can enhance pulmonary nodule images and vascular images respectively. The enhancement of the lung nodule and blood vessel image at differents cales is achieved using different values of σ in Gauss's functions. The blood vessel image is removed from the enhanced pulmonary nodule imagine. The 3D lung images of suspected pulmonary nodules are obtained.

1. Nodule and vessel model

In the 3D images of the lung, the nodules and blood vessels exhibit spherical and tubular characteristics, respectively, in the morphology, and the 3D images of the nodules and blood vessels can be considered a sphere and a cylinder, respectively. Thus, three ideal models are constructed in 3D space that represent a point, line, and plane [18]. The expression is shown in formula (3), and the ideal sphere model and cylinder model are shown in Fig. 5.

\(\left\{\begin{array}{l} \text { dot }: d(x, y, z)=\exp \left\{-\frac{x^{2}+y^{2}+z^{2}}{2 \sigma^{2}}\right\} \\ \text { line: } l(x, y, z)=\exp \left\{-\frac{x^{2}+y^{2}}{2 \sigma^{2}}\right\} \\ \text { plane: } p(x, y, z)=\exp \left\{-\frac{x^{2}}{2 \sigma^{2}}\right\} \end{array}\right.\)       (3)

Fig. 5. Ideal sphere model and cylinder model

Although realization of this ideal model in real images is difficult, we can obtain an approximate model by Gauss filtering of a certain size. The dot model can be understood as asphere that consists of circles of different gray spherical surface. The gray value of the sphere decreases from the center to the surrounding area in the Gauss distribution. For the line model, which represents a cylindrical body, the gray value of each point in the cylindrical cross-section declines from the center to the surrounding area in the Gauss distribution. For the plane model, which represents the thickness of the cross section, the gray values decreases from the center to the periphery in the Gauss distribution. The gray levels of these three models have diminishing characteristics, however, their decreasing directions differ, Thus, their edges are relatively blurred, which is consistent with the edge features of the nodules and vessels.

A Hessian matrix is a square matrix that comprises the second derivative of the multivariate function, which describes the local curvature of the function. A Hessian matrix is a method for extracting the image feature direction using a high- order differential [19]. In the ideal spherical model, the expression of the Hessian matrix H corresponds to every voxel V(x, y, z) (4). By calculation, H is a three order symmetric matrix, and the values of six mixed partial derivatives \(f_{x y}, f_{x z}, f_{y x}, f_{y z}, f_{z x}, f_{z y}\) are zero. The formulas for \(f_{xx}, f_{yy}\) , and \(f_{zz}\) are shown in (5).

\(H=\nabla^{2} f=\left(\begin{array}{lll} f_{x x} & f_{x y} & f_{x z} \\ f_{y x} & f_{y y} & f_{y z} \\ f_{z x} & f_{z y} & f_{z z} \end{array}\right)\)        (4)

\(\left\{\begin{array}{l} f_{x x}=-\frac{1}{\sigma^{2}} \exp \left(-\frac{x^{2}+y^{2}+z^{2}}{2 \sigma^{2}}\right)\left(1-\frac{x^{2}}{\sigma^{2}}\right)=-\frac{f}{\sigma^{2}}\left(1-\frac{x^{2}}{\sigma^{2}}\right) \\ f_{y y}=-\frac{1}{\sigma^{2}} \exp \left(-\frac{x^{2}+y^{2}+z^{2}}{2 \sigma^{2}}\right)\left(1-\frac{y^{2}}{\sigma^{2}}\right)=-\frac{f}{\sigma^{2}}\left(1-\frac{y^{2}}{\sigma^{2}}\right) \\ f_{z z}=-\frac{1}{\sigma^{2}} \exp \left(-\frac{x^{2}+y^{2}+z^{2}}{2 \sigma^{2}}\right)\left(1-\frac{z^{2}}{\sigma^{2}}\right)=-\frac{f}{\sigma^{2}}\left(1-\frac{z^{2}}{\sigma^{2}}\right) \end{array}\right.\)       (5)

In addition, each of the Hessian matrix can be decomposed by formula (6), where λ₁, λ₂, and  λ₃ are eigenvalues of the matrix, and|λ1| >= |λ2| >= |λ3|, the corresponding eigenvector are e₁, e₂ and e₃, respectively.

\( H=\left[\begin{array}{lll}\mathbf{e}_{1} & \mathbf{e}_{2} & \mathbf{e}_{3}\end{array}\right]\left[\begin{array}{ccc}\lambda_{1} & 0 & 0 \\ 0 & \lambda_{2} & 0 \\ 0 & 0 & \lambda_{3}\end{array}\right]\left[\begin{array}{c}\mathbf{e}_{1}^{\mathrm{T}} \\ \mathbf{e}_{2}^{\mathrm{T}} \\ \mathbf{e}_{3}^{\mathrm{T}}\end{array}\right]=\lambda_{1} \mathbf{e}_{1} \mathbf{e}_{1}^{\mathrm{T}}+\lambda_{2} \mathbf{e}_{2} \mathbf{e}_{2}^{\mathrm{T}}+\lambda_{3} \mathbf{e}_{3} \mathbf{e}_{3}^{\mathrm{T}}\)       (6)

Since the Hessian matrix in equation (4) is a diagonal matrix, the eigenvalues of this matrix are diagonal elements. The formula are shown in (7). The eigenvectors that correspond to the eigenvalues are pairwise orthogonal. The direction of the eigenvector corresponds to the main direction of the three-dimensional ellipsoidal axis, and the size of the eigenvalues corresponds to the length of each axis. Both the direction and the length reflect the shape of the object and the size of the object. The eigenvalues of the spherical model and the cylindrical model are shown in Fig. 6.

\(\left\{\begin{array}{l} \lambda_{1}=-\frac{f}{\sigma^{2}}\left(1-\frac{x^{2}}{\sigma^{2}}\right) \\ \lambda_{2}=-\frac{f}{\sigma^{2}}\left(1-\frac{y^{2}}{\sigma^{2}}\right) \\ \lambda_{3}=-\frac{f}{\sigma^{2}}\left(1-\frac{z^{2}}{\sigma^{2}}\right) \end{array}\right.\)       (7)

Fig. 6. Eigenvalues of the spherical model and the cylindrical model

By formula (7) and Fig. 6, we find that the eigenvalues of the interior point within the object depend on the coordinate of the point, that is, the eigenvalues of a point are related to its distance from the point to the center of the object. We note that the distance of a point from the center is nearer, and the absolute values of eigenvalues are larger, when the point is the center of the sphere. The absolute values of the 3 eigenvalue values reach maximum and equal values, which are shown in formula (8). Similarly, the distance of a point from the center line of a cylinder is nearer, and the absolute values of eigenvalues are larger, when the point is on the center line. The absolute values of its 2 eigenvalues reach maximum and equal values, and the third eigenvalue is zero.

\(\lambda_{1}=\lambda_{2}=\lambda_{3}=-\frac{f}{\sigma^{2}}\)        (8)

For the different ideal models, the corresponding relations of their eigenvalues are summarized in Table 1.

Table 1. The relationship of the eigenvalues for different ideal models

2. Enhanced filter for nodules and vessels

According to the eigenvalue characteristics of the Hessian matrix that corresponds to threedifferent ideal models, Sato designed two enhancement functions to enhance the nodule image and vascular image. S_blob [20] and S_line [21] are defined as (9) and (10) respectively.

\(S_{b l o b}=\left\{\begin{array}{l} \left|\lambda_{3}\right|\left(\frac{\lambda_{2}}{\lambda_{3}}\right)^{\gamma}\left(\frac{\lambda_{1}}{\lambda_{2}}\right)^{\gamma} \\ 0, \text { otherwise } \end{array} \quad\left(\left|\lambda_{1}\right| \geq\left|\lambda_{2}\right| \geq\left|\lambda_{3}\right|\right)\right.\)       (9)

\(S_{\text {line}}=\left\{\begin{array}{l} \left|\lambda_{1}\right|\left(\frac{\lambda_{2}}{\lambda_{1}}\right)^{\gamma}\left(1+\frac{\lambda_{3}}{\lambda_{2}}\right)^{\gamma} \quad\left(\left|\lambda_{1}\right| \geq\left|\lambda_{2}\right| \geq\left|\lambda_{3}\right|, \lambda_{1}<0, \lambda_{2}<0, \lambda_{3} \rightarrow 0\right) \\ 0, \text { otherwise } \end{array}\right.\)        (10)

When the value of γis 1, S_blob = |λ₃|，S_line = |λ₂|+|λ₃|.

A certain relationship between the eigenvalues and its coordinate values exists. Determining whether the point belongs to the nodules that only rely on λ₃ is difficult. Based on the Sato’s method, Li designed the products magnitude and likelihood to enhance the nodule and vessel images [18]. The enhancement functions Z_dot and Z_line are defined as formula (11) and formula (12), respectively. 

\(z_{\text {dot}}=\left\{\begin{array}{l} g_{\text {dot}} \times k_{\text {dot}}=\left|\lambda_{3}\right| \times \frac{\left|\lambda_{3}\right|}{\left|\lambda_{1}\right|} \quad\left(\left|\lambda_{1}\right| \geq\left|\lambda_{2}\right| \geq\left|\lambda_{3}\right|, \lambda_{1}<0, \lambda_{2}<0, \lambda_{3}<0\right) \\ 0, \text { otherwise } \end{array}\right.\)       (11)

\(z_{\text {line}}=\left\{\begin{array}{l} g_{\text {line}} \times k_{\text {line}}=\left|\lambda_{2}\right| \times \frac{\left|\lambda_{2}\right|-\left|\lambda_{3}\right|}{\left|\lambda_{1}\right|} \quad\left(\left|\lambda_{1}\right| \geq\left|\lambda_{2}\right| \geq\left|\lambda_{3}\right|, \lambda_{1}<0, \lambda_{2}<0, \lambda_{3} \rightarrow 0\right) \\ 0, \text { otherwise } \end{array}\right.\)       (12)

The experimental results indicate that the enhancement function proposed by Li is effective. However, we note that the eigenvalue λ₂ does not play a role in the enhancement function in the formula (11). When the eigenvalues λ₁ and λ₃ of two points are equal in the image, the output of the enhancement filter is the same even if the eigenvalues λ₂ differ. To improve the filter effect, we consider the relationships among the eigenvalues and improve the enhancement function proposed by Li. The ratio of eigenvalues are the input of the enhancement function, and the new nodule and vessel enhancement functions are presented. We define function (13).

\(\Im(m, n)=\frac{2|m||n|}{m^{2}+n^{2}}=\frac{2\left|\frac{n}{m}\right|}{1+\left|\frac{n}{m}\right|^{2}},|m| \geq|n|\)       (13)

The values of the input (n/m) and output\(\mathfrak{I}(m, n)\) of the function are [0, 1] and monotonically increasing. When the values of m and n are closer, the output of the function is larger, when the values of m and n are equal, the maximum output is 1. To consider the relationships of the eigenvalues, we define the new nodule enhancement function E_ball as formula (14).

\(E_{\text {ball}}=\Im\left(\lambda_{1}, \lambda_{2}\right) \times \Im\left(\lambda_{2}, \lambda_{3}\right),\left(\left|\lambda_{1}\right| \geq\left|\lambda_{2}\right| \geq\left|\lambda_{3}\right|\right)\)       (14)

The final nodule enhancement function E_ball as formula (15) is obtained by simplification.  

\(E_{\text {ball }}=\left\{\begin{array}{l} \frac{4 \lambda_{3} / \lambda_{1}}{1+\left(\frac{\lambda_{2}}{\lambda_{1}}\right)^{2}+\left(\frac{\lambda_{3}}{\lambda_{2}}\right)^{2}+\left(\frac{\lambda_{3}}{\lambda_{1}}\right)^{2}}\left(\lambda_{1}|\geq| \lambda_{2}|\geq| \lambda_{3} |, \lambda_{1}<0, \lambda_{2}<0, \lambda_{3}<0\right) \\ 0, \text { otherwise } \end{array}\right.\)       (15)

Similarly, the new vessel enhancement function is define as formula (16).

\(E_{line}=\exp \left(\lambda_{3}\right) \times \Im\left(\lambda_{1}, \lambda_{2}\right),\left(\left|\lambda_{1}\right| \geq\left|\lambda_{2}\right| \geq\left|\lambda_{3}\right|\right)\)       (16)

The final vessel enhancement function E_line as formula (17) is obtained by simplification.

\(E_{\text {line}}=\left\{\begin{array}{l} \exp \left(\lambda_{3}\right) \times \frac{2 \lambda_{2} / \lambda_{1}}{1+\left(\frac{\lambda_{2}}{\lambda_{1}}\right)^{2}}\left(\left|\lambda_{1}\right| \geq\left|\lambda_{2}\right| \geq\left|\lambda_{3}\right|, \lambda_{1}<0, \lambda_{2}<0, \lambda_{3} \rightarrow 0\right) \\ 0, \text { otherwise } \end{array}\right.\)       (17)

We use λ₂/λ₁ and λ₃/λ₂ as the coordinate input to draw the output curve of the nodule enhancement function E_ball , as shown in Fig. 7. The curve indicates that λ₂/λ₁ and λ₃/λ₂ have the same contribution value to the value of E_ball, and the range of Eball is [0, 1] and monotonically increasing. Therefore, the value of E_ball of a voxel can express the possibility that the voxel belongs to the nodule. The greater the value of E_ball is, the greater the possibility is that the voxel belongs to the nodule. Conversely, the smaller E_ball is, the smaller the possibility of the nodule is. Thus, the nodule image can be enhanced by the enhancement function E_ball.

Fig. 7. Output curve of nodule enhancement function Eball

Similarly, we use λ₂/λ₁ and λ₃ as the coordinate input to construct the output curve of the vessel enhancement function E_line, as shown in Fig. 8. As shown in the curve, the E_line value monotonically increases at [0, 1]. Therefore, the E_line value of a voxel can express the possibility that the voxel belongs to the blood vessel. the greater the E_line value is, the greater the possibility is that the voxel belongs to the blood vessel; Conversely, the smaller the E_line value is, the smaller the possibility of the blood vessel is. Therefore, the image of the blood vessel can be enhanced by the enhancement function E_line. We note that the weights of the two inputs in the vascular enhancement function Eline, are not equal in their response, and hope to improve this result in future studies.

Fig. 8. Output curve of blood vessel enhancement function Eline

3. Multi-scale computation of Gauss function

In lung CT images, the size of nodules is uncertain, a substantial amount of image noise exists, and the computational process of the second-partial-derivative has a strong sensitivity to image noise. If an enhancement filter is directly applied to an image, poor results are likely. To effectively detect nodules of different sizes, a multi-scale filtering method based on a Gaussfunction is adopted in this paper. The Gauss function is used to simultaneously convolve the image, remove noise from the image and smooth the image, which enhances the image on multi-scale[22-29].

For a Gaussian function with a standard deviation of σ, a minimum 95% area of the Gauss function is included in the range [-4σ, 4σ]. When designing the scale σ of the enhancement function, the σ of the function should be d/4 when the nodule diameter is d. Before the nodule is enhanced, the image should be convolved by the Gaussian function with the standard deviation σ, which not only eliminates the image noise but also improves the enhancement effect of the filter.

Assume that the diameter of the object to be enhanced is [d₀, d₁], to enhance all objects in this range, the range of the Gauss filter can be set to [d₀/4, d₁/4]. Select N different σ values, and apply convolution operations on images to smooth. The larger N is, the higher the detection accuracy of nodules, and the complexity of the calculation are. The smaller N is, the lower the detection accuracy of nodules is and the simpler the calculation is. Therefore, the selection of N has a significant influence on the accuracy of the nodule detection results. In this paper, the diameter of the candidate nodules is [3mm, 30mm], and the number of scales N is five. When N is determined, the size of each scale is calculated as formula (18).

\(\left\{\begin{array}{l} \sigma_{1}=\frac{d_{0}}{4}, \sigma_{2}=r \sigma_{1}, \ldots, \sigma_{N}=r^{N-1} \sigma_{1}=\frac{d_{1}}{4} \\ r=\left(\frac{d_{1}}{d_{0}}\right)^{\frac{1}{N-1}} \end{array}\right.\)       (18)

The specific steps of using a multi-scale filter to detect suspected nodules are listed as follows:

Step 1: Determine the range [d₀, d₁] of the object to be enhanced and the scale number N of the filter. Use formula (12) to calculate the size of each scale σ_N.

Step 2: For each scale σ_N, execute Steps 3-6.

Step 3: Smooth the image with Gauss filtering, and execute Steps 4-6 for each voxel.

Step 4: Construct a Hessian matrix, compute the three eigenvalues λ₁, λ₂, and λ₃, and use formulas (15) and (17) to calculate the nodule enhancement function E_ball and the vessel enhancement function E_line , respectively, of each voxel respectively.

Step 5: Nodule and vessel enhancement filter are used to enhance the voxels.

Step 6: End the cycle for voxels and scales.

Step 7: For each voxel, the final output is the maximum value computed at different scales, and the nodule enhancement image and the vessel enhanced image are obtained.

Step 8: 3D Lung images of suspected pulmonary nodules are obtained by removing the overlapped part between the nodule image and the vessel image from the nodule image.

Step 9: Output the final nodule detection result.

2.2.3 Feature extraction and pulmonary-nodule classification

In Section 2.2.2, nodules and a large number of vascular crossings have been enhanced. We observe that removing the overlap part from the nodule image can eliminate the blood vessels. However, some false positive nodules, such as vascular intersection points will remain. Many studies have demonstrated that the false positive nodules of a pulmonary dection system primarily include vessels, vascular cross and vascular bending, which is known as vascular-type false-positive nodules [30]. To improve the accuracy of pulmonary nodules detection, we need to examine feature extraction of suspected pulmonary nodules.

In a CAD system, Feature extraction is the core problem of lung nodule detection and diagnosis, which determines the performance of subsequent classification [31]. An acceptable feature requires a certain degree of image translation, rotation, and scaling invariance. Therefore, characteristics such as area, volume, diameter, roundness, sphericity, compactness and global shape have been extensively applied to describe lung nodules and remove false positive nodules. In addition, feature factors based on a histogram have rotation, translation and scale invariance, which are important feature extraction techniques in image classification and computer vision. A histogram of oriented gradients (HOG) and scale invariant feature transformation (SIFT) are commonly used.

To better describe the characteristics of pulmonary nodules, we propose the feature descriptors Surface normal orientation angles histogram (SNOAH) for lung nodules. First, use an 18- neighborhood 3D region growing method to extract all candidate nodules. Second, use a Hessian matrix to extract the surface normal vector histogram of each voxel for the final nodal classification. This method does not depend on the segmentation results of lung nodules and has a certain robustness.

1. Nodules extraction

In Section 2.2.2, the final output of each voxel is the maximum value of the enhancement function for all scales. The larger the value of the nodular enhancement function is, the greater the probability is that it will be attached to pulmonary nodules. Therefore, the coordinate of candidate nodules can be obtained by analyzing the value of the nodule enhancement function for voxel points. In addition, we use a 5-dimension feature vector [E, g, x, y, z]^T for each voxel, where E represents the value of the nodule enhanced function E_ball, g is the gray values, x, y are the point coordinates, and z is the CT serial number. The Euclidean distance between the eigenvectors of two voxels is used to measure the similarity among pixels.

Region growing [32] is a classic image segmentation method. The algorithm can usually segment the connected region with the same features, and provide excellent boundary information and segmentation results. In this paper, we propose an18-neighborhood 3D region growing method for the extraction of each nodule image as shown in Fig. 9. First, the threshold T_E is set. When the nodal enhancement function value E_ball of the voxel is larger than the threshold T_E, these voxels are selected as the nodule initial seed. Second, for all seed points, calculate the Euclidean distance D of the feature vectors between the 18 neighborhood voxels and the seed point. If D is less than the threshold T_d, the voxel point will be expanded into the seed region as the seed point to a new cycle until the size of the seed region does not change. The specific steps to extract pulmonary nodules are listed as follows.

Step 1: Set the threshold T_E of the nodule enhancement function E_ball, and obtain the voxel seed point set {Seed₁, Seed_2, …, Seed_i}.

Step 2: Extend the 18-neighborhood region of each seed point.

Step 3: Calculate the Euclidean distance D between each voxel point in the 18-neighborhood domain and the existing seed point.

Step 4: If D < T_d, the voxel will be included in the seed region and used as the new seed point.

Step 5: Repeat Steps 2-4 until the size of the seed region does not changed.

Step 6: Output detected 3D images of all suspected nodules.

Fig. 9. 18-neighborhood 3D region growing method. The black point is the seed point and the white point is the 18-neighborhood point

2. Surface normal orientation angles histogram

The lung nodule feature descriptor Surface normal orientation angles histogram (SNOAH) can reflect the probability distribution of the normal vector direction angle of a pulmonary nodule surface. Fig. 10 shows the normal vector direction of each voxel surface of the nodules and vessels. The normal direction of the nodule is uniformly radial, without the main direction, whereas the normal direction of the blood vessel is consistent with the direction of blood flow in the blood vessel. Therefore, the nodule and vessel can be distinguished by the normal orientation distribution of voxels. In this paper, the SNOAH feature extraction procedure is described as follows. Input candidate voxel coordinates for solving eigenvalues and eigenvectors by Hessian matrix decomposition. Calculate the vector direction angle for each voxel, and perform a statistical analysis to obtain the voxel SNOAH.

Fig. 10. Normal vector direction of the nodule (a) and the vessel (b)

Sectio 2.2.2 we obtain formula (6) by decomposing the Hessian matrix and obtain formula (19) from formula (6).

\(H=\left(\lambda_{1}-\lambda_{2}\right) \mathbf{e}_{1} \mathbf{e}_{1}^{\mathrm{T}}+\left(\lambda_{2}-\lambda_{3}\right)\left(\mathbf{e}_{1} \mathbf{e}_{1}^{\mathrm{T}}+\mathbf{e}_{2} \mathbf{e}_{2}^{\mathrm{T}}\right)+\lambda_{3}\left(\mathbf{e}_{1} \mathbf{e}_{1}^{\mathrm{T}}+\mathbf{e}_{2} \mathbf{e}_{2}^{\mathrm{T}}+\mathbf{e}_{3} \mathbf{e}_{3}^{\mathrm{T}}\right)\)         (19)

where \(e_1e_1^T, \ e_1e_1^T+e_2e_2^T, \ e_1e_1^T+e_2e_2^T+e_3e_3^T\) represent the plane tensor, cylindrical tensorand spherical tensor, respectively. The coefficients \(\lambda_{1}-\lambda_{2}, \lambda_{2}-\lambda_{3}\) and \(\lambda_{3}\) represent the corresponding weights. For these three types, three types of features expression are definded as local surface (surface-ness), curve (curve-ness) and local point or ball (point-ness).

surface-ness：direction is the same with e1, weight is λ₁ - λ₂;

curve-ness：direction is the same with e3, weight is λ₂ - λ₃;

point-ness：no direction, weight is λ₃.

From the definition, we determine that the direction of surface-ness is the same for e₁, which corresponds to λ₁. Therefore, we can determine the direction of surface-ness by e₁. In the space coordinate system, any point in space can be represented as (r, θ, φ), where r denotes the radius, θ denotes the elevation angle and φ is the azimuth angle, as shown in Fig. 11. θ and φ can be calculated according to formula (20).

\(\left\{\begin{array}{l} \theta=\arccos \left(e_{1}^{(z)}\right) \\ \varphi=\arctan \left(\frac{e_{1}^{(y)}}{e_{1}^{(x)}}\right) \end{array}\right.\)       (20)

where \(e_{1}^{(x)}, e_{1}^{(y)}, \text { and } e_{1}^{(z)}\) represent the component in the x, y and z directions, respectively, of the surface normal vector e₁. The range of the elevation angle θ is [0, 180], and the range of the azimuth angle φ is [0, 360]. A two-dimensional statistical histogram is obtained by the statistics of the elevation angle θ and the azimuth φ. In the statistical surface normal vectorelevation angle θ and the azimuth angle φ, θ and φ are equally divided into m parts. Each part of θ and φ covers an angle range of (180/m) degrees and (360/m) degrees, respectively, that is, the dimension of the final extracted feature vector is 2m.

Fig. 11. Diagram of elevation angle θ and azimuth angle φ

3. Support Vector Machine (SVM)-Based Classification

A SVM is one of the best single classifiers for the diagnosis of benign and malignant pulmonary nodules, it has a significant effect on solving small samples, and high-dimensional and nonlinear data. In this paper, a sample set of N suspected pulmonary nodules is constructed as \(\boldsymbol{X}=\left\{\left(\boldsymbol{x}_{i}, \boldsymbol{y}_{i}\right)\right\}, \boldsymbol{x}_{i} \in R^{n}, \quad i=1,2, \ldots N, \boldsymbol{y}_{i} \in\{1,-1\}\). The purpose of the SVM training sample is to minimize the error function, that is, to solve the optimization problem of formula (21).

\( \text {Minimize}\left\{\frac{1}{2} w^{T} w+\mathrm{C} \sum_\limits{i=1}^{N} \xi_{i}\right\} \)     

\(Subject\ to \\ \left\{\begin{array}{c} f\left(x_{i}\right)=\sum_\limits{i=1}^{N} w^{T} \kappa\left(x_{i}\right)+b \geq 1-\xi_{i} \\ \xi_{i} \geq 0 \end{array} \quad \forall i \in[1, \mathrm{N}]\right. \)        (21)

where f(x_i) is a decision function, and κ is a function mapping, that can map x_i into a high-dimensional feature space. The purpose of a SVM is to obtain a hyper-plane (w, b) to satisfy the classification requirements, and extend the sample point distance from the plane as far as possible. C is a penalty factor that is greater than zero.

Introducing the kernel mapping SVM can solve the linearly non-separable problem. In this manner, the linearly non-separable problem in the original sample space becomes linearly separable in a high-dimensional sample space. The kernel function is defined as K(x_i , x)≡ κ(x_i)^Tκ(x). In this paper, we use the polynomial kernel K_p(x_i , x), the RBF kernel K_r(x_i , x), and the mixed kernel K_mix(x_i , x) as kernel functions, their formulas are shown as formula (22) , formula (23) and formula (24), respectively, where β is a mixed weighting factor. In addition, we use k-folding cross-validation to train and test the classifier.

\(K_{p}\left(x_{i}, x\right)=\left\|x_{i} \cdot x+1\right\|^{p}\)       (22)

\(K_{r}\left(x_{i}, x\right)=\exp \left(\frac{-\left\|x_{i}-x\right\|^{2}}{2 \sigma^{2}}\right)\)       (23)

\(K_{m i x}\left(x_{i}, x\right)=\beta K_{p}\left(x_{i}, x\right)+(1-\beta) K_{r}\left(x_{i}, x\right)\)       (24)

3. Results

The experimental platform of this algorithm is Visual Studio 2010 and MATLAB 2012b. The PC processor is Intel Core i7-3770, the main frequency is 3.40 GHz and the memory is 8 GB. We establish some necessary parameters to ensure the accuracy and effectiveness of our method. The values of these parameters are listed in Table 2.

Table 2. Parameter values settings

3.1 Lung segmentation

Using the adaptive iterative threshold and the region growing method, we segment all lung images sequences. Fig. 12 shows the processes of lung segmentation. Column (a) contain 5 original lung images from the top to the bottom of the lung. The initial threshold T0 is set to 500 Hu, and the optimal threshold T_opt is obtained by the iterative threshold method for lung image binaryzation. The binary images are shown in column (b). ROI images are extracted as shown in column (c). Lung parenchyma images are coarsely segmented using the region growing methods shown in column (d). We conduct trachea removal the results are listed in column (e). Morphological processing is used to generate lung parenchyma masks to separate lung parenchyma results, as shown in column (f) and column (g).

Fig. 12. Lung parenchyma segmentation in of CT image sequences  

3.2 Suspected nodule detection

The detection method proposed in our paper is used to detect suspected pulmonary nodules onlung images. First, two multi-scale enhancement filters are constructed to enhance similar spheres pulmonary nodule images and cylindrical vascular images. In this paper, the diameter of the nodule image to be enhanced ranges between 3 mm and 30 mm. We use 5 different filter scales, the size of each scale is 0.75mm, 1.33mm, 2.37mm, 4.21mm and 7.5mm, respectively. The detection results of the suspected pulmonary nodules for 2 cases are shown in Fig. 13 and Fig. 14.

Fig. 13. Comparison of image enhancement result for Data1 

Fig. 14. Comparison of image enhancement result for Data2

3.3 Feature extraction

The SNOAH features extracted from the nodules, vessels and pleura are shown in Fig. 15, Fig. 16 and Fig. 17, respectively. Figure (b) and (c) display the statistics histograms of elevation angle and azimuth angle. As shown in figure (c), the azimuth angles of nodules do not have a distinct main direction, they are evenly distributed from 0-360 degrees, they have a distinct main direction and they are centered at 180 and 360 degrees. Conversely, the azimuth angles of pleural are only distributed at several special angles. The feature descriptor SNOAH proposed in this paper can distinguish the shape of objects. We have performed angle histogram statistics for several types of special images. The result of the pleural traction image is shown in Fig. 18 and the vascular intersection image is shown in Fig. 19.

Fig. 15. SNOAH features of nodules (a)Nodule. (b) Elevation. (c) Azimuth 

Fig. 16. SNOAH features of vessel (a)Vessel. (b) Elevation. (c) Azimuth

Fig. 17. SNOAH features of pleural (a)Pleural. (b) Elevation. (c) Azimuth 

Fig. 18. SNOAH features of pleural traction (a) Pleural traction. (b) Elevation. (c) Azimuth 

Fig. 19. SNOAH features of vascular intersection (a) Vascular intersection. (b) Elevation. (c) Azimuth 

3.4 Classification with SVM

In our experiences, the polynomial kernel SVM_P, the radial basis kernel SVM_R and the hybrid kernel SVM_M are used to classify the features of pulmonary nodules. Each classifier was trained and tested using K-fold cross-verification. The experimental results werecompared by sensitivity (SE), specificity (SPE), accuracy (ACC) and receiver operating characteristic (ROC) curves. The formulations of SE, SPE and ACC are shown as (25), (26) and (27), respectively.

\(S E=\frac{T P}{T P+F N}\)       (25)

\(S P E=1-\frac{F P}{F P+T N} \)       (26)

\(A C C=\frac{T P+T N}{T P+F P+T N+F N}\)       (27)

where TP, FN, FP and TN represent true positive, false negative, false positive and truenegative, respectively. In addition, the experimental results of three classifiers with different k values are shown in Table 3, Table 4 and Table 5, and the corresponding ROC is shown in Fig. 20, Fig. 21 and Fig. 22, respectively.

Table 3. k-foled cross-validation results for SVM_P-based nodule detection

Table 4. k-foled cross-validation results for SVM_R-based nodule detection

Table 5. k-foled cross-validation results for SVM_M-based nodule detection

Fig. 20. ROC curve of SVM_P

Fig. 21. ROC curve of SVM_R

Fig. 22. ROC curve of SVM_M

The best result for each classifier is shown in Table 6, where the value of k is 6. Table 6 indicate that the SVM_M obtained the best classification performance with a sensitivity of 96.88% and a false positive of 0.93, which outperformed the SVM_P and SVM_R. The classification sensitivities of SVM_P and SVM_R were also 91.49% and 93.61%, which verified the effectiveness of our features.

Table 6. Best performance of each classifier, where the value of k is 6

Table 7 shows the comparison results between the existing CAD method and our method. The sensitivity of our method is 96.88%, which is lower than the sensitivity of 97.66% in one study [19] and the sensitivity of 100% in another study [42]. However, the false positive is 0.93, which is lower than the sensitivity of other methods. Thus, our method can effectively reduce false positive results and maintain high detection sensitivity and accuracy.

Table 7. Performance comparison of reported methods

4. Conclusion

Our work indicates that the proposed method is more accurate with lower false positive rates than the rates of traditional methods. In our method, an adaptive threshold iteration method and region growing algorithm are used for segmentation of lung parenchyma image sequences. Nodular and vascular enhancement filters are used to enhance the nodular images and vascular images, respectively, in the 3D lung. The blood vessel images in the enhanced pulmonary nodule images are removed, and the suspected nodular images are obtained, which significantly reduces the vascular type false-positive nodules. The SNOAH features of suspected pulmonary nodule images are extracted and classified by a SVM classifier with different kernel functions, which can remove the false positive nodules from suspected nodules by extracting effective features. We compare the results of our method with the results of existing CAD methods for detecting lung nodules. The experimental results indicate that the proposed method can achieve 0.9688 sensitivity and 0.9467 specificity. The accuracy of pulmonary nodule detection is 0.9511, and the average number of false positive nodules in each case is 0.93, which is significantly better than other detection methods. The results demonstrates that the proposed lung nodule feature descriptor SNOAH is effective for distinguishing nodules and blood vessels.