 Research
 Open Access
 Published:
3D vasculature segmentation using localized hybrid levelset method
BioMedical Engineering OnLinevolume 13, Article number: 169 (2014)
Abstract
Background
Intensity inhomogeneity occurs in many medical images, especially in vessel images. Overcoming the difficulty due to image inhomogeneity is crucial for the segmentation of vessel image.
Methods
This paper proposes a localized hybrid levelset method for the segmentation of 3D vessel image. The proposed method integrates both local region information and boundary information for vessel segmentation, which is essential for the accurate extraction of tiny vessel structures. The local intensity information is firstly embedded into a regionbased contour model, and then incorporated into the levelset formulation of the geodesic active contour model. Compared with the preset global threshold based method, the use of automatically calculated local thresholds enables the extraction of the local image information, which is essential for the segmentation of vessel images.
Results
Experiments carried out on the segmentation of 3D vessel images demonstrate the strengths of using locally specified dynamic thresholds in our levelset method. Furthermore, both qualitative comparison and quantitative validations have been performed to evaluate the effectiveness of our proposed model.
Conclusions
Experimental results and validations demonstrate that our proposed model can achieve more promising segmentation results than the original hybrid method does.
Background
The reconstruction of blood vessel plays an important role for many clinical tasks such as diagnosis of vessel diseases, surgery planning and blood flow simulation [1]. The first task to vessel reconstruction is to identify the vessel region from initial vessel images. This is usually achieved by using the technique of segmentation, which plays an important role in medical image processing. Although many segmentation techniques have been proposed for various image segmentation problems, blood vessel image segmentation still remains a challenging task [2], due to the high complexity of vessel branching and thinning geometry as well as the decrease in image contrast from the root of the vessel to its thin branches [1].
Among all the segmentation techniques, the deformable model (i.e. active contour) [3] is one of the most popular approaches during the past two decades. The basic idea of deformable model method is to evolve the initial curve/surface to the boundaries of target objects driven by the combination of internal forces determined by the geometry of the evolving curve and the external forces induced from the image [4]. Depending on the way the underlying curved contour is represented, the deformable active contour can be implemented either parametrically or implicitly [5]. In parametric models, active contours are explicitly represented in parametric forms in a Lagrangian framework [6]. A multitude of powerful methods based on parametrical deformable model have been proposed for the segmentation of medical image [7]. Klein et al. [8] presented a method for the extraction of vessel boundaries using deformable surface models represented by Bspline functions. However, the parameterized deformable models are difficult to handle the complex topological changes for the whole vasculatures during evolution, as reparameterization is a tough problem to be solved for parametric model [9]. On the other hand, active contours can be implicitly represented as levelset [10] function, which is embedded in higher dimensional spaces and evolves according to the Partial Differential Equation (PDE) in an Eulerian framework [6]. Compared with parametric model, levelset model is easier to perform shape operations, and adapt freely into complex topologies of objects. Therefore, various levelset based active models have been developed and applied to image segmentation. Some of the wellknown models are the geodesic active contour model [11] that utilizes image gradient to stop evolving contours on the object boundaries, and ChanVese model [12] that solves the leakage problem by using region information of the target boundary for segmentation. There are also several localized active models [13, 14] introducing the local image fitting energy to extract the local image information been presented to improve the performance of the active contours. Furthermore, some techniques combining edge and region information [4, 15, 16] have been developed and applied to segmentation of medical imaging.
In the area of 3D vessel image segmentation, Lorigo et al. [17] proposed the ’CURVES’ algorithm to extract vessels by using a geodesic active contour based on the codimension two level set method [18]. Nain et al. [1] proposed a regionbased geometric deformable model for MRA segmentation, which is driven by the combination of image statistics and shape information. Yan and Kassim [19] presented a segmentation scheme for the extraction of vasculature from MRA images by using capillary active contour. In [20], the vessel contours are segmented from the enhanced dataset by means of a level set evolution, which relies on the image intensity statistics to either expand or shrink the evolving contour.In [21], a shape functional is used to regularize the ‘CURVES’ and flux maximizing flow functional [22], which is claimed to be a promising tool to improve the efficiency of both techniques. Law and Chung [23] presented a vessel extraction approach based on a weighted local variance based edge detection scheme, which is effective for the segmentation of vessel image with low contrast edges. Recently, Wang and Jiang presented a nonparametric shape constrained algorithm for segmentation of coronary arteries in computed tomography images within the framework of active contours [24]. However, as most of these active contour models only take the edge information or region information of the image into account, they may have difficulty in correctly extracting tiny vessels from images, due to the low image contrast of the thin vessel branches.
In this paper, a localized hybrid levelset method for the segmentation of 3D vessel image is proposed. Our technique is inspired by Zhang et al.’s hybrid levelset model [4]. The hybrid model integrates both region and boundary information for the segmentation of medical images, in which the boundary information can help to detect the precise location of the target object and the region information can help to prevent the boundary leakage problem. However, the original hybrid model utilized a preset global threshold indicating the lower bound of target object to specify the term concerning regional information, which is not quite appropriate, especially for medical images with intensity inhomogeneity. In fact, intensity inhomogeneity occurs in many medical images [13], especially in vessel images, as the vessel branching is highly complex and the image contrast decreases from the root of the vessel to its thin branches [1]. On the other hand, our proposed localized hybrid levelset model utilizes the automatically calculated dynamic thresholds to indicate the lower bound of target object to specify the term concerning regional information instead of using a preset global threshold, which is crucial for the segmentation of vessel images [15]. In our proposed technique, the local intensity information is firstly embedded into a regionbased contour model, and then incorporated into the levelset formulation of the geodesic active contour model. Compared with the global threshold based method, the use of locally specified dynamic thresholds enables the extraction of the local image information, which is essential for the segmentation of vessel images. Experimental results on 3D medical datasets and validations are presented to demonstrate the strengths of our localized hybrid levelsetmethod.
Methods
Theory
In the proposed localized hybrid model, the local intensity information is firstly embedded into a regionbased contour model, and then incorporated into the levelset formulation of the geodesic active contour model. Compared with original hybrid levelset model [4], our proposed model calculates the locally specified dynamic thresholds to indicate the lower bound of target object, instead of using a preset global threshold to specify the term concerning regional information.Consider a given image I:Ω→R, in which Ω⊂R^{n} (n=2 or 3) is the image domain, and let ϕ:Ω→R be a contour in the image domain. The proposed energy functional to be minimised is defined as:
where u∈Ω, and α and β are the preset weights to balance the two terms. g:[0,∞]→R^{+} is the decreasing function, such as g(h)=1/(1+c h^{2}) with c controlling the scope, and H(ϕ) is the Heaviside function defined as:
μ(u)is the automatically calculated localized thresholds indicating the lower bound of target object. The definition of μ(u)is as below:
where k∈[0.5,1] is an adjustment coefficient, and K _{ σ } is the Gaussian kernel function with a localization property, such as ${K}_{\sigma}\left(\mathbf{u}\right)=\frac{1}{2\pi {\sigma}^{2}}{e}^{{\left\mathbf{u}\right}^{2}/2{\sigma}^{2}}$ with a scale parameter σ. The coefficient k is effective for preventing the active contours stopping evolving inside the target areas before reaching the boundary, which makes our proposed model more robust for the segmentation of 3D medical dataset. The first term of functional (1) defines the region term in our localized hybrid levelset model. The second term is equal to the geodesic active contour functional, which is the same as the second term of the original hybrid model. In addition to keeping the advantages of original hybrid model, the localized hybrid model is able to extract more accurate local image information by utilizing the locally specified dynamic thresholds μ(u) to indicate the lower bound of target object, which is essential for segmenting images with intensity inhomogeneity.
The associated PDE can be derived from the gradient descent flow [13] applied to functional (4)
Implementation
The original Heaviside function is discontinuous and is unable to achieve smooth transition at the boundary. In practice, this problem is generally solved using a kind of smooth Heaviside function, which is usually called smooth step function or smooth unit step function. Depending on the problem of applications, different smooth step functions have been introduced to approximate smoothly the Heaviside function [25]. More recently, a piecewise polynomial smooth unit step functions is developed and used for geometric shape design, which can be considered as a generalization of the Heaviside unit step function [26].
In our implementation, we adopt the smooth function H _{ ε }(s) introduced in [12]
and the corresponding Dirac function δ _{ ε }(s) is defined as
Then, H in Equation (1) is replaced with H _{ ε }, and δ in Equation (4) is replaced with δ _{ ε }. In addition, in order to avoid problems of developing shocks, sharp or flat shapes during evolution, it is essential to reinitialize ϕ to be a signed distance function (SDF) [27]. And if ϕ is a SDF, then ∇ϕ=1. Therefore, Equation (4) can be written as
In the proposed model, the temporal partial derivative ∂ ϕ / ∂ t is approximated by the forward difference scheme. Then, the approximation of Equation (4) can be expressed as [6]
where ϕ^{(k+1)} and ϕ^{(k)} denote the embedding function ϕ in the (k+1)^{th} and (k)^{th} iterations respectively, and Δ t is the predefined time step. The difference equation (8) can also be rewritten as the following iteration
Finally, the principal steps in the original hybrid model [4] is adopted to update from ϕ^{(k)} to ϕ^{(k+1)}:

1.
Reinitialize ϕ ^{(k)} to be SDF;

2.
Calculate the localized thresholds μ(u) automatically according to Equation 3;

3.
Compute the local region energy functional δ _{ ε }(ϕ ^{(k)})(I−μ(u));

4.
Update ϕ ^{(k)} to ϕ ^{(k)} ^{′} using ϕ ^{(k)} ^{′}=ϕ ^{(k)}+Δ t α δ _{ ε }(ϕ ^{(k)})(I−μ(u));

5.
Reinitialize ϕ ^{(k)} ^{′} to be SDF;

6.
Update ϕ ^{(k)} ^{′} to obtain ϕ ^{(k+1)} via ϕ ^{(k+1)}=ϕ ^{(k)} ^{′}+Δ t β δ _{ ε }(ϕ ^{(k)})d i v(g(∇I)∇ϕ ^{(k)}).
Results and discussion
Results of our localized hybrid model
The proposed localized hybrid model has been applied to more than ten 3D medical datasets for the segmentation of vasculatures. The datasets are supplied by Intelligent Bioinformatics Systems Division, Institute of Automation, the Chinese Academy of Sciences, in the format of DICOM (Digital Imaging and Communications in Medicine). In all these experiments, the same parameters are used, namely, Δ t=4.0, α=0.01, β=0.5, σ=3.0, ε=1.0, and k=0.9. In the following, we present some typical segmentation results as well as the Maximum Intensity Projections (MIPs) of the corresponding original datasets. MIPs are very useful to evaluate current 3D angiography images since the overall shapes and paths of the vessel can be clearly visualized [28].
The first example is the segmentation of cerebral vasculatures for 3D MRA images with a resolution of 352×448×114 and spacing of 0.49 m m×0.49 m m×0.80 m m. Figure 1 shows the segmentation results for some 2D slices of the dataset. As mentioned before, Magnetic Resonance (MR) images are typically intensity inhomogeneous. As can be seen from Figure 1, the intensity of the middle region is quite lower than that of other regions of the vessel. Owning to the utilization of locally specified dynamic thresholds, the proposed method can successfully depict the vessels contours (indicated by blue curves) in these images. Figure 2 shows the MIP of this dataset, and Figure 3 presents the segmentation result of whole complex cerebral vasculature successfully extracted from the 3D MRA dataset.
The second example is the segmentation of the common iliac artery for 3D MRA images with a resolution of 384×384×72 and spacing of 1.0 m m×1.0 m m×1.5 m m. The MIP of this dataset is shown in Figure 4, while the segmentation result using our localized hybrid technique is shown in Figure 5. Although the raw dataset has great intensity inhomogeneity, the whole iliac artery is still successfully extracted.
Another example is the segmentation of peripheral artery for the 3D CTA images with a resolution of 512×512×240 and spacing of 0.83 m m×0.83 m m×1.0 m m. As can been seen from the MIP of this dataset (Figure 6), the visualization of vessel structures is greatly distracted by bone areas. However, our proposed model can still successfully extract the whole peripheral artery structures from the raw images despite of the distraction of bone structures (see Figure 7).
Comparison with the original hybrid levelset model
In this section, we investigate why the original hybrid model may fail in the segmentation of real medical data with intensity inhomogeneity whereas our localized model works. According to our prior knowledge, the intensity (i.e. gray level) of most vessels in the first dataset (i.e. cerebral MRA images) is larger than 200. However, the intensity of some small vessels and their branches can be varied, for some of them can even be lower than 100. Therefore, for the original model, it is quite difficult to predefine the global threshold μ indicating the lower bound of gray level of target vessel. If the global threshold μ were set to be 200, it would be unable to extract the small vasculatures with intensity smaller than 200; if μ were set to be 100, the segmentation result could include some irrelevant objects, which would lead to inaccurate extraction of vasculatures.
For the purpose of comparison, we use the same set of parameter values (please refer to the beginning of this section) for the common parameters of the original model and our localized model. In addition, we also set the same seed contours for both models. As can be seen from Figure 8(a), it is very clear that our method has the ability to extract small branches with satisfactory results. Figure 8(b) shows the segmented result of the original model with μ=200. In this case, it is unable to extract small branches with an intensity less than 200, which just provides us with the unconnected vessel tree. On the other hand, when μ is set as 100, the original hybrid model can produce a more connected vessel tree, but some of the individual branches become much thicker than expected, which can be seen clearly in Figure 8(c). In other words, when the value of μ is set too low in the original hybrid model, it may also lead to inaccurate segmentation of vasculatures.
Moreover, the original hybrid model may be still unable to extract thin vessels in low contrast situations even though μ is set to be very low. As can be seen from Figure 9,(a) is a detailed look at the MIP of the 3D dataset, and the rectangular region in (b) indicates the successful extraction of thin vessels with low contrast using our localized model. Figure 9(c) shows the segmentation result using the original hybrid model with μ=100, which still cannot extract thin vessels.
In a word, this qualitative comparison indicates that our localized hybrid model can achieve a more satisfactory segmentation result of the vascular structures than the original model with different μ s.
Validations
In this study, a test dataset as shown in Figure 10 with the resolution of 100×70×106 voxels is used to evaluate the performance of the proposed localized hybrid model. In this dataset, an actual vessel structure reconstructed from the real MRA images is chosen as the known vessel pattern for validation, which is shown in Figure 10(a). The maximum intensity value in the original dataset is 500, and the intensity value of vessel boundary is 250. Then the original dataset is corrupted by several zeromean Gaussian noise levels. Furthermore, the intensity value of part of the vessel branches are also modified to simulate the situation of intensity inhomogeneity, which occurs in many medical images [13]. Figure 10(b)(d) demonstrate the MIP images of the dataset corrupted by different Gaussian noises with variances of 0.001, 0.01 and 0.015, and parts of the vassel branches look like disconnected due to the effect of intensity inhomogeneity.In these experiments, for the purpose of comparison, we set the same seed contours for both the original and our localized models. Figure 11 shows the the segmentation results obtained from noisy dataset with variance of 0.001. Ellipse areas of panel (b) demonstrate that the original model is unable to extract the connected vessel branches due to the effect of intensity inhomogeneity, even though the global threshold indicating the lower bound of intensity value of target vessel is set very low. On the other hand, our proposed model can solve this problem by calculating the local threshold automatically according to the local image information, which can extract the connected vessel branches shown in ellipse areas of panel (a).
Besides the visual inspection, we calculate the segmentation errors using the Dice similarity coefficient [29] for quantitative validation:
where A and B are the target and obtained segmentation sets, and n is the voxel number. The smaller value of SE means the more accuracy of the segmentation result.
Table 1 illustrates the segmentation errors of ChanVese (CV) model [12], RegionScalable Fitting Energy(RSFE) based method [14], Nonparametric Shape Prior Constrained Active Contour (NSPCAC) Model [24], Original Hybrid Model [4], and the proposed Localized Hybrid Model for different Gaussian noise variance levels applied on the test dataset. The other four methods that we compared with are typical: CV is a classical region based active contours model; RSFE is a wellknown segmentation method that utilizes the local image information; NSPCAC is a recently developed method for vessel image segmentation; The Original Hybrid Model is the most relevant method to the proposed Localized Hybrid Model. As have been shown in the table, the SE values of our proposed method are smaller than that of CV model, RSFE based method, NSPCAC model, Original Hybrid Model for each different variance levels. That is, by keeping the advantages of hybrid levelset model as well as using locally specified dynamic thresholds, our method can achieve more accurate segmentation results than other methods.
Conclusions
Vessel image segmentation is an essential step for many clinical tasks such as diagnosis of vessel diseases, surgery planning and blood flow simulation. However, due to the high complexity of vessel branching and thinning geometry, blood vessel image segmentation still remains a challenging task. In this paper, a localized hybrid levelset model for vessel image segmentation is proposed. Our technique is inspired by Zhang et al.’s hybrid levelset model [4]. However, instead of using a preset global threshold to specify the term concerning regional information, the proposed method compute the local thresholds automatically to indicate the lower bound of target object, which is crucial for the segmentation of vessel images. Compared with the global threshold based method, the use of locally specified dynamic thresholds enables the extraction of the local image information. Our proposed method has been successfully applied to over ten 3D medical datasets for the segmentation of vasculatures. Some experimental results are presented to demonstrate the strengths of our localized hybrid levelset method. In addition, both the qualitative and the quantitative validations have been performed to indicate that our proposed model can achieve more accurate segmentation results than original hybrid model and CV model do.
References
 1.
Nain D, Yezzi A, Turk G: Vessel segmentation using a shape driven flow. Lect Notes Comput Sci 2004, 3216: 51–59. 10.1007/9783540301356_7
 2.
Lesage D, Angelini ED, Bloch I, FunkaLea G: A review of 3d vessel lumen segmentation techniques: Models, features and extraction schemes. Med Image Anal 2009, 13: 819–845. 10.1016/j.media.2009.07.011
 3.
Kass M, Witkin A, Terzopoulos D: Snakes: active contour models. Int J Comput Vis 1988, 1: 321–331. 10.1007/BF00133570
 4.
Zhang Y, Matuszewski BJ, Shark LK, Moore CJ: Medical image segmentation using new hybrid levelset method. In Proceedings of Fifth International Conference BioMedical Visualization: Information Visualization in Medical and Biomedical Informatics. London, UK: IEEE 2008, 71–76.
 5.
A Review on the Current Segmentation Algorithms for Medical Images http://www.researchgate.net/publication/221116197_A_Review_on_the_Current_Segmentation_Algorithms_for_Medical_Images
 6.
Li C, Xu C, Gui C, Fox MD: Level set evolution without reinitialization: a new variational formulation. In Proceedings of IEEE International Conference on Computer Vision and Pattern Recognition (CVPR). San Diego. USA: IEEE 2005, 430–436.
 7.
McInerney T, Terzopoulos D: Deformable models in medical image analysis: a survey. Med Image Anal 1996,1(2):91–108. 10.1016/S13618415(96)800077
 8.
Klein AK, Lee F, Amini AA: Quantitative coronary angiography with deformable spline models. IEEE Trans Med Imaging 1997, 16: 468–482. 10.1109/42.640737
 9.
Antiga L: Patientspecific modelling of geometry and blood flow in large arteries. PhD thesis. 2002.
 10.
Osher S, Sethian JA: Fronts propagating with curvature dependent speed: algorithms based on hamiltonjacobi formulations. J Comput Phys 1988, 79: 12–49. 10.1016/00219991(88)900022
 11.
Caselles V, Kimmel R, Sapiro G: Geodesic active contours. Int J Comput Vis 1997, 22: 61–79. 10.1023/A:1007979827043
 12.
Chan T, Vese L: Active contours without edges. IEEE Trans Image Processs 2001,10(2):266–277. 10.1109/83.902291
 13.
Li C, Kao C, Gore J, Ding Z: Implicit active contours driven by local binary fitting energy. In Proceedings of 2007 IEEE Conf. Computer Vision and Pattern Recognition. Minneapolis, USA: IEEE 2007, 1–7.
 14.
Li C, Kao CY, Gore JC, Ding Z: Minimization of RegionScalable Fitting Energy for Image Segmentation. IEEE Trans Image Process 2008,17(10):1940–1949.
 15.
Hong Q, Wang B: Segmentation of vessel images using a local hybrid levelset method. In Proceedings of 6th International Congress on Image and Signal Processing. Hangzhou, China: IEEE 2013, 631–635.
 16.
Li B, Chui C, Chang S, Ong S: A new unified level set method for semiautomatic liver tumor segmentation on contrastenhanced ct images. Expert Syst Appl 2012, 39: 9661–9668. 10.1016/j.eswa.2012.02.095
 17.
Lorigo LM, Faugeras O, Grimson WEL, Keriven R, Kikinis R, Nabavi A, Westin CF: Curves: Curve evolution for vessel segmentation. Med Image Anal 2001, 5: 195–206. 10.1016/S13618415(01)000408
 18.
Ambrosio L, Soner H. M: Level set approach to mean curvature flow in arbitrary codimension. J Differential Geometry 1996, 43: 693–737.
 19.
Yan P, Kassim AA: Segmentation of volumetric mra images by using capillary active contour. Med Image Anal 2006,10(3):317–329. 10.1016/j.media.2005.12.002
 20.
Wu X, Luboz V, Krissian K, Cotin S, Dawson S: Segmentation and reconstruction of vascular structures for 3d realtime simulation. Med Image Anal 2010,15(1):22–34.
 21.
Gooya A, Liao H, Matsumiya K, Masamune K, Masutani Y, Dohi T: A variational method for geometric regularization of vascular segmentation in medical images. IEEE Trans Image Process 2008,17(8):1295–1312.
 22.
Vasilevskiy A, Siddiqi K: Flux maximizing geometricflows. IEEE Trans Pattern Anal Mach Intell 2002,24(12):1565–1578. 10.1109/TPAMI.2002.1114849
 23.
Law MWK, Chung ACS: Weighted local variancebased edge detection and its application to vascular segmentation in magnetic resonance angiography. IEEE Trans Med Imaging 2007,26(9):1224–1241.
 24.
Wang Y, Jiang H: A nonparametric shape prior constrained active contour model for segmentation of coronaries in cta images. Comput Math Methods Med 2014, 2014: 1–11.
 25.
Li Q, Griffiths JG, Ward J: Constructive implicit fitting. Comput Aided Geom Des 2006,23(1):17–44. 10.1016/j.cagd.2005.04.011
 26.
Li Q, Tian J: Partial shapepreserving splines. Comput Aided Des 2011, 43: 394–409. 10.1016/j.cad.2011.01.007
 27.
Osher S, Fedkiw R: Level Set Methods and Dynamic Implicit Surfaces. New York: Springer; 2002.
 28.
Orkisz MM, Bresson C, Magnin IE, Champin O, Douek PC: Improved vessel visualization in mr angiography by nonlinear anisotropic filtering. Magn Reson Med 1997, 37: 914–919. 10.1002/mrm.1910370617
 29.
Zijdenbos A, Dawant B, Margolin R, Palmer A: Morphometric analysis of white matter lesions in mr images: Method and validation. IEEE Trans Med Imag 1994,13(4):716–724. 10.1109/42.363096
Acknowledgements
This work is supported by the Fundamental Research Funds for the Central Universities (Grant No. 2013121030), grant of Scientific and Technology Emphasis Project of Fujian Province (2011H0031), projects from National Natural Science Foundation of China (No.61100106, No.61174161, No.61402389), and MOE (Ministry of Education in China) Project of Humanities and Social Sciences (Grant No. 11YJC870027). The authors would like to thank the reviewers for their constructive comments.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
QH designed the study, implemented the algorithm, performed the experiments and was responsible for the data analysis. QL contributed to the discussion and the suggestion through this project and, along with QH drafted the manuscript. BW, YL, JY, KL, and QW contributed to the discussion and suggestion through this project. All authors have read and approved the final manuscript.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
About this article
Received
Accepted
Published
DOI
Keywords
 Segmentation
 Vessel image
 Levelset
 Intensity inhomogeneity