936

IEEE JOURNAL OF BIOMEDICAL AND HEALTH INFORMATICS, VOL. 17, NO. 5, SEPTEMBER 2013

Segmentation of the Thoracic Aorta in Noncontrast Cardiac CT Images Olga C. Avila-Montes, Uday Kurkure, Member, IEEE, Ryo Nakazato, Daniel S. Berman, Damini Dey, and Ioannis A. Kakadiaris, Senior Member, IEEE

Abstract—Studies have shown that aortic calcification is associated with cardiovascular disease. In this study, a method for localization, centerline extraction, and segmentation of the thoracic aorta in noncontrast cardiac-computed tomography (CT) images, toward the detection of aortic calcification, is presented. The localization of the right coronary artery ostium slice is formulated as a regression problem whose input variables are obtained from simple intensity features computed from a pyramid representation of the slice. The localization, centerline extraction, and segmentation of the aorta are formulated as optimal path detection problems. Dynamic programming is applied in the Hough space for localizing key center points in the aorta which guide the centerline tracing using a fast marching-based minimal path extraction framework. The input volume is then resampled into a stack of 2-D cross-sectional planes orthogonal to the obtained centerline. Dynamic programming is again applied for the segmentation of the aorta in each slice of the resampled volume. The obtained segmentation is finally mapped back to its original volume space. The performance of the proposed method was assessed on cardiac noncontrast CT scans and promising results were obtained. Index Terms—Aorta, image segmentation, noncontrast computed tomography, regression, slice localization.

I. INTRODUCTION ORONARY heart disease is the primary cause of morbidity and mortality in the United States and around the globe. According to the National Heart, Lung and Blood Institute, the cardiovascular diseases, including atherosclerosis, were responsible for about one million deaths in 2008 in the United States [1]. Studies have found that the coronary artery calcium is a distinct marker for the presence of atherosclerosis [2]. Such calcifications can also occur in extra-coronary structures including the ascending aorta, the aortic arch, the descending aorta, and the abdominal aorta. Several studies have shown that the

C

Manuscript received July 17, 2012; revised December 6, 2012; accepted February 24, 2013. Date of publication June 18, 2013; date of current version August 30, 2013. This work was supported in part by the National Institutes of Health under Grant 1R21EB006829-01A2 and in part by UH Eckhard Pfeiffer Endowment Fund. Any opinions, findings, conclusions or recommendations expressed in this material are of the authors and may not reflect the views of the NIH or UH. O. C. Avila-Montes, U. Kurkure, and I. A. Kakadiaris are with the Computational Biomedicine Lab, Deptartment of Computer Science, University of Houston, Houston, TX 77035 USA (e-mail: [email protected]; [email protected]; [email protected]). R. Nakazato, D. S. Berman, and D. Dey are with the Department of Imaging, Cedars-Sinai Medical Center, Los Angeles, CA 90048 USA (e-mail: [email protected]; [email protected]; [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/JBHI.2013.2269292

calcification of the thoracic aorta, aortic arch, and aortic valve is related to an increased probability of being at risk of cardiovascular disease [3]–[6]. Recent research has also shown a close relation between thoracic aorta calcium, aortic valve calcium, and advanced coronary artery calcium [7]–[10] that suggest the presence of systemic atherosclerosis. Prior to including aortic calcium measurements in cardiovascular risk assessment, longitudinal studies need to be conducted to assess the incremental value of thoracic aorta calcium in addition to coronary artery calcium. The aorta and coronary artery calcifications can be quantified by noncontrast computed tomography (CT) imaging [11]. In addition, the thoracic aorta calcifications can be measured from standard cardiac calcium scoring CT scans, without requiring any additional scanning, which is especially advantageous for retrospective studies. Currently, aortic calcium scoring is not routinely performed in clinical practice because the available calcium scoring tools require manual annotation of the calcifications, which is a laborintensive and time consuming process. One way to automate the process of aortic calcium scoring is to determine the aortic boundaries and extract calcifications within these boundaries automatically or with minimal user interaction. In this paper, we present computational methods that are developed for the automated segmentation of the thoracic aorta, which is comprised of the ascending aorta, aortic arch, and descending aorta in noncontrast CT imaging data. The proposed method consists of the following steps (see Fig. 1): 1) A robust regression model to predict the right coronary artery (RCA) ostium slice in a cardiac CT volume is constructed by cascading multiple regressors, where each regression function refines the estimates of the previous one. Given a new volume, the RCA ostium slice location is predicted by the constructed cascaded regressor [12] [see Fig. 1(a)]. 2) The RCA ostium slice is used to obtain regions of interest (ROIs) for the ascending and descending aorta, and key waypoints of the aorta are determined using a combination of Hough transform and dynamic programming [13] [see Fig. 1(b)]. 3) The points obtained in the previous step are used to guide the extraction of an approximated aorta centerline by employing a fast marching-based minimal path extraction algorithm. The centerline is then traversed in order to construct a transformed volume composed of 2-D crosssectional planes orthogonal to the centerline [see Fig. 1(c)]. 4) The aorta boundary is determined using dynamic programming in each slice of the resampled volume [14]. Finally, the resulting segmentation is transformed back to the CT volume space [see Fig. 1(d)]. The parts of this study have appeared in [12]–[14]. In this paper, we present the complete framework to segment the whole aorta, i.e., the ascending aorta, the descending

2168-2194 © 2013 IEEE

AVILA-MONTES et al.: SEGMENTATION OF THE THORACIC AORTA IN NONCONTRAST CARDIAC CT IMAGES

Fig. 1. Overview of the proposed thoracic aorta segmentation method: (a) right coronary artery ostium slice localization [see Section III-A, for a 2-D visualization of this slice, please refer to Fig. 3(a)], (b) detection of key ROIs and aorta center-points (see Section III-B1), (c) aorta centerline and crosssectional plane extraction (see Section III-B2), and (d) resulting thoracic aorta segmentation (see Section III-B3). 3-D isosurface of the manually annotated aorta (in yellow) added for reference.

aorta, and the aortic arch, from noncontrast cardiac CT scans. We demonstrate the performance of the proposed framework through extensive experiments on a large dataset. Our contributions are the following: 1) we introduce a novel framework for estimating the position of a slice of interest within a CT scan. Unlike other methods that require a similarity function to retrieve a slice of interest, we used regression techniques to directly estimate the position of the target slice without the need for a similarity criterion. 2) We present a robust and unified method for the automated segmentation of the thoracic aorta comprising the ascending aorta, the aortic arch, and the descending aorta in the noncontrast CT data. 3) We propose an entropy-based cost function for aortic boundary detection in order to overcome the effect of lack of strong contrast and noise for which a gradient-based method may fail. The rest of the paper is organized as follows. In Section II, related study is reviewed. In Section. III, a detailed description of the proposed methods is provided. The performance of the method is evaluated in Section IV by validating the results on clinical data. Section V provides the discussion, and finally Section VI concludes the paper. II. RELATED WORK Considerable research has been conducted toward thoracic aorta segmentation in magnetic resonance (MR), magnetic resonance angiography (MRA), and computed tomography angiography (CTA) images. Behrens et al. [15] proposed an aorta segmentation method for MRA images in which the tubular structures are modeled as generalized cylinders. The aorta is detected and tracked using an extended randomized Hough transform and a Kalman filtering method, respectively. However, the method provides only a coarse segmentation and additional processing will be required to obtained accurate aorta boundaries. In addition, it requires a starting point, an approximate direction, and an approximate radius as input. Adame et al. [16] proposed

937

to detect aorta using the Hough transform and obtain boundary contours using a minimal cost approach in MR images. Their method was applied only for the descending aorta segmentation and it relies on strong gradient information present in the MR images to detect the boundaries. Worz and Rohr [17] proposed a model-based vessel segmentation method that fits a 3-D cylindrical parametric intensity model using a Kalman filter in MRA and CTA images. However, intensity profile-based methods have difficulties in segmenting vessels in regions characterized by low contrast and high noise. In addition, their method requires a starting point and an approximate orientation as initialization. The deformable model-based approaches rely on high contrast edges to determine the boundaries. Rueckert et al. [18] proposed to track aorta using a medialness function and the final segmentation is obtained through an energy minimizing deformable model defined on a Markov-random field framework. Lorigo et al. [19] proposed the curves method to segment blood vessels in MRA images by extending geodesic active contours and minimal surfaces methods to model the vessels as 3-D curves. Kovacs et al. [20] proposed to segment thoracic aorta in CTA images using an elastic mass-spring deformable model to refine the initially detected contour through Hough transform. Most of the methods rely on the presence of high contrast between lumen and its surroundings in the MR, MRA, and CTA images to obtain sufficiently large intensity gradients at the vessel boundary only. The lack of contrast between blood pool regions, muscle walls and fat in the noncontrast CT imaging renders the aorta segmentation as a challenging task. In addition, the methods employing Hough transform detect the circle of interest by selecting the point of global maxima in the Hough space. The point of global maxima does not always correspond to the aorta in noncontrast CT data due to the lack of contrast and the presence of noise. In model-based approaches, Zheng et al. [21] presented a part-based aorta segmentation method for intraoperative C-arm CT images. Marginal Space Learning is employed to train a detector for the aortic root and arch, while the ascending and descending aorta are tracked using a 2-D circle detector. The aorta boundary is then refined iteratively using a learning-based 3-D boundary detector. Their method was evaluated only on the contrast-enhanced data and the scans with poor contrast were not included in quantitative evaluation, though they claim that their method can segment the aorta in those scans too. Escabert et al. [22] presented a model-based method for segmenting the heart and the attached great vessels in cardiac CT angiography images. The model is comprised of the heart chambers and vessel trunks, and the great vessels are represented by the concatenation of short tubular segments. The aorta is segmented by using parametric and deformable mesh adaptation techniques to match the tubular segments to the image. The method was developed and evaluated on contrast-enhanced data. Their method performed well for aorta segmentation because the aorta was well contrasted. However, it performed poorly on other vessels that were poorly contrasted. An alternative approach for aorta segmentation is through registration. Isgum et al. [23] applied a novel atlas-based segmentation technique to 29 cases of no-contrast CT thoracic scans.

938

IEEE JOURNAL OF BIOMEDICAL AND HEALTH INFORMATICS, VOL. 17, NO. 5, SEPTEMBER 2013

The final segmentation is obtained by fusing the results of registration of multiple atlases to a target image. A novelty in their method is its capability of local evaluation of registration success, and therefore directing decision fusion of the propagated labels during registration. They also apply sequential forward selection, from statistical pattern recognition theory to select the optimal subset of atlas images. Their method overcomes the problem of bias in a single atlas-based methods. However, the presence of low contrast boundaries pose a problem in the detection of accurate boundaries in their method which is common to all atlas-based segmentation methods. Also, the low contrast problem is magnified in the lower resolution images (e.g., calcium scoring CT scans) due to the blurring of boundaries. This poses a significant challenge in detecting the aortic arch boundaries because the resolution of the cross-sectional plane of the aortic arch is up to five times lower than the resolution of the cross-sectional plane of the ascending and descending aorta. Recently, Taeprasartsit and Higgins [24] proposed a method that employs centerline models of the aorta, model fitting, and reverse Euclidean distance transform to obtain the aorta segmentation. Their method was tested for both the contrast enhanced and noncontrast multi-detector CT scans by detecting the presence of contrast and using different parameter settings. A combination of centerline models were used to obtain a point at the beginning of the ascending aorta, which is used to initialize the tracing of the aorta centerline. If the automated centerline extraction fails, the authors proposed a semiautomated method that requires the user to input three landmarks to guide the centerline extraction at the ascending aorta, the descending aorta, and the aortic arch. However, it is not clearly reported in how many cases the automated algorithm failed, requiring the input of three anatomical landmarks for the successful segmentation of the aorta using their semiautomated method.

predicted for the ith slice. To obtain R1 , the label yi is defined by E(τα , τi ), where E(·) is a measure of the relative error, which is the difference between the location of the target slice τα and the initial location of the ith slice. In subsequent levels, l > 1, the label yi is defined by E(τα , Rl−1 (xi )), where Rl−1 (xi ) is the resulting prediction from the previous regressor Rl−1 . The process is repeated until either the regressor Rl is incapable of reducing the prediction error compared to the previous level or when a predefined number of L levels is reached. 2) Testing Phase: In the testing phase, the input images undergo the same process of ROI detection and alignment. Then, the cropped images are mapped into the low-dimensional subspace which was learned in the training stage. The cascaded regressor is evaluated at each slice by computing τ L = τ l−1 + Rl (xτ l −1 ) from l = 1, . . . , L, where xτ l −1 is the feature vector computed from the slice at position τ l−1 , and Rl (xτ l −1 ) is the resulting prediction from the regressor at level l. Finally, the location of the target slice is estimated by computing the probability density function of the resulting predictions and selecting the location corresponding to the maximum density. Fig. 1(a) depicts a 3-D isosurface visualization of the manually annotated aorta and the RCA ostium slice. B. Aorta Segmentation In this section, the methods for the detection, centerline extraction, and segmentation of the aorta are presented. The outline of the method is presented in Algorithm 1.

III. METHODS A. Slice Localization In this section, we present a regression-based method to predict the location of a slice of interest (in this case, the RCA ostium slice) in a given noncontrast CT volume. 1) Learning Phase: In the training stage, with 2-D cardiac CT slices as the input to the method, a region of interest (ROI) is detected in order to include only relevant information. First, the inner thoracic cavity is extracted using a graph-based method [25]. Then, the slices are properly rotated to be aligned and the ROI is further refined by thresholding and connectedcomponent analysis. In order to capture all the information from the ROI and also maintain the spatial ordering of the features with some tolerance to intersubject variations, a spatial pyramid representation [26] with three levels is used to aggregate local statistical features (mean, standard deviation, skewness, kurtosis, and entropy) as a slice descriptor. Then, a low-dimensional data representation of the location manifold is learned, so that it is sufficiently discriminative. Finally, a robust cascaded regressor R = (R1 , . . . , RL ) is constructed for learning the embedding function of the dataset D = {(x1 , y1 ), . . . , (xD , yD )}, where D is the total number of slices in the dataset, xi is the extracted feature vector for the ith slice, and yi is the label to be

1) Detection of the Aorta Location: In spite of heart dynamics, the aorta maintains its global tubular shape with minor local deformations. Since the ascending and descending aorta run mainly vertically, their appearance in axial slices approximates a circular shape which is extracted using the Hough transform. A dynamic programming-based optimal path search is performed on the Hough accumulators of subsequent axial slices to obtain a series of optimal best-fit circles for the aorta. Fig. 1(b) depicts the detected aorta key points (black spheres) with a 3-D isosurface of the manually annotated aorta added for reference. The boxes depict the ROIs in red, green, and blue for the ascending aorta, descending aorta and lower descending aorta (LDA), respectively. Note that that we could have used a single ROI for the descending aorta that includes both of its ROIs. However, the computation of the Hough space for such a large ROI is computationally expensive and memory intensive.

AVILA-MONTES et al.: SEGMENTATION OF THE THORACIC AORTA IN NONCONTRAST CARDIAC CT IMAGES

We found through experimentation that we can instead use two small ROIs at the start and end locations of the descending aorta without compromising the results because the descending aorta runs almost vertically across the axial plane. a) Computation of hough accumulator: First, in order to accentuate the tissue of interest, such as blood and muscles, a nonlinear gray-scale modification [27] is performed using ˆ y) = M (I(x, y)), where piecewise linear mapping models, I(x, ˆ y) are the input pixel value and the output I(x, y) and I(x, pixel value at coordinates (x, y), respectively, and M (·) is the mapping equation ⎧ for I(x, y) < −30 ⎪ ⎨ −30 M (I(x, y)) = 100 for I(x, y) > 100 . (1) ⎪ ⎩ I(x, y) otherwise. Then, anisotropic diffusion [28] is used to smooth the ROIs. Anisotropic diffusion has the advantage of preserving the edges while smoothing regions and is governed by: ∂ Iˆ ˆ ˆ ∂ t = div(g(∇I)∇I), where div is the divergence operator, ˆ is the gradient magnitude, g(·) is the so-called diffusion ∇I coefficient, and t denotes the evolution of time. We used the diffusion function that favors wide regions over smaller regions, and is defined by ⎛  2 ⎞−1 ∇Iˆ ˆ = ⎝1 + ⎠ g(∇I) κ where κ controls the conduction as a function of gradient. Finally, the Hough accumulator, of size U × V (U and V are the x- and y-dimensions of the input image, respectively), is computed by the convolution operation H = E ∗ k, where H is the Hough accumulator, E is the Canny [29] edge response, and k is a circular kernel. The most probable circle corresponding to the aorta appears as a maxima in the Hough accumulator. However, in practice, the global maximum point does not always correspond to the aorta due to the presence of other structures similar to the aorta. To overcome this problem, the Hough accumulator is treated as a medialness feature space (center coordinates and radius are the dimensions of the space) and a dynamic programming-based method is used to select the correct Hough circle. b) Localization using dynamic programming: The goal is to select a single point in the Hough accumulator of each axial image that corresponds to the aorta. In terms of dynamic programming, the problem is reformulated as determining an optimal combination of points from the Hough accumulators of each axial image that corresponds to the centerline of the aorta. This is achieved by minimizing a cost function. If there are S axial images in consideration, the solution of dynamic programming is a set PS of S points (p1 , p2 , . . . , pS ) ∈ PS , where each point pi is represented by its x-coordinate, y-coordinate, radius, Hough value, and z-coordinate (x, y, r, v, z). The set PˆS is the optimal combination of points PS that has the minimum cost for the cost function C(PS ), given by: C(p1 , p2 , . . . , pS ) =

S The cost associated with each point is expressed i=1 cH (pi ). as: cH (pi ) = K j =1 hj (pi )(i = 1, . . . , S), where the cost terms hj (pi )(j = 1, . . . , K) characterize the relationship of pi with

939

its adjacent points, and K corresponds to the number of characterizations used. The following features are used to construct the feature-based cost function. Horizontal offset: h1 (pi ) = (xi − xi−1 )2 + (yi − yi−1 )2 . This term is associated with the change in the horizontal position of adjacent circles. It imposes 3-D verticality of the medial axis of the aorta by constraining the amount of horizontal shift allowed for aortic circles in adjacent axial slices. Smoothness: h2 (pi ) = |ri − ri−1 |. This term is associated with the change in radius of adjacent circles. It imposes a smooth tubular shape by restricting the amount of change in size of aortic circles allowed in adjacent axial slices. Hough values: h3 (pi ) = (1 + vi )−1 . This term is associated with the number of votes in the Hough accumulator. It forces the selection of circle points which have the maximum edge response. A dynamic programming-based approach, such as the one described previously, allows us to avoid an exhaustive search for the optimal combination of circle points, and exploits the 3-D tubular shape characteristic of the aorta. To further reduce the computational cost, the search is limited to Q points in each Hough accumulator which have the maximum Hough value. This reduces the search space from U × V × W × S points to Q × S points. 2) Detection of the Aorta Centerline: Anatomically, the ascending aorta, the aortic arch, and the descending aorta resemble a candy cane shape. Therefore, the goal of this step is to obtain a transformed 3-D volume, where the aorta is observed as a straightened tubular structure. Specifically, the steps for centerline extraction are the following. a) Construction of the likelihood image: The speed image used to detect the centerline of the aorta should be able to describe the likelihood of each voxel being closer or further away from the aorta centerline. To this end, the image is binarized by defining a threshold of τ HU to eliminate all voxels belonging to other tissues such as fat and air. Then, in order to remove small holes, a closing operation is performed using a ball structuring element with radius ρ. In addition, anatomically, the aorta extends upward from the base of the left ventricle of the heart, and then, arches like a candy cane running downward, lying on the anterior surface of the vertebral column [30]. This prior information about the aorta shape is very valuable and can be used to guide the centerline extraction. To exploit this prior information, one approach can be to use a mean shape of the aortic arch to guide the minimal path computation. However, we observed that the aortic arch exhibits considerable shape variability and therefore, the mean shape may not be sufficient to account for such variations. Another approach can be to perform image registration prior to mean shape computation but that increases the computational cost significantly. In addition, if the mean shape is far away from the real aortic arch, it can misguide the minimal path computation. Therefore, we opted to include the shape prior by removing an ellipsoidal region that conforms to the shape of aorta and let the data drive the computation of the minimal path without enforcing any artificial shape constraints. Toward that, we define an ellipsoidal region [see Fig. 2(a)] that masks the region below the aortic arch and between the ascending and descending aorta to steer the centerline toward the arch and avoid any leakage

940

IEEE JOURNAL OF BIOMEDICAL AND HEALTH INFORMATICS, VOL. 17, NO. 5, SEPTEMBER 2013

Fig. 2. Depiction of the aorta centerline extraction: (a) 3-D isosurface visualization of the ellipsoid, and (b) extracted aorta centerline (black line) with spheres that depict the center-points used for guidance. 3-D isosurface of the manually annotated aorta added for reference.

through the heart chambers. The ellipsoid is positioned at the midpoint of the line segment connecting the detected ascending and the descending aorta center points at the RCA ostium slice. This line segment also defines the z-axis of the ellipsoid. The x-axis of the ellipsoid is computed as the perpendicular to the z-axis in the axial slice plane. Finally the y-axis of the ellipsoid is computed as the orthogonal axis to the x- and z-axes of the ellipsoid. The size parameters of the ellipsoid are determined based on a learned model relating the distance between the center points of the ascending and the descending aorta, and the distance between the center points of the aortic arch and the RCA ostium slice. The values of all the voxels belonging to the obtained ellipsoid were set to zero in the binary image. Then, the Euclidean distance map is computed for the remaining foreground voxels using the algorithm developed by Danielsson [31]. In the distance map image, each voxel indicates its distance to the nearest background voxel. All the voxels that have distance values greater than γ mm are removed from the distance map. This threshold is applied to the likelihood image to avoid the leakage of the extracted centerline into other organs/tissues which have a similar intensity range as that of the aorta. Another option is to replace the threshold with a decreasing function of distance beyond a fixed distance value. However, its usefulness in avoiding the leakage is questionable and has to be evaluated through further experimentation. The resulting likelihood image is finally scaled from sl to sh . The process of obtaining the likelihood image and the inclusion of the aorta shape as a prior is depicted in Fig. 3. b) Centerline extraction: To obtain the approximated centerline of the aorta, the fast marching-based minimal path extraction framework proposed by Mueller et al. [32] was used. Given the start and end points, the minimal path is extracted in the following manner: 1) a fast marching arrival function is computed by expanding a wave front along a speed function from the end point until the front reaches the start point and 2) the minimal path is traced by back propagating from the start point toward the end point. In this study, the start point is the center point of the ascending aorta at the RCA ostium slice, while the end point is the center point of the descending aorta at the last slice of the scan. We include all ROIs in the final extraction of the centerline because the previously detected center points in those ROIs can be used as waypoints to guide the full centerline extraction. This way, we can obtain a smooth continuous

Fig. 3. Depiction of the construction of the likelihood image: (a) original CT image, (b) threshold mask of tissue of interest after closing operation, (c) mask after inclusion of ellipsoid, and (d) Euclidean distance map with valid distance values (within the aorta radius range).

centerline, and also correct for the errors, if any, in the previous step. An advantage of using this fast marching-based method is that it also allows for the incorporation of waypoints. The waypoints are the points that are located on/near the centerline between the start and end points of the centerline. These waypoints guide the minimal path to pass through the centerline/medial axis by steering the path toward themselves. The modified method of extracting the path using the waypoints is as follows: 1) the front is propagated from the first waypoint until it reaches the start point; and 2) the minimal path is traced from the start point toward the first waypoint. Then, the process is repeated for each remaining waypoint, with the previous waypoints acting as the start points, until the path reaches the end point. The resultant path is formed by the concatenation of each path segment. In Mueller’s method, a parameter called termination value is used to determine how close the path will get to each of the guide points. While tracing the path from one point to another, the tracing stops when the path has reached an arrival value less than the termination value. This parameter affects both the path segment reaching a waypoint, and the path segment reaching the end point. In the case of the aorta centerline extraction, this behavior is disadvantageous because, ideally, the waypoints should serve as soft guidance points so the path will tend toward them but may not reach them (i.e., using a large termination value), but at the same time, the path should terminate very close to the end point (i.e., using a small termination value). Toward this goal, we modified the method and added a new parameter called end termination value, so that the termination value now regulates how close the path will get to the waypoints, and the end termination value regulates how close the path will get to the end point.

AVILA-MONTES et al.: SEGMENTATION OF THE THORACIC AORTA IN NONCONTRAST CARDIAC CT IMAGES

The speed image obtained as described in Section III-B2a, and the guidance points belonging to the ascending, descending, and LDA-ROIs computed as described in Section III-B1, were used in order to guide the minimal path extraction. The back-propagation is achieved by using the gradient descent optimizer. Fig. 2(b) depicts the detected aorta center points and the approximated aorta centerline obtained by the fast marching minimal path extraction method. c) Volume transformation: To facilitate the aorta segmentation, the original volume is resampled into a stack of 2-D slices where the aorta resembles the shape of a circle. In other words, the arched structure of the aorta is straightened. In order to obtain multiple centered planes orthogonal to the centerline path [see Fig. 1(c)], the following needs to be provided for each plane: 1) a center point and 2) a local coordinate system to orient the cross-sectional plane. The approximated aorta centerline extracted, as described in Section III-B2b, is resampled to obtain equally spaced points along the curve. Each resampled point is then used as the center point for a crosssectional plane. The local coordinate system is composed by local x-, y- and z-directions. The local y-direction is defined as the normal of the plane that best fits the resampled centerline points, which is obtained by using the Random Sample Consensus (RANSAC) algorithm [33]. The RANSAC algorithm starts by obtaining the parameters of the plane fitting three randomly selected centerline points. Then, it detects all the points located within a given distance from the plane and classifies them as inliers. These two procedures are repeated until the number of inliers converges to a maximum, and the final plane parameters are obtained. The local z-direction is defined as the direction vector that describes the line segment between two contiguous center points. Once the y- and z-directions are defined, the local x-direction is simply the vector cross-product of the two directions. Cross-sectional images are obtained by resampling the input volume along the x- and y-directions of the local coordinate system and within the neighborhood of the center point, rendering in an image of a certain size. The obtained 2-D slices are then stacked in order to construct the 3-D volume. 3) Segmentation of the Aorta: The goal of this step is to determine the boundary contours of the aorta in the transformed volume computed as described in Section III-B2c and map them to the original image space. Using a polar coordinate system, the problem of finding the aortic boundary reduces to a horizontal boundary detection which can be efficiently computed using the dynamic programming method. An approach similar to the one described in Section III-B1 is used, with the difference being in the manner in which the problem is formulated. In this formulation, an optimal boundary contour is detected in a transformed polar coordinate system. In addition, this formulation can be extended to detect a closed contour that deviates away from the circular shape by modifying the coordinate system of the search space, if an initial contour is available. One of the dimensions of this modified coordinate system corresponds to the parametric representation of the initial contour and the other corresponds to the normals of the contour. The search space can be restricted by controlling the length along the normals on either side of the contour. However, if the boundary has a high curvature segment, the normals on

941

the boundary may cross each other and this might result in the folding of the search space depending on the length of the normals. To overcome the aforementioned limitation, a novel method is used to compute the narrow band search space and the coordinate system. The narrow band is obtained by limiting the space at a certain distance from the contour through a distance transform. The inner and outer boundaries of the narrow band are parametrized with an equal number of points, which form one of the dimensions of the search space. A correspondence is sought between the inner and outer boundary points in order to obtain connecting line segments. These line segments are also parametrized with an equal number of points, which form the other dimension of the search space. This formulation allows the contour to deviate from a circular shape and does not cause the folding of the search space. The contour can be refined further following an iterative procedure. The details of the steps involved in the segmentation of the aorta are described next. a) Construction of the cost function: The optimal path detection problem can be solved using dynamic programming. However, its success depends largely on the cost function used. The cost terms used to construct the cost function should aid in describing the boundary of interest. Given an M × N image, the cost of a path is defined as the sum of the local energies along the path, and is given by N  (cF (vi ) + cO (vi−1 , vi )) C(v1 , v2 , . . . , vN ) = cF (v1 ) + i=2

where {v1 , . . . , vi , . . . , vN } are the N vertices of the path BN , and each vertex is represented by its x-coordinate, y-coordinate, ˆN is the and intensity value (x, y, q). The desired boundary B optimal polygonal line BN that has the global minimum cost C(BN ). The cost term cF (vi ) exploits the features that characterize the boundary and the term cO (vi−1 , vi ) imposes the spatial continuity thereby yielding a smooth boundary. These two terms

˜ are expressed as, cF (vi ) = K j =1 wj fj (vi )(i = 1, . . . , N ), and cO (vi−1 , vi ) = wK +1 (yi − yi−1 )2 (i = 2, . . . , N ), where the local cost terms f˜j (vi ) (j = 1, . . . , K) are transformations of the image features fj (vi ), the term (yi − yi−1 )2 is related to the vertical distance between vertices vi−1 and vi , K corresponds to the number of features used, and w represents the weights for the respective cost components. The feature-based cost function has the following features: ax(f 1 )−f 1 (v i ) Intensity: f˜1 (vi ) = mmax(f 1 )−m in(f 1 ) where ⎧ −190 for qi < −30 ⎪ ⎪ ⎪ ⎨ q + 60 for 40 < q < 100 i i f1 (vi ) = ⎪ > 100 100 for q i ⎪ ⎪ ⎩ otherwise. qi This term is computed by performing a nonlinear gray-scale modification in the input image, in order to accentuate the tissue of interest such as blood, then its inverse is obtained and is finally scaled from 0 to 1. It forces the boundary to follow the homogeneous path that contains the pixels with high blood intensity values.

942

IEEE JOURNAL OF BIOMEDICAL AND HEALTH INFORMATICS, VOL. 17, NO. 5, SEPTEMBER 2013

Gradient: f˜2 (vi ) = where

m ax(f 2 )−f 2 (v i ) m ax(f 2 )−m in(f 2 )

⎡

0

⎢ f2 (vi ) = h ∗ f1 (vi ), h = ⎣ 0 0

1 0

0

⎤

⎥ 0⎦

−1 0

and ∗ indicates convolution with 3 × 3 neighborhood of the vertex vi . The gradient cost component term is computed by performing a convolution operation between the filter h and the neighborhood of the vertex vi , then, it is inverted and scaled from 0 to 1. This term is responsible for moving the boundary toward the points that have a high gradient value in a direction perpendicular to the boundary. ax(f 3 )−f 3 (v i ) where f3 (vi ) = Local Entropy: f˜3 (vi ) = mmax(f 3 )−m in(f 3 )

− p log(p), and p is the empirical probability of the occurrence of a gray level within a 3 × 3 neighborhood of the vertex vi . The local entropy measures the randomness or the degree of disorder of an image. In homogenous regions, the changes in the gray levels are small, and thus, the local entropy value is low. However, in regions with gray level transitions (e.g., edge and line features) the local entropy is high. Since the local entropy feature is able to highlight such ROIs, the entropy-based cost component term forces the optimal path to follow the boundary of the aorta. 4 (v i )−m in(f 4 ) where f4 (vi ) is the disContinuity: f˜4 (vi ) = mfax(f 4 )−m in(f 4 ) tance from the vertex vi to the nearest vertex of the detected boundary in the previous slice. This cost term is added only when the aorta boundary contour is available from the previous slice. This cost component is obtained by computing the distance from the previous contour, and then it is scaled from 0 to 1. This term provides a 3-D surface continuity factor guiding the initialization of the subsequent slices and benefiting the low contrast slices from their neighboring slices. b) Initial boundary detection: First, a polar coordinate system is constructed using the pair of (ˆ x, yˆ) coordinates and radius rˆ as center and radius estimates, respectively. Then, an optimal path is detected using the cost function described previously through dynamic programming in the constructed polar coordinate system to obtain an initial contour. c) Boundary refinement: The computed boundary is further refined by recomputing the optimal path in a new coordinate system constructed by considering a narrow band around the current estimation of the aorta contour in both directions (inward and outward) as a local search space. Thus, the contour is updated only within this band through dynamic programming, and the narrow band is also updated in each iteration as the contour changes. The narrow band is bounded on either side by two curves (Bi and Bo ), which are at an Euclidean distance δ from the contour. The distance δ determines the capture range of the contour. If it is too large, the contour may get attracted to undesired features. If it is too small, the contour may not be able to capture the real aorta boundary. The inner and the outer boundaries of the narrow band are defined as circular linked lists of vertices by sampling them into N points. To determine the correspondence between the points of inner and outer contours, both contours’ lists are rotated such that the start positions of both lists have the shortest possible distance. Each line segment

Fig. 4. Depiction of the segmented aorta mesh in: (a) the transformed 3-D volume and (b) the original volume space.

connecting the corresponding pair of points from the curves Bi and Bo is discretized into M points. The new coordinate system is formed by rearranging the N boundary points from the inner and the outer contours, and the M points from each line segment connecting the N points, yielding an image grid of size M × N . The boundary estimation is refined by applying this step iteratively for T number of iterations. d) Transformation to original volume space: Then, a triangular mesh of the aorta mask extracted in the transformed 3-D volume is obtained [see Fig. 4(a)], and each point of the mesh is mapped back to its coordinates in the original volume [see Fig. 4(b)]. Finally, the binarized volume from the triangular mesh in the original volume space is computed to obtain the final segmentation. IV. RESULTS A. Experimental Data We used noncontrast, cardiac CT scans from 45 subjects to evaluate our proposed method. The images were collected at the Department of Imaging, Cedars-Sinai Medical Center, as part of the early identification of subclinical atherosclerosis using noninvasive imaging research (EISNER) registry. The noncontrast CT scans were acquired using an electron beam CT (EBCT) scanner, and the acquisition protocol has been described in detail by Dey et al. [34]. Each scan consisted of a stack of 49 to 79 axial slices with a matrix size of 512 × 512, a resolution of 0.68 mm in the x- and y- axes, and 3 mm in the z-axis. B. Aorta Segmentation The proposed aorta segmentation method was developed using MATLAB and Insight Segmentation and Registration Toolkit (ITK) libraries [35]. This study also benefited from submissions to the Insight Journal such as those by Mueller [36], [37], Yaniv [38], and Tustison [39]. Accuracy measure: In order to evaluate the performance of the proposed aorta segmentation method, three accuracy measures are defined. 1) Point Deviation: distance (in mm) between annotated start/end centerline point (xα , yα , zα ) and x, yˆ, zˆ) obtained by the method  the points (ˆ (xα − x ˆ)2 + (yα − yˆ)2 + (zα − zˆ)2 . 2) Point Set Deviation: Hausdorff distance (in mm) between a set of annotated points Pα , and the method’s output point set Pˆ . The Hausdorff distance is used as

AVILA-MONTES et al.: SEGMENTATION OF THE THORACIC AORTA IN NONCONTRAST CARDIAC CT IMAGES

a proximity measure between two sets of points, and is described as: H(Pα , Pˆ ) = max(h(Pα , Pˆ ), h(Pˆ , Pα )), where h(Pα , Pˆ ) = maxp α ∈P α minp∈ ˆ), and d(·) ˆ Pˆ d(pα , p is the Euclidean distance. In other words, the Hausdorff distance is the longest distance to travel between the closest points of Pα and Pˆ . 3) Volume Overlap: Dice similarity coefficient (DSC) and Jaccard similarity coefficient (JSC) between the manual ˆ segmentation method. Both (Sα ) and the proposed (S) similarity measures range from 0 to 1, with 1 indicating complete overlap. The ground truth available is the manual aorta segmentation Sα and the center points (xα , yα , zα ) for the ascending aorta (at the RCA ostium slice) and the descending aorta (at the last slice of the volume). The annotations were performed by graduate students trained in the process of annotation. Since the descending and ascending aorta run mainly vertically on the CT scan, the center of mass could be used to describe their centerlines, but the same does not hold true for the aortic arch. Therefore, using the center of mass at each slice may not provide an accurate portrayal of the full aorta centerline. To overcome this issue we computed a 3-D distance map from the annotations and obtained the centerline using the fast-marching minimal path method. Thus, the ground-truth centerline Pα is derived from the manual annotations as follows: 1) compute the distance map from the manually segmented edges and 2) find the minimal path between the annotated start and end center points by using the minimal path extraction framework (as in Section III-B2b). To this end, the termination value was set to 1.0, and the optimizer parameters learning rate and number of iterations are set to 0.5 and 500, respectively. Parameter selection: The parameters of the method were obtained by statistical analysis and by performing different experiments on a set of 15 randomly selected volumes. In the slice localization method, we computed a descriptor of length 105 from the three pyramid levels (5 + 20 + 80) and it was reduced to length 40 using PCA. The conduction parameter κ for anisotropic diffusion was set to 100. The intensity threshold parameter ν to compute the distance map was set to 10 HU. The radius ρ of the structuring element to perform morphological operations was set to 1. The distance threshold parameter γ for the likelihood image was set to 22.5 mm because a study performed by Garcier et al. [40] determined that the maximum radius of a normal thoracic aorta is 22.5 mm. The scaling parameters for the likelihood image sl and sh were set to 0.1 and 1, respectively. We observed that only couple of iterations were required to obtain a refined boundary contour. Therefore, the value for the number of iterations T for the aorta boundary contour refinement was set to three. The transformed volume was obtained by computing cross-sectional images of size 100 × 100 pixels. For the aorta position detection method described in Section III-B1, the following parameters need to be provided: 1) start- and end-slice for processing (or the number of slices to be processed S and the end slice); 2) number of points with maximum Hough value taken into account for dynamic programming Q; and 3) slice processing order, which can be from

943

the top slice to the bottom slice of the volume, or from the bottom slice to the top slice of the volume. For the purpose of the aorta segmentation method proposed in this study, there are three regions of interest where the aorta center points are obtained: the ascending aorta (AA), whose end slice is defined by the RCA ostium slice, the descending aorta (DA), whose end slice is also defined by the RCA ostium slice, and the LDA, whose end slice is defined as the last slice of the scan. The order for processing the AA region is set from the bottom to the top slice, and for the DA and LDA regions is set from the top to the bottom slice. We experimentally determined the number of Hough-maxima points by setting the number of slices to be processed to 17 for the AA and DA regions, and 21 for the LDA region. The aorta center points were detected using Hough-maxima points in a range from 5 to 100, and the number leading to the highest accuracy, which is measured in terms of Point Deviation and Point Set Deviation, was selected. Table I provides the accuracy measures of these experiments. After obtaining the optimal values for the number of Houghmaxima points for each region (30 for the AA and DA regions, and 20 for the LDA region), we determined the number of slices that should be processed to obtain the center points. The intent is to extract the maximum number of aorta center points, without compromising the accuracy. We processed a range of slices from 5 to 17 for the AA and DA regions, and from 5 to 21 slices for the LDA region to detect the center points. The results of this experiment are shown in Table II. From the results, we selected best values of 17 for the AA region, and 13 for the DA and LDA regions as the number of slices to be processed. For the detection of the aorta centerline described in Section III-B2, the required parameters are: for Section III-B2a, the ellipsoid radii along the 1) x-, 2) y-, and 3) z-axes; for Section III-B2b, the 4) learning rate, 5) number of iterations for the gradient descent optimizer, 6) termination value for the way-points, and 7) termination value for the end point; and for Section III-B2c the 8) maximal distance from plane, and 9) desired probability for no outliers. The distance d between the center points of the ascending and descending aorta at the RCA ostium slice was used to determine the ellipsoid radii along the x-, y-, and z-axes. The ellipsoid radii along the x- and z-axes are determined by the equation rxz = m · d + c − 22.5, where 22.5 is the maximum radius (in mm) of a normal thoracic aorta [40], and m and c were obtained by the statistical analysis of the linear relationship between d and the distance from the RCA ostium slice to a center point in the aortic arch (m = 0.2807 and c = 62.69). The ellipsoid radius along the y-axis was determined by the equation ry = (p · d)/2, where p is a scale factor selected empirically (p = 0.7). The settings for the learning rate and the number of iterations were determined based on the paired comparison matrix in Table III. The learning rate was determined in the range from 0.1 to 1, whereas the number of iterations was determined in the range from 250 to 2500. The accuracy for each pair of parameters was assessed and the pair with highest accuracy was selected. The learning rate and the number of iterations parameters were determined to be 0.1 and 2,250, respectively. Note that for larger learning rates, the number of iterations does not

944

IEEE JOURNAL OF BIOMEDICAL AND HEALTH INFORMATICS, VOL. 17, NO. 5, SEPTEMBER 2013

TABLE I ACCURACY MEASURES FOR THE SELECTION OF THE NUMBER OF MAXIMUM HOUGH POINTS

TABLE II ACCURACY MEASURES FOR THE SELECTION OF THE NUMBER OF SLICES TO BE PROCESSED

TABLE III POINT SET DEVIATION (MM) FOR THE SELECTION OF THE LEARNING RATE AND NUMBER OF ITERATIONS

TABLE IV ACCURACY MEASURES FOR THE SELECTION OF THE TERMINATION VALUE FOR THE WAYPOINTS

affect the accuracy and for smaller learning rates, more number of iterations are preferred. The termination value provides the rigidity/flexibility to the way the minimal path travels near the center points. A smaller termination value forces the centerline in the proximity of the points irrespective of the likelihood values. Thus, a smaller value is suitable for the end point. However, for the waypoints, we want to be more flexible allowing the centerline to follow the likelihood image but with some constraints imposed by

the waypoints. Thus, a higher termination value is more suitable for the waypoints. The termination value parameter was evaluated in the range from 5 to 100. The resulting accuracy measures are provided in Table IV. Note that the deviation from the centerline increases if the termination value for the waypoints is very small or very large. Also, a more relaxed termination value yields similar accuracy measurements. The termination value for the waypoints was selected as 30. The termination value of the end point was set to 1.0. The parameters

AVILA-MONTES et al.: SEGMENTATION OF THE THORACIC AORTA IN NONCONTRAST CARDIAC CT IMAGES

945

TABLE VI DESCRIPTIVE STATISTICS OF THE POINT SET DEVIATION MEASURE FOR THE CENTERLINE EXTRACTED BY THE AORTA CENTERLINE DETECTION METHOD, AS WELL AS THE DICE VOLUME OVERLAP AND JACCARD VOLUME OVERLAP MEASURES FOR THE SEGMENTED VOLUME OBTAINED BY THE AORTA SEGMENTATION METHOD

Fig. 5. Depiction of the influence of the likelihood image (high likelihood voxels are represented as red and low likelihood voxels are represented as green) and the guidance points on the centerline extraction. 3-D visualization of the likelihood image (a) without using the ellipsoid mask, and (b) using the ellipsoid mask. Notice the change in the likelihood image at the area pointed by the arrow. Depiction of the centerline extracted (c) using just the start- and endpoints and the likelihood image without using the ellipsoid mask, (d) using the start-, end- and way-points and the likelihood image without using the ellipsoid mask, (e) using just the start- and end-points and the likelihood image using the ellipsoid mask, and (f) using the start-, end-, and way-points and the likelihood image using the ellipsoid mask. 3-D isosurface of the manually annotated aorta added for reference. TABLE V DESCRIPTIVE STATISTICS OF THE POINT DEVIATION MEASURE FOR THE START AND END POINTS EXTRACTED BY THE AORTA POSITION DETECTION METHOD, AS WELL AS THE POINT SET DEVIATION MEASURE FOR THE WAYPOINTS

Fig. 6. (a) and (b) Box and whisker plots for the aorta position detection, (c) the aorta centerline extraction, and (d) the aorta segmentation methods: the boxes depict the 25th and 75th percentiles, the line within the box is the median, the lines projecting out of the box contain the adjacent values which are not more than 1.5 times the width of the box, and all remaining points are outliers.

maximal distance from the plane and the desired probability for no outliers, were empirically selected to be 0.5 and 0.999, respectively. The influence of the likelihood image and the waypoints on the centerline extraction is depicted in Fig. 5. It can be observed that the centerline extraction is successful when both the likelihood image with the prior information about the aorta shape and the waypoints obtained from the three ROIs of the aorta are used together. Validation: The aorta location detection method was assessed in terms of point deviation and point set deviation for the AA and LDA regions, and in terms of point set deviation for the DA region; the centerline extraction method was assessed in terms of point set deviation; and the segmentation method was assessed in terms of volume overlap. Tables V and VI (column a) provide descriptive statistics of these measurements. The accuracy for

the ascending aorta is lower than that for the descending aorta. This can be attributed to the presence similar structures (e.g., pulmonary artery and heart chambers) in the neighborhood of the ascending aorta. Therefore, the detection of the ascending aorta is more challenging than the detection of the descending aorta. Box and whisker plots (see Fig. 6) depict the distribution of the error measurements. Note that there more outliers in the detection of the ascending aorta. Additionally, Fig. 7 depicts a one-to-one comparison of the processed data at each step of the aorta segmentation method. The volume overlap measures are considerably lower for the subjects that have higher errors in the previous steps. Note that these results were obtained with the annotated RCA ostium slice.

946

IEEE JOURNAL OF BIOMEDICAL AND HEALTH INFORMATICS, VOL. 17, NO. 5, SEPTEMBER 2013

Fig. 7. Accuracy measures for each subject in the dataset at each step of the aorta segmentation method: (a) start and end point deviation for AA and LDA regions, (b) point set deviation of the waypoints for the three processed regions (AA, DA, and LAD), (c) point set deviation of the aorta centerline, and (d) volume overlap of the segmented aorta.

C. Complete Framework For the validation of the complete framework for aorta segmentation proposed in this paper, which includes the automated slice localization method to find the RCA ostium slice, the dataset was split in two groups: 1) a learning set of 15 randomly selected volumes, resulting in a total of 894 CT slices and 2) a testing set with the remaining 30 volumes, with a total of 1807 CT slices. The assessment of the method was performed in terms of the accuracy measures mentioned in Section IV-B. Descriptive statistics of these measurements are provided in Tables V and VI (column b). It can be observed that the errors in the automated slice localization affects the accuracy of the detection of the ascending aorta more than that of the descending aorta. It also reduces the accuracy of the extracted centerline. However, the overlap measures for the segmented aorta are comparable with the results with the manual slice localization. Box and whisker plots were generated to show the

distribution of the error measurements (see Fig. 8). A one-toone comparison of the processed data at each step of the aorta segmentation framework is depicted in Fig. 9. It depicts how errors at each stage affect the following stages and the overall performance of the method. While slice localization errors do not affect the end point detected in the LDA region, they do affect negatively the detection of the start point in region AA. This error can also propagate to the localization of the waypoints in the region AA. Nonetheless, it is important to point out that the centerline extraction method is able, in many cases, to cope with these errors and still provide an acceptable centerline path. Tolerance to such errors is quite advantageous, because it provides robustness to our method and allows for a successful aorta segmentation, rendering on a DSC of 0.85 ± 0.11, compared to 0.89 ± 0.06 obtained by manually annotating the RCA ostium slice. Fig. 10 shows the final segmentation of the thoracic aorta for a few example CT volumes.

AVILA-MONTES et al.: SEGMENTATION OF THE THORACIC AORTA IN NONCONTRAST CARDIAC CT IMAGES

947

landmarks for the successful segmentation of the aorta using their semiautomated method. An advantage of the method proposed in this study over these previously proposed methods, is that, besides the location of the RCA ostium slice which can be provided manually, the proposed method does not require any kind of additional input data such as the multiple aorta centerline models in the work by Taeprasartsit et al. [24], or the manually labeled images in the work by Isgum et al. [23], which can take about 60 min per scan to obtain.

V. DISCUSSION

Fig. 8. (a) and (b) Box and whisker plots for the aorta position detection, (c) the aorta centerline extraction, and (d) and the aorta segmentation methods using the RCA ostium slice obtained from the slice localization method: the boxes depict the 25th and 75th percentiles, the line within the box is the median, the lines projecting out of the box contain the adjacent values which are not more than 1.5 times the width of the box, and all remaining points are outliers.

D. Comparison With Related Work It is important to note that a direct comparison with previously proposed methods for aorta segmentation is difficult because the input datasets are different. The need to perform multiple registrations in the method proposed by Isgum et al. [23] is a disadvantage, because registering an atlas to a target volume is a very time-consuming task: a single registration typically took 15 min. In comparison, the running time of the method proposed in this paper for obtaining the aorta segmentation of one scan was approximately 9 min. The authors reported the performance of their method in terms of the JSC, obtaining an average of 0.78 ± 0.04, compared to the 0.80 ± 0.09 obtained by the aorta segmentation proposed in this paper when using the annotated RCA ostium slice, and the 0.75 ± 0.14 when using the automated RCA ostium slice. Taeprasartsit et al. [24] evaluated their method in 40 datasets, 17 of which were contrast enhanced scans (about 43% of the data). The JSC was used as a measure of performance, and the evaluation was performed using all possible combinations among eight models (255 model combinations in total). The best reported average JSC was 0.92 ± 0.03, but this value cannot be used for comparison with the aorta segmentation method proposed in this paper because 1) the authors report their performance in a mixed set with contrast and noncontrast data and 2) it is not clearly reported in how many cases the automated algorithm failed, requiring the input of three anatomical

The noncontrast cardiac CT scans that are routinely acquired for calcium scoring suffer mainly from low contrast and diffused boundaries between different anatomical regions. This problem poses significant challenges toward the segmentation of the aorta because the pulmonary artery and the heart chambers are located in the neighborhood. This problem is further magnified along z-axis because the resolution in z-axis is almost five times lower than in the other axes. This poses a challenge to segment the aortic arch when its cross-sectional plane deviates away from the x–y plane. To overcome this problem of low contrast and diffused boundaries in the lower resolution noncontrast CT data, we proposed an entropy-based cost function to determine the aorta boundaries. Though the segmentation results in the ascending aorta region are promising, the segmentation in the aortic arch region suffers from leakage in some cases. To cope with this problem, a stronger shape constraint could be used to restrict the segmentation contour within acceptable shape space. Also, the thoracic aorta has coronary arteries as branches and supraaortic branches. Therefore, we rely on the continuity (curvature) constraint to overcome the effect of branches. We have performed extensive experiments to determine the optimal set of parameters within a wide range of parameter values. Our method is quite robust to small changes in the parameters. For the detection of the aorta location the accuracy measures are very similar to each other regardless of the parameters selected as can be observed from the Tables I and II. In case of the centerline extraction, the accuracy is adversely affected by the selection of a small learning rate, a small number of iterations and a small termination value. But as can be observed in the Table III, for larger learning rates, the number of iterations does not affect the accuracy anymore, and a more relaxed termination value yields similar accuracy measurements. Moreover, it is important to address the problems arising due to pathologies. In cardiac noncontrast CT scans, the calcifications are defined as the voxels that have HU intensity values higher than 130 HU. We performed a nonlinear gray-scale transformation before computing the intensity and the gradient features for segmentation to reduce the effect of dense calcified plaques. This way the calcifications are included within the boundary of the segmented aorta. Our method can handle the aneurysms as it involves segmenting a dilated aorta, which may require to use a higher distance threshold in the likelihood map. However, our method cannot detect the presence of a dissected aorta and therefore, cannot segment all the parts of the dissected aorta.

948

IEEE JOURNAL OF BIOMEDICAL AND HEALTH INFORMATICS, VOL. 17, NO. 5, SEPTEMBER 2013

Fig. 9. Accuracy measures for each subject in the dataset at each step of the complete aorta segmentation framework: (a) slice localization error, (b) start and end point deviation for AA and LDA regions, (c) point set deviation of the waypoints for the three processed regions (AA, DA, and LAD), (d) point set deviation of the aorta centerline, and (e) volume overlap of the segmented aorta.

VI. CONCLUSION In this paper, we have proposed novel methods for the thoracic aorta segmentation in the noncontrast CT imaging data. Our method is fully automated and is able to segment the whole aorta, including the ascending and descending aorta and the aortic arch, from the cardiac CT scans. Our method achieved promising results on a challenging dataset of noncontrast, cardiac CT scans from 45 subjects. Nevertheless, the slice localization method can be further improved and extended. There are two primary directions we wish to investigate. The first direction is to use features from the neighboring slices and/or 3-D descriptors to learn the regression functions. The second direction is to extend the method to preserve the ordering of the slices in a volume, which is critical when multiple slices need to be detected (e.g., basal and apical slices for left ventricle detection in MR images).

Fig. 10.

Depiction of the segmented thoracic aorta in for 4 CT scan examples.

AVILA-MONTES et al.: SEGMENTATION OF THE THORACIC AORTA IN NONCONTRAST CARDIAC CT IMAGES

REFERENCES [1] NHLBI Fact Book, National Heart, Lung, and Blood Institute, Bethesda, MD, USA, 2011. [2] A. Lembcke, “Coronary artery calcifications: A critical assessment of imaging techniques,” Blood Purific., vol. 25, no. 1, pp. 115–119, 2007. [3] C. Iribarren, S. Sidney, B. Sternfeld, and W. Browner, “Calcification of the aortic arch: Risk factors and association with coronary heart disease, stroke, and peripheral vascular disease,” J. Amer. Med. Assoc., vol. 283, no. 21, pp. 2810–2815, 2000. [4] N. Rodondi, B. Taylor, D. Bauer, L. Lui, M. Vogt, H. Fink, W. Browner, S. Cummings, and K. Ensrud, “Association between aortic calcification and total and cardiovascular mortality in older women,” J. Internal Med., vol. 261, no. 3, pp. 238–244, 2007. [5] J. Witteman, F. Kok, J. van Saase, and H. Valkenburg, “Aortic calcification as a predictor of cardiovascular mortality,” Lancet, vol. 2, no. 8516, pp. 1120–1122, 1986. [6] C. Otto, B. Lind, D. Kitzman, B. Gersh, and D. Siscovick, “Association of aortic-valve sclerosis with cardiovascular mortality and morbidity in the elderly,” N. Engl. J. Med., vol. 341, no. 3, pp. 142–147, 1999. [7] I. Litovchik, R. Krakover, A. Blatt, and Z. Vered, “Coronary and aortic calcification: Is the relationship important?” Israel Med. Assoc. J., vol. 9, no. 4, pp. 328–330, 2007. [8] N. Wong, M. Sciammarella, Y. Arad, R. Miranda-Peats, D. Polk, R. Hachamovich, J. Friedman, S. Hayes, A. Daniell, and D. Berman, “Relation of thoracic aortic and aortic valve calcium to coronary artery calcium and risk assessment,” Amer. J. Cardiol., vol. 92, no. 8, pp. 951– 955, 2003. [9] Y. Adler, M. Motro, A. Tenenbaum, D. Tanne, E. Fisman, I. Wiser, B. Hovav, D. Stolero, and J. Shemesh, “Aortic valve calcium on spiral computed tomography is associated with calcification of the thoracic aorta in hypertensive patients,” Amer. J. Cardiol., vol. 89, no. 5, pp. 632–635, 2002. [10] Y. Adler, J. Shemesh, A. Tenenbaum, B. Hovav, E. Z. Fisman, and M. Motro, “Aortic valve calcium on spiral computed tomography (dual slice mode) is associated with advanced coronary calcium in hypertensive patients,” Coronary Artery Disease, vol. 13, no. 4, pp. 209–213, 2002. [11] P. Raggi, T. Callister, B. Cooil, Z. He, N. Lippolis, D. Russo, A. Zelinger, and J. Mahmarian, “Identification of patients at increased risk of first unheralded acute myocardial infarction by electron-beam computed tomography,” Circulation, vol. 101, no. 8, pp. 850–855, 2000. [12] O. Avila-Montes, U. Kurkure, R. Nakazato, D. Berman, D. Dey, and I. Kakadiaris, “Cascaded regression for CT slice localization,” in Proc. IEEE Int. Symp. Biomed. Imag. Nano, Macro, Chicago, IL, Mar./Apr. 2011, pp. 1881–1884. [13] U. Kurkure, O. Avila-Montes, and I. Kakadiaris, “Automated segmentation of thoracic aorta in non-contrast CT images,” in Proc. IEEE Int. Symp. Biomed. Imag. Nano, Macro, Paris, France, May 2008, pp. 29–32. [14] O. Avila-Montes, U. Kurkure, and I. Kakadiaris, “Aorta segmentation in non-contrast cardiac CT images using an entropy-based cost function,” in Proc. Soc. Photograph. Instrum. Eng. Med. Imag. Conf., San Diego, CA, USA, Feb. 2010, vol. 7623, p. 76233J. [15] T. Behrens, K. Rohr, and H. Stiehl, “Robust segmentation of tubular structures in 3-D medical images by parametric object detection and tracking,” IEEE Trans. Syst. Man Cybern., vol. 33, no. 4, pp. 554–561, Aug. 2003. [16] I. Adame, R. van der Geest, D. Bluemke, J. Lima, J. Reiber, and B. Lelieveldt, “Automatic vessel wall contour detection and quantification of wall thickness in in vivo MR images of the human aorta,” J. Magn. Res. Imag., vol. 24, no. 3, pp. 595–602, 2006. [17] S. W¨orz and K. Rohr, “Segmentation and quantification of human vessels using a 3-D cylindrical intensity model,” IEEE Trans. Med. Imag., vol. 16, no. 8, pp. 1994–2004, Aug. 2007. [18] D. Rueckert, P. Burger, S. Forbat, R. Mohiaddin, and G. Yang, “Automatic tracking of the aorta in cardiovascular MR images using deformable models,” IEEE Trans. Med. Imag., vol. 16, no. 5, pp. 581–590, Oct. 1997. [19] L. Lorigo, O. Faugeras, W. Grimson, R. Keriven, R. Kikinis, A. Nabavi, and C. Westin, “CURVES: Curve evolution for vessel segmentation,” Med. Imag. Anal., vol. 5, no. 3, pp. 195–206, 2001. [20] T. K´ovacs, P. Cattin, H. Alkadhi, S. Wildermuth, and G. Szekely, “Automatic segmentation of the vessel lumen from 3D CTA images of aortic dissection,” in Proc. Bildverarbeitung Med., 2006, pp. 161–165.

949

[21] Y. Zheng, M. John, R. Liao, J. Boese, U. Kirschstein, B. Georgescu, S. Zhou, J. Kempfert, T. Walther, G. Brockmann, and D. Comaniciu, “Automatic aorta segmentation and valve landmark detection in c-arm ct: Application to aortic valve implantation,” in Proc. Med. Imag. Comput. Comput.-Assist. Intervent., 2010, vol. 6361, pp. 476–483. [22] O. Ecabert, J. Peters, M. J. Walker, T. Ivanc, C. Lorenz, J. von Berg, J. Lessick, M. Vembar, and J. Weese, “Segmentation of the heart and great vessels in ct images using a model-based adaptation framework,” Med. Imag. Anal., vol. 15, no. 6, pp. 863–876, 2011. [23] I. Isgum, M. Staring, A. Rutten, M. Prokop, M. Viergever, and B. van Ginneken, “Multi-atlas-based segmentation with local decision fusion– application to cardiac and aortic segmentation in CT scans,” IEEE Trans. Med. Imag., vol. 28, no. 7, pp. 1000–1010, Jul. 2009. [24] P. Taeprasartsit and W. Higgins, “Robust method for extracting the pulmonary vascular trees from 3-D MDCT images,” in Proc. SPIE, 2011, vol. 7962, p. 796237. [25] D. Chittajallu, P. Balan, and I. Kakadiaris, “Automatic delineation of the inner thoracic region in non-contrast CT data,” in Proc. Int. Conf. IEEE Eng. Med. Biol. Soc., Minneapolis, MN, USA, Sep. 2009, pp. 3569–3572. [26] S. Lazebnik, C. Schmid, and J. Ponce, “Beyond bags of features: Spatial pyramid matching for recognizing natural scene categories,” in Proc. IEEE Comput. Soc. Conf. Comput. Vis. Patt. Recognit., Washington, DC, USA, 2006, vol. 2, pp. 2169–2178. [27] S. E. Umbaugh, Computer Imaging: Digital Image Analysis and Processing. Boca Raton, FL, USA: CRC Press, 2005. [28] P. Perona and J. Malik, “Scale-space and edge detection using anisotropic diffusion,” IEEE Trans. Patt. Anal. Mach. Intell., vol. 12, no. 7, pp. 629– 639, Jul. 1990. [29] J. Canny, “A computational approach to edge detection,” IEEE Trans. Patt. Anal. Mach. Intell., vol. 8, no. 6, pp. 679–698, Nov. 1986. [30] J. Devarajan and B. Subramaniam, “Applied anatomy of the aorta,” in Proc. Anesthesia Perioperat. Care Aortic Surg., Springer, New York, 2011, pp. 1–15. [31] P.-E. Danielsson, “Euclidean distance mapping,” Comput. Graph. Imag. Process., vol. 14, no. 3, pp. 227–248, 1980. [32] D. Mueller and A. Maeder, “Robust semi-automated path extraction for visualising stenosis of the coronary arteries,” Comput. Med. Imag. Graph., vol. 32, no. 6, pp. 463–475, 2008. [33] M. Fischler and R. Bolles, “Random sample consensus: A paradigm for model fitting with applications to image analysis and automated cartography,” Commun. Assoc. Comput. Mach., vol. 24, no. 6, pp. 381–395, 1981. [34] D. Dey, N. Wong, B. Tamarappoo, R. Nakazato, H. Gransar, V. Cheng, A. Ramesh, I. Kakadiaris, G. Germano, P. Slomka, and D. Berman, “Computer-aided non-contrast CT-based quantification of pericardial and thoracic fat and their associations with coronary calcium and metabolic syndrome,” Atherosclerosis, vol. 209, no. 1, pp. 136–141, 2010. [35] Kitware. (2010, May). The insight segmentation and registration toolkit. [Online]. Available: http://www.itk.org [36] D. Mueller. (2008, Jan.). Fast marching minimal path extraction in ITK. [Online]. Available: http://www.insight-journal.org/browse/publication/ 213 [37] D. Mueller. (2007, Jul.). Lookat transform initializer and oblique section image filter. [Online]. Available: http://www.insight-journal.org/browse/ publication/168 [38] Z. Yaniv. (2010, Jul.). Random sample consensus (RANSAC) algorithm, a generic implementation. [Online]. Available: http://www.insight-journal. org/browse/publication/769 [39] N. Tustison and J. Gee. (2009, Jul.). Introducing dice, jaccard, and other label overlap measures to ITK. [Online]. Available: http://www. insight-journal.org/browse/publication/707 [40] J. M. Garcier, V. Petitcolin, M. Filaire, R. Mofid, K. Azarnouch, A. Ravel, G. Vanneuville, and L. Boyer, “Normal diameter of the thoracic aorta in adults: A magnetic resonance imaging study,” Surg. Radiol. Anatomy, vol. 25, no. 3–4, pp. 322–329, 2003.

Authors’ photographs and biographies not available at the time of publication.