This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/JBHI.2014.2332955, IEEE Journal of Biomedical and Health Informatics
> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT)
REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) < Following this approach, Kubo et al. used the VanderBurg linear feature detector and morphological operators to find fissures in 25 lungs, 12 of which contained pulmonary nodules [10]. Zhang et al. combined a pulmonary atlas, a ridgeness operator and fuzzy reasoning to identify fissures in 7 normal and 5 diseased lungs [11]. Using a similar ridgeness operator to enhance the fissures, Wang et al. employed Bayesian theory to identify fissures in 10 lungs with pulmonary nodules [12]. Their algorithm yielded high segmentation accuracy with manually delineated fissures for initialization. To eliminate manual intervention, we developed an algorithm using wavelet transform to locate fissures in 9 lungs with pulmonary nodules [13]. In general, the 2D fissure approach works well for identifying disrupted and deformed fissures in individual CT images; however, it does not always produce complete segmentation of lung lobes, due to the lack of continuity between fissures in adjacent CT images. 2) 3D Surface Approach: Another approach of segmenting lung lobes is by finding fissure surfaces using 3D anatomic structures within the lungs. With a Kuhnigk et al. used a segmented vascular tree and an interactive watershed transform (IWT) to find fissure surfaces in 24 lungs [14]. Ukil et al. extended this algorithm by automating the IWT using segmented vascular and bronchial trees [15]. They then used a graph search to refine fissure surfaces in 12 normal and 17 emphysema lungs. Their algorithm worked well for lungs with emphysema; however, one quarter (25%) of their lungs required manual correction [15]. Following the same approach, Zhou et al. used Voronoi division to locate fissure surfaces in 7 normal lungs [16]. Similarly, Saita et al. employed a sheet-emphasis filter to extract fissure surfaces in 16 normal lungs and 4 lungs with unknown diseases [17]. Using a slightly different strategy, Van Rikxxort et al. first enhanced fissures in individual CT images using a set of supervised filters [18], and then utilized a segmented bronchial tree to find fissure surfaces in 20 unspecified normal/abnormal lungs [19]. Lastly, Pu et al. used a computational geometry method to first locate fissure patches [20] and then used implicit radial basis functions to extend these patches into fissure surfaces in 65 lungs with chronic obstructive pulmonary disease [21]. While the 3D surface approach yields complete segmentation of the lung lobes, it often requires difficult segmentation of anatomic structures such as the bronchial and vascular trees [22-23]. More importantly, current 3D surface approaches have difficulties in handling disrupted and deformed fissures, which are often caused by cancer lesions. To the best of our knowledge, no existing methods have been tested under all three constraints of: (1) automation, (2) dealing with multiple diseases, and (3) producing a complete segmentation in the 3D space, as listed in Table I. To meet these constraints, we developed a new algorithm employing a hybrid 2D/3D approach by first finding fissures in 2D CT images using a texture analysis method, followed by a surface fitting technique based on ridge regression to model actual fissure surfaces in 3D. Compared to our previous algorithms in [13] and [24], the novelties of this algorithm include the following two aspects: fully automatic segmentation of 3D major (oblique) and minor (horizontal) fissures using ridge regression
2
TABLE I STATE-OF-THE-ART ALGORITHMS FOR SEGMENTING ALL LUNG LOBES Evaluated on Segmentation Authors Automation Multiple in 3D Diseases 2D Fissure Approach Kubo et al. [10] Zhang et al. [11] Wang et al. [12] Wei et al. [13]
YES NO NO NO
NO YES NO NO
NO NO NO NO
NO NO NO NO NO NO
YES YES YES YES YES YES
3D Surface Approach Kuhnigk et al. [14] Ukil et al. [15] Zhou et al. [16] Saita et al. [17] Van Rikxoort et al. [18, 19] Pu et al. [20, 21]
NO NO YES YES YES YES
in diseased lungs which often have disrupted, deformed and/or incomplete fissures; and an increase of segmentation accuracy in about 10 percent. We evaluated our algorithm using lungs with a variety of different diseases: emphysema, small cell lung cancer, large cell lung cancer and bronchiectasis. To present our algorithm, we organized this paper as follows: Section II describes our methodology. Section III and IV provides the results and discussion, respectively. Lastly, Section V gives the conclusion and future work. II. METHODOLOGY Our algorithm (developed in Matlab) uses a pipelined process to segment diseased lung lobes. The algorithm employs four stages: A) preprocessing, B) 2D oblique fissure recognition, C) 2D horizontal fissure recognition and D) 3D lung lobe segmentation. Fig. 1 depicts the flow diagram of this algorithm. The following context presents these four stages and evaluation method in five subsections, respectively. A. Preprocessing In this stage, we preprocessed a stack of CT images using the Wiener filter and lung segmentation methods that we developed previously [13]. The Wiener filter reduces the noise inherited in isotropic CT images, while the lung segmentation method decreases computation time by removing area outside of the lungs. Fig. 2(a) and (b) illustrate example cross-sectional and sagittal images (see subsection C) from a preprocessed CT stack, respectively. B. 2D Oblique Fissure Recognition To identify the oblique fissures in a stack of cross-sectional CT images, we used a method based on our preliminary algorithm [24]. In short, this method uses three steps to recognize the oblique fissures in pathological lungs: (1) texture analysis, (2) fissure region analysis, and (3) fissure identification. In the first step, we performed 2D texture analysis (instead of 3D to reduce computation) by dividing all images in the preprocessed CT stack into overlapping
2168-2194 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/JBHI.2014.2332955, IEEE Journal of Biomedical and Health Informatics
> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <
X
CT Image Stacks
3
Y
A) Preprocessing
Wiener Filter
Lung Segmentation
Y B) Oblique Fissure Recognition
Texture Analysis using Neural Network
C) Horizontal Fissure Recognition
Fissure Identification
(a)
Stack Rotation Left Lung Removal
Oblique Fissure Region Analysis
Z
Lower Lobe Removal Horizontal Fissure Region Analysis
Z
(b) X
Fig 2. Examples of (a): cross-sectional images from a preprocessed CT stack and (b): sagittal images from the same CT stack.
(a)
(b)
(c)
(d)
Fissure Identification
D) Lung Lobe Segmentation
Fissure Verification
Surface Modeling
Fig. 1. Flowchart of the lobe segmentation algorithm. 100 Y (Pixels)
region-of-interests (ROIs) with dimensions of 10 × 10 pixels. The selection of 10 pixels corresponds to the width of typical fissure regions. We used gray level run length texture features [25] and a feed-forward back-propagation neural network (NN) [26] to classify these ROIs as either a fissure or non-fissure ROI. For NN training, we employed a set of predefined ROIs in the lungs. As a result, the first step of texture analysis yields a stack of binary images containing classified fissure and non-fissure ROIs. Fig. 3(a) and (b) illustrates an original CT image and the matching binary image yielded by the first step, respectively. To find regions containing the fissures (fissure regions) in the second step of fissure region analysis, we performed a stack projection followed by a line fitting. The stack projection transforms the binary stack into two accumulation images (a left and a right accumulation image) showing the propagation of fissure regions in the y-z plane (see Fig. 2). The two accumulation images correspond to the left and right lungs. For each accumulation image, the line fitting calculates a straight line (propagation line) to extract the locations of the fissure regions for each image in the binary stack. The propagation line calculated from the right accumulation image is also used for the horizontal fissure recognition (see subsection C). Fig 3(c) and (d) presents an example right accumulation image and the corresponding propagation line calculated from the line fitting, respectively. Fig. 3(e) depicts fissure regions (white pixels) found within the binary image in Fig. 3(b) using the calculated propagation line. In the last step of fissure identification, we combined the found fissure regions in the binary stack and the images in the preprocessed CT stack to identify the actual fissures. Fig. 3(f) gives example fissures found within the fissure regions shown
80 60 40 20 0
(e)
0
20 40 60 80 100 X (Pixels)
(f)
Fig. 3. Example images from oblique fissure recognition: (a) original CT image; (b) binary image containing classified fissure (white pixels) or non-fissure ROIs (black pixels); (c) right accumulation image; (d) line fitting of pixels on the accumulation image shown in (c); (e) found fissure regions (white pixels); (f) partial fissures found within fissure regions shown in (e).
in Fig. 3(e). As illustrated in Fig. 3(f), fissures found by the algorithm do not extend to edges of the lungs due to diseases and incomplete fissures (addressed in subsection D). More detailed description of the above method can be found in [24]. C. 2D Horizontal Fissure Recognition Our previous work [24] did not identify the horizontal fissure that is more challenging to segment than the oblique fissures. Compared to the oblique fissures, the horizontal
2168-2194 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/JBHI.2014.2332955, IEEE Journal of Biomedical and Health Informatics
> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT)
REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <
(a)
denotes the ith cross-sectional CT image that the fissure pixel, pi, resides in. Using this set of n validated voxels, {vi(xi, yi, zi), i=1,2…n}, we modeled a fissure surface in the form of:
(b)
z f ( x, y ). (c)
Fig. 5. Horizontal fissure region analysis: (a) binary image showing the largest region (white pixels) and the selected neighboring regions (gray pixels), as the horizontal fissure region; (b) horizontal fissure region found using bounding box of the selected regions in (a); and (c) horizontal fissure found within the fissure region using dynamic programming.
5) Fissure Identification: We used the bounding box of all selected regions in the binary images to define the fissure region in the sagittal CT images. Fig. 5(b) shows an example fissure region found from Fig. 5(a). Within the fissure region, we performed dynamic programming using Eq. (2) to find the optimal path extending from the left to the right of the fissure region. Fig. 5(c) illustrates a fissure found from the fissure region in Fig. 5(b) using dynamic programming. D. 3D Lung Lobe Segmentation The oblique and horizontal fissure recognition uses 2D texture analysis and fissure anatomy to find the fissures. Although this approach yields good approximations of the fissures in individual 2D CT images, it does not produce complete segmentation 3D. This is due to the fact that many fissures appear incomplete in 2D CT images. As well, lung diseases may disrupt and deform existing fissures causing discontinuities between fissures in adjacent CT images. Fissures appear naturally as 3D surfaces separating adjacent lung lobes. Therefore, the most logical method for producing a realistic segmentation of the lung lobes is to model a 3D surface (fissure surface) for each of the LOF, ROF, and RHF, using fissures found in individual 2D CT images. This hybrid 2D/3D approach also provides the benefit of being able to fill in incomplete and disrupted fissures automatically. To ensure only the correct fissures are used to model the fissure surface, we first performed fissure verification by comparing fissures in adjacent CT images. For each fissure on a cross-sectional CT image, we computed its accuracy as the percentage of its pixels satisfying the following condition: ( y1 y 2 ) 2 ( x1 x2 ) 2 2 pixels ,
5
(3)
where (x1, y1) and (x2, y2) are the coordinates of a fissure pixel and its closest counterpart in the previous adjacent CT image, respectively. The 2 pixel criterion in Eq. (3) ensures that fissures in neighboring CT images are connected. We chose all fissures with an accuracy of >80% (based on experimentation), in 5 or more consecutive CT images, as the correct fissures. We selected a range of 5 or more consecutive images to prevent false verification due to a single incorrectly found fissure. Once the fissures are verified, we represented each fissure pixel pi(xi, yi), as a voxel (a point in 3D), vi(xi, yi, zi), where zi
(4)
Many 3D modeling methods exist [29, 30]. One method is the minimization of energy functions, as used by Pu et al. [21]. While this method produced smooth approximations of fissure surfaces, it is computationally intensive (Pu et al. had to down sample from 200,000 points to 3,000 points during their surface modeling). For our surface modeling, we used a method based on ridge regression [31]. Ridge regression introduces regularization to the standard least squares (LS) formulation. This method carries several advantages over common fissure modeling techniques: (1) provides control over the smoothness of the fissure in the 3D plane; (2) reduces the sensitivity to erroneously found fissures; (3) is computationally and memory efficient. To derive ridge regression for our purpose of fissure modeling, we start by using linear regression to model the relationship between the z and x-y coordinates of each fissure voxel with the equation: zˆi w0 w1 xi w2 yi e ,
(5)
where zˆi is the estimated value of zi, while e represents the modeling error. Similarly, for n voxels, we can model the same relationship using a set of n equations, expressed in a matrix notation: z Aw e ,
(6)
which has the expanded form of: z1 1 x1 z 1 x 2 2 z n 1 xn
y1 e1 w0 y2 e2 w1 . w2 yn en
(7)
Thus, Eq. (6) represents a linear regression model, where w denotes the parameters of our fissure surface. The simplest method for solving this linear regression model would be the least squares (LS) method. The LS method estimates w, by minimizing the sum of squares of the errors:
w min min Aw z
2
.
(8)
Despite the simplicity of the LS method, it carries two drawbacks: (1) local inconsistency with the continuity of input data; (2) oversensitivity of small errors in the input causing a few outliers to skew the estimation or prediction model. Ridge regression alleviates these two drawbacks of the LS method by using regularization. Given a linear regression model shown in Eq. (6), ridge regression estimates the parameters, w, by performing the following minimization:
w min min Aw z
2
w
2
,
(9)
where Γ is the regularization term. The implementation of ridge regression is as follows: (1) transform the LS formulation
2168-2194 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/JBHI.2014.2332955, IEEE Journal of Biomedical and Health Informatics
> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) < to obtain uncorrelated components, (2) apply regularization to the model to achieve a desired property and (3) restore the original coordinate system by reversing the original transform. In the first step of ridge regression, we transformed the (x, y) coordinates of the voxels into a linear index aij by using the equation: 1， if j yi xi Dx aij , otherwise 0， i [1, n], j [1, m], m Dx D y ,
(10)
where Dx and Dy are the pixel dimensions of a cross-sectional CT image in the x- and y-axes, respectively. Thus, the ridge regression formulation becomes: z1 a11 z a 2 21 z n an1
a12 a1m w1 e1 a 22 a2 m w2 e2 . a n 2 anm wm en
(11)
In the second step of ridge regression, we applied regularization. Regularization provides several advantages for fissure surface modeling. One, it allows the modeling of a non-linear fissure surface using a simple linear model. Two, it allows control over the smoothness of the surface. Three, it reduces sensitivity to outliers in the data, i.e. erroneously found fissures. In order to produce a smooth fissure surface, we used Laplacian functions (along the two independent axes, x and y) as our regularization term. That is, 2 f ( x, y ) 2 2 x , f ( x, y ) 2 y
(12)
where α is a scale factor that controls the effect of regularization. The minimization of the Laplacian functions ensures minimal changes in the gradient of our fissure surface. This allows us to control the smoothness of the fissure surface by varying the scale factor, α. We selected an optimal α value of 1 for our evaluation (see Section III, regarding the effect of α on the accuracy of segmentation). In the last step of ridge regression, we computed our fissure surface in the original x-y coordinates, z(xi, yi), by reversing the transformation of Eq. (10). We reversed the transformation by using the equation: z ( xi , y j ) w min (i jDx ),
(13)
where wmin is the result of Eq. (9) while i, j, and Dx are as defined in Eq. (10). E. Evaluation Method For evaluation, we used isotropic CT image stacks from 24 anonymous patients (following Canadian tri-council ethics guidelines) at the Foothills Medical Center, Calgary, AB, Canada. A Siemens Sensation 16 multi-slice CT scanner acquired these image stacks using the standard Foothills
Patient 1–6 7–9 10 – 13 14 – 18 19 – 20 21 22 23 24
6
TABLE II PATIENT PATHOLOGY Types of Lung Disease Mild Emphysema Moderate Emphysema Large Cell Lung Cancer (LCLC) Small Cell Lung Cancer (SCLC) Bronchiectasis Elevated Left Hemidiaphragm Hilar Mass Neoplastic No Diagnosis
Medical Center Thoracic Abdominal 3D Research protocol and iodine contrast agent. All images have a resolution of 512 × 512 pixels with a thickness of 0.6 mm. Of the 24 patients, 23 had lung diseases and one had no diagnosis. Table II lists the types of diseases associated with these patients. To evaluate our algorithm, we employed a two-fold cross validation method [26] by splitting the datasets from the 24 patients into two groups. For each group, we randomly selected half of the datasets containing each disease. This ensures a good representation of all diseases for each group. For diseases where only one dataset is present (Patients 21-24), we randomly placed the dataset into one of the two groups. Of the two groups, we used one group for training and the remaining group for testing. We swapped the two groups to test all 24 datasets. Finally, we repeated the above procedure 10 times and averaged the result. For performance evaluation, we compared the output of our algorithm to that of manual segmentation. While our algorithm yields complete segmentation of the lung lobes, it is extremely time consuming for a radiologist to completely trace all fissures in 24 patients with an average of 300 CT images per patient. To make the evaluation process more feasible, we used manually placed anchor points to mimic manual segmentation. Under the guidance of a radiologist, we manually placed anchor points (at about 5 pixels apart) on fissures for every 5th image of a CT stack, which is approximately the thickness of the CT images that surgeons/radiologists read in clinical settings. Using this manual segmentation, we performed the following analysis: 1) One-way Multivariate Analysis of Variance (MANOVA): One-way MANOVA examines statistically a significant difference between two groups of measurements [32]. In our case, the two groups of measurements are the fissure pixel coordinates, (xi, yi, zi), produced by our algorithm and manual segmentation. By computing a corresponding pair of Wilk’s lambda (λ) and p-value, one-way MANOVA tests the null hypothesis that the two groups of pixel coordinates are different. If the p-value is greater than 0.05, the null hypothesis is rejected, meaning that no significant difference exists between the pixel coordinates produced by the two methods of segmentation. 2) Agreement Analysis: MANOVA is designated to test a significant difference between two groups of measurements. If the test fails, the one-way MANOVA does not indicate whether these two groups of measurements are in agreement. Thus, we
2168-2194 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/JBHI.2014.2332955, IEEE Journal of Biomedical and Health Informatics
> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT)
REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) < TABLE III ACCURACY FOR SEGMENTING THE LOF Pat. # 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 Avg.
MANOVA λ
p
1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
0.99 0.98 1.00 1.00 1.00 1.00 0.98 0.98 0.83 1.00 1.00 0.52
1.00 1.00 1.00 1.00 1.00 0.44 1.00 1.00 1.00 1.00 1.00 -
0.96 1.00 0.90 0.99 1.00 0.00 0.99 1.00 0.97 0.98 0.99 -
Mean Diff. (mm)
RMS (mm)
Std. Dev. (mm)
TABLE V ACCURACY FOR SEGMENTING THE RHF Max (mm)
0.68 1.08 0.84 5.50 1.56 2.28 1.66 13.05 0.82 1.49 1.25 9.49 1.15 1.89 1.51 9.24 1.01 1.54 1.16 8.61 1.28 1.81 1.28 7.68 1.90 2.70 1.92 12.95 3.44 5.14 3.82 26.25 1.60 2.45 1.86 15.40 0.89 1.41 1.10 8.42 1.05 2.00 1.71 12.83 1.99 2.84 2.03 11.80 Surgically Removed 0.76 1.02 0.68 5.13 0.96 1.25 0.81 5.27 3.89 6.03 4.61 24.19 1.51 2.22 1.63 9.71 1.34 1.72 1.08 7.32 N/A (p < 0.05) 1.26 1.78 1.26 7.12 1.77 2.49 1.76 9.24 1.72 2.21 1.40 9.33 0.75 1.29 1.05 8.88 1.44 1.98 1.37 7.68 1.49 2.21 1.63 10.69
±3 mm (%)
Pat. #
97.28 83.16 90.09 89.67 93.82 90.69 80.20 61.94 84.93 96.07 90.92 73.10
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 Avg.
96.99 97.80 58.32 87.31 93.77 89.89 81.94 85.31 95.11 86.76 86.59
TABLE IV ACCURACY FOR SEGMENTING THE ROF
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 Avg.
MANOVA λ 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 -
MANOVA λ
p
1.00 1.00 0.98 1.00 1.00 1.00 1.00 1.00 0.98 1.00 1.00 0.16 1.00 0.95 1.00 1.00 0.96 0.11 1.00 0.99 1.00 1.00 1.00 0.98 -
1.00 1.00 0.82 1.00 0.88 0.81 0.91 0.94 0.11 0.52 1.00 0 0.91 0.63 0.94 1.00 0.13 0 0.98 0.54 1.00 0.95 0.98 0.63 -
Mean Diff. (mm)
RMS (mm)
0.89 0.89 2.60 0.53 2.50 1.47 1.26 1.23 2.73 0.69 0.52
1.49 1.26 3.62 0.86 3.83 2.17 1.88 1.67 3.91 0.98 0.90
2.57 3.55 2.11 1.24 3.03 2.44 1.60 1.94 0.65 0.50 1.91 1.66
Std. Dev. (mm)
Max (mm)
1.20 5.01 0.89 3.29 2.58 8.38 0.68 3.50 2.91 18.87 1.61 9.08 1.40 7.25 1.13 5.47 2.81 13.17 0.69 3.29 0.74 4.02 N/A (p < 0.05) 4.01 3.09 10.95 3.92 1.70 6.81 3.03 2.19 8.44 1.73 1.22 4.95 4.87 3.84 15.48 N/A (p < 0.05) 2.90 1.57 5.47 2.43 1.83 9.04 2.58 1.70 8.92 1.00 0.76 3.98 0.81 0.64 3.29 3.14 2.51 8.16 2.38 1.69 7.50
±3 mm (%) 98.15 98.00 71.43 98.00 75.31 89.25 87.71 90.74 65.00 99.54 97.59 74.14 38.89 72.22 96.54 69.81 63.74 82.33 76.92 99.35 99.35 73.68 82.62
4
p
Mean Diff. (mm)
RMS (mm)
Std. Dev. (mm)
Max (mm)
±3 mm (%)
0.99 0.99 0.97 0.99 1.00 1.00 0.98 0.58 0.98 0.99 1.00 0.99 0.99 1.00 1.00 0.96 0.87 1.00 1.00 0.97 0.94 0.68 1.00 0.99 -
1.16 2.01 1.28 0.65 0.88 1.83 1.13 2.45 2.34 0.97 0.93 0.81 1.70 0.67 1.26 1.43 2.71 1.50 1.26 1.94 2.21 3.67 0.63 1.13 1.52
1.90 3.55 2.48 1.12 1.70 2.72 1.68 4.03 3.25 1.55 1.50 1.23 2.62 1.41 2.12 2.80 4.57 2.16 1.84 3.35 2.92 7.12 1.10 1.66 2.51
1.50 2.94 2.13 0.91 1.46 2.01 1.24 3.20 2.26 1.20 1.17 0.93 2.00 1.24 1.70 2.41 3.68 1.55 1.34 2.73 1.91 6.11 0.90 1.22 1.99
8.88 16.74 18.41 7.25 11.01 13.16 6.57 20.06 14.27 8.92 6.18 12.02 9.91 9.02 10.20 23.00 22.16 12.19 7.98 14.42 12.67 33.05 7.25 6.63 13.00
89.46 80.14 87.55 97.02 90.76 79.78 90.13 76.43 72.62 90.78 93.27 99.05 78.66 90.30 84.81 85.12 66.90 86.26 91.57 78.55 72.25 67.11 96.46 90.16 84.80
segmentation of the right horizontal fissure despite the fuzziness caused by the lesion. Table III, IV, and V present the evaluation result for segmenting the LOF, ROF, and RHF, respectively, in all 24 patients. Within a confidence interval of 95% (p > 0.05), our one-way MANOVA indicated that our algorithm was able to
RMS (mm)
Pat. #
8
3.5 3 2.5 2 -2 10
-1
0
10 10 10 α, the smoothing coefficient
1
10
2
Fig. 9. Average RMS error (of all patients) for different values of the smoothing coefficient, α. [The star indicates the selected α value].
successfully segment (no significant difference found between manual and automatic segmentation) the LOF, ROF, and RHF in 23, 24, and 22 of the 24 patients, respectively. Within these patients, our algorithm yielded root mean square (RMS) errors of 2.21 ± 1.21 mm, 2.51 ± 1.36 mm, and 2.38 ± 1.27 mm for segmenting the LOF, ROF, and RHF, respectively. The corresponding mean differences for the LOF, ROF, and RHF are 1.49 mm, 1.52 mm, and 1.66 mm with maximum distances (worst case errors) of 10.69 mm, 13.00 mm, and 7.50 mm. Using our ±3 mm percentile measure, our algorithm produced accuracies of 86.59%, 84.80%, and 82.62% for segmenting the LOF, ROF, and RHF, respectively. Fig. 9 plots the average RMS error (of all patients) for different values of the smoothing coefficient, α. We selected an α value of 1 to provide the best balance between fissure smoothness and segmentation accuracy. Using this α value, our
2168-2194 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/JBHI.2014.2332955, IEEE Journal of Biomedical and Health Informatics
> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT)
REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) < 2.21 ± 1.21 mm, 2.51 ± 1.36 mm, and 2.38 ± 1.27 mm for segmenting the LOF, ROF, and RHF, respectively. The average accuracies for segmenting the LOF, ROF, and RHF are 86.6%, 84.8%, and 82.62%, correspondingly. These results indicate the feasibility of developing an automatic algorithm for complete segmentation of diseased lung lobes. Future work includes: further evaluation of our algorithm using more diseased lungs and enhance the performance of the algorithm by porting it from Matlab into C/C++. REFERENCES [1] [2] [3]
[4]
[5] [6] [7] [8] [9] [10]
[11] [12] [13] [14] [15] [16]
[17]
[18]
American Cancer Society, Cancer Facts and Figures 2008. Atlanta:American Cancer Society, 2008. Y. Hu and R. A. Malthaner, “The feasibility of three-dimensional displays of the thorax for preoperative planning in the surgical treatment of lung cancer,” Europ. Cardio-Thoracic Surg., vol. 31, pp. 506-511, 2007. J.-M. Kuhnigk, V. Dicken, S. Zidowitz, L. Bornemann, B. Kuemmerlen, S. Krass, H.-O. Peitgen, S. Yuval, H.-H. Fend, W. S. Rau, and T. Achenbach, “Newtools for computer assistance in thoracic CT. Part 1. Functional analysis of lungs, lung lobes, and bronchopulmonary segments,” RadioGraphics, vol. 25, pp. 525–536, 2005. B.M.Hemminger, P. L. Molina,T. M. Egan, F.C.Detterbeck,K. E.Muller, C. S. Coffrey, and J. K. Lee, “Assessment of real-time 3D visualization for cardiothoracic diagnostic evaluation and surgery planning,” J. Digit. Imag., vol. 18, pp. 145–153, 2005. W. R. Webb, N. L. Muller, and D. P. Naidich, High-resolution CT of the lung, 3rd ed., Philadelphia, PA: Lippincott, Williams and Wilkins, 2001. K. Hayashi, A. Aziz, K. Ashizawa, H. Hayashi, K. Nagaoki, and H. Otsuji, “Radiographic and CT appearances of the major fissures,” Radiographics, vol. 21, pp. 861-874, 2001. B. N. Raasch, E. W. Carsky, E. J. Lane, J. P. O’Callaghan, and E. R. Heitzman, “Radiographic anatomy of the interlobar fissures: a study of 100 specimens,” AJR, vol. 138, pp. 647-554, 1982. M. Kass, A. Witkin, and D. Terzopoulos, “Snakes: active contour models,” Int. J. Comp. Vis., vol. v1, pp. 321-331, 1987. R. Malladi, J. A. Sethian, and B. C. Vemuri, “Shape modeling with front propagation: a level set approach,” IEEE Trans. Pattern Anal. Mach. Intell.,vol. 17, pp. 158-175, 1995. M. Kubo, N. Niki, S. Nakagawa, K. Eguchi, M. Kaneko, N. Moriyama, H. Omatsu, R. Kakinuma, and N. Yamaguchi, “Extraction algorithm of pulmonary fissures from thin-section CT images based on linear feature detector method,” IEEE Trans. Nucl. Sci., vol. 46, pp. 2128–2133, 1999. L. Zhang, E. A. Hoffman, and J. M. Reinhardt, “Atlas-driven lung lobe segmentation in volumetric x-ray CT images,” IEEE Trans. Med. Imag., vol. 25, pp. 1-16, 2006. J. Wang, M. Betke, and J. P. Ko, “Segmentation of pulmonary fissures on diagnostic CT – preliminary experience,” Proc. Int. Conf. Diagn. Imag. Anal., 2002, pp.107-112. Q. Wei, Y. Hu, G. Gelfand, and J. H. Macgregor, “Segmentation of lung lobes in high-resolution isotropic CT images,” IEEE Trans. Biomed. Eng., vol. 56, pp. 1383-1393, 2009. J.-M. Kuhnigk, H. K. Hahn, et.al, “Lung lobe segmentation by anatomy-guided 3-D watershed transform,” Proc. SPIE (Med. Imag.), vol. 5032, 2003, pp. 1482-1490. S. Ukil and J. M. Reinhardt, “Anatomy-guided lung lobe segmentation in x-ray CT images,” IEEE Trans. Med. Imag., vol. 28, pp. 202-214, 2009. X. Zhou, T. Hayashi, T. Hara, H. Fujita, R. Yokoyama, T. Kiryu, and H. Hoshi, “Automatic segmentation and recognition of anatomical lung structures from high-resolution chest ct images,” Compu. Med. Imag. Graphics, vol. 30, pp. 299-313, 2006. S. Saita, M. Kubo, Y. Kawata, N. Niki, H. Ohmatsu, and N. Moriyama, “An algorithm for the extraction of pulmonary fissures from low-dose multislice CT images,” Systems and Computers in Japan, vol. 37, pp. 63-76, 2006. E. M. van Rikxxort and B. V. Ginneken, “Supervised enhancement filters: application to fissure detection in chest CT scan,” IEEE Trans. Med. Imag., vol. 27, pp. 1-10, 2008.
10
[19] E. M. van Rikxoort, M. Prokop, B. Hoop, M. A. Viergever, J. P. W. Pluim, and B. van Ginneken, “Automatic Segmentation of Pulmonary Lobes Robust Against Incomplete Fissures,” IEEE Trans. Med. Imag., vol. 29, pp. 1286-1296, 2010. [20] J. Pu, J. K. Leader, B. Zheng, F. Knollmann, C. Fuhrman, F. C. Sciurba, and D. Gur, “A computational geometry approach to automated pulmonary fissure segmentation in CT examinations,” IEEE Trans. Med. Imag., vol. 28, pp. 710-719, 2009. [21] J. Pu, B. Zheng, J. K. Leader, C. Fuhrman, F. Knollmann, A. Klym, and D. Gur, “Pulmonary lobe segmentation in CT examinations using implicit surface fitting,” IEEE Trans. Med. Imag., vol. 28, pp. 1986-1996, 2009. [22] D. Aykac, E. A. Hoffman, G. McLennan, and J. M. Reinhardt “Segmentation and analysis of the human airway tree from three-dimensional X-ray CT imags,” IEEE Trans. Med. Imag., vol. 22, pp. 940-950, 2003. [23] J. Tschirren, E. A. Hoffman, G. McLennan, and M. Sonka, “Intrathoracic airway trees: segmentation and airway morphology analysis from low-dose CT scans,” IEEE Trans. Med. Imag., vol. 24, pp. 1529-1539, 2005. [24] Q. Wei, Y. Hu, J. H. MacGregor and G. Gelfand, “Automatic recognition of major fissures in human lungs,” Int. J. CARS, vol. 7, pp. 111-123, 2012. [25] M. Galloway, “Texture analysis using gray level run lengths,” Comput. Graphics Image Proc., vol. 4, pp. 172-199, 1974. [26] J. A. Anderson, An Introduction to Neural Networks. Cambridge, MA: MIT Press, 1995. [27] M. Sonka, V. Hlavac, and R. Boyle, Image Processing, Analysis, and Machine Vision, 3rd ed., Toronto, ON: Thomson Learning, 2008. [28] A. Aziz, K. Ashizawa, K. Nagaoki, and K. Hayashi, “High resolution CT anatomy of the pulmonary fissures,” J. Thoracic Imag., vol. 19, pp. 186–191, 2004. [29] A. Helzer, M. Barzohar, and D. Mala, “Stable fitting of 2D curves and 3D surfaces by implicit polynomials,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 26, pp. 1283-1294, 2004. [30] V. Weiss, L. Andor, G. Renner, and T. Varady, “Advanced surface fitting techniques,” Computer Aided Geometric Design, vol. 19, pp. 19-42, 2002. [31] T. Tasdizen, J.-P. Tarel and D. B. Cooper, “Improving the stability of algebraic curves for applications,” IEEE Trans. Image Process., vol. 9, pp. 405-416, 2000. [32] B. G. Tabachick and L. S. Fidell, Experimental Designs Using ANOVA, Pacific Grove, CA: Buxbury Press, 2006. [33] J. M. Bland and D. G. Altman, “Statistical methods for assessing agreement between two methods of clinical measurement,” Lancert, vol. i, pp. 307-310, 1986.
2168-2194 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT)
REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) < Following this approach, Kubo et al. used the VanderBurg linear feature detector and morphological operators to find fissures in 25 lungs, 12 of which contained pulmonary nodules [10]. Zhang et al. combined a pulmonary atlas, a ridgeness operator and fuzzy reasoning to identify fissures in 7 normal and 5 diseased lungs [11]. Using a similar ridgeness operator to enhance the fissures, Wang et al. employed Bayesian theory to identify fissures in 10 lungs with pulmonary nodules [12]. Their algorithm yielded high segmentation accuracy with manually delineated fissures for initialization. To eliminate manual intervention, we developed an algorithm using wavelet transform to locate fissures in 9 lungs with pulmonary nodules [13]. In general, the 2D fissure approach works well for identifying disrupted and deformed fissures in individual CT images; however, it does not always produce complete segmentation of lung lobes, due to the lack of continuity between fissures in adjacent CT images. 2) 3D Surface Approach: Another approach of segmenting lung lobes is by finding fissure surfaces using 3D anatomic structures within the lungs. With a Kuhnigk et al. used a segmented vascular tree and an interactive watershed transform (IWT) to find fissure surfaces in 24 lungs [14]. Ukil et al. extended this algorithm by automating the IWT using segmented vascular and bronchial trees [15]. They then used a graph search to refine fissure surfaces in 12 normal and 17 emphysema lungs. Their algorithm worked well for lungs with emphysema; however, one quarter (25%) of their lungs required manual correction [15]. Following the same approach, Zhou et al. used Voronoi division to locate fissure surfaces in 7 normal lungs [16]. Similarly, Saita et al. employed a sheet-emphasis filter to extract fissure surfaces in 16 normal lungs and 4 lungs with unknown diseases [17]. Using a slightly different strategy, Van Rikxxort et al. first enhanced fissures in individual CT images using a set of supervised filters [18], and then utilized a segmented bronchial tree to find fissure surfaces in 20 unspecified normal/abnormal lungs [19]. Lastly, Pu et al. used a computational geometry method to first locate fissure patches [20] and then used implicit radial basis functions to extend these patches into fissure surfaces in 65 lungs with chronic obstructive pulmonary disease [21]. While the 3D surface approach yields complete segmentation of the lung lobes, it often requires difficult segmentation of anatomic structures such as the bronchial and vascular trees [22-23]. More importantly, current 3D surface approaches have difficulties in handling disrupted and deformed fissures, which are often caused by cancer lesions. To the best of our knowledge, no existing methods have been tested under all three constraints of: (1) automation, (2) dealing with multiple diseases, and (3) producing a complete segmentation in the 3D space, as listed in Table I. To meet these constraints, we developed a new algorithm employing a hybrid 2D/3D approach by first finding fissures in 2D CT images using a texture analysis method, followed by a surface fitting technique based on ridge regression to model actual fissure surfaces in 3D. Compared to our previous algorithms in [13] and [24], the novelties of this algorithm include the following two aspects: fully automatic segmentation of 3D major (oblique) and minor (horizontal) fissures using ridge regression
2
TABLE I STATE-OF-THE-ART ALGORITHMS FOR SEGMENTING ALL LUNG LOBES Evaluated on Segmentation Authors Automation Multiple in 3D Diseases 2D Fissure Approach Kubo et al. [10] Zhang et al. [11] Wang et al. [12] Wei et al. [13]
YES NO NO NO
NO YES NO NO
NO NO NO NO
NO NO NO NO NO NO
YES YES YES YES YES YES
3D Surface Approach Kuhnigk et al. [14] Ukil et al. [15] Zhou et al. [16] Saita et al. [17] Van Rikxoort et al. [18, 19] Pu et al. [20, 21]
NO NO YES YES YES YES
in diseased lungs which often have disrupted, deformed and/or incomplete fissures; and an increase of segmentation accuracy in about 10 percent. We evaluated our algorithm using lungs with a variety of different diseases: emphysema, small cell lung cancer, large cell lung cancer and bronchiectasis. To present our algorithm, we organized this paper as follows: Section II describes our methodology. Section III and IV provides the results and discussion, respectively. Lastly, Section V gives the conclusion and future work. II. METHODOLOGY Our algorithm (developed in Matlab) uses a pipelined process to segment diseased lung lobes. The algorithm employs four stages: A) preprocessing, B) 2D oblique fissure recognition, C) 2D horizontal fissure recognition and D) 3D lung lobe segmentation. Fig. 1 depicts the flow diagram of this algorithm. The following context presents these four stages and evaluation method in five subsections, respectively. A. Preprocessing In this stage, we preprocessed a stack of CT images using the Wiener filter and lung segmentation methods that we developed previously [13]. The Wiener filter reduces the noise inherited in isotropic CT images, while the lung segmentation method decreases computation time by removing area outside of the lungs. Fig. 2(a) and (b) illustrate example cross-sectional and sagittal images (see subsection C) from a preprocessed CT stack, respectively. B. 2D Oblique Fissure Recognition To identify the oblique fissures in a stack of cross-sectional CT images, we used a method based on our preliminary algorithm [24]. In short, this method uses three steps to recognize the oblique fissures in pathological lungs: (1) texture analysis, (2) fissure region analysis, and (3) fissure identification. In the first step, we performed 2D texture analysis (instead of 3D to reduce computation) by dividing all images in the preprocessed CT stack into overlapping
2168-2194 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/JBHI.2014.2332955, IEEE Journal of Biomedical and Health Informatics
> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <
X
CT Image Stacks
3
Y
A) Preprocessing
Wiener Filter
Lung Segmentation
Y B) Oblique Fissure Recognition
Texture Analysis using Neural Network
C) Horizontal Fissure Recognition
Fissure Identification
(a)
Stack Rotation Left Lung Removal
Oblique Fissure Region Analysis
Z
Lower Lobe Removal Horizontal Fissure Region Analysis
Z
(b) X
Fig 2. Examples of (a): cross-sectional images from a preprocessed CT stack and (b): sagittal images from the same CT stack.
(a)
(b)
(c)
(d)
Fissure Identification
D) Lung Lobe Segmentation
Fissure Verification
Surface Modeling
Fig. 1. Flowchart of the lobe segmentation algorithm. 100 Y (Pixels)
region-of-interests (ROIs) with dimensions of 10 × 10 pixels. The selection of 10 pixels corresponds to the width of typical fissure regions. We used gray level run length texture features [25] and a feed-forward back-propagation neural network (NN) [26] to classify these ROIs as either a fissure or non-fissure ROI. For NN training, we employed a set of predefined ROIs in the lungs. As a result, the first step of texture analysis yields a stack of binary images containing classified fissure and non-fissure ROIs. Fig. 3(a) and (b) illustrates an original CT image and the matching binary image yielded by the first step, respectively. To find regions containing the fissures (fissure regions) in the second step of fissure region analysis, we performed a stack projection followed by a line fitting. The stack projection transforms the binary stack into two accumulation images (a left and a right accumulation image) showing the propagation of fissure regions in the y-z plane (see Fig. 2). The two accumulation images correspond to the left and right lungs. For each accumulation image, the line fitting calculates a straight line (propagation line) to extract the locations of the fissure regions for each image in the binary stack. The propagation line calculated from the right accumulation image is also used for the horizontal fissure recognition (see subsection C). Fig 3(c) and (d) presents an example right accumulation image and the corresponding propagation line calculated from the line fitting, respectively. Fig. 3(e) depicts fissure regions (white pixels) found within the binary image in Fig. 3(b) using the calculated propagation line. In the last step of fissure identification, we combined the found fissure regions in the binary stack and the images in the preprocessed CT stack to identify the actual fissures. Fig. 3(f) gives example fissures found within the fissure regions shown
80 60 40 20 0
(e)
0
20 40 60 80 100 X (Pixels)
(f)
Fig. 3. Example images from oblique fissure recognition: (a) original CT image; (b) binary image containing classified fissure (white pixels) or non-fissure ROIs (black pixels); (c) right accumulation image; (d) line fitting of pixels on the accumulation image shown in (c); (e) found fissure regions (white pixels); (f) partial fissures found within fissure regions shown in (e).
in Fig. 3(e). As illustrated in Fig. 3(f), fissures found by the algorithm do not extend to edges of the lungs due to diseases and incomplete fissures (addressed in subsection D). More detailed description of the above method can be found in [24]. C. 2D Horizontal Fissure Recognition Our previous work [24] did not identify the horizontal fissure that is more challenging to segment than the oblique fissures. Compared to the oblique fissures, the horizontal
2168-2194 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/JBHI.2014.2332955, IEEE Journal of Biomedical and Health Informatics
> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT)
REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <
(a)
denotes the ith cross-sectional CT image that the fissure pixel, pi, resides in. Using this set of n validated voxels, {vi(xi, yi, zi), i=1,2…n}, we modeled a fissure surface in the form of:
(b)
z f ( x, y ). (c)
Fig. 5. Horizontal fissure region analysis: (a) binary image showing the largest region (white pixels) and the selected neighboring regions (gray pixels), as the horizontal fissure region; (b) horizontal fissure region found using bounding box of the selected regions in (a); and (c) horizontal fissure found within the fissure region using dynamic programming.
5) Fissure Identification: We used the bounding box of all selected regions in the binary images to define the fissure region in the sagittal CT images. Fig. 5(b) shows an example fissure region found from Fig. 5(a). Within the fissure region, we performed dynamic programming using Eq. (2) to find the optimal path extending from the left to the right of the fissure region. Fig. 5(c) illustrates a fissure found from the fissure region in Fig. 5(b) using dynamic programming. D. 3D Lung Lobe Segmentation The oblique and horizontal fissure recognition uses 2D texture analysis and fissure anatomy to find the fissures. Although this approach yields good approximations of the fissures in individual 2D CT images, it does not produce complete segmentation 3D. This is due to the fact that many fissures appear incomplete in 2D CT images. As well, lung diseases may disrupt and deform existing fissures causing discontinuities between fissures in adjacent CT images. Fissures appear naturally as 3D surfaces separating adjacent lung lobes. Therefore, the most logical method for producing a realistic segmentation of the lung lobes is to model a 3D surface (fissure surface) for each of the LOF, ROF, and RHF, using fissures found in individual 2D CT images. This hybrid 2D/3D approach also provides the benefit of being able to fill in incomplete and disrupted fissures automatically. To ensure only the correct fissures are used to model the fissure surface, we first performed fissure verification by comparing fissures in adjacent CT images. For each fissure on a cross-sectional CT image, we computed its accuracy as the percentage of its pixels satisfying the following condition: ( y1 y 2 ) 2 ( x1 x2 ) 2 2 pixels ,
5
(3)
where (x1, y1) and (x2, y2) are the coordinates of a fissure pixel and its closest counterpart in the previous adjacent CT image, respectively. The 2 pixel criterion in Eq. (3) ensures that fissures in neighboring CT images are connected. We chose all fissures with an accuracy of >80% (based on experimentation), in 5 or more consecutive CT images, as the correct fissures. We selected a range of 5 or more consecutive images to prevent false verification due to a single incorrectly found fissure. Once the fissures are verified, we represented each fissure pixel pi(xi, yi), as a voxel (a point in 3D), vi(xi, yi, zi), where zi
(4)
Many 3D modeling methods exist [29, 30]. One method is the minimization of energy functions, as used by Pu et al. [21]. While this method produced smooth approximations of fissure surfaces, it is computationally intensive (Pu et al. had to down sample from 200,000 points to 3,000 points during their surface modeling). For our surface modeling, we used a method based on ridge regression [31]. Ridge regression introduces regularization to the standard least squares (LS) formulation. This method carries several advantages over common fissure modeling techniques: (1) provides control over the smoothness of the fissure in the 3D plane; (2) reduces the sensitivity to erroneously found fissures; (3) is computationally and memory efficient. To derive ridge regression for our purpose of fissure modeling, we start by using linear regression to model the relationship between the z and x-y coordinates of each fissure voxel with the equation: zˆi w0 w1 xi w2 yi e ,
(5)
where zˆi is the estimated value of zi, while e represents the modeling error. Similarly, for n voxels, we can model the same relationship using a set of n equations, expressed in a matrix notation: z Aw e ,
(6)
which has the expanded form of: z1 1 x1 z 1 x 2 2 z n 1 xn
y1 e1 w0 y2 e2 w1 . w2 yn en
(7)
Thus, Eq. (6) represents a linear regression model, where w denotes the parameters of our fissure surface. The simplest method for solving this linear regression model would be the least squares (LS) method. The LS method estimates w, by minimizing the sum of squares of the errors:
w min min Aw z
2
.
(8)
Despite the simplicity of the LS method, it carries two drawbacks: (1) local inconsistency with the continuity of input data; (2) oversensitivity of small errors in the input causing a few outliers to skew the estimation or prediction model. Ridge regression alleviates these two drawbacks of the LS method by using regularization. Given a linear regression model shown in Eq. (6), ridge regression estimates the parameters, w, by performing the following minimization:
w min min Aw z
2
w
2
,
(9)
where Γ is the regularization term. The implementation of ridge regression is as follows: (1) transform the LS formulation
2168-2194 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/JBHI.2014.2332955, IEEE Journal of Biomedical and Health Informatics
> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) < to obtain uncorrelated components, (2) apply regularization to the model to achieve a desired property and (3) restore the original coordinate system by reversing the original transform. In the first step of ridge regression, we transformed the (x, y) coordinates of the voxels into a linear index aij by using the equation: 1， if j yi xi Dx aij , otherwise 0， i [1, n], j [1, m], m Dx D y ,
(10)
where Dx and Dy are the pixel dimensions of a cross-sectional CT image in the x- and y-axes, respectively. Thus, the ridge regression formulation becomes: z1 a11 z a 2 21 z n an1
a12 a1m w1 e1 a 22 a2 m w2 e2 . a n 2 anm wm en
(11)
In the second step of ridge regression, we applied regularization. Regularization provides several advantages for fissure surface modeling. One, it allows the modeling of a non-linear fissure surface using a simple linear model. Two, it allows control over the smoothness of the surface. Three, it reduces sensitivity to outliers in the data, i.e. erroneously found fissures. In order to produce a smooth fissure surface, we used Laplacian functions (along the two independent axes, x and y) as our regularization term. That is, 2 f ( x, y ) 2 2 x , f ( x, y ) 2 y
(12)
where α is a scale factor that controls the effect of regularization. The minimization of the Laplacian functions ensures minimal changes in the gradient of our fissure surface. This allows us to control the smoothness of the fissure surface by varying the scale factor, α. We selected an optimal α value of 1 for our evaluation (see Section III, regarding the effect of α on the accuracy of segmentation). In the last step of ridge regression, we computed our fissure surface in the original x-y coordinates, z(xi, yi), by reversing the transformation of Eq. (10). We reversed the transformation by using the equation: z ( xi , y j ) w min (i jDx ),
(13)
where wmin is the result of Eq. (9) while i, j, and Dx are as defined in Eq. (10). E. Evaluation Method For evaluation, we used isotropic CT image stacks from 24 anonymous patients (following Canadian tri-council ethics guidelines) at the Foothills Medical Center, Calgary, AB, Canada. A Siemens Sensation 16 multi-slice CT scanner acquired these image stacks using the standard Foothills
Patient 1–6 7–9 10 – 13 14 – 18 19 – 20 21 22 23 24
6
TABLE II PATIENT PATHOLOGY Types of Lung Disease Mild Emphysema Moderate Emphysema Large Cell Lung Cancer (LCLC) Small Cell Lung Cancer (SCLC) Bronchiectasis Elevated Left Hemidiaphragm Hilar Mass Neoplastic No Diagnosis
Medical Center Thoracic Abdominal 3D Research protocol and iodine contrast agent. All images have a resolution of 512 × 512 pixels with a thickness of 0.6 mm. Of the 24 patients, 23 had lung diseases and one had no diagnosis. Table II lists the types of diseases associated with these patients. To evaluate our algorithm, we employed a two-fold cross validation method [26] by splitting the datasets from the 24 patients into two groups. For each group, we randomly selected half of the datasets containing each disease. This ensures a good representation of all diseases for each group. For diseases where only one dataset is present (Patients 21-24), we randomly placed the dataset into one of the two groups. Of the two groups, we used one group for training and the remaining group for testing. We swapped the two groups to test all 24 datasets. Finally, we repeated the above procedure 10 times and averaged the result. For performance evaluation, we compared the output of our algorithm to that of manual segmentation. While our algorithm yields complete segmentation of the lung lobes, it is extremely time consuming for a radiologist to completely trace all fissures in 24 patients with an average of 300 CT images per patient. To make the evaluation process more feasible, we used manually placed anchor points to mimic manual segmentation. Under the guidance of a radiologist, we manually placed anchor points (at about 5 pixels apart) on fissures for every 5th image of a CT stack, which is approximately the thickness of the CT images that surgeons/radiologists read in clinical settings. Using this manual segmentation, we performed the following analysis: 1) One-way Multivariate Analysis of Variance (MANOVA): One-way MANOVA examines statistically a significant difference between two groups of measurements [32]. In our case, the two groups of measurements are the fissure pixel coordinates, (xi, yi, zi), produced by our algorithm and manual segmentation. By computing a corresponding pair of Wilk’s lambda (λ) and p-value, one-way MANOVA tests the null hypothesis that the two groups of pixel coordinates are different. If the p-value is greater than 0.05, the null hypothesis is rejected, meaning that no significant difference exists between the pixel coordinates produced by the two methods of segmentation. 2) Agreement Analysis: MANOVA is designated to test a significant difference between two groups of measurements. If the test fails, the one-way MANOVA does not indicate whether these two groups of measurements are in agreement. Thus, we
2168-2194 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/JBHI.2014.2332955, IEEE Journal of Biomedical and Health Informatics
> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT)
REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) < TABLE III ACCURACY FOR SEGMENTING THE LOF Pat. # 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 Avg.
MANOVA λ
p
1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
0.99 0.98 1.00 1.00 1.00 1.00 0.98 0.98 0.83 1.00 1.00 0.52
1.00 1.00 1.00 1.00 1.00 0.44 1.00 1.00 1.00 1.00 1.00 -
0.96 1.00 0.90 0.99 1.00 0.00 0.99 1.00 0.97 0.98 0.99 -
Mean Diff. (mm)
RMS (mm)
Std. Dev. (mm)
TABLE V ACCURACY FOR SEGMENTING THE RHF Max (mm)
0.68 1.08 0.84 5.50 1.56 2.28 1.66 13.05 0.82 1.49 1.25 9.49 1.15 1.89 1.51 9.24 1.01 1.54 1.16 8.61 1.28 1.81 1.28 7.68 1.90 2.70 1.92 12.95 3.44 5.14 3.82 26.25 1.60 2.45 1.86 15.40 0.89 1.41 1.10 8.42 1.05 2.00 1.71 12.83 1.99 2.84 2.03 11.80 Surgically Removed 0.76 1.02 0.68 5.13 0.96 1.25 0.81 5.27 3.89 6.03 4.61 24.19 1.51 2.22 1.63 9.71 1.34 1.72 1.08 7.32 N/A (p < 0.05) 1.26 1.78 1.26 7.12 1.77 2.49 1.76 9.24 1.72 2.21 1.40 9.33 0.75 1.29 1.05 8.88 1.44 1.98 1.37 7.68 1.49 2.21 1.63 10.69
±3 mm (%)
Pat. #
97.28 83.16 90.09 89.67 93.82 90.69 80.20 61.94 84.93 96.07 90.92 73.10
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 Avg.
96.99 97.80 58.32 87.31 93.77 89.89 81.94 85.31 95.11 86.76 86.59
TABLE IV ACCURACY FOR SEGMENTING THE ROF
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 Avg.
MANOVA λ 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 -
MANOVA λ
p
1.00 1.00 0.98 1.00 1.00 1.00 1.00 1.00 0.98 1.00 1.00 0.16 1.00 0.95 1.00 1.00 0.96 0.11 1.00 0.99 1.00 1.00 1.00 0.98 -
1.00 1.00 0.82 1.00 0.88 0.81 0.91 0.94 0.11 0.52 1.00 0 0.91 0.63 0.94 1.00 0.13 0 0.98 0.54 1.00 0.95 0.98 0.63 -
Mean Diff. (mm)
RMS (mm)
0.89 0.89 2.60 0.53 2.50 1.47 1.26 1.23 2.73 0.69 0.52
1.49 1.26 3.62 0.86 3.83 2.17 1.88 1.67 3.91 0.98 0.90
2.57 3.55 2.11 1.24 3.03 2.44 1.60 1.94 0.65 0.50 1.91 1.66
Std. Dev. (mm)
Max (mm)
1.20 5.01 0.89 3.29 2.58 8.38 0.68 3.50 2.91 18.87 1.61 9.08 1.40 7.25 1.13 5.47 2.81 13.17 0.69 3.29 0.74 4.02 N/A (p < 0.05) 4.01 3.09 10.95 3.92 1.70 6.81 3.03 2.19 8.44 1.73 1.22 4.95 4.87 3.84 15.48 N/A (p < 0.05) 2.90 1.57 5.47 2.43 1.83 9.04 2.58 1.70 8.92 1.00 0.76 3.98 0.81 0.64 3.29 3.14 2.51 8.16 2.38 1.69 7.50
±3 mm (%) 98.15 98.00 71.43 98.00 75.31 89.25 87.71 90.74 65.00 99.54 97.59 74.14 38.89 72.22 96.54 69.81 63.74 82.33 76.92 99.35 99.35 73.68 82.62
4
p
Mean Diff. (mm)
RMS (mm)
Std. Dev. (mm)
Max (mm)
±3 mm (%)
0.99 0.99 0.97 0.99 1.00 1.00 0.98 0.58 0.98 0.99 1.00 0.99 0.99 1.00 1.00 0.96 0.87 1.00 1.00 0.97 0.94 0.68 1.00 0.99 -
1.16 2.01 1.28 0.65 0.88 1.83 1.13 2.45 2.34 0.97 0.93 0.81 1.70 0.67 1.26 1.43 2.71 1.50 1.26 1.94 2.21 3.67 0.63 1.13 1.52
1.90 3.55 2.48 1.12 1.70 2.72 1.68 4.03 3.25 1.55 1.50 1.23 2.62 1.41 2.12 2.80 4.57 2.16 1.84 3.35 2.92 7.12 1.10 1.66 2.51
1.50 2.94 2.13 0.91 1.46 2.01 1.24 3.20 2.26 1.20 1.17 0.93 2.00 1.24 1.70 2.41 3.68 1.55 1.34 2.73 1.91 6.11 0.90 1.22 1.99
8.88 16.74 18.41 7.25 11.01 13.16 6.57 20.06 14.27 8.92 6.18 12.02 9.91 9.02 10.20 23.00 22.16 12.19 7.98 14.42 12.67 33.05 7.25 6.63 13.00
89.46 80.14 87.55 97.02 90.76 79.78 90.13 76.43 72.62 90.78 93.27 99.05 78.66 90.30 84.81 85.12 66.90 86.26 91.57 78.55 72.25 67.11 96.46 90.16 84.80
segmentation of the right horizontal fissure despite the fuzziness caused by the lesion. Table III, IV, and V present the evaluation result for segmenting the LOF, ROF, and RHF, respectively, in all 24 patients. Within a confidence interval of 95% (p > 0.05), our one-way MANOVA indicated that our algorithm was able to
RMS (mm)
Pat. #
8
3.5 3 2.5 2 -2 10
-1
0
10 10 10 α, the smoothing coefficient
1
10
2
Fig. 9. Average RMS error (of all patients) for different values of the smoothing coefficient, α. [The star indicates the selected α value].
successfully segment (no significant difference found between manual and automatic segmentation) the LOF, ROF, and RHF in 23, 24, and 22 of the 24 patients, respectively. Within these patients, our algorithm yielded root mean square (RMS) errors of 2.21 ± 1.21 mm, 2.51 ± 1.36 mm, and 2.38 ± 1.27 mm for segmenting the LOF, ROF, and RHF, respectively. The corresponding mean differences for the LOF, ROF, and RHF are 1.49 mm, 1.52 mm, and 1.66 mm with maximum distances (worst case errors) of 10.69 mm, 13.00 mm, and 7.50 mm. Using our ±3 mm percentile measure, our algorithm produced accuracies of 86.59%, 84.80%, and 82.62% for segmenting the LOF, ROF, and RHF, respectively. Fig. 9 plots the average RMS error (of all patients) for different values of the smoothing coefficient, α. We selected an α value of 1 to provide the best balance between fissure smoothness and segmentation accuracy. Using this α value, our
2168-2194 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/JBHI.2014.2332955, IEEE Journal of Biomedical and Health Informatics
> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT)
REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) < 2.21 ± 1.21 mm, 2.51 ± 1.36 mm, and 2.38 ± 1.27 mm for segmenting the LOF, ROF, and RHF, respectively. The average accuracies for segmenting the LOF, ROF, and RHF are 86.6%, 84.8%, and 82.62%, correspondingly. These results indicate the feasibility of developing an automatic algorithm for complete segmentation of diseased lung lobes. Future work includes: further evaluation of our algorithm using more diseased lungs and enhance the performance of the algorithm by porting it from Matlab into C/C++. REFERENCES [1] [2] [3]
[4]
[5] [6] [7] [8] [9] [10]
[11] [12] [13] [14] [15] [16]
[17]
[18]
American Cancer Society, Cancer Facts and Figures 2008. Atlanta:American Cancer Society, 2008. Y. Hu and R. A. Malthaner, “The feasibility of three-dimensional displays of the thorax for preoperative planning in the surgical treatment of lung cancer,” Europ. Cardio-Thoracic Surg., vol. 31, pp. 506-511, 2007. J.-M. Kuhnigk, V. Dicken, S. Zidowitz, L. Bornemann, B. Kuemmerlen, S. Krass, H.-O. Peitgen, S. Yuval, H.-H. Fend, W. S. Rau, and T. Achenbach, “Newtools for computer assistance in thoracic CT. Part 1. Functional analysis of lungs, lung lobes, and bronchopulmonary segments,” RadioGraphics, vol. 25, pp. 525–536, 2005. B.M.Hemminger, P. L. Molina,T. M. Egan, F.C.Detterbeck,K. E.Muller, C. S. Coffrey, and J. K. Lee, “Assessment of real-time 3D visualization for cardiothoracic diagnostic evaluation and surgery planning,” J. Digit. Imag., vol. 18, pp. 145–153, 2005. W. R. Webb, N. L. Muller, and D. P. Naidich, High-resolution CT of the lung, 3rd ed., Philadelphia, PA: Lippincott, Williams and Wilkins, 2001. K. Hayashi, A. Aziz, K. Ashizawa, H. Hayashi, K. Nagaoki, and H. Otsuji, “Radiographic and CT appearances of the major fissures,” Radiographics, vol. 21, pp. 861-874, 2001. B. N. Raasch, E. W. Carsky, E. J. Lane, J. P. O’Callaghan, and E. R. Heitzman, “Radiographic anatomy of the interlobar fissures: a study of 100 specimens,” AJR, vol. 138, pp. 647-554, 1982. M. Kass, A. Witkin, and D. Terzopoulos, “Snakes: active contour models,” Int. J. Comp. Vis., vol. v1, pp. 321-331, 1987. R. Malladi, J. A. Sethian, and B. C. Vemuri, “Shape modeling with front propagation: a level set approach,” IEEE Trans. Pattern Anal. Mach. Intell.,vol. 17, pp. 158-175, 1995. M. Kubo, N. Niki, S. Nakagawa, K. Eguchi, M. Kaneko, N. Moriyama, H. Omatsu, R. Kakinuma, and N. Yamaguchi, “Extraction algorithm of pulmonary fissures from thin-section CT images based on linear feature detector method,” IEEE Trans. Nucl. Sci., vol. 46, pp. 2128–2133, 1999. L. Zhang, E. A. Hoffman, and J. M. Reinhardt, “Atlas-driven lung lobe segmentation in volumetric x-ray CT images,” IEEE Trans. Med. Imag., vol. 25, pp. 1-16, 2006. J. Wang, M. Betke, and J. P. Ko, “Segmentation of pulmonary fissures on diagnostic CT – preliminary experience,” Proc. Int. Conf. Diagn. Imag. Anal., 2002, pp.107-112. Q. Wei, Y. Hu, G. Gelfand, and J. H. Macgregor, “Segmentation of lung lobes in high-resolution isotropic CT images,” IEEE Trans. Biomed. Eng., vol. 56, pp. 1383-1393, 2009. J.-M. Kuhnigk, H. K. Hahn, et.al, “Lung lobe segmentation by anatomy-guided 3-D watershed transform,” Proc. SPIE (Med. Imag.), vol. 5032, 2003, pp. 1482-1490. S. Ukil and J. M. Reinhardt, “Anatomy-guided lung lobe segmentation in x-ray CT images,” IEEE Trans. Med. Imag., vol. 28, pp. 202-214, 2009. X. Zhou, T. Hayashi, T. Hara, H. Fujita, R. Yokoyama, T. Kiryu, and H. Hoshi, “Automatic segmentation and recognition of anatomical lung structures from high-resolution chest ct images,” Compu. Med. Imag. Graphics, vol. 30, pp. 299-313, 2006. S. Saita, M. Kubo, Y. Kawata, N. Niki, H. Ohmatsu, and N. Moriyama, “An algorithm for the extraction of pulmonary fissures from low-dose multislice CT images,” Systems and Computers in Japan, vol. 37, pp. 63-76, 2006. E. M. van Rikxxort and B. V. Ginneken, “Supervised enhancement filters: application to fissure detection in chest CT scan,” IEEE Trans. Med. Imag., vol. 27, pp. 1-10, 2008.
10
[19] E. M. van Rikxoort, M. Prokop, B. Hoop, M. A. Viergever, J. P. W. Pluim, and B. van Ginneken, “Automatic Segmentation of Pulmonary Lobes Robust Against Incomplete Fissures,” IEEE Trans. Med. Imag., vol. 29, pp. 1286-1296, 2010. [20] J. Pu, J. K. Leader, B. Zheng, F. Knollmann, C. Fuhrman, F. C. Sciurba, and D. Gur, “A computational geometry approach to automated pulmonary fissure segmentation in CT examinations,” IEEE Trans. Med. Imag., vol. 28, pp. 710-719, 2009. [21] J. Pu, B. Zheng, J. K. Leader, C. Fuhrman, F. Knollmann, A. Klym, and D. Gur, “Pulmonary lobe segmentation in CT examinations using implicit surface fitting,” IEEE Trans. Med. Imag., vol. 28, pp. 1986-1996, 2009. [22] D. Aykac, E. A. Hoffman, G. McLennan, and J. M. Reinhardt “Segmentation and analysis of the human airway tree from three-dimensional X-ray CT imags,” IEEE Trans. Med. Imag., vol. 22, pp. 940-950, 2003. [23] J. Tschirren, E. A. Hoffman, G. McLennan, and M. Sonka, “Intrathoracic airway trees: segmentation and airway morphology analysis from low-dose CT scans,” IEEE Trans. Med. Imag., vol. 24, pp. 1529-1539, 2005. [24] Q. Wei, Y. Hu, J. H. MacGregor and G. Gelfand, “Automatic recognition of major fissures in human lungs,” Int. J. CARS, vol. 7, pp. 111-123, 2012. [25] M. Galloway, “Texture analysis using gray level run lengths,” Comput. Graphics Image Proc., vol. 4, pp. 172-199, 1974. [26] J. A. Anderson, An Introduction to Neural Networks. Cambridge, MA: MIT Press, 1995. [27] M. Sonka, V. Hlavac, and R. Boyle, Image Processing, Analysis, and Machine Vision, 3rd ed., Toronto, ON: Thomson Learning, 2008. [28] A. Aziz, K. Ashizawa, K. Nagaoki, and K. Hayashi, “High resolution CT anatomy of the pulmonary fissures,” J. Thoracic Imag., vol. 19, pp. 186–191, 2004. [29] A. Helzer, M. Barzohar, and D. Mala, “Stable fitting of 2D curves and 3D surfaces by implicit polynomials,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 26, pp. 1283-1294, 2004. [30] V. Weiss, L. Andor, G. Renner, and T. Varady, “Advanced surface fitting techniques,” Computer Aided Geometric Design, vol. 19, pp. 19-42, 2002. [31] T. Tasdizen, J.-P. Tarel and D. B. Cooper, “Improving the stability of algebraic curves for applications,” IEEE Trans. Image Process., vol. 9, pp. 405-416, 2000. [32] B. G. Tabachick and L. S. Fidell, Experimental Designs Using ANOVA, Pacific Grove, CA: Buxbury Press, 2006. [33] J. M. Bland and D. G. Altman, “Statistical methods for assessing agreement between two methods of clinical measurement,” Lancert, vol. i, pp. 307-310, 1986.
2168-2194 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
