Accurate vessel segmentation with constrained B-snake.

2440

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 24, NO. 8, AUGUST 2015

Accurate Vessel Segmentation With Constrained B-Snake Yuanzhi Cheng, Member, IEEE, Xin Hu, Ji Wang, Yadong Wang, Member, IEEE, and Shinichi Tamura, Fellow, IEEE Abstract— We describean active contour framework with accurate shape and size constraints on the vessel cross-sectional planes to produce the vessel segmentation. It starts with a multiscale vessel axis tracing in a 3D computed tomography (CT) data, followed by vessel boundary delineation on the cross-sectional planes derived from the extracted axis. The vessel boundary surface is deformed under constrained movements on the cross sections and is voxelized to produce the final vascular segmentation. The novelty of this paper lies in the accurate contour point detection of thin vessels based on the CT scanning model, in the efficient implementation of missing contour points in the problematic regions and in the active contour model with accurate shape and size constraints. The main advantage of our framework is that it avoids disconnected and incomplete segmentation of the vessels in the problematic regions that contain touching vessels (vessels in close proximity to each other), diseased portions (pathologic structure attached to a vessel), and thin vessels. It is particularly suitable for accurate segmentation of thin and low contrast vessels. Our method is evaluated and demonstrated on CT data sets from our partner site, and its results are compared with three related methods. Our method is also tested on two publicly available databases and its results are compared with the recently published method. The applicability of the proposed method to some challenging clinical problems, the segmentation of the vessels in the problematic regions, is demonstrated with good results on both quantitative and qualitative experimentations; our segmentation algorithm can delineate vessel boundaries that have level of variability similar to those obtained manually. Index Terms— Vessel segmentation, constrained B-snake, computed tomography angiography (CTA), vessel-axis-tracing, shape and size constraints.

I. I NTRODUCTION

A

UTOMATIC vessel segmentation in threedimensional (3D) medical CT images has become increasingly important for diagnosis and quantification of

Manuscript received March 15, 2014; revised September 24, 2014 and December 26, 2014; accepted January 10, 2015. Date of publication April 7, 2015; date of current version April 29, 2015. This work was supported in part by the Fundamental Research Funds for the Central Universities under Grant HIT. IBRSEM. 201328 and in part by the Scientific Research Project through the Heilongjiang Provincial Education Department under Grant 12541164. The associate editor coordinating the review of this manuscript and approving it for publication was Prof. Jose M. Bioucas-Dias. Y. Cheng, X. Hu, and Y. Wang are with the School of Computer Science and Technology, Harbin Institute of Technology, Harbin 150001, China (e-mail: [email protected]; [email protected]; [email protected]). J. Wang is with the Department of Computer Science, Chiba Institute of Technology, Narashino 275-0016, Japan (e-mail: [email protected]). S. Tamura is with the Center for Advanced Medical Engineering and Informatics, Osaka University, Suita 565-0871, Japan (e-mail: [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/TIP.2015.2417683

vascular disease, for vascular surgery planning, and for patient-specific flow simulations. In the past decade, different approaches for vessel segmentation have been developed, and several good reviews have been published in [1] and [2]. Previous work on vessel segmentation can be roughly classified into three categories [2]. The first approach is based on features extracted from image content. Multiscale analysis of second derivatives is widely used for enhancement or detection of curvilinear structures in two-dimensional (2D) and 3D medical images [3], [4]. They share a common approach: the images are convolved with 3D Gaussian filters at multiple scales and the eigenvalues of the Hessian matrix at each voxel are analyzed in terms of a response function to determine the shape of the local structures in the image. The eigenvalues for the voxels that correspond to a linear structure would be different from those that correspond to a planar structure, blob, noise, or no structure. The response of the enhancement filter reaches its maximum when the scale of the filter matches the size of the local structures. The local structures can then be extracted using the local maxima [3]. Although these have proved to be useful for enhancement of line features, the final output of such procedures is not a direct segmentation of the input image. Skeleton-based methods [5] are algorithms in which subsequent 2D slices of vessels are resolved using tubular shape priors for ridge detection. Recently, another tracking methodology was proposed by Tyrrell et al. [6] using 3D cylindroidal superellipsoids and local regional statistics to extract topological information from microvasculature networks. These methods are shown to be robust against noise, however, their explicit parametrical shape priors are too exclusive, in case of complex vessel boundaries. The second approach is to use the minimal path techniques for vessel extraction. Minimal path techniques are commonly employed in interactive frameworks, requiring the definition of start and end points for each target vessel. Some works have proposed the definition of termination criteria to automatically stop the path propagation and relax the need for end points. Such criteria can for instance be devised through heuristic thresholds on the underlying features [7]. Classical vessel-dedicated minimal path algorithms focus solely on extracting the centerline curve. Notable exceptions are the works from [8]–[10], which share the core idea of incorporating an additional dimension corresponding to the vessel radius. The spatial position of the centerline and its associated radius are then optimized conjointly. The final approach is to apply vessel model for vessel segmentation. Model-based approach can be further

1057-7149 © 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

CHENG et al.: ACCURATE VESSEL SEGMENTATION WITH CONSTRAINED B-SNAKE

categorized into: 1) centerline + cross-section models (generalized cylinders), and 2) deformable models. In the first subcategory, the vessel surface is modeled with respect to a centerline curve, following generalized cylinder methods. A first approach in this vast family of models is to define vessel walls as 2D cross-sectional contours sweeping along the centerline. Classical cross-sectional models include circles, ellipses, star patterns [11] or parametric curves [12]. Cross-section models help reducing the complexity of the extraction procedure, under the hypothesis that the centerline estimation is reliable. Generalized cylinder models differ in the kind of centerlines and cross-sections they accept and the way they maintain the 3D coherence of the object. The straight homogeneous generalized cylinder model [13] is limited to linear centerlines and cross-sections varying as scaled versions of a reference shape. Of particular interest, the spatial intensity distribution is modeled as radial variations in the corresponding cross-sections, introducing the notion of intensity profiles. Mathematically appealing profiles include parabolic profiles and 3D Gaussian line profiles [14], for which saliency criteria and estimation bounds were analytically derived. In practice, vessels (especially large ones) often do not have rounded but rather plateau-like profiles, which are not adequately described by the aforementioned models. The simplest model accounting for that observation is the bar-like profile [15]–[17], which corresponds to a local cylinder with homogeneous intensity. The second subcategory is based on deformable models that evolve an interface through different forces: external forces, derived from the image, and internal force, constraining the contour geometry and its regularity. The first deformable model for image segmentation, known as the snakes model was introduced in [18]. Since its publication [18], deformable models have become one of the most active and successful research areas in medical image segmentation. Frangi et al. [19] proposed to reconstruct 2D vessel boundaries or 3D vessel walls using deformable surface models represented by B-spline surfaces. However, it is not possible to employ parameterized geometric models to effectively deal with whole vessel trees, as the models would be required to change topology during evolution. Other papers employing geometric deformable models can be found in [20] and [21]. Later, the geodesic active contour model was proposed by Casseles et al. in [22] as a geometric-variational alternative for snakes. The idea, similar to the snake model, is a minimization of a functional that integrates over an edge indicator function along a contour. However, the arbitrary parametrization in the snake is replaced with the curve’s arclength. The edge indicator function obtains low values in image locations where the gradient is high. Apparently, the gradient magnitude edge indicator was not enough for capturing thin structures. In [23] Vasilevskiy and Siddiqi used maximization of the inner product between a vector field and the surface normal in order to construct an evolution that is used for segmentation of thin structures. At the other end, Holtzman-Gazit et al. [24] integrate the better qualities of the above geometric methods in order to segment thin structures in volumetric medical images.

2441

Nevertheless, each of the existing techniques in the literature has limitations when used on challenging images (cases). The main challenges may be outlined as follows. i) There is a wide range of vessel width, from less than a voxel to 15 voxels wide. Thin vessels often have the lowest contrast. These methods may produce an incomplete segmentation of the vessels. Moreover, they can not detect the thin vessels (in diameter < 2.5 voxels) accurately. ii) Highly curved and touching vessels are more difficult to handle. In such situations, a false detection may occur. iii) The presence of pathologic structure attached to a vessel are particularly challenging for automatic vessel extraction. The existing methods cannot be used straightforwardly to achieve satisfactory segmentation result in the aforementioned cases, since challenges involved are different. Overcoming the above limitations is exactly what we are concerned with in this paper. In this paper, we take a very different approach to segment the vessels of interest. First, we extract vessel axes from a gray-scale 3D angiogram without any form of segmentation on the cross-sectional planes. Our method can produce continuous vessel axes. With the extracted axis, we then segment the vessels on the cross-sections that are derived from the axis. Such decoupling of segmentation and vessel-axis-tracing avoids incomplete segmentation of the vessel at problematic regions. Our approach is evaluated and demonstrated on CT data sets from our partner site. These data sets include: 1) acrylic rod phantom, 2) liver arteries, 3) carotid arteries, and 4) pulmonary vessel trees. The results of our work are compared with those from three closely related methods [20], [21], [24]. These methods will be discussed in detail in the following section. Our method is also tested on two publicly available databases and its results are compared with the state-of-the-art method [9]. II. L ITERATURE R EVIEW In this section, we review the other closely related approaches on active contour models or snake, and explain how differ from ours. An approach, proposed by Hernández-Hoyos et al. [20], is to segment the vessel in the cross-sectional planes based on the extracted vessel axis. In [20], the first step towards this purpose is the extraction of the vessel central axis using an expansible vessel skeleton. Contours are then detected in the cross-sectional plane locally orthogonal to the centerline with 2D active contours. The advantage of this approach resides in the simplicity of the model initialization and in the flexibility of its evolution. The algorithm starts from a single point inside the vascular segment to extract, and a simulated propagation of this point along the vessel defines its skeleton. In similar way, the active model for the contour detection starts from the single point determined by the axis position in the cross-section. Recently, Li et al. [21] presented a new region-based active contour model for segmentation of tubular structure.

2442


Lis method draws upon intensity information in local regions at a controllable scale. This has the advantage that the method is able to segment image with intensity inhomogeneity, and has desirable performance for thin vessels with weak object boundaries. Holtzman-Gazi et al. [24] explore a hierarchical approach to segment thin vessels in volumetric data. The method is based on geometric active surfaces that evolve according to geometric partial differential equations until they stop at the boundaries of the objects. The main advantage is that the method is able to extract two adjacent objects with similar gray values such as the intracranial blood vessels and the adjacent bones. In [26], Lee et al. use the 2D sampling to offer a speed advantage over the 3D operations for vascular reconstruction. This approach satisfies these two criteria: 1) automatic process and sub-voxel accuracy. It combines automatic seeding and tracking of vessels with radius detection based on active contours. First, an extrapolation step is made along the current tangential direction of the vessel and the cross-section of the vessel is extracted therein. Correction can be further refined with two distinct methods, namely calculating the intensity-weighted centroid of the pixels to improve the estimate of the orthogonal plane and fitting a 2D active contours to perform an accurate segmentation of the cross-section. Apparently, several problems remain when using the aforementioned methods. First, gradient vector flow is directed toward the vascular boundary. If the initial snake does not contain the vessel axis (i.e., center of the circle), the external force field will force the snake to collapse on one side of the vascular boundary. Second, snake-based tracking is to use the position of the captured object from the previous plane as the initial snake for the subsequent plane. So, if one uses snake for the external force in object tracking with this initialization, then the maximum displacement from one plane to the next plane is less than the radius of the vessel. Beyond this maximum displacement the snake will fail to capture the object. Third, snake is not able to detect the desired contour at problematic regions that contain diseased portions (pathologic structure attached to a vessel), touching vessels (vessels in close proximity to each other) and thin vessels with low contrast (both vessel and background appear with similar gray values). One remedy is to incorporate priors which can be integrated and used for known structures. Prior shape information typically represented as a statistical shape model (SSM) [10], is useful for robust segmentation. All techniques for learning shape models presented in the previous section assume that the training shapes are given as landmark vectors with inherent point correspondence. In practice, the training shapes are usually extracted from expert segmentations. Given a set of aligned corresponding training data, the SSM is built by learning a principal subspace spanned by the eigenvectors of the covariance matrix of the data via principal component analysis (PCA). For a given corresponding test data that is aligned to the model coordinate system, it is constrained by projecting into the learned principal subspace and imposing bounds on its the parameter vector, and is then transformed

back to its original coordinate system. Yet, another option is to impose the constraints into an active contour model to steer the surface evolution. For example, Lai and Chin [27] introduce shape constraints into active contour evolution. The energy functional is designed for shape constraint, and the constraint is derived through the energy minimization principle for the active contours in terms of geometric primitives from the basic principles of the calculus of variations. The evolution of the curve is computed as a weighted sum of a shape force and the active contour force. The latter manner is applied in the present study. The main differences from the traditional methods and our novel contributions may be enumerated as follows: i) Unlike traditional centerline-based method [9] that are applied to vascular segmentation by incorporating the spatial position of the centerline and its associated radius. Thus, these methods heavily depend on an accurate vessel centerline. In contrast, our algorithm does not need an accurate vessel centerline. Vessel centerline is used to provide the capturing path for the active contour model. Active contours can move and capture 2D vessel boundary on the cross-sections along the centerline path. Correction can be refined with the constrained active contours performing an accurate segmentation of the cross-section. ii) Active contour model with accurate shape and size constraints can move and capture vessel boundary on the cross-sections along the vessel axis path rather than freely in a 2D plane or 3D space. This helps to detect thin and low contrast vessels while avoiding false detection of vessels in the problematic region. iii) The boundary detection method based on the CT scanning model opens the road for overcoming the inherent limits on segmentation accuracy arising from finite resolution. This allows us to obtain the accurate boundary of the thin vessel. iv) An interpolate method can handle missing contour points at problematic region. v) The B-spline incorporates the constrained snake so that fewer sample points are require to implement contour evolution, thus, the evolution efficiency is improved. III. D ETECTION AND T RACKING OF V ESSEL A XES In this section, vessel regions are enhanced and extracted by applying multiscale filtering method. At this stage, the delineation of extracted regions is not assumed to be very accurate, and the detected regions are used as mask regions for vessel axis points (centerpoints) detection in the subsequent stages. Subsequently, centerpoints defined as intensity ridge points are detected at multiple scales by using our previous method [15]. Finally, tracking and connecting of centerpoints detected at multiple scales are performed in scale space to reconstruct continuous centerlines and their branching structures. A. Vessel Enhancement On CT angiography (CTA) images, the vessel intensity may vary within a relatively wide range due to blood flow rate and vessel dimensions. Vessels with various diameters on CTA


images can be regarded as 3D line-structures. Therefore, the problem of vessel enhancement may be cast as a task of multiple scale line analysis. x ) = [Ixi x j ]3×3 (1 ≤ i, j ≤ 3) be the Hessian Let ∇ 2 I ( matrix of a 3D function I (x, y, z) about a point x = [x y z] = [x 1 x 2 x 3 ], and λ( x ) = [λ1 ( x ) λ2 ( x ) λ3 ( x )] with λ1 ( x) > λ2 ( x ) > λ3 ( x ) be the eigenvalues of H with corresponding x ), e2 ( x ), e3 ( x ), respectively. Using eigenvectors given by e1 ( the matrix of the eigenvectors, we have ⎡ ⎡ ⎤T ⎤ x) x) e3 ( e3 ( x ) ⎦ ∇ 2 I ( x) ⎦ x )E = ⎣ e2 ( x ) ⎣ e2 ( E T ∇ 2 I ( e1 ( e1 ( x) x) ⎛ ⎞ λ3 ( x) 0 0 λ2 ( x) 0 ⎠ =⎝ 0 0 0 λ1 ( x) ⎛ ⎞ Ix x 0 0 Iy y 0 ⎠ =⎝ 0 (1) 0 0 Iz z Equation (1) states that the matrix H describes the second-order structure of local intensity variations around each point of I ( x ), and the local second-order features of I ( x ) can x ), λ2 ( x ), and λ3 ( x) be directly obtained by the eigenvalues λ1 ( which are equivalent to second partial derivatives of I ( x) with respect to a rotated coordinate system x = [x y z ]T with x = E x . The local cross-sectional profile of an ideal 3D line model is then considered to have a 2D Gaussian shape x 2 + y 2 ) ex p(− 2σr2

(2)

where the standard deviation σr is related to the width of the line, [x y z ]T = R[x − cx y − c y z − cz ]T , where R is a 3 × 3 rotation matrix, and [cx c y cz ]T is the center of the profile. According to equation (1) and equation (2), the absolute values of λ2 ( x ) and λ3 ( x ) associated with a point on a line-structure should be large. Therefore, we can utilize the following condition to determine if a point is on a line-structure or not

x )| |λ2 ( x )| ≤ |λ3 ( x )| |λ1 ( (3) λ2 ( x ) < 0 and λ3 ( x) < 0 x ) and λ2 ( x) On the other hand, the absolute values of λ1 ( associated with a point on a plate-structure should be small, but x ) associated with such a point should the absolute value of λ3 ( be large. We can utilize the following condition to eliminate such plate-structures.

x )| ≤ |λ2 ( x )| |λ3 ( x )| λ1 ( (4) x) < 0 λ3 ( We use normalized second partial derivatives of the Gaussian function with standard deviation σ f (corresponding to width of the filter) to determine the elements of the Hessian matrix x; σ f ) = G xi x j (

∂2 1 x 2 + y2 + z2 ex p(− ) 3/2 (2π) σ f ∂ x i ∂ x j 2σ 2f

(5)

2443

where x i or x j denotes x, y or z. The Hessian matrix H can be obtained by H ( x ; σ f ) = [G xi x j ( x ; σ f ) ∗ I ( x )]3×3 1 ≤ i, j ≤ 3

(6)

where ∗ denotes the convolution operator. The eigenvalues of H ( x ; σ f ) are also given by λ( x; σ f ) = x ; σ f ) λ2 ( x ; σ f ) λ3 ( x ; σ f )] with λ1 ( x; σ f ) > [λ1 ( λ2 ( x ; σ f ) > λ3 ( x ; σ f ). We define L( x ; σ f ) to determine if λ( x ; σ f ) satisfies equation (3) ⎧ 2 ⎨e(−α(L (x ;σ f )−1) ) , i f λ2 < 0, λ3 < 0, L( x; σ f ) = (7) |λ1 | |λ2 | ≤ |λ3 |. ⎩ 0, other wi se. where α > 0 and 2λ1 − λ2 − λ3 L ( x; σ f ) = ( x; σ f ) 2 λ21 + λ22 + λ23 − λ1 λ2 − λ1 λ3 − λ2 λ3 As stated in equation (4) for a point x on a platestructure, the corresponding |λ3 ( x ; σ f )| should be great and |λ2 ( x ; σ f )| should be small. On the other hand, both x ; σ f )| and |λ3 ( x ; σ f )| corresponding to a point x on |λ2 ( a line-structure, should be greater according to equation (3). Therefore, we can use L( x ; σ f ) given in equation (7) and x ; σ f ) to measure lineness and generate the line-structure λ2 ( enhanced image for a given σ f from Sline (I ; σ f ) = −L( x ; σ f )λ2 ( x; σ f )

(8)

where I can be regarded as the intensity itself for a point x. Enhancement of vessels can be achieved by a nonlinear combination of multiple scale filtering equation (7). That is Mline (I ; σ1 , n) = max Sline {I ; σi } 1≤i≤n

(9)

where σ f = σi = σ1 · 2(i−1)/2 , σ1 = 1 and i = 1, 2, 3, 4. A detailed description for justifying the choice of using the max operator in equation (9) has been published in [3]. B. Detection of Vessel Axis Points Assuming that e1 , e2 and e3 are eigenvectors corresponding to λ1 , λ2 and λ3 , respectively. For the vessel structure, e1 corresponding to λ1 is expected to give its tangential direction, both |λ2 | and |λ3 |, directional second derivatives orthogonal to e1 , should be large on its medial axis. We assume that the tangential direction is given by e1 at the voxel x0 along the medial axis. The 2D intensity function, C(u, v), on the cross-sectional plane of I ( x ; σ ) orthogonal to e1 , should have its peak on the medial axis. The second-order approximation of C(u, v) is given by T T u u 1 u 2 C(u, v) = I ( x0 ; σ ) + ∇C0 + 2 ∇ C0 v v v (10) where ue2 + ve3 = x − x0 , ∇C0 = [∇ I · e2 ∇ I · e3 ]T (∇ I is the gradient vector, that is, ∇ I ( x 0 ; σ )), and λ2 0 ∇ 2 C0 = 0 λ3

2444


C(u, v) should have its peak on the medial axis of the line. The peak is located at the position satisfying

find the local maximum, which we regard as the starting axis point. We denote it by ( x 0 , s0 ), that is (x 0 , y0 , z 0 , s0 ). In this procedure, the starting point is searched from the largest scale. This means that the connecting process prioritizes the extraction of distinct large-diameter vessel segments. This coarse-to-fine strategy stabilizes the extraction of the whole vessel structure. Given the starting point ( x 0 , s0 ), axis points are connected toward the opposite two directions o( x 0 , s0 ) and −o( x 0, s0 ) from the starting point. Among adjacent axis points which can be continuously connected, the axis point having the maximum filter response r ( x , s) is connected sequentially. In the following, a method for connecting of the axis points only toward direction o( x 0 , s0 ) described since connecting of them toward direction −o( x 0, s0 ) is essentially the same. We assume that axis point ( x i , si ) has already been connected and the sign of o( x i , si ) has been changed so as to x i − xi−1 ) > 0. Firstly, axis point candidates satisfy o( x i , si ) · ( satisfying the continuity conditions of axial direction and position are selected among adjacent voxels in the direction of travel. Secondly, the axis point having the maximum filter response r ( x , s) is finally selected among the selected axis point candidates as ( x i+1 , si+1 ). The criteria for selecting the axis point to be connected next, ( x i+1 , si+1 ), are formally described as follows. Let Vad j denote the adjacent voxels in the direction of travel. Vad j is described as

∂ ∂ C(u, v) = 0 and C(u, v) = 0 (11) ∂u ∂v By solving equation (11), we have the offset vector, 2 q = [q x q y qz ]T = le2 + te3 , where l = − ∇ λI ·e 2 ∇ I ·e3 and t = − λ3 . For the medial axis to exist at the voxel x0 , the peak of C(u, v) needs to be located in the territory of voxel x0 . Thus, the medial axis is detected only if qx ≤ 1/2, q y ≤ 1/2 and qz ≤ 1/2. We assume that σ is discretize into n levels, and denoted as σs , where s = 1, 2, . . . , n and σ1 < σ2 < · · · < σn . The partial second-order derivatives and Hessian matrix computed at each scale are analyzed to judge whether the vessel axis passes through each voxel and determine the subvoxel position and orientation of the axis point. Let B( x , s) denote four-dimensional (4D) array of binary volume data computed at different scale σs , which is defined as ⎧ ⎨1, I f axi s poi nt exi sts i nsi de the voxel p. B( x , s) = (12) ⎩ 0, Other wi se. The subvoxel position ( x , s) and direction o( x , s) of the 3 are assigned to axis point, and filter response r ( x , s) = λ2 +λ 2 the voxel satisfying B( x , s) = 1, where |o( x , s)| = 1. C. Tracking and Connecting of Vessel Axis Points The connecting algorithm repeats the process consisting of selecting the starting axis point ( x 0 , s0 ) and connecting axis points from the starting axis point ( x 0 , s0 ) toward the opposite two directions o( x , s) and −o( x , s) until no new axis points satisfy the criteria for connecting. In the following, criteria for starting axis point selection and subsequent connecting are described. The starting axis point is searched in the ( x , s)-space, that is, the (x, y, z, s)-space. Here, we refer to the initial position for searching “starting axis point” in the ( x , s)-space as “initial position” to distinguish “starting axis point” from it. The initial position is firstly selected among unconnected axis points, which have not been connected to any axis yet, in the (x, y, z)-subspace at the largest scale, and then hill-climbing is performed to find an unconnected axis point having the local maximum of filter response in the ( x , s)-space. Let σsinit = σn be the initial scale for initial position determination, that is, sinit = n. The detailed description of the starting axis point search method is as follows: a) The initial position xinit = (x init , yinit , z init ) having the largest filter response, that is, the maximum point in the r ( x , sinit ), is selected. If no unconnected axis point is found in the r ( x , sinit ), the initial scale is moved to the next largest scale. That is, scale level sinit is set as sinit = sinit − 1, and selection of the initial position is performed again. When sinit < 1, the whole multiscale connecting process terminates since that means no unconnected axis point remains. b) Starting from initial position (x init , yinit , z init , sinit ) hill-climbing is performed in the (x, y, z, s)-space to

x , s)|B( x , s) = 1, o( x i , si )( x i − xi−1 ) > 0, Vad j = {( x i , si ), ( x , s)) ≤ R} (13) di stmax (( where R is the radius defining an adjacent area of ( x i , si ), and di stmax is the discrete distance in n-dimensional array defined as = max{|a1 − b1 |, |a2 − b2 |, . . . |an − bn |} a , b) di stmax ( (14) which is essentially n-dimensional extension of eight neighbor distance of 2D image. Let Vcont denote the axis point candidates satisfying the continuity conditions. Vcont is described as x , s)|( x , s) ∈ Vad j , |o( x i , si ) · o( x , s)| > T1 , Vcont = {( x , s), (( x i , si ), o( x i , si ))} < T2 } di st point −line {( (15) are threshold values, and where T1 and T2 di st point −line {( x , s), (( x i , si ), o( x i , si ))} denotes the distance between 3D point ( x , s) and the 3D straight line defined x i , si ). by ( x i , si ) and o( In our implementation, we firstly use R = 1 to determine Vad j , and then the radius defining the adjacency is changed to R = 2 is used only when Vcont is an empty set using R = 1. If Vcont is still an empty set using R = 2, the connecting process terminates. The other connecting process is performed toward the opposite direction. After it terminates, the starting point for the next connecting process is selected.


IV. V ESSEL S EGMENTATION Active contours or snakes are parametric or nonparametric, closed, or open curves that can move on the image plane and capture an object boundary. Snake used here is parametric, closed curves that can move on the vessel cross-section planes and capture vessel contour. Constraints for general shape have been introduced for active contour evolution in a statistical framework [27], [28]. In this paper, we introduce position, shape and size constraints on the cross-sections and incorporate the constraint into an active contour model. In addition, the Gradient Vector Flow (GVF) [29] force is combined with the B-spline snake. B-spline snakes have several characteristics which make them well suitable for describing vessel boundary as well as snake evolution: 1) The B-spline implicitly incorporates contour smoothness and avoids the adhoc tension and rigidity parameters of the traditional GVF snake. 2) Fewer sample points are required to implement contour evolution for the B-spline, and thus, the evolution efficiency is improved. A. Parametric Snake Model 1) The B-Snake: Once a vessel axis is extracted in a 3D space and regularized, we define a steady 3D vector field on a regular grid that is identical (in terms of sample spacing and dimensions) to the original 3D angiogram. The vectors are the gradient of the Euclidean distance transform (EDT) of the vessel axis. With this dedicated steady vector field, we employ a commonly used technique in flow visualization, called stream surfaces [30], to construct the cross-sections. After the cross-sections are generated, the next step is to delineate the vessel boundary on those surfaces. In general, the vessel boundary to be outlined is round on the cross-section. Therefore a circular (closed and round) contour drawn on the cross-section around the axis point is sufficient for the boundary delineation. The traditional snake can detect the desired boundaries under the influence of internal forces from the curve itself and external forces from the image data. However, the evolution of traditional snake is time-consuming as we know. A more economical realization of snake can be reached by using far fewer state variables by B-Splines. The B-Splines are piecewise polynomial functions that provide local approximations to contours using a small number of parameters (control points), and what more, internal forces are not required in deforming B-Snake since the representation of B-Spline has guaranteed smoothness by hard implicit constraints. The B-Splines can be constructed of any order (with the complexity of the possible conformation increasing with spline order), but the third order (cubic) is most often encountered. The cubic B-Splines exhibit a favorable tradeoff between their shape complexity for describing natural curves and computational burden required to solve them. Therefore, we are concerned with the implementation of cubic B-Spline in this paper. Here we define a closed uniform cubic B-Spline with a given set of n + 1 control points, qi = [u i v i ]T , i = 0, 1, . . . , n. The B-Spline consists of n + 1 connected curve segments

2445

gi (s) = [u i (s) v i (s)]T . Each curve segment is a linear combination of four cubic polynomials by the parameter s, s ∈ [0 1]. That is ⎤ ⎡ qi mod(n+1) ⎢ q(i+1)mod(n+1) ⎥ ⎥ (16) gi (s) = M R (s) ⎢ ⎣ q(i+2)mod(n+1) ⎦ q(i+3)mod(n+1) where M R (s) is the matrix of the B-Spline bases, ⎡ 1 1 1 − − ⎢ 6 2 2 ⎢ ⎢ 1 1 ⎢ −1 3 ⎢ 2 2 2 ⎢ M R (s) = s s s 1 ⎢ 1 1 ⎢− 0 ⎢ 2 2 ⎢ ⎣ 1 2 1 6 3 6

given by ⎤ 1 6⎥ ⎥ ⎥ 0⎥ ⎥ ⎥ ⎥ 0⎥ ⎥ ⎥ ⎦ 0

Then the cubic B-Spline is defined as follows: C(s) = gi (s), 0 ≤ s ≤ 1

(17)

i

Since the internal forces are not required in deforming B-Snake, the B-Snake energy is calculated by the sum of the external forces E ext along the B-Spline curve C(s) as follows: 1 E ext (C(s))ds (18) E snake = 0

2) The Proposed Snake: We introduce position, shape and size constraints, and incorporate the constraints into B-snake. The total energy of the constrained snake is represented as E snake = λ1 E ext + λ2 E posit ion + λ3 E shape−size + λ4 E sampling

(19)

where E posit ion , E shape−size and E sampling are respectively, the position, shape, size, and sampling constraints and are defined in Section IV-C. The nonnegative terms {λ1 , λ2 , λ3 , λ4 } give the relative strengths of the respective energy components and are selected empirically. Here, in order to increase the capture range of the snake, the GVF [29] is selected as the definition of external forces for the B-Snake tracking. B. Detecting Contour Point and Dealing With Missing Contour Points In this section, we introduce a model-based method for producing the accurate contour points. An interpolate method is subsequently described for handling the missing contour points at problematic region. 1) Contour Point Detection on the Vessel Cross-Section: Zero-crossings of the second-order derivative along the gradient direction were introduced by Haralik [31] and then used by Canny [32] as 2D edge detectors. Due to the intrinsic limitations in the finite spatial resolution, for small vessels (diameter< 2.5 voxel), the conventional edge methods can not detect the boundary of the vessel accurately. If the acquisition process is accurately modeled, it is possible to find the

2446


Fig. 2. Missing contour point at problematic region. The model-based method cannot locate a contour point (arrow) along ray at problematic region.

Fig. 1. Estimating the contour points with the new model-based method. (a) 2D vessel cross-section plane orthogonal to the vessel axis. (b) 1D profile of the ideal intensity distribution along one ray from the axis point. (c) Graph shows the modeled intensity profile and actual intensity profile along ray after applying the Levenberg-Marquardt algorithm. The fit procedure yields the radial distance d from the vessel center-point.

accurate boundaries of the smaller vessels. Therefore, a new method based on the scanning model is proposed for this application. We use the scanning model to predict the shape of the gray level profiles along rays cast from the centroid of the vessel cross-section and use a maximum-likelihood method to the vessel radius and boundary point. Fig. 1 shows a simple 2D model of a single slice of an ideal vessel cut orthogonal to the vessel axis. Assuming that the vessel cross-section is not perfectly round, 6 rays are spaced every 60° from the axis point on the cross-section and then 6 boundary positions are estimated along the rays. The vessel center is at the origin, and 6 radial distance di , i = 1, 2, . . . , 6, from the vessel center to the vessel boundary point are estimated. Estimate of radial distance, di , is given in Appendix A. Here, the extracted centerline is not assumed to be very accurate, so, we determine 6 different boundary positions. After determining radial distance di , the boundary positions (contour points) of the vessel on the cross-section can be expressed as bi = [di cos θi di sin θi ]

(20)

in which θi = 0° + i ∗ 60°, i = 0, 1, . . . , 5. Similarly, for each of the cross-sections, six rays are spaced every 60° from the axis point. Vessel shape and size on the cross-section are represented by its axis point and the six points of intersection between the individual rays and the contour. A rough delineation of each vessel cross-section (i.e., discrete samples obtained from the model-based method) is used to obtain position and shape constraints in Section IV-C. The contour points, bi , are therefore called constraint points. Note that for larger vessels (in diameter> 2.5 voxels), zero-crossings method is used to detect the contour points of the vessels on the cross-section. This is because zero-crossings method can detect the larger vessels accurately [15]. 2) Dealing With Missing Contour Points at Problematic Region: One final important issue for locating the contour points is what the model-based method should do when it is unable to locate a contour point along ray (i.e., determine the

radius) from a given centroid. This might happen, for example, there are touching vessels and pathologic structure attached to a vessel in those problematic regions. In such situation, the radius at the problematic region is at least two times as large as the mean of determining the radii on the same cross-section [Fig. 2]. As shown in Fig. 2, mean of 5 determined radii is used for locating the edge point at problematic region. Another case where this typically happens is the regions, where the vessels are so low contrast or surrounded by noise, which does not exhibit a well difference from the peripheral tissue. Since the algorithm does not have any other means of locating the contour points, the only viable solution to this problem is to interpolate or extrapolate the contour point from neighboring contour points. The algorithm can be described as follows. First of all, the radius of the vessel is extracted for each axis point. After this, if there is a gap in some vessels due to the low contrast or the presence of noise, i.e., if the radius of the vessel is undefined at some axis point, but there are some points in front and behind the current axis point that have a defined radius, the radius for the current axis point is obtained by linear interpolation. This can be formalized as follows. Let m be the index of the last point and n be the index of the next point with a defined radius, respectively. Let a be the length of the axis from m to the current point k and b be the total axis length from m to n. Then the radius of the current point k is given by b−a a dm + dn (21) b b This scheme can easily be extended to the case where either m or n are undefined, i.e., the vessel radius (i.e. from axis point to contour point) is undefined at either end of the vessel. dk =

C. Various Constraints We observed that the vessel boundaries to be segmented on the cross-sections are mostly round (appear approximately circle on the cross-sections), compared with those on either axial, coronal or sagittal image planes. This is reasonable since the vessel lumen is often circular under physiological perfusion pressures, even when atherosclerosis is present in the artery wall [31]. Therefore, by segmenting vessels on the cross-sections, we can reduce the comparatively complex 3D segmentation problem into a series of simpler 2D vessel boundary delineation problems. Active contours can move and capture 2D vessel boundary on the cross-sections with position, shape and size constraints rather than freely in a 3D space. This helps to avoid shrinkage of the model and


bunching of the boundary points at regions that two adjacent objects have the similar gray values. 1) Position Constraint: The B-snake evolves by driving the movements for control points, however, the snake energy transmit to each control point is summed by the points along the snake curve. Therefore, we can structure constraint energies along the snake curve to constraint the movements of control points. As described above, the vessel boundary appears approximately circle on the cross-section. The centroid of this circle can be estimated by the obtained constraint points (contour points) bi , i = 1, 2, . . . , 6. Then, the evolving B-snake center will be constrained to this centroid. A detailed description is given in Appendix B. 2) Shape and Size Constraint: Based on the above description, a circle model can be structured by the constraint points bi with estimated centroid Cc and average radius. Then the evolving snake is forced to this model. Consequently, the snake curve points are forced to maintain a circle shape and transmit the constraint energy to each control point. A detailed description is given in Appendix B. 3) Sampling Constraint: During the traditional snake evolution, some portion of the snake will be stretched while the other portion will be shortened. So, compression as well as rarefaction of the snake points occurs during evolution. In the B-snake evolution, these actions also happen as compression and rarefaction of the control points. This gives a call for the insertion of control point at the place with maximal energy force to keep the uniform space [35]. Unfortunately, the insertion process represents a significant computational expense. Thus, we use an implicit sampling technique to eliminate the need for the explicit inserting of the control points without additional computation. The resultant addition to the B-snake energy functional is referred to as the sampling constraint. Because the evolution of B-snake is discrete, we introduce the sampling constraint energy in the discrete domain. The continuous parameter s used so far to denote the snake point position C(s) is indexed by j ∈ 0, 1, . . . , m − 1 with m being the total number of snake curve points. So, we have a discrete snake curve as C(s j ) = [u j v j ]. A detailed description is given in Appendix B. D. B-Snake Evolution and Constraints Incorporation So far, we have done appropriate preparations: we have an active contour model combined with the B-spline; we have different constraints containing position, shape and size, and sampling. The thing remaining is to: 1) implement the B-snake evolution and 2) incorporate the constraints into the B-snake. 1) Implementation of the B-Snake Evolution: Once we discretize snake point position C(s) as C(s j ) = [u j v j ], equation (17) can be given as a matrix form Q = [M T M]−1 M T C

(22)

where M is the matrix of the B-Spline bases [36], Q is the matrix of control points, Q = [q0 , q1 , . . . , qn ]T . To avoid expensive computation we employ standard steepest descent method to minimize equation (18). The

2447

steepest descent technique for evolving conventional B-snake gives the following equations: ∇C = (−

∂ ∂ (λ1 E ext ), − (λ1 E ext )) ∂u ∂v

(23)

∂ (λ1 E ext ) To facilitate description we denote that E u = − ∂u ∂ and E v = − ∂v (λ1 E ext ). Combine equation (22) and equation (23), the B-snake evolution can be described as

∇ Q = [M T M]−1 M T (E u , E v )

(24)

2) Constraints Incorporation: We also use the calculus of variations to minimize the constraint energies, and arrive at several Euler equations corresponding to different constraints. All the Euler equations derived from calculus of variations are linear. This property makes it very easy to incorporate all the constraints stated so far in the B-snake evolution. Then we will introduce the Euler equations of the constraints and incorporate the constraints into the B-snake. A detailed description is given in Appendix C. E. Voxelization of B-Snake Model Once the B-snake model converges, the vessel contours are generated by sweeping along the vessel axis. As we described in Section III-A, e1 can be expected as the normal direction of the cross-sectional plane, while e2 and e3 are the binormal vectors orthogonal to each other. For each cross-section plane, the obtained contour points C(s j ) = [u j v j ], j = 1, 2, . . . , m, can be described below: C j = pc + r j cosθ j e2 + r j si nθ j e3

(25)

where C j = [x j y j z j ] is the coordinates of contour points, [r j θ j ] represents the location of the contour point C(s j ) in polar coordinates of the local cross-sectional plane, pc = [x c yc z c ] is the axis point coordinate. With equation (25) the contours C(s j ) which outcome from the B-snake model delineate the boundary of a vessel along its axis in 3D space. Subsequently, a 3D vessel boundary surface model can then be constructed from the extracted contours by joining the vertices of adjacent ones (triangulation). Due to the same in the number of vertices on each extracted contour, it is easy to join the adjacent vertices to construct a 3D vessel boundary surface model. As a final step, the surface model is voxelized into a binary volume with voxel size and dimensions identical to the original angiogram. V. E XPERIMENTAL R ESULTS Our validation experiments are divided into two stages. First, our method is tested and demonstrated on CT datasets from our partner site (The 2nd affiliated Hospital of Harbin Medical University). These data sets include: 1) acrylic rod phantom, 2) liver arteries, 3) carotid arteries, and 4) pulmonary vessel trees. The results of our work are compared with those from other related methods [20], [21], [24]. In the second stage, our approach is tested on two publicly available databases and its results are compared with the recently published method [9]. 2011 VascuSynth Sample (3 datasets) and March 2013 VascuSynth Sample (10 datasets)

2448


Fig. 3.

Physical appearance of the acrylic rod phantom.

are available on line. We use the 10 datasets on the March 2013 to evaluate the performance of our algorithm; moreover, 3 datasets on the 2011 are used in the noise experiments to compare with the Cetin’s method [9]. We also demonstrate the quantitative performance of our centerline extraction step on 32 CT cardiac angiography data sets used in [9] and comparison to Cetin’s method. 32 CT data sets can be obtained from Rotterdam Coronary Artery Evaluation Framework [38], which provides eight training datasets with ground truths and 24 test datasets. A. Illustration and Validation of the Algorithm on CT Data Sets From Our Partner Site In this section, we first give the parameters of the algorithm. Subsequently, the algorithm is evaluated on CT data sets from our partner site and compared with other related methods. In the images, the gray value of each voxel is adjusted as follows: a gray value (GV) of 0 corresponds to −1024 Hounsfield units (HU) and 1024 GV corresponds to 0 HU (i.e, GV( x ) = HU( x ) + 1024). 1) Estimation of Parameters: The parameters of the algorithm are estimated in the training sets (two Liver artery datasets, one carotid artery dataset and two lung vessel datasets) from several experiments as follows. (i) The parameters for connecting the vessel centerpoints: We perform experiments on the acrylic rod phantom, liver artery and lung vessel tree to estimate the T1 and T2 used in equation (15) for connecting vessel centerpoints. In the acrylic rod phantom, the curve acrylic rods (shown in [Fig. 3]) are used to perform the estimation of parameters. T1 and T2 should be maintained within 0.1–1.0 voxels. The most accurate centerline result is obtained for the T1 = 0.9 and T2 = 0.5. Also, the radius R is selected 1 or 2 for this study based on visual assessment. (ii) The parameters for detecting the accurate contour points on the vessel cross-section: Acrylic rod phantom with the known radius is imaged using CT scanner for estimating scanner point spread function σ ps f in equation (28). A detailed description for estimating σ ps f has been published in our previous study [15]. The parameter K used in equation (28) may be estimated from the gray level at the centroid of the vessel cross-section, I (0), using the approximation that 2 2 2 e −di /2σ ps f . Moreover, the initial value I (0) ≈ 2π K σ ps f for is determined from the gray level of the peripheral 2 . The values tissue, I (t), t > di , using I (t) ≈ 2π K ασ ps f

for K and α are estimated along 120 rays (N = 120), the average of the results is regarded as the final K and α. These parameters were estimated for different types of images as follows: · Phantom dataset: σ ps f = 0.416, K = 3000 and α = 0.06. · Liver artery dataset: σ ps f = 0.502, K = 2000 and α = 0.35. · Lung vessel dataset: σ ps f = 0.502, K = 2200 and α = 0.05. (iii) The parameters for the proposed snake: The positive parameters λ1 , λ2 , λ3 and λ4 used in equation (19) give the relative strengths of the respective energy components and are chosen empirically. Since usually λ1 and λ2 should be maintained within 0.2 − 2.0, and λ1 and λ2 should be smaller than λ3 and λ4 , we use λ1 = λ2 = {0.1, 0.2, . . . , n − 0.1, n (n = 2.1)}, and λ3 = λ4 = {3, 4, . . . , m − 1 (m = 16)}. These parameters were modified through experiments on a small subset from different types of images [Abdominal CTA, Carotid CTA, lung CT, etc.], but for a certain type of images, we used the same set of parameters. These parameters were estimated for different types of images as follows: · Phantom dataset: λ1 = 0.2, λ2 = 2, λ3 = 4, λ4 = 10. · Liver artery dataset: λ1 = 0.5, λ2 = 2, λ3 = 5, λ4 = 12. · Lung vessel dataset: λ1 = 0.3, λ2 = 1.8, λ3 = 5, λ4 = 7. · Carotid artery dataset: λ1 = 0.4, λ2 = 2, λ3 = 5, λ4 = 9. 2) Accuracy Analysis Using an Acrylic Rod Phantom: Segmentation accuracy was determined using a Japanese Spine Phantom (JSP), an accepted standard in vessel densitometry [Fig. 3]. The JSP phantom was acquired using a Toshiba Aquilion CT scanner. Isotropic voxel imaging was performed with voxel size of 0.5 mm3 . It is composed of the straight acrylic rods and the curve straight acrylic rods. The diameters of the straight acrylic rods were 0.5, 0.75, 1.0, 1.5, 2.0 and 3.0 mm while those of the curve acrylic rods were 0.5, 0.75, 1.0, 1.5, and 2.0 mm. The straight rods were trimmed into seven sets and the angles of their axes relative to the x y-plane of the CT coordinate system were from 0° to 90° at an angle increment of 15°. We investigated segmentation accuracy for the radii of the acrylic rods. The horizontal and vertical acrylic rods were used for accuracy analysis. Starting from the center determined from a circle manually positioned around the rays at an increment of 30° was used to measure the true diameter parameters (0.5, 0.75, 1.0, 1.5, 2.0 and 3.0 mm). We compare our method with the Hernández-Hoyos [20], Li et al. [21] and Holtzman-Gazit et al. [24] methods. Fig. 4 indicates the mean estimation error and standard error of the estimation error for the Hernández-Hoyos’s, Li’s, Holtzman-Gazit’s and the proposed methods. For four methods, 120 rays were cast from the centroid of each of


2449

TABLE I JACCARD S IMILARITY C OEFFICIENT (JSC) VALUES C ALCULATED A MONG E VERY PAIR OF THE F OUR S EGMENTATIONS−O BTAINED F ROM E ACH OF C OMPUTATIONAL M ETHODS (H ERNÁNDEZ -H OYOS ’ S , L I ’ S , H OLTZMAN -G AZIT ’ S , O UR M ETHODS ) AND S LICE -B Y-S LICE S EGMENTATION A PPROACH B Y T WO H UMAN R ATERS − ON THE F OUR A RTERY V ESSEL S YSTEMS (A LL VALUES G IVEN IN %)

Fig. 4. The mean measurement error and standard error of the measurement error for the Hernández-Hoyos’s, Li’s, Holtzman-Gazit’s and the proposed methods.

Fig. 5. Slices from the arterial phase of a CT angiography (CTA) data set. The main problem in segmenting the liver arteries is that they are small or very poorly contrasted in the CTA data. Original image slice (left). Corresponding enhanced image (right).

the acrylic rods with 6 different known diameters. Estimates of the vessel wall locations were made independently for each ray. The graphs show the proposed method achieved the highest segmentation accuracy compared with the other three segmentation approaches. Differences in method measurement accuracy were assessed using paired t-tests. The t-test was applied to the measurements made on each ray for each of the 6 parameters. The results show that these methods are significantly different with the proposed approach showing greater accuracy. The cases for which the methods are significantly different are marked with an asterisk in Fig. 4. 3) Liver Artery Vessel Segmentation: In liver surgery planning, the liver arteries play an important role in determining optimal resections and in predicting blood supply to remain liver tissue. The main problem in segmenting the liver arteries is that they are small and very poorly contrasted in the CT data, see Fig. 5. Three abdominal CT angiography data sets were acquired on a Philips Brilliance 64 CT scanner. The acquired image volumes are typically of the size 512 × 512 × 400 voxels with a voxel size of about 0.68 × 0.68 × 0.6 mm. The liver regions in the images were roughly hand-segmented by a radiology specialist and used as a mask. Firstly, quantitative analysis on the segmentations produced with computational methods is conducted on the basis of Jaccard similarity coefficient (JSC) between the obtained segmentations and the truth. JSC is defined as the ratio of the size of the intersection volume to the size of the union volume of the two given segmentations [39]. Due to the fact that the truth segmentations of those clinical data sets do not exist, we have asked two human raters to compose segmentations by marking boundaries in slices. The JSC values are calculated among every pair of the three segmentations (obtained from our algorithm and the two human raters) per artery vessel system.

Fig. 6. Comparison of the three segmentation methods for three different liver artery systems. The first row, the second row and third row illustrate the three different liver arteries, respectively. First column shows our segmentation, the second column shows the Holtzman-Gazit’s segmentation, and the third column shows Li’s segmentation. White rectangles highlight the ability of the different methods to extract the end branches of the target liver arteries. As can be easily observed by naked eye, the small vessels appear to be segmented slightly thinner in our method than in the Holtzman-Gazit’s and Li’s methods.

Jaccard similarity coefficient (JSC) values are calculated among every pair of the three segmentations (obtained from each of computational methods and the two human raters) per artery vessel system. The results are tabulated in Table I. It is evident that our framework is capable of producing vessel segmentations that have level of variability similar to those obtained from the human raters. This study suggests that our algorithm can be good enough to replace the slice-by-slice manual segmentation approach, which indeed is very tedious and time consuming, to delineating vessels. We also show the sample visualizations of the results obtained from Li’s, Holtzman-Gazit’s and our methods [Fig. 6]. It is observed that, in comparison with our method, Holtzman-Gazit’s and Li’s methods gave some disconnected segmentation of the vessels in the regions that contain the low-contrast vessels and, more importantly, Li’s method produced an incomplete segmentation (i.e., some

2450

Fig. 7. Sample slices and carotid artery vessel boundary surface of the sub-volume from one clinical CT angiography data set. (a) One slice of the carotid artery. (b) Result of our method on (a). (c) Result of the Holtzman-Gazit’s method on (a). (d) Result of the Li’s method on (a). (e) Another slice of the carotid artery. (f) Result of our method on (e). (g) Result of the Holtzman-Gazit’s method on (e). (h) Result of the Li’s method on (e). (i)-(k): A portion of 3D carotid artery vessel obtained from our method, Holtzman-Gazit’s method and Li’s method, respectively. In figure, rectangles indicate the ability of the different methods to detect the highly curved and touching vessels.

Fig. 8. Original slices of part of the thorax CT. Rectangles indicate the presences of a nodule attached to a vessel, while arrows indicate the vessels with small caliber.

branches can not be detected) of the vessels in regions that contain the small branches. 4) Carotid Artery Vessel Segmentation: Several challenges of vessel extraction in the carotid artery images are illustrated by the images shown in Fig. 7. Highly curved [Fig. 7(a)] and touching [Fig. 7(e)] vessels are more difficult to handle. In such situations, a false detection may occur. We have applied the three methods to sub-volume of one carotid data set. The size of image volume is 256 × 256 × 256 voxels with isotropic voxel in 0.33 mm. As shown in Fig. 7 (rectangles), our method successfully stop the snake at the desired object boundary, while the Holtzman-Gazit et al. and Li et al. methods fail in doing this. Note that the zero-crossing method is used for obtaining the contour points on the vessel cross-section. The JSC values are calculated among every pair of the three segmentations (obtained from each of computational methods and the two human raters). The results are tabulated in Table I. Compared with the other three methods, our method achieved the highest segmentation accuracy. 5) Lung Vessel Segmentation: We also compare our algorithm with Holtzman-Gazit’s and Li’s methods for the extraction of the lung vessel tree. The data are acquired by


Fig. 9. Segmentation results of part of the thorax CT using three different segmentation methods. The first row and the second row illustrate the two different data sets on Fig. 8, respectively. First column shows our segmentation, the second column shows the Holtzman-Gazit’s segmentation, and the third column shows Li’s segmentation. In figure, rectangles indicate the presences of a nodule attached to a vessel; arrows highlight the ability of the different methods to detect the small vessel with low contrast; circles indicate the spatial discontinuity of Holtzman-Gazit’s method and incomplete segmentation of Li’s method.

the same Philips Brilliance 64 CT. The whole volume in this case had a size of 512×512×500 voxels with an in-slice pixel resolution of 0.68 × 0.68 mm2 and a 0.67 mm slice thickness. A visual inspection was performed to evaluate those vessels with pathologic structure (rectangles) and small caliber (arrows), as shown in Fig. 8. It can be seen that many of the small vessels not located by the other two methods were extracted by our method (arrows in Fig. 9). These small vessels were subsegmental peripheral vessels and they were also estimated to be smaller than 2 voxels in diameter. Although the Holtzman-Gazit’s method can detect some small vessels, spatial discontinuity is encountered (circles in Fig. 9). As shown in Fig. 8 (rectangles), the existences of a nodule (pathologic structure) attached to a vessel are particularly challenging for automatic vessel extraction. In Fig. 9 (rectangles), our method distinguishes vessels and nodules to reconstruct vessels while suppressing nodules and other noise. However, the segmentation result obtained with the Holtzman-Gazit’s method goes astray in regions that contain the nodules. Li’s method gives the similar results. B. Validation on Publicly Available Database In the second stage, our approach is tested on two publicly available databases and its results are compared with the recently published methods [9], [38]. 1) Vascular Tree Segmentation on Synthetic Datasets: Synthetic vascular trees are generated by using the VascuSynth Software [37]. The vascular trees of multiple complexities were generated by changing the number of terminal nodes from 5 to 1000 within a volume of 100 × 100 × 100 voxels


2451

TABLE II S EGMENTATION R ESULTS OF O UR M ETHOD ON 10 S YNTHETIC VASCULAR I MAGES (A LL VALUES G IVEN IN %)

TABLE III C OMPARISON OF THE C ETIN ’ S M ETHOD W ITH O UR M ETHOD IN THE

Fig. 10. Extracted vessel tree from the synthetic vascular dataset 1 in the 2011 VascuSynth Sample Data. Left: Extracted centerline. Right: Segmented surface. The first row shows the extracted results without Gaussian noise. The second row shows the extracted results with Gaussian noise of σ 2 = 20.

P RESENCE OF T WO L EVELS OF G AUSSIAN N OISE :

σ 2 = 20, σ 2 = 60 (A LL VALUES G IVEN IN %)

of 0.6 mm. 2011 VascuSynth Sample (3 datasets) and March 2013 VascuSynth Sample (10 datasets) are available on-line. We use the 10 datasets on the March 2013 to evaluate the performance of our algorithm; moreover, 3 datasets on the 2011 are used in the noise experiments to compare with the Cetin’s method [9]. We use four different quantitative measures for the synthetic R ), False Negatives validation as True Positives (TP = N BN∩N R N R −N B ∩N R B ∩N R ), False Positives (FP = N B −N ) and (FN = NR NR N B ∩N R Overlap Measure (OM = 2 N B +N R ) between the obtained vessel segmentation and the ground truth vessel segmentation. N R is the number of reference ground truth voxels, and N B is the number of voxels detected by our algorithm in synthetic data. Table II tabulates average TP, FN, FP, and OM rates on the 10 datasets. We also added Gaussian noise with two different variances to the 2011 VascuSynth Sample Data, and compared with Cetin’s method. Table III shows the comparison of the Cetin’s method with our method. The performance of our algorithm is slightly better than the Cetin’s method with the presence of low level Gaussian noise, and the performance of the two methods is close when the level of noise increases. Fig. 10 shows the visual results of our algorithm on the dataset 1 in the 2011 VascuSynth Sample Data.

2) Vessel Centerline Extraction on Datasets of Rotterdam Coronary Artery Evaluation Framework: We demonstrate the quantitative performance of our centerline extraction step on Rotterdam Coronary Artery Evaluation Framework [38], which facilitates testing over worldwide standard datasets and comparison to existing state-of-the-art algorithms in the literature. They provide 8 training datasets with ground truths and 24 test datasets. In our study, method that calculates the different constraints requires the accurate estimation value of vessel radius di (i = 0, 1, . . . , 5) to be known, which are obtained from equation (27) and (28). For this reason, the scanner pointspread function σ ps f is required to be estimated beforehand. To refine the estimate of σ ps f , it is necessary to perform the model-based optimization with the phantom images (shown in Fig. 3) and known tube radius. Thus, the constrained B-snake model cannot be used to perform the vessel segmentation. We only performed quantitative evaluation of centerline extraction step described in Section III over the Rotterdam Coronary Artery datasets, and comparison with Cetin’s method. The proposed centerline extraction algorithm was applied to all 32 dataset and the results were submitted to the organizers for extracting the evaluation parameters. The evaluation framework consists of three categories that describe the amount of user-interaction required: Automatic extraction method, methods with minimum user interaction (semi-automatic) and interactive extraction methods. Separate rankings are made by comparing with 23 centerline extraction algorithms, which contain 9 interactive methods, 5 semi-automatic methods, and 9 automatic methods. Among different evaluation measures, the overlap measure is with the most important for our algorithm. There are three overlap measures: Overlap (OV), Overlap until f ir sterr or (OF) and Overlap wi th > 1.5 mm vessel (OT). Our algorithm obtained a middle score with 92.6% OV, 53.3% OF, 93.7% OT, as shown in Table IV. Finally, the table shows statistics for the ’rank’ among the 23 methods. Our algorithm

2452


TABLE IV OVERLAP S CORES FOR THE E XTRACTED V ESSEL C ENTERLINES (C HALLENGE 1: AUTOMATIC M ETHODS , C HALLENGE 2: S EMI -AUTOMATIC M ETHODS , AND C HALLENGE 2: I NTERACTIVE M ETHODS )

ranks second in the automatic category. The processing time for the centerline extraction step is 1 min and 36 sec. Note that the core idea of Cetin’s method is to extract the centerline of the vessels by a new minimal path technique. This idea was applied to vascular segmentation by incorporating the spatial position of the centerline and its associated radius. Thus, Cetin’s method heavily depends on an accurate vessel centerline. In contrast, our algorithm does not need an accurate vessel centerline. Vessel centerline is used to provide the capturing path for the active contour model. Active contours can move and capture 2D vessel boundary on the cross-sections along the centerline path. Correction can be refined with the constrained active contours performing an accurate segmentation of the cross-section. VI. C ONCLUSION AND D ISCUSSION We have presented a new method to extract vessel axes and delineate vessel boundaries in a decoupled fashion. This method avoids disconnected and incomplete segmentation of the vessel in the problematic regions. It is particularly suitable for accurate segmentation of thin and low contrast vessels. Our method starts with the extraction of vessel axis from the original angiogram, followed is the shape and size constrained snake-based vessel boundary delineation on cross-sections defined by the extracted axis. The design philosophy of this method is to promote a more intuitive vessel boundary delineation on the vessel cross-sections, where the boundaries are mostly round, compared with those on either axial, coronal or sagittal image planes. The vessel boundary surface is then deformed under constrained movements on the cross-sections and is voxelized to produce the final vascular segmentation. A comparison has been made with two state-of-the-art vessel segmentation methods (active contour models). Quantitative results on phantom data indicate that our method is more

accurate than these classical models for small or thin vessels. Furthermore, experimental results on clinical data have shown that the vessels extracted by our method are continuous and complete in the problematic regions as compared with these classical models; our segmentation algorithm with cross-sections can delineate vessel boundaries that have level of variability similar to those obtained manually. In the Holtzman-Gazit’s and Li’s methods, for thin and low contrast vessels, the evolving surface splits, leaving two parts of the same blood vessel as disconnected components (Fig.9). To overcome this problem, we first exploit a multiscale vessel tracking method for obtaining the complete and continuous vessel axis. With the extracted axis, the contour points are then interpolated from neighboring contour point for handling missing contour points. Such decoupling of segmentation and vessel-axis-tracing avoids disconnected and incomplete segmentation of the vessel of interest. Once the vessel axis has been determined, the vessel boundary may be initialized. This can be accomplished by computing a radius function as the distance between the vessel axis and the geodesic path. Since the geodesic path lies somewhere close to the vessel boundary, its distance to the vessel axis provides an approximation of the vessel radius at every point along the axis. Our approach has some limitations. The method performs the vessel segmentation according to the assumption of approximately circular shape on the cross-sectional plane for vessels. Although the method holds for healthy vessels, it is inappropriate for some vessel disease, such as aneurysm or stenosis. Vessel disease affects the morphology and the appearance of vessels in medical images, which often breaks theses assumption. In such cases, our method will sometimes lead to inaccurate segmentation for the irregular shape.


Further research will focus on extending the method to irregular cross sections for stenosis detection and quantification.

Here, radial distance di from the vessel axis points is the number of samples in the model and the number of pixels on the cross-sectional plane. The vessel is surrounded by tissue of constant density α (0 ≤ α < 1). Assuming that Pr o(u, v) represent the 2D density profile of ideal vessel model [Fig. 1]. The scanning process is modeled as the linear convolution of the density function Pr o(u, v) with the scanner point spread function (PSF) Gau(u, v) [21] and written as ∞ ∞ Pr o(ξ, w)Gau(u − ξ, v − w)dξ dw f (u, v) = K −∞ −∞

(26) where K is a constant scale factor. Assuming that Gau(u, v) can be modeled as a symmetric 2D Gaussian with standard 2 2 2 deviation σ ps f : Gau(u, v) = e−(u +v )/2σ ps f . Scanner σ ps f can be obtained from our previous study [15]. Without loss of generality, consider the single ray running from u = 0 along the u axis in the positive u direction. For this single ray, equation (26) can be written as 2π di 2 2 2 e−(u +r −2ru)/2σ ps f r dr dθ f (u, 0) = K 0

2π

+ Kα

0

∞

e−(u

2 +r 2 −2ru)/2σ 2 ps f

r dr dθ

+ 2π K α

∞

e−(u

di

2 +r 2 )/2σ 2 ps f

2 l0 (r u/σ ps f )r dr

A. Position Constraint The deviations of the B-snake center from this centroid is penalized and, hence, the position constraint energy function E pos is E pos (C(s)) 1 1 1 1 = ( u(s)ds − u c )2 + ( v(s)ds − v c )2 (30) 2 0 2 0 where C(s) = [u(s) v(s)] is the B-Snake curve, Cc = [u c v c ] is the centroid results from the obtained constraint points bi . The centroid Cc can be estimated as follows: 1 1 di cosθi , v c = di si nθi (31) uc = 6 6 i

i

where [di cosθi di si nθi ] is the coordinate of the constraint point bi on cross-section. B. Shape and Size Constraint

(28)

The constraint energy can be obtained by summing along the snake curve. Let E shape−size specify the energy term for shape and size constraint, E shape−size is expressed as 1 1 (u(s) − u c − Rcos(2πs))2 ds E shape−size (C(s)) = 2 0 1 1 + (v(s) − v c − Rsi n(2πs))2 ds 2 0 1 R= || pi − Cc || (32) 6

where l0 (u) is the modified Bessel function of order zero [34]. Equation (28) allows us to predict the shape of the signal intensity profile at any contour location along the ray for a given vessel geometry. On the vessel cross-section, a nonlinear optimization technique is used to match the actual CT intensity profile with the predicted profile model. Using the model of CT image process, we can obtain the one-dimensional (1D) profile of modeled intensity, f (u, 0), along ray [Fig. 1(a)]. Similarly, the 1D profile of actual intensity, I (u), is derived from CT image along ray. Reconstruction of the 1D profile of actual intensity, I (u), is performed at sub-voxel resolution by using a bilinear interpolation. Let u j be the j th voxel position along ray in the actual image. The observed profile is sampled at M discrete points along ray in the actual CT image, di are estimated by finding the values of di minimizing M−1 E(di ) = (I (u j ) − f (u j , 0))2 (29) j =0

Active contours can move and capture 2D vessel boundary on the cross-sections with position, shape and size constraints rather than freely in a 3D space. This helps to avoid shrinkage of the model and bunching of the boundary points at regions that two adjacent objects have similar gray values. Different constraints are given below.

(27)

di

Simplifying, we obtain di 2 −(u 2 +r 2 )/2σ ps f l (r u/σ 2 )r dr e f (u, 0) = 2π K 0 ps f 0

We employ the Levenberg-Marquardt algorithm to solve this non-linear least square problem. A PPENDIX B VARIOUS C ONSTRAINTS

A PPENDIX A D ETERMINING R ADIAL D ISTANCE di

0

2453

i

where Cc = [u c v c ] is the estimated centroid as defined in equation (31), R is the average radius of constraint points bi . The constraint energy term penalizes the deviation of the B-snake from the molded circle, which has a mean radius R and a center Cc . C. Sampling Constraint The sampling constraint energy makes a curve point maintain equal distance from its immediate left and right neighboring points on the snake curve. For the application at hand, the curve points maintain a circle shape, the j th and ( j − 1)th contour points follow the relationships as 2π j 2π( j − 1) ) − Rcos( ) u j − u j −1 = Rcos( m m 2π j 2π( j − 1) v j − v j −1 = Rsi n( ) − Rsi n( ) (33) m m

2454


where R is the average radius of the constraint points, as already defined in equation (32). Also, we note C(sm ) = C(s0 ), that is u m = u 0 and v m = v 0 . We now introduce the following energy function,

F. Sampling Constraint E sampling We incorporate the sampling energy term in snake evolution equation. The sampling energy equation (34) can be written in the following matrix-vector form

1 = [(u j − u j −1 − d uj )2 + (v j − v j −1 − d vj )2 ] 2 m

E sampling

E sampling =

j =1

(34) where d uj and d vj are shorthand for the right-hand sides of equation (33) and are written as 2π j 2π( j − 1) ) + Rcos( ) m m 2π j 2π( j − 1) d vj = Rsi n( ) + Rsi n( ) m m

d uj = Rcos(

(35)

A PPENDIX C C ONSTRAINTS I NCORPORATION Here, we introduce the Euler equations of the constraints and incorporate the constraints into the B-snake. D. Position Constraint E pos The energy functional expressed in equation (30) is indeed amenable to the analysis in a continuous framework given the assumptions of the continuity on the curve, and all the integrands defined in equation (30) are continuous and integrable. To obtain the snake position that minimizes the constraint energy equation (30), we use the calculus of variations and arrive at the following Euler equation ∇ E pos = (

∂ ∂ E pos , E pos ) = 0 ∂u ∂v

or, equivalently, as 1 u(s)ds − u c = 0 and 0

1

(36)

As we proposed, the snake curve C(s) was discretized as C(si ) = [u i v i ]. We can write them out in the matrix-vector form and compute the gradient for the position constraint energy

[u, . . . , u]T ,

where u = [u 1 , u 2 , . . . , u m ]T and v = [v 1 , v 2 , . . . , v m ]T are the vectors of discrete snake curve point locations du = [d1u , d2u , . . . , dmu ]T , dv = [d1v , d2v , . . . , dmv ]T , d uj and d vj are defined in are defined in equation (35) and H , G are mby-m matrices as follows: ⎡ ⎤ 2 −1 −1 ⎢ −1 2 ⎥ ⎢ ⎥ ⎢ ⎥ . . . ⎢ ⎥ .. .. .. H =⎢ ⎥ ⎢ ⎥ 2 −1 ⎦ ⎣ −1 −1 2 ⎤ ⎡ 1 −1 −1 ⎥ ⎢ −1 1 ⎥ ⎢ ⎥ ⎢ . . . ⎥ ⎢ .. .. .. G=⎢ ⎥ ⎥ ⎢ 1 −1 ⎦ ⎣ −1 −1 1 The energy E sampling is in the quadratic form, so one can now easily obtain the gradient of the energy functional as follows: ∂ ∂ E sampling ) ∇ E sampling = ( E sampling , ∂u ∂v = (H u − Gdu, H v − Gdv) (40) Since there is no internal force, we just add the vectors λ5 Gdu and λ5 Gdv to E u and E v , respectively.

v(s)ds − v c = 0

0

∇ E pos = (u − uc , v − vc )

(37)

T v = [v,. . . , v]T , uc = [u c, . . . , uc ] , m 1 1 m u = m j =1 u j , v = m j =1 v j .

where u = vc = [v c , . . . , v c ]T and We incorporate the position constraint into the B-snake evolution equation (24) by adding λ2 (u c − u) and λ2 (v c − v) each element of E u and E v , respectively.

ACKNOWLEDGMENT The authors would like to thank the anonymous reviewers for their valuable comments and help suggestions that greatly improved the paper’s quality. R EFERENCES [1] [2] [3]

E. Shape and Size Constraint E shape−size For the shape and size constraints energy equation (32), we obtained two Euler equations

[4] [5]

u(s) − u c − Rcos(2πs) = 0 v(s) − v c − Rsi n(2πs) = 0

1 T 1 u H u + v T H v − uT Gu − v T Gv 2 2 1 1 (39) + (du)T (du) + (dv)T (dv) 2 2

(38)

Similarly, λ3 (u c + Rcos(2π j/m)) and λ3 (v c + Rsi n(2π j/m)) are added to the j th element of E u and E v .

[6]

C. Kirbas and F. Quek, “A review of vessel extraction techniques and algorithms,” ACM Comput. Surv., vol. 36, no. 2, pp. 81–121, Jun. 2004. D. Lesage, E. D. Angelini, I. Bloch, and G. Funka-Lea, “A review of 3D vessel lumen segmentation techniques: Models, features and extraction schemes,” Med. Imag. Anal., vol. 13, no. 6, pp. 819–845, Dec. 2009. Y. Sato et al., “Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images,” Med. Image Anal., vol. 2, no. 2, pp. 143–168, Dec. 1998. J. Chen and A. A. Amini, “Quantifying 3D vascular structures in MRA images using hybrid PDE and geometric deformable models,” IEEE Trans. Med. Imag., vol. 23, no. 10, pp. 1251–1262, Oct. 2004. Y. Fridman, S. M. Pizer, S. Aylward, and E. Bullitt, “Segmenting 3D branching tubular structures using cores,” in Proc. Med. Image Comput. Comput.-Assist. Intervent. (MICCAI), Montreal, QC, Canada, Nov. 2003, pp. 570–577. J. A. Tyrrell et al., “Robust 3D modeling of vasculature imagery using superellipsoids,” IEEE Trans. Med. Imag., vol. 26, no. 2, pp. 223–237, Feb. 2007.


[7] [8] [9] [10] [11]

[12] [13] [14] [15]

[16] [17] [18] [19]

[20]

[21] [22] [23] [24] [25] [26] [27] [28] [29]

M. A. Gülsün and H. Tek, “Robust vessel tree modeling,” in Medical Image Computing and Computer-Assisted Intervention—MICCAI. Berlin, Germany: Springer-Verlag, 2008, pp. 602–611. H. Li and A. Yezzi, “Vessels as 4D curves: Global minimal 4D paths to extract 3D tubular surfaces and centerlines,” IEEE Trans. Med. Imag., vol. 26, no. 9, pp. 1213–1223, Sep. 2007. S. Cetin, A. Demir, A. Yezzi, M. Degertekin, and G. Unal, “Vessel tractography using an intensity based tensor model with branch detection,” IEEE Trans. Med. Imag., vol. 32, no. 2, pp. 348–363, Feb. 2013. T. F. Cootes, C. J. Taylor, D. H. Cooper, and J. Graham, “Active shape models—Their training and application,” Comput. Vis. Image Understand., vol. 61, no. 1, pp. 38–59, Jan. 1995. J. P. Williams, J. K. Johnstone, and L. B. Wolff, “Rational discrete generalized cylinders and their application to shape recovery in medical images,” in Proc. IEEE Comput. Soc. Conf. Comput. Vis. Pattern Recognit., Jun. 1997, pp. 387–392. R. Li and S. Ourselin, “Accurate curvilinear modeling for precise measurements of tubular structures,” in Proc. Int. Conf. Digit. Image Comput., Techn. Appl., 2003, pp. 243–252. J. Ponce, D. Chelberg, and W. B. Mann, “Invariant properties of straight homogeneous generalized cylinders and their contours,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 11, no. 9, pp. 951–966, Sep. 1989. K. Rohr and S. Worz, “High-precision localization and quantification of 3D tubular structures,” in Proc. 3rd IEEE Int. Symp. Biomed. Imag., Nano Macro, Apr. 2006, pp. 1160–1163. Y. Sato, S. Yamamoto, and S. Tamura, “Accurate quantification of small-diameter tubular structures in isotropic CT volume data based on multiscale line filter responses,” in Proc. 7th Int. Conf. Med. Image Comput. Comput.-Assist. Intervent. (MICCAI), 2004, pp. 508–515. J. M. Reinhardt, N. D’Souza, and E. A. Hoffman, “Accurate measurement of intrathoracic airways,” IEEE Trans. Med. Imag., vol. 16, no. 6, pp. 820–827, Dec. 1997. C. Steger, “An unbiased detector of curvilinear structures,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 20, no. 2, pp. 113–125, Feb. 1998. M. Kass, A. Witkin, and D. Terzopoulos, “Snakes: Active contour models,” Int. J. Comput. Vis., vol. 1, no. 4, pp. 321–331, Jan. 1988. A. F. Frangi, W. J. Niessen, R. M. Hoogeveen, T. van Walsum, and M. A. Viergever, “Model-based quantitation of 3D magnetic resonance angiographic images,” IEEE Trans. Med. Imag., vol. 18, no. 10, pp. 946–956, Oct. 1999. M. Hernández-Hoyos, “Segmentation anisotrope 3D pour la quantification en imagerie vasculaire par rTsonance magnTtique (in French),” Ph.D. dissertation, École Doctorale des Sciences de l’IngTnieur de Lyon, CREATIS UMR-CNRS 5515, Lyon, France, 2002. C. Li, C.-Y. Kao, J. C. Gore, and Z. Ding, “Minimization of regionscalable fitting energy for image segmentation,” IEEE Trans. Image Process., vol. 17, no. 10, pp. 1940–1949, Oct. 2008. V. Caselles, R. Kimmel, and G. Sapiro, “Geodesic active contours,” Int. J. Comput. Vis., vol. 22, no. 1, pp. 61–79, Feb. 1997. A. Vasilevskiy and K. Siddiqi, “Flux maximizing geometric flows,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 24, no. 12, pp. 1565–1578, Dec. 2002. M. Holtzman-Gazit, R. Kimmel, N. Peled, and D. Goldsher, “Segmentation of thin structures in volumetric medical images,” IEEE Trans. Image Process., vol. 15, no. 2, pp. 354–363, Feb. 2006. J. Tang and S. T. Acton, “Vessel boundary tracking for intravital microscopy via multiscale gradient vector flow snakes,” IEEE Trans. Biomed. Eng., vol. 51, no. 2, pp. 316–324, Feb. 2004. J. Lee, P. Beighley, E. Ritman, and N. Smith, “Automatic segmentation of 3D micro-CT coronary vascular images,” Med. Image Anal., vol. 11, no. 6, pp. 630–647, Dec. 2007. K. F. Lai and R. T. Chin, “Deformable contours: Modeling and extraction,” in Proc. IEEE Comput. Soc. Conf. Comput. Vis. Pattern Recognit., Jun. 1994, pp. 601–608. N. Ray, S. T. Acton, and K. Ley, “Tracking leukocytes in vivo with shape and size constrained active contours,” IEEE Trans. Image Process., vol. 21, no. 10, pp. 1222–1235, Oct. 2002. C. Xu and J. L. Prince, “Snakes, shapes, and gradient vector flow,” IEEE Trans. Image Process., vol. 7, no. 3, pp. 359–369, Mar. 1998.

2455

[30] W. Schroeder, K. Martin, and B. Lorensen, Visualization Toolkit: An Object-Oriented Approach to 3D Graphics, 4th ed. New York, NY, USA: Kitware, Dec. 2006. [31] R. M. Haralick, “Digital step edges from zero crossing of second directional derivatives,” IEEE Trans. Pattern Anal. Mach. Intell., vol. PAMI-6, no. 1, pp. 58–68, Jan. 1984. [32] J. Canny, “A computational approach to edge detection,” IEEE Trans. Pattern Anal. Mach. Intell., vol. PAMI-8, no. 6, pp. 679–698, Nov. 1986. [33] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions. New York, NY, USA: Dover, 1972. [34] M. J. Davies, “A macro and micro view of coronary vascular insult in ischemic heart disease,” Circulation, vol. 82, no. 3, pp. 38–46, 1990. [35] Y. Wang and E. K. Teoh, “Object contour extraction using adaptive B-snake model,” J. Math. Imag. Vis., vol. 24, no. 3, pp. 295–306, May 2006. [36] Y. Wang, E. K. Teoh, and D. Shen, “Lane detection and tracking using B-snake,” Image Vis. Comput., vol. 22, no. 4, pp. 269–280, Apr. 2004. [37] G. Hamarneh and P. Jassi, “VascuSynth: Simulating vascular trees for generating volumetric image data with ground-truth segmentation and tree analysis,” Comput. Med. Imag. Graph., vol. 34, no. 8, pp. 605–616, Dec. 2010. [38] M. Schaap et al., “Standardized evaluation methodology and reference database for evaluating coronary artery centerline extraction algorithms,” Med. Image Anal., vol. 13, no. 5, pp. 701–714, Oct. 2009. [39] K. Van Leemput, F. Maes, F. Bello, D. Vandermeulen, A. Colchester, and P. Suetens, “Automated segmentation of MS lesions from multichannel MR images,” in Medical Image Computing and ComputerAssisted Intervention (Lecture Notes in Computer Science), vol. 1679. Berlin, Germany: Springer-Verlag, 1999, pp. 11–21.

Yuanzhi Cheng is currently a Professor with the School of Computer Science and Technology, Harbin Institute of Technology. He has authored over 30 papers in scientific journals. His research interests are concentrated on pattern recognition, image processing, and computer-assisted surgical system.

Xin Hu is currently pursuing the degree with the Department of Computer Science, Harbin Institute of Technology, China.

Ji Wang is currently pursuing the degree with the Department of Computer Science, Chiba Institute of Technology, Japan.

Yadong Wang is currently a Professor with the School of Computer Science and Technology, Harbin Institute of Technology. He has authored over 60 papers in scientific journals. His research interests include pattern recognition, machine learning, and bioinformatics.

Shinichi Tamura (F’09) is currently a Professor with the Center for Advanced Medical Engineering and Informatics, Osaka University. He has authored over 260 papers in scientific journals. His current research activities include medical image analysis and its applications. He received several awards from journals, including Pattern Recognition and Investigative Radiology.

NUCLEI SEGMENTATION VIA SPARSITY CONSTRAINED CONVOLUTIONAL REGRESSION.

CONSTRAINED SPECTRAL CLUSTERING FOR IMAGE SEGMENTATION.

Iterative Vessel Segmentation of Fundus Images.

Robust vessel segmentation in fundus images.

Unsupervised Retinal Vessel Segmentation Using Combined Filters.

3-D Lung Segmentation by Incremental Constrained Nonnegative Matrix Factorization.

Shape-Constrained Multi-Atlas Segmentation of Spleen in CT.

Left ventricle segmentation by dynamic shape constrained random walks.

A spatial shape constrained clustering method for mammographic mass segmentation.

Accurate de novo design of hyperstable constrained peptides.

Blood Vessel Segmentation of Fundus Images by Major Vessel Extraction and Subimage Classification.

Robust Retinal Blood Vessel Segmentation Based on Reinforcement Local Descriptions.

A vessel active contour model for vascular segmentation.

Fast and accurate circle detection using gradient-direction-based segmentation.

Adaptively constrained convex optimization for accurate fiber orientation estimation with high order spherical harmonics.

Accurate cell segmentation in microscopy images using membrane patterns.

Robust, accurate and fast automatic segmentation of the spinal cord.

Medially constrained deformable modeling for segmentation of branching medial structures: Application to aortic valve segmentation and morphometry.

Correction: Retinal Vessel Segmentation: an Efficient Graph Cut Approach with Retinex and Local Phase.

Active Contours Driven by Multi-Feature Gaussian Distribution Fitting Energy with Application to Vessel Segmentation.

Leveraging Multiscale Hessian-Based Enhancement With a Novel Exudate Inpainting Technique for Retinal Vessel Segmentation.

Joint 3-D vessel segmentation and centerline extraction using oblique Hough forests with steerable filters.

Retinal vessel segmentation: an efficient graph cut approach with retinex and local phase.

Accurate white matter lesion segmentation by k nearest neighbor classification with tissue type priors (kNN-TTPs).