A “loop” shape descriptor and its application to automated segmentation of airways from CT scans Jiantao Pu, Chenwang Jin, Nan Yu, Yongqiang Qian, Xiaohua Wang, Xin Meng, and Youmin Guo Citation: Medical Physics 42, 3076 (2015); doi: 10.1118/1.4921139 View online: http://dx.doi.org/10.1118/1.4921139 View Table of Contents: http://scitation.aip.org/content/aapm/journal/medphys/42/6?ver=pdfcov Published by the American Association of Physicists in Medicine Articles you may be interested in Lung texture in serial thoracic CT scans: Assessment of change introduced by image registrationa) Med. Phys. 39, 4679 (2012); 10.1118/1.4730505 Automated segmentation of a motion mask to preserve sliding motion in deformable registration of thoracic CT Med. Phys. 39, 1006 (2012); 10.1118/1.3679009 A modified gradient correlation filter for image segmentation: Application to airway and bowel Med. Phys. 36, 480 (2009); 10.1118/1.3056461 Automated detection and quantitative assessment of pulmonary airways depicted on CT images Med. Phys. 34, 2844 (2007); 10.1118/1.2742777 Automatic segmentation of phase-correlated CT scans through nonrigid image registration using geometrically regularized free-form deformation Med. Phys. 34, 3054 (2007); 10.1118/1.2740467

A “loop” shape descriptor and its application to automated segmentation of airways from CT scans Jiantao Pu Department of Radiology, First Affiliated Hospital of Medical College, Xi’an Jiaotong University, Shaanxi 710061, People’s Republic of China, and Departments of Radiology and Bioengineering, University of Pittsburgh, Pittsburgh, Pennsylvania 15213

Chenwang Jin,a) Nan Yu, and Yongqiang Qian Department of Radiology, First Affiliated Hospital of Medical College, Xi’an Jiaotong University, Shaanxi 710061, People’s Republic of China

Xiaohua Wang Third Affiliated Hospital, Peking University, Beijing, People’s Republic of China, 100029

Xin Meng Department of Radiology, University of Pittsburgh, Pittsburgh, Pennsylvania 15213

Youmin Guo Department of Radiology, First Affiliated Hospital of Medical College, Xi’an Jiaotong University, Shaanxi 710061, People’s Republic of China

(Received 22 October 2014; revised 2 April 2015; accepted for publication 4 May 2015; published 28 May 2015) Purpose: A novel shape descriptor is presented to aid an automated identification of the airways depicted on computed tomography (CT) images. Methods: Instead of simplifying the tubular characteristic of the airways as an ideal mathematical cylindrical or circular shape, the proposed “loop” shape descriptor exploits the fact that the cross sections of any tubular structure (regardless of its regularity) always appear as a loop. In implementation, the authors first reconstruct the anatomical structures in volumetric CT as a three-dimensional surface model using the classical marching cubes algorithm. Then, the loop descriptor is applied to locate the airways with a concave loop cross section. To deal with the variation of the airway walls in density as depicted on CT images, a multiple threshold strategy is proposed. A publicly available chest CT database consisting of 20 CT scans, which was designed specifically for evaluating an airway segmentation algorithm, was used for quantitative performance assessment. Measures, including length, branch count, and generations, were computed under the aid of a skeletonization operation. Results: For the test dataset, the airway length ranged from 64.6 to 429.8 cm, the generation ranged from 7 to 11, and the branch number ranged from 48 to 312. These results were comparable to the performance of the state-of-the-art algorithms validated on the same dataset. Conclusions: The authors’ quantitative experiment demonstrated the feasibility and reliability of the developed shape descriptor in identifying lung airways. C 2015 American Association of Physicists in Medicine. [http://dx.doi.org/10.1118/1.4921139] Key words: airways, tubular structure, shape descriptor, loop, vessels, CT 1. INTRODUCTION Airways constitute the primary conductive structure of the human respiratory system for delivering oxygen and carrying away carbon dioxide. Morphological characteristics of airways, or variations thereof, may have a direct impact on airflow, thereby alternating pulmonary function1–3 and possibly the progression of diseases (e.g., COPD or asthma).4–6 Due to its high temporal and spatial resolutions, computed tomography (CT) is widely used to noninvasively investigate in vivo airway structures and characteristics. However, a single CT examination typically contains a large number of images. It is extremely difficult and very time-consuming to manually trace the airways. Sonka et al. reported7 that manual segmentation of the airway tree in a single CT examination (slice 3076

Med. Phys. 42 (6), June 2015

thickness in their study was 3.0 mm) required ∼7 h of analysis. In addition, manual operations may introduce large interand/or intraoperator variability as well as possible biases.8,9 At the same time, large variations in performance exist among human experts.10 Therefore, it is desirable to develop fully (or at a minimum semi) automated computerized schemes for efficient, consistent, and accurate analysis of the airways as depicted on CT exams. Among available methods for airway tree segmentation, a three-dimensional (3D) region growing procedure is often used as a preprocessing step.11 The underlying motivation is to leverage the relatively high contrast between the airway lumen and the airway wall. However, despite its simplicity and efficiency, a purely region growing based operation frequently leads to leakage into the lung parenchyma (i.e., a sudden

0094-2405/2015/42(6)/3076/9/$30.00

© 2015 Am. Assoc. Phys. Med.

3076

3077

Pu et al.: A loop shape descriptor

explosion) because of partial volume effects and/or image noise (artifacts), in particular, in the presence of diseases (e.g., emphysema). To alleviate the leakage issue and meanwhile identify more airways, a number of solutions have been developed, including (1) morphological based methods,12–15 (2) knowledge or rule based methods,7,16–18 (3) template matching based methods,19–23 (4) machine learning classifiers based methods,24–26 and (5) shape analysis based methods.27–32 A relatively detailed description of these approaches can be found in a review article.33 Despite the intensive efforts, available algorithms still miss a large fraction of small airways. Whereas abnormalities, such as obstruction, frequently occur in peripheral regions, and small airways always constitute a region of interest for investigating various lung diseases.34 Therefore, it is highly important to develop algorithms that are capable of identifying small peripheral airways. In this study, we described a novel “loop” shape descriptor to automatically identify the airways depicted in a volumetric CT examination. Its most distinctive characteristic is the way of exploiting the “tubular” characteristic of the airways, where whether a point located on a tubular shape is determined by whether a loop can be formed around this point. The performance of the developed approach was quantitatively assessed on a publicly available dataset consisting of 20 chest CT examinations acquired on different protocols (e.g., dose, scanners, and reconstruction kernels),35 which were specifically collected for validating airway tree segmentation algorithms. 2. METHODS AND MATERIALS 2.A. Scheme overview

The developed airway segmentation scheme consists of three main steps as illustrated by the flowchart in Fig. 1. First, given a volumetric chest CT examination, the lung volume regions are identified using a thresholding operation by taking advantage of the relatively high contrast between lung parenchyma and surrounding structures [Fig. 2(b)]. This step limits the image analysis within the lung regions and meanwhile improves the efficiency in memory and time. Second, the marching cubes algorithm (MCA)36–38 is used to model the 3D lung anatomical structures as a triangle mesh surface, where airways, vessels, and other diseased or noisy regions (e.g., emphysema) are modeled [Figs. 2(c) and 2(d)]. Third, a loop shape descriptor is used to identify the airways on the basis of the surface model. To handle the variation of the airway walls depicted on CT in density, which directly affects the MCA-based surface modeling, the second step and the third step are performed repeatedly at multiple isovalue thresholds. The identified airways at these thresholds are combined using a Boolean union operation. The implementations of these steps are described separately in Secs. 2.B–2.F.

3077

2.B. Lung volume segmentation

Given the low density of lung parenchyma and airway lumen depicted on CT, a thresholding operation was applied to identify the lung regions, where the threshold was simply set at −400 HU. This threshold is low enough to guarantee the inclusion of lung parenchyma and airways. The example in Fig. 2(b) showed the segmentation result after the application of the described thresholding procedure to the CT scan in Fig. 2(a). This procedure aims to limit the analysis within the lung volume, and it not only improves the computational efficiency in both memory and time but also avoids possible incorrect identification of false positives (e.g., esophagus). 2.C. Lung anatomical structure modeling

Application of the MCA to the scalar field of the identified lung volume will result in a 3D isosurface model of the lung structures as shown in Figs. 2(c) and 2(d). In this study, the front surface whose normal vector points to the outside of an object (i.e., from high intensity to low intensity) is displayed in red, while the back faces whose normal vectors point to the inside of an object (i.e., from low intensity to high intensity) are displayed in green. The directions of the surface normal vectors are used to differentiate convex surfaces from concave surfaces [Fig. 3(b)]. It is notable that the high resolution characteristic of CT typically leads to a huge number of triangles. For example, given a CT exam with a slice thickness of 0.625 mm, modeling the segmented lung volume will result in around 20 000 000 triangles and modeling the entire chest CT scan will result in around 30 000 000 triangles. Hence, a simple lung volume segmentation is recommended to reduce the computational costs in memory and computation. 2.D. A loop shape descriptor 2.D.1. The definition of loop

As illustrated by the examples in Fig. 3(b), a loop shape is defined here as having a path that is circular or curved over on itself, where the end is connected to the beginning. According to this definition, even a structure that has a circular shape but without a closed path, as shown in Fig. 3(c), will be classified as one with nonloop shape. Previously, the tubular structures (e.g., airways or vessels) are widely treated ideally as a circleor cylinderlike shape [Fig. 3(a)], such as the one in Fig. 3(a). Although an ideal mathematical shape assumption or simplification makes it relatively straightforward to use some wellestablished mathematical methods (e.g., curvature analysis33) for detection purpose, there is no biological structure in reality that consistently appears as a circle or a cylinder in a strict

F. 1. Schematic flowchart of the airway segmentation algorithm based on a loop shape descriptor. Medical Physics, Vol. 42, No. 6, June 2015

3078

Pu et al.: A loop shape descriptor

3078

F. 2. Lung volume segmentation and lung anatomy modeling: (a) a CT examination, (b) the identified lung volume in overlay, (c) the surface model of the lung anatomy in the cubic subvolume indicated by the box in (a), and (d) triangle meshes of the local regions indicated by the box in (c).

sense. It is difficult, even impossible to represent the airways consistently and faithfully using an ideal mathematical model, thereby making it challenging to identify the airway trees on the basis of the circular or cylindrical shape assumption. Therefore, it is desirable to have a novel shape descriptor to reliably characterize the tubelike shape in a relatively “loose” manner. As demonstrated by the examples in Fig. 3(b), there exist significant local differences among them in shape, but all of them meet the definition of a loop globally as proposed here. Our motivation here is to explore the common global shape characteristic, but ignore the local variation. In particular, a loop can be classified into two categories, namely, concave loop and convex loop, in terms of the relationship between the normal vector at a point and the vector from this point to the center of the loop. As compared with vessels, which typically appear as high-density regions as compared with surrounding structures, the airway lumen regions have relatively low densities. In geometric shape, the airways are often regarded as having concave tubular shape and the vessels are often regarded as having convex tubular shape (Fig. 4).

2.D.2. Loop shape descriptor and its implementation

On the basis of the loop definition, we propose a loop shape descriptor in this study to characterize the shape of tubular structures and apply this descriptor to identify airway trees with concave tubular shape. Specifically, given a surface model (e.g., the example in Fig. 5), the loop shape descriptor for a point pi is defined as whether there exist a loop across

this point and whether this loop has limited size (e.g., the largest diameter smaller than a predefined size). This can be implemented by computing the intersection (a curve) between the plane perpendicular to the principal direction T i min and the surface model, and then checking whether the intersection curve is closed and meets the size constraint. The size constraint is necessary in practice, because a specific type of tubular biological structure generally has a limited size, namely, a size within a certain range. For example, despite its variations across individual subjects, the largest diameter of the airways (e.g., at the trachea) is less than 40 mm. Whereas the surface model in our study is formed by a number discrete triangle, the proposed loop descriptor can be implemented using a triangle-based region growing operation, which consists of the following steps: Step(1): given a point pi on a surface model and its principal curvatures, pi is pushed into a stack S. Step(2): pop an element qi from S. Step(3): check the neighbor vertices of qi that share the same triangle with qi and add the ones that have never been visited and meet a predefined constraint to stack S. As illustrated by the example in Fig. 5, the constraint is that the vertex has a distance to pi along the minimal principal direction of pi less than a small value d (e.g., d = 2 mm). Due to the discrete characteristic of the triangle mesh based surface representation, the band thickness d should be larger than the size of a triangle in the lung surface model but small enough

F. 3. An ideal cylindrical shape is shown in (a), and several looplike shapes are shown in (b). The shapes listed in (c) are not regarded as loops because their ends are not connected to their beginnings. The arrows indicate the normal directions, which could be used to differentiate convex surfaces [e.g., the examples in the top row in (b)] from concave surfaces [e.g., the examples in the bottom row in (b)]. Medical Physics, Vol. 42, No. 6, June 2015

3079

Pu et al.: A loop shape descriptor

3079

F. 4. Illustration of the airways with concave shape and the vessels with convex shape. The shadow regions indicate the high-density tissues. The arrows indicate the normal directions of the surfaces from high intensity to low intensity.

to assure that the resulting ringlike triangle mesh has a small thickness. The triangles generated by the MCA are formed by the vertices on the edges of the cubes defined by the eight neighboring voxels. Hence, the size of a triangle is similar to that of a cube (i.e., the size of a voxel). Step(4): repeat Step (2) and Step (3) until there is no element in the stack. The vertices that have been pushed into the stack will be recorded. Their combination will form a ringlike structure (e.g., the band with a height of 2d in Fig. 5). Step(5): determine whether the ring appears as a concave or convex surface. This can be accomplished by

testing the sign of the dot product between (c −r j ) and the normal vector n j of r i , where r i is a point on the ringlike structure and c is the centroid of the ringlike structure. If the sign of (c − r j ) · n j is negative, then the ringlike structure is convex; otherwise, the ringlike structure is concave (Fig. 4). Step(6): mark the points forming the ring as visited pointsof-interest if the identified ring is determined as expected tubular structure; otherwise, only the seed point is marked as visited but not points-ofinterest. Step(7): repeat Step (1)–Step (6) until all the points on the surface model are marked as visited.

F. 5. Illustration of the loop across a point p i . T i min is the minimum principal direction along which the curvature takes the minimum value, T i max is the maximum principal direction along which the curvature takes the maximum value, and d is band thickness (i.e., the thickness or height of a loop).

Given our specific aim, namely, the identification of an airway tree, additional operations can be taken to improve the efficiency of the above procedures. First, it is unnecessary to check every point on the lung surface model in Step (1) but only the points that meet two criteria: (a) the maximum curvature is larger than −0.05 (unit: 1/mm), corresponding to a radius of 20 mm and (b) the absolute value of the minimum curvature is larger than 0.2 (unit: 1/mm), corresponding to a radius of 5 mm. The first criterion excludes the planar regions and convex regions as region growing seeds. Second, whereas the maximum diameter of the airways (i.e., trachea) is typically less than 25 mm, step (4) can be terminated when there is a point on the ringlike structure that has a distance larger than 30 mm to pi . In order to compute the principal curvatures and the principal directions of each vertex on the triangle surface, we used the finite-difference approach proposed by Rusinkiewicz.39 As compared with previous approaches,39–41 this method offers reasonable accuracy with a linear computational complexity in space and time [i.e., O(N), where N is the number of vertices]. As demonstrated by the example in Fig. 6(a), the curvature-based filtering operation discarded most nontubular surfaces and/or convex surfaces (e.g., vessels). At

Medical Physics, Vol. 42, No. 6, June 2015

3080

Pu et al.: A loop shape descriptor

3080

F. 6. (a) shows the results after the application of the curvature filtering operation to the surface model obtained at an isovalue of −800 HU, (b) shows the airway tree identified at an isovalue of −800 HU after the application of the loop shape descriptor, where the vertices in (a) were used as the seeds, (c) shows the airway tree identified at an isovalue of −900 HU after the application of the loop shape descriptor, and [(d)–(f)] are the local enlargements of the regions indicated by the boxes in [(a)–(c)], respectively.

the same time, because of the second derivative characteristic of the curvature, a certain number of nonairway structures are incorrectly kept and some airway structures were incorrectly discarded. Regardless, the remaining vertices are used as the seeds for the loop shape descriptor. The results in Figs. 6(b) and 6(c) demonstrated that the loop shape descriptor picked up the regions that are located in a loop and discarded those that did not meet the definition of the loop in a reliable manner.

combined. The combination is a Boolean union operation that combines multiple sets into a single set. Each set is comprised of voxels forming the airways identified at a specific isovalue (threshold). The multiple isovalues are selected on the basis of the density spectrum of the airway walls on CT scans, which should be within [−700, −930 HU]. In our implementation, the multiple isovalues that we used are −700, −800, −900, and −930 HU.

2.E. The impact of the isovalue for surface modeling on airway identification

2.F. Performance assessment

When modeling the lung anatomical structures using the MCA, an isovalue or threshold is needed because the basic idea of the MCA is to treat neighboring voxels as the corners of a cube and identify the isosurface across these corners. As the example in Fig. 6(e) shows, there were some holes that should be part of the airways but are not identified as loops as expected. This issue is caused by the fact that the missed regions are fused with surrounding tissues at the given isovalue; hence, there are no closed paths curved over on themselves around these points. However, if the isovalue is set at a lower value (e.g., −900 HU), this fusion disappears and the airways at this region are successfully identified. Given the variety of CT scanning protocols and manufacturers in reality, which may significantly affect the image density and image noise/artifact levels, it is extremely difficult to adaptively determine a threshold that could achieve the optimal performance. Therefore, we proposed a multilevel thresholding solution, where multiple isovalues are used and the final segmentation results obtained at these isovalues are

We assessed the performance of the developed loop shape descriptor in airway identification using a publicly available dataset (http://image.diku.dk/exact/index.php),35 which was collected specifically for evaluating the performance of an airway tree segmentation scheme. This dataset consists of 20 chest CT examinations acquired under different protocols (e.g., dose, scanner, and reconstruction kernels). Fifteen different airway segmentation algorithms were evaluated on this dataset and their results were reported in Ref. 35. Performing an assessment on this dataset allows a direct comparison of the performance of the developed scheme with available approaches. A detailed description of these CT examinations and the experiments can be found elsewhere.35 The measures, including the airway length (excluding the trachea), the branch number, and the generation number, were computed. In this study, we computed the same measures under the help of an airway skeletonization scheme presented in Ref. 42. This skeletonization scheme could identify both the skeletons and the branching points of any 3D object, thereby making the

Medical Physics, Vol. 42, No. 6, June 2015

3081

Pu et al.: A loop shape descriptor

3081

F. 7. Illustration of the performance of the airway skeletonization scheme (Ref. 42).

computation of the proposed measures relatively easy. When computing the generations, we simply assume that the airways appear as a bifurcating tree. An example in Fig. 7 shows the performance of this skeletonization scheme. 3. EXPERIMENTS AND RESULTS After the application of the scheme to the above mentioned database, the three performance measures, namely, the total tree length, the generation number, and the branch number, were listed in Table I. Among the testing cases, the airway length ranged from 64.6 to 429.8 cm, the generation ranged from 7 to 11, and the branch number ranged from 48 to 312. In particular, identification of the airway tree in a typical CT examination consisting of ∼300 image slices takes approximately 8 min for a desktop PC (AMD Athlon™ 64 × 2 Dual 2.11 GHz central processor and 4 GB RAM). To enable a direct visual inspection of the performance, we listed the screen-shots of the identified airway trees in the Appendix, Fig. 10.

T I. Summary of airway segmentation performance levels for the 20 test cases in EXACT’09 (http://image.diku.dk/exact/).

Case

Total tree length (cm)

Number of generations

Number of branches

21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

114.6 264.8 215.6 210.9 270.1 92.3 100 125.7 157 177.9 220.8 234.6 163.5 340.6 339.1 364.7 212.1 64.6 360.7 429.8

9 10 11 8 9 5 7 8 9 9 9 11 10 11 11 11 8 7 11 10

78 173 138 118 152 48 55 76 101 114 125 125 100 190 230 201 105 41 214 312

Medical Physics, Vol. 42, No. 6, June 2015

4. DISCUSSION In this study, we proposed a novel scheme termed loop descriptor to identify airways on CT images and validated its performance using a publicly available database that was acquired under significant different protocols. The most distinctive characteristics of the loop shape descriptor are its global perspective of viewing a tubular structure and its geometric implementation. The loose definition of the tubular structure as proposed in this study makes it possible to consider the variety of tubular structures such as airways in a transparent manner. The geometric characteristic makes it straightforward and simple to identify the loop at a given point by computing the intersection between a plane and a 3D surface model, which can be performed using a triangle-based flooding operation (Sec. 2.D). Also, unlike many previous approaches,12–31 the developed method does not perform any path tracing procedure. In addition, this scheme is general in methodology and could be used to identify vascular structures. As a way to demonstrate the generalizability of this descriptor, we present some examples as shown in Fig. 8 after its application to the segmentation of cardiac vessels, lung vessels, and brain vessels. The proposed approach involved three parameters (Sec. 2.D), namely, the size constraint, the band thickness (i.e., 2d in Fig. 5) of the loop, and the isovalue for geometric surface modeling, which are independent of each other. Properly setting the size constraint could avoid the false positive detection. Only the regions with a dimension less than the predefined size will be identified. For example, we defined the size limit at 10 mm for pulmonary vascular tree segmentation [Fig. 8(a)]; thus, the aorta will not be detected as part of the lung vessels. Mathematically, a loop does not have thickness; however, given the discrete triangle-based representation of a surface model, the band thickness considers the size of a triangle and assures that a triangle-based loop could be formed by performing a triangle-based flooding operation. The thickness of a loop can be set at a relatively small value, but larger than the size of a triangle. In this study, we set this parameter conservatively at 3 mm because the voxel size is typically less than 1.5 mm. The third parameter, namely, the isovalue, was used for the MCA to model the meaningful structures embedded in the volumetric images. The selection of the isovalue depends on the density range of the structure-of-interest on images. Since the density of the airway walls varies significantly on CT, it is impossible to use a single isovalue to fully model the airway structures. As a solution, we proposed to use multiple isovalues to capture

3082

Pu et al.: A loop shape descriptor

3082

F. 8. Examples illustrating the performance of the loop shape descriptors in identifying different vascular structures, including lung vessels, cardiac arteries, and brain vessels.

the variation of the airways in density. In our experiment, we simply chose several thresholds on the basis of the density range of the airway walls (i.e., −930 to −700 HU). In theory, more isovalues will lead to the identification of more airways; however, this could result in more computational cost in time. When validating the performance of the developed shape descriptor in airway identification, we used a publicly available database (EXACT) consisting of 20 test cases acquired under different protocols. The variety of these cases in image quality and resolution not only makes this database a valuable source for validating an airway segmentation scheme but also enables a straightforward and relatively fair comparison with available approaches in performance. Previously, Lo et al.35 assessed the performance of the airway tree segmentation schemes developed by 15 teams on this dataset. As compared with these algorithms, our scheme demonstrated an overall better performance in terms of airway length, branch number, and generations as well as in terms of visual assessment (Appendix, Fig. 10).

Similar to our previously developed algorithm,32 this scheme could also be regarded as a computational geometry solution, including the utilization of the well-known MCA and the principal curvatures. However, there are distinctive differences in methodology (or concept) between the two methods. As a differential geometric approach, the previous method took advantage of both principal curvatures and principal directions in differentiating airways from other structures and at the same time used a “puzzle game” to discard nonairway regions and pick up airway regions. When there exists obvious image noise/artifact, as the example in Fig. 9(a) shows, the airway tree may be identified incompletely at a specific isovalue. In contrast, a novel shape descriptor termed loop shape descriptor was proposed in this study. Its loose definition of tubular structures enables a reliable identification of a relatively “complete” airway tree at a specific isovalue. As demonstrated by the examples in Fig. 9, the loop shape descriptor is capable of identifying the airways at the branching

F. 9. The airway trees identified at the isovalue of −850 HU using the loop shape descriptor (bottom row) and the differential geometry approach in Ref. 32 (top row), respectively. (B) and (E) show the local enlargement of the branching regions, and (C) and (F) show the local enlargement of the trachea. Medical Physics, Vol. 42, No. 6, June 2015

3083

Pu et al.: A loop shape descriptor

regions and those with high curvatures. In particular, the new method has a relatively higher performance than the previous one in both computational efficiency and small airway delineation. Even though, we are aware that additional efforts are needed to significantly improve the computational efficiency. Considering that the efficiency issue is primarily caused by the large number of triangles generated by the MCA and the utilization of multiple isovalues, potential optimization efforts could be taken by using a triangle decimation procedure and/or parallel programming techniques. It is notable that the triangle decimation procedure should be capable of maintaining the tiny details of the surfaces, because we are particularly interested in identifying small airways. Given the fact that the airways are identified at multiple isovalues, the multicore processor in modern computers could be leveraged to improve the efficiency by parallelly running the algorithm at different isovalues.

3083

5. CONCLUSION A novel loop shape descriptor was described in this study to identify airway from CT scans. We tested its feasibility and performance on a publicly available database that consists of 20 chest CT images. Our quantitative assessment demonstrated that this shape descriptor had a number of advantages, such as generalizability, simplicity, and insensitivity to the existence of image noise or artifacts. ACKNOWLEDGMENTS This work is supported in part by a Grant No. R01 HL096613 from National Institutes of Health, a Grant No. 201402013 from National Health and Family Planning Commission of the People’s Republic of China, and a Grant No. 2012KTCL03-07 from Shaanxi Science and Technology Innovation Program.

APPENDIX: THE SCREENSHOTS OF THE SEGMENTATION RESULTS AFTER THE APPLICATION OF THE “LOOP” SHAPE DESCRIPTOR TO THE EXACT’09 DATASET

F. 10. Segmentation results on the publically available 20 cases. Medical Physics, Vol. 42, No. 6, June 2015

3084

Pu et al.: A loop shape descriptor

a)Author

to whom correspondence should be addressed. Electronic mail: [email protected] 1Y.Nakano, S.Muro, H.Sakai, T.Hirai, K.Chin, M.Tsukino, K.Nishimura, H. Itoh, P. D. Paré, J. C. Hogg, and M. Mishima, “Computed tomographic measurementsofairwaydimensionsandemphysemainsmokers.Correlationwith lung function,” Am. J. Respir. Crit. Care Med. 162(3), 1102–1108 (2000). 2H. Arakawa, K. Fujimoto, Y. Fukushima, and Y. Kaji, “Thin-section CT imaging that correlates with pulmonary function tests in obstructive airway disease,” Eur. J. Radiol. 80(2), e157–e163 (2011). 3J. C. Hogg, F. Chu, S. Utokaparch, R. Woods, W. M. Elliott, L. Buzatu, R. M. Cherniack, R. M. Rogers, F. C. Sciurba, H. O. Coxson, and P. D. Paré, “The nature of small-airway obstruction in chronic obstructive pulmonary disease,” N. Engl. J. Med. 350(26), 2645–2653 (2004). 4W. J. Kim, E. K. Silverman, E. Hoffman, G. J. Criner, Z. Mosenifar, F. C. Sciurba, B. J. Make, V. Carey, R. S. Estepar, A. Diaz, J. J. Reilly, F. J. Martinez, and G. R. Washko, “CT metrics of airway disease and emphysema in severe COPD,” Chest 136(2), 396–404 (2009). 5A. L. James and S. Wenzel, “Clinical relevance of airway remodelling in airway diseases,” Eur. Respir. J. 30(1), 134–155 (2007). 6D. O. Wilson, J. L. Weissfeld, A. Balkan, J. G. Schragin, C. R. Fuhrman, S. N. Fisher, J. Wilson, J. K. Leader, J. M. Siegfried, S. D. Shapiro, and F. C. Sciurba, “Association of radiographic emphysema and airflow obstruction with lung cancer,” Am. J. Respir. Crit. Care Med. 178(7), 738–744 (2008). 7M. Sonka, W. Park, and E. A. Hoffman, “Rule-based detection of intrathoracic airway trees,” IEEE Trans. Med. Imaging 15(3), 314–326 (1996). 8R. Millioni, S. Sbrignadello, A. Tura, E. Iori, E. Murphy, and P. Tessari, “The inter- and intra-operator variability in manual spot segmentation and its effect on spot quantitation in two-dimensional electrophoresis analysis,” Electrophoresis 31(10), 1739–1742 (2010). 9H. Fujita, Y. Uchiyama, T. Nakagawa, D. Fukuoka, Y. Hatanaka, T. Hara, G. N. Lee, Y. Hayashi, Y. Ikedo, X. Gao, and X. Zhou, “Computer-aided diagnosis: The emerging of three CAD systems induced by Japanese health care needs,” Comput. Methods Programs Biomed. 92(3), 238–248 (2008). 10X. Artaechevarria, D. Pérez-Martín, M. Ceresa, G. de Biurrun, D. Blanco, L. M. Montuenga, B. van Ginneken, C. Ortiz-de-Solorzano, and A. MuñozBarrutia, “Airway segmentation and analysis for the study of mouse models of lung disease using micro-CT,” Phys. Med. Biol. 54(22), 7009–7024 (2009). 11J. Rosell and P. Cabras, “A three-stage method for the 3D reconstruction of the tracheobronchial tree from CT scans,” Comput. Med. Imaging Graphics 37(7-8), 430–437 (2013). 12D. Aykac, E. A. Hoffman, G. McLennan, and J. M. Reinhardt, “Segmentation and analysis of the human airway tree from three-dimensional x-ray CT images,” IEEE Trans. Med. Imaging 22(8), 940–950 (2003). 13A. P. Kiraly, W. E. Higgins, G. McLennan, E. A. Hoffman, and J. M. Reinhardt, “Three-dimensional human airway segmentation methods for clinical virtual bronchoscopy,” Acad. Radiol. 9, 1153–1168 (2002). 14A. Fabija´ nska, “Two-pass region growing algorithm for segmenting airway tree from MDCT chest scans,” Comput. Med. Imaging Graphics 33(7), 537–546 (2009). 15C. I. Fetita, F. Prêteux, C. Beigelman-Aubry, and P. Grenier, “Pulmonary airways: 3-D reconstruction from multislice CT and clinical investigation,” IEEE Trans. Med. Imaging 23(11), 1353–1364 (2004). 16W. Park, E. A. Hoffman, and M. Sonka, “Segmentation of intrathoracic airway trees: A fuzzy logic approach,” IEEE Trans. Med. Imaging 17(4), 489–497 (1998). 17L. Fan and C. W. Chen, “Reconstruction of airway tree based on topology and morphological operations,” Proc. SPIE 3978, 46–57 (2000). 18T. Kitasaka, K. Mori, J. Hasegawa, and J. Toriwaki, “A method for extraction of bronchus regions from 3D chest x-ray CT images by analyzing structural features of the bronchus,” in Medical Image Computing and ComputerAssisted Intervention - MICCAI 2003, Lecture Notes in Computer Science Vol. 2879 (Springer Berlin Heidelberg, Montréal, Canada, 2003), pp. 603–610. 19D. Bartz, D. Mayer, J. Fischer, S. Ley, A. del Rio, S. Thust, C. P. Heussel, H. U. Kauczor, and W. Strasser, “Hybrid segmentation and exploration of the human lungs,” in IEEE Visualization (IEEE, Washington, DC, 2003), pp. 177–184. 20D. Mayer, D. Bartz, J. Fischer, S. Ley, A. del Rio, S. Thust, H. U. Kauczor, and C. P. Heussel, “Hybrid segmentation and virtual bronchoscopy based on CT images,” Acad. Radiol. 11(5), 551–565 (2004). 21J. Tschirren, E. A. Hoffman, G. McLennan, and M. Sonka, “Intrathoracic airway trees: Segmentation and airway morphology analysis from low-dose Medical Physics, Vol. 42, No. 6, June 2015

3084 CT scans,” IEEE Trans. Med. Imaging 24(12), 1529–1539 (2005). Bülow, C. Lorenz, and S. Renisch, “A general framework for tree segmentation and reconstruction from medical volume data,” in Proceedings of the Medical Image Computing and Computer-Assisted Intervention— MICCAI 2004, Vol. 3216 (Springer, Saint-Malo, France, 2004), pp. 533–540. 23W. B. van Ginneken and E. M. van Rikxoort, “Robust segmentation and anatomical labeling of the airway tree from thoracic CT scans,” in Medical Image Computing and Computer-Assisted Intervention – MICCAI 2008 (Springer Berlin Heidelberg, New York, NY, 2008), pp. 219–226. 24P. Lo and M. de Bruijne, “Voxel classification based airway tree segmentation,” Proc. SPIE 6914, 69141K (2008). 25P. Lo, J. Sporring, J. J. Pedersen, and M. de Bruijne, “Airway tree extraction with locally optimal paths,” Med. Image Comput. Comput.-Assisted Intervention 12(2), 51–58 (2009). 26P. Lo, J. Sporring, H. Ashraf, J. J. Pedersen, and M. de Bruijne, “Vesselguided airway tree segmentation: A voxel classification approach,” Med. Image Anal. 14(4), 527–538 (2010). 27M. W. Graham, J. D. Gibbs, D. C. Cornish, and W. E. Higgins, “Robust 3-D airway tree segmentation for image-guided peripheral bronchoscopy,” IEEE Trans. Med. Imaging 29(4), 982–997 (2010). 28Q. Li, S. Sone, and K. Doi, “Selective enhancement filters for nodules, vessels, and airway walls in two- and three-dimensional CT scans,” Med. Phys. 30(8), 2040–2051 (2003). 29Y. Sato, S. Nakajima, N. Shiraga, H. Atsumi, S. Yoshida, T. Koller, G. Gerig, and R. Kikinis, “Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images,” Med. Image Anal. 2(2), 143–168 (1998). 30K. Krissiana, G. Malandaina, N. Ayachea, R. Vaillantb, and Y. Troussetb, “Model-based detection of tubular structures in 3D images,” Comput. Vis. Image Understanding 80, 130–171 (2000). 31C. Bauer and H. A. Bischof, “Novel approach for detection of tubular objects and its application to medical image analysis,” in Pattern Recognition, Lecture Notes in Computer Science Vol. 5096 (Springer Berlin Heidelberg, Munich, Germany, 2008), pp. 163–172. 32J. Pu, C. Fuhrman, W. F. Good, F. C. Sciurba, and D. Gur, “A differential geometric approach to automated segmentation of human airway tree,” IEEE Trans. Med. Imaging 30(2), 266–278 (2011). 33J. Pu, S. Gu, S. Liu, S. Zhu, D. Wilson, J. M. Siegfried, and D. Gur, “CT based computerized identification and analysis of human airways: A review,” Med. Phys. 39(5), 2603–2616 (2012). 34J. Petersen, M. Nielsen, P. Lo, Z. Saghir, A. Dirksen, and M. de Bruijne, “Optimal graph based segmentation using flow lines with application to airway wall segmentation,” Inf. Process. Med. Imaging 22, 49–60 (2011). 35P. Lo, B. van Ginneken, J. M. Reinhardt, T. Yavarna, P. A. de Jong, B. Irving, C. Fetita, M. Ortner, R. Pinho, J. Sijbers, M. Feuerstein, A. Fabija´nska, C. Bauer, R. Beichel, C. S. Mendoza, R. Wiemker, J. Lee, A. P. Reeves, S. Born, O. Weinheimer, E. M. van Rikxoort, J. Tschirren, K. Mori, B. Odry, D. P. Naidich, I. Hartmann, E. A. Hoffman, M. Prokop, J. H. Pedersen, and M. de Bruijne, “Extraction of airways from CT (EXACT’09),” IEEE Trans. Med. Imaging 31(11), 2093–2107 (2012). 36W. E. Lorensen and H. E. Cline, “Marching cubes: A high resolution threedimensional surface construction algorithm,” ACM SIGGRAPH Comput. Graphics 21(4), 163–169 (1987). 37L. Custodioa, T. Etieneb, S. Pescoa, and C. Silvac, “Practical considerations on marching cubes 33 topological correctness,” Comput. Graphics 37(7), 840–850 (2013). 38X. Chen and F. Schmitt, “Intrinsic surface properties from sur-face triangulation,” in Proceedings of European Conference on Computer Vision (Springer-Verlag, London, UK, 1992), pp. 739–743. 39S. Rusinkiewicz, “Estimating curvatures and their derivatives on triangle meshes,” in Proceedings of 3D Data Processing, Visualization, and Transmission (3DPVT) (IEEE, Washington, DC, 2004), pp. 486–493. 40J. Goldfeather and V. Interrante, “A novel cubic-order algorithm for approximating principal direction vectors,” ACM Trans. Graphics 23(1), 45–63 (2004). 41M. Meyer, M. Desbrun, P. Schroder, and A. H. Barr, “Discrete differential geometry operators for triangulated 2-manifolds,” in Visualization and Mathematics III, edited by H. C. Hege and K. Polthier (Springer-Verlag, Heidelberg, Germany, 2003), pp. 35–57. 42N. D. Cornea, D. Silver, X. Yuan, and R. Balasubramanian, “Computing hierarchical curve-skeletons of 3D objects,” Visual Comput. 21(11), 945–955 (2005). 22T.