An improved level set method for vertebra CT image segmentation
© Huang et al.; licensee BioMed Central Ltd. 2013
Received: 5 January 2013
Accepted: 22 May 2013
Published: 28 May 2013
Clinical diagnosis and therapy for the lumbar disc herniation requires accurate vertebra segmentation. The complex anatomical structure and the degenerative deformations of the vertebrae makes its segmentation challenging.
An improved level set method, namely edge- and region-based level set method (ERBLS), is proposed for vertebra CT images segmentation. By considering the gradient information and local region characteristics of images, the proposed model can efficiently segment images with intensity inhomogeneity and blurry or discontinuous boundaries. To reduce the dependency on manual initialization in many active contour models and for an automatic segmentation, a simple initialization method for the level set function is built, which utilizes the Otsu threshold. In addition, the need of the costly re-initialization procedure is completely eliminated.
Experimental results on both synthetic and real images demonstrated that the proposed ERBLS model is very robust and efficient. Compared with the well-known local binary fitting (LBF) model, our method is much more computationally efficient and much less sensitive to the initial contour. The proposed method has also applied to 56 patient data sets and produced very promising results.
An improved level set method suitable for vertebra CT images segmentation is proposed. It has the flexibility of segmenting the vertebra CT images with blurry or discontinuous edges, internal inhomogeneity and no need of re-initialization.
Lumber disc herniation is an important cause of lower back pains. Clinical diagnosis and therapy for the lumbar disc herniation requires the knowledge of the stress and strain throughout the lumbar region . The finite element method based on medical images is able to analyze the biomedical characteristic of lumbar in the compression. We are sure that accurate 2D vertebra segmentation will help us reconstruct 3D vertebra geometric model because 3D vertebra segmentation modeling is fundamentally performed based on a set of axial slices. The understanding of geometrical information about the normal anatomy and the degenerative bony deformations of the spine necessitates vertebra CT image segmentation for the clinical diagnosis and the preoperative planning of spinal diseases.
There are several proposed approaches in the literature for vertebra segmentation. Statistical shape models (SSMs) [2, 3], which generated mean shapes using their own shape parameters by Fourier and wavelet descriptors, used shape constraints to overcome ambiguous boundary information. Active shape models (ASMs)  was a kind of SSMs that iteratively searched a boundary while maintaining shape constraints. Although SSMs and ASMs could overcome an ambiguous boundary problem, they could not converge into an unusual shape or represent small variations in a boundary. Active appearance models (AAMs)  which combined appearance information and shape constraints, could provide better robust results than ASMs in many medical segmentation applications. However, its application to vertebra segmentation was difficult because the texture patterns of vertebra bodies are different among patients. A deformable spine model  using landmarks exploited shape information and gray-level inhomogeneities using necklace and string models. The necklace model captured variations in vertebra structures while the string model represented spinal curvatures. However the deformable spine model could be trapped into a local minimum and failed to segment abnormal vertebra. Yao J  segmented a vertebra by fitting a four-part vertebra model, but the segmentation could not separate the vertebra region into composing vertebra bones, where a spinous process belonging to the upper vertebra exists with a transverse process pertaining to the current vertebra. Hong S  proposed the concept of localized priors which guided the level set to avoid leakage and local minimum at the places where most necessary, then segmented the completed individual vertebras from the complex neighboring structures. Multiple level set methods [9, 10] were used to extract only vertebra bodies but not to segment spinous parts. Kim Y  presented a fully automatic vertebra segmentation method using 3D deformable fences (3DDF) for 3D CT images, which extracted 2D curve with a deformable model that utilized 3D valley information and was expanded to form a 3D surface. However, it was not robust to segment the vertebral images with weak valley information occurring in abnormal cases. Klinder T  first used various kinds of models, such as shape, gradient, and appearance information, and applied 3D deformable model approach to segment the vertebra CT images. Although they achieved very competitive identification rates for vertebrae, their algorithm depends heavily on spatial registration of the deformable model, which is computationally very expensive. Interactive tools for spine segmentation  were developed to achieve more accurate results. Although the interactive method provides protocols for segmentation, it still required a laborious manual process. Poay et al.  focused on 3D segmentation firstly introducing willmore flow into the level set method (WFLS). The framework incorporated prior shape knowledge through the KDE and local geometrical features by introducing Willmore flow into the level set segmentation and obtained good 3D segmentation results of normal spinal vertebra images.
The shape of the vertebra exhibits complicated topological characteristics. The boundaries in vertebra CT images are ambiguous and discontinuous, while the intensity in vertebra CT images is highly inhomogeneous. The complex shape and inhomogeneous intensity in the vertebra CT images makes its segmentation challenging. In this paper, we developed an improved level set model to achieve a 2D vertebra extraction method. By introducing the edge detection function (edf) and region detection function (rdf) into the proposed model, the images with ambiguous or discontinuous boundary and intensity inhomogenity can be effectively segmented. At the same time, we automatically initialize the level set function by Otsu threshold, thus roughly obtain the regions of interest and multiple initial curves. The curves evolve stably and quickly according to the evolution equation, with its zero level set curves converged to the exact boundary of regions of interest. The algorithm can obtain the accurate segmentation results not only when the internal intensity of vertebra CT images is inhomogeneous, but also when the boundary in CT images is ambiguous or discontinuous. Besides, the algorithm needs no costly re-initialization because of the regularization term, which improves the segmentation speed greatly.
This paper is organized as follows. We briefly review some vertebra segmentation methods and well-known level set methods in "Background". Our edge- and region- based level set (ERBLS) model is presented in "Methods". In "Results", our proposed model is validated by some experiments on synthetic and real images. In "Discussion", we discussed our proposed method and compared our segmentation results with those of 3DDF method  and WFLS method . Finally, some conclusive remarks are included in "Conclusion".
The related methods
Level set method was developed by Osher and Sethian in 1988 , which was an effective method of contour evolution. It utilizes dynamic variational boundaries for image segmentation and can be categorized into two types: edge-based models  and region-based models .
where G σ * I stands for the convolution of the image I with a smoothing Gaussian kernel G σ . The range of edf (I) is between 0 and 1. This edge detector has low values close to 0 at the object boundary, and high value closes to 1 at homogenous background.
where δ(φ) denotes the Dirac function. The parameters μ, λ and ν control the individual contributions of these terms.
In essence, the term η(g,φ) attracts φ towards the variational boundary, which is similar to the standard level set method. The penalty term P(φ) eliminates the computationally expensive re-initialization for signed distance functions. This modification leads to a fast level set algorithm for image segmentation. However, the edge-based level set method only uses the edge detecting function to stop evolving curves, which results in the active contours leaking out the ideal contours when the edges are ambiguous.
To solve the boundary leaking problem, Zhang et al.  proposed a region-based active contour model with a region-based signed pressure force (SPF) function which can efficiently stop the contours at weak or blurred edges. This model only uses the image statistical information of the entire region inside and outside the contour, and is unable to successfully segment images with intensity inhomogeneity.
The LBF model is able to obtain desirable segmentation sometimes in the presence of intensity inhomogeneity. At the same time, the time-consuming re-initialization is avoided. However, the computational cost is still very high, which is pointed out by Zhang et al. . In addition, the LBF model is sensitive to initialization to some extent , which limits its practical applications. Recently, Liu et al.  proposed LRCV model, which have similar capability of handling intensity inhomogeneity as LBF model.
where g σ (x − y)is a Gaussian kernel function with an averaging filter of k × k size and can be considered as the weight assigned to each intensity I(y) at y. Due to the location property of the kernel function g σ (x-y), the contribution of the intensity I(y) to c1(x) and c2(x) decrease and approach to zero as the point y goes away from the center point x. Therefore, c1(x) and c2(x) are determined by the intensities of the points in the neighborhood of the point x. Then the region detection function (rdf) is also dominated by the intensities of the points in the neighborhood of the point x.
where β and γ are fixed constants.
It is obvious that the above equation has the merits of both edge-based models and region-based models. When the contour is far away from object boundaries, the force from the local region intensity information is dominant and has a certain capture range. When the contour is close to the object boundaries, the force from the gradient information becomes dominant, which attracts the contours and finally stop the contours at the object boundaries. The technique of using local region information can improve the robustness to the initialization of contours. When the boundary is blurred or discontinuous, the interference from the local intensity force is able to result in contours’ stopping at the real object boundary. Furthermore, due to the region stopping function making use of local region information, the ERBLS model is able to provide desirable segmentation results even in the presence of the images with intensity inhomogeneity. Besides, our method introduces a new penalizing energy to the regularization term, therefore the computational cost is heavily decreased.
In the proposed ERBLS model, the main computational cost comes from computing c1(x) and c2(x) in Equation (10). At the first sight, there are four convolutions to compute c1(x) and c2(x). It can be noticed that the expression can be rewritten to the combination of the four convolutions: ∫ Ω g σ (x − y)dy, ∫ Ω g σ (x − y)H(ϕ(y))dy, ∫ Ω g σ (x − y)I(y)dy and ∫ Ω g σ (x − y)(I(y)H(ϕ(y)))dy. Because the two convolutions ∫ Ω g σ (x − y)dy and ∫ Ω g σ (x − y)I(y)dy can be computed only once before the iterations, the terms ∫ Ω g σ (x − y)dy and ∫ Ω g σ (x − y)I(y)dy do not depend on the evolution of level set functionφ. Therefore there are only two convolutions ∫ Ω g σ (x − y)H(ϕ(y))dy and ∫ Ω g σ (x − y)(I(y)H(ϕ(y)))dy to be computed at each iteration. In comparison, there are at least four convolutions in the LBF model . As a result, the computational cost of the ERBLS model is about half that of the LBF model for each iteration.
The region detecting function (rdf(I(x)))
Automatic initialization by Otsu threshold
With the above procedures, the initialization of level set function by Otsu threshold can be completely automatic without any human interfaces. The segmentation result can then be taken as the initial contours for the evolution of the ERBLS model.
The clinical image data set was acquired by the Department of Neurosurgery of Beijing Xuanwu Hospital Affiliated to Capital Medical University in China in compliance with the Helsinki Declaration approved by Guojun Zhang. The data set consists of 56 CT images of intervertebral disc protrusion images of patients aged 18 to 66. The patients are carefully selected by radiologists to form a representative group. These images are acquired from 64-detector row Siemens CT System. The in-plane resolution for these images is 1 mm with slice thickness of 1.5 mm. Original images have fixed sizes of 512×512 and the total number of vertebrae is 293. The ground truths are delineated by clinical experts. Our algorithm is implemented in Matlab 7.0 on 2.79-GHz Intel Pentium IV PC. Unless otherwise specified, we used the following parameters in our model: σ=3.0, ϵ=1.0, β=5.0, γ=2.0, time increment Δt=1.0.
Iteration number and processing time for the LBF model and proposed ERBLS model in segmenting the images in Figure 3
CPU time (s)
CPU time (s)
Iteration number and processing time for the LBF model and the proposed ERBLS model in segmenting the images in Figure 4
CPU time (s)
CPU time (s)
In order to further evaluate our segmentation algorithm, we reconstructed 3D vertebra images based on 2D segmentation results by using our proposed method. The proposed method has been applied to 56 patient data sets and the segmentation results are compared with those of 3D deformable fences method (3DDF)  and introducing willmore flow into level set segmentation (WFLS) .
where and represent the volumes of segmented object and the ground-truth respectively. The measurement (varies from 0 to 1) indicates the correspondence between two volumes, i.e., 0 indicates the two volumes do not overlap and 1 shows they are perfectly matched.
where A and B are the boundaries of two different segmented volumes, respectively. It measures the distance between the farthest point of a set to the nearest point of the other. The measurement (varies from 0 to ∞ theoretically) represents the difference between two closed surfaces, e.g., 0 shows that both volumes share exactly the same boundaries, and larger HD values indicates larger distances between the boundaries. In summary, a high DSC and a low HD are desirable for good segmentation.
Average DSC and HD (mm) with standard deviation for segmentation of vertebra CT images using our ERBLS method, 3DDF method and WFLS method
We have described an edge- and region-based level set method for accurate segmentation of vertebra CT images. The ERBLS model can efficiently segment the images with intensity inhomogenity and blurry or discontinuous boundaries by employing the image gradient information and the local image information. Meanwhile, the level set function is automatically initialized by Otsu threshold, which segmentation result is taken as the initial contours of the EBRLS model. Experimental results on both synthetic and real images demonstrated that the proposed ERBLS model is very robust and efficient. Compared with the well-known local binary fitting (LBF) model, the ERBLS model is not only much more computationally efficient and but also much less sensitive to the initial contours. The proposed method has also applied to 56 patient data sets and has produced very promising results.
About the Author—Juying Huang, PhD, is working at the medical physics education and research work. Her current research interests are medical image processing, CT image reconstruction, medical system modelling and computing.
About the Author—Fengzeng Jian, PhD, Chief physician, is working at clinical and research work in the spine and skull base surgery.
About the Author—Hao Wu, PhD, associate Chief physician, is mainly working at surgical operation in the spinal cord disease.
About the Author—Haiyun Li received his PhD degree in the Dept. of Biomedical engineering, Sun Yat-Sen University in 1997. In 2000/6-2002/10, he was Postdoctoral Fellow in vision and image processing Lab, National University of Singapore, Singapore. In 2008/10-2009/10, he was the visiting scholar, Henry H. Wheeler Jr. Brain Imaging Center, Helen Wills Neuroscience Institute, University of California, Berkeley. USA. Currently, he is a Professor in School of Biomedical Engineering, Capital Medical University, Beijing, China.His current research interests include Computer simulation and Medical Image Computing, MRI/fMRI, Medical system modelling and computing.
The subject is jointly sponsored by the National Natural Science Foundation of China (Grant No.30670576), Beijing Natural Science Foundation (Grant No. 4122018) and National Natural Science Foundation of China (Grant No. 81271519).
- Li H, Wang Z: Intervertebral disc biomechanical analysis using the finite element modeling based on medical images. Comput Med Imaging Graph 2006, 30: 363–370. 10.1016/j.compmedimag.2006.09.004View ArticleGoogle Scholar
- Robin S, Skalli W, Lavaste F: Influence of geometrical factors on the behavior of lumbar spine segments: a finite element analysis. Eur Spine J 1994, 3(2):84–90. 10.1007/BF02221445View ArticleGoogle Scholar
- Benameur S, Mignotte M, Parent S, Labelle H, Skalli W, DeGuise J: 3D/2D registration and segmentation of scoliotic vertebrae using statistical models. Comput Med Imag Graph 2003, 27: 321–327.View ArticleGoogle Scholar
- Aubin CE, Dansereau J, Petit Y: Three-dimensional measurement of wedged scoliotic vertebrae and intervertebral disks. Eur Spine J 1998, 7: 59–65. 10.1007/s005860050029View ArticleGoogle Scholar
- Roberts MG, Cootes TF, Adams JE: Automatic segmentation of lumbar vertebrae on digital radiographs using linked active appearance models. Proc Med Image Underst Anal 2006, 2: 120–124.Google Scholar
- Ghebreab S, Smeulders AWM: Combining strings and necklaces for interactive three-dimensional segmentation of spinal images using an integral deformable spine model. IEEE Trans Biomed Eng 2004, 51: 1821–1829. 10.1109/TBME.2004.831540View ArticleGoogle Scholar
- Yao J, O’Connor S, Summers RM: Automated spinal column extraction and partitioning. IEEE 2006 international symposium on biomedical imaging (isbi 2006) proceedings 2006, 9054328: 390–393.Google Scholar
- Hong S, Andrew L: In Localized Priors for the Precise Segmentation of Individual Vertebras from CT Volume Data. Edited by: Metaxas D. 2008, 367–375. MICCAIGoogle Scholar
- Tan S, Yao J, Ward MM, et al.: Level set based vertebra segmentation algorithm for the evaluation of Ankylosing Spondylitis. Proc SPIE Med Imaging 2006, 6144: 58–67.Google Scholar
- Tan S, Yao J, Ward MM, et al.: 3D multi-scale level set segmentation of vertebra. Biomedical imaging: from Nano to Macro, ISBI 2007. 4th IEEE Int Symp 2007, 9528684: 896–899.Google Scholar
- Kim Y, Kim D: A fully automatic vertebra segmentation method using 3D deformable fences. Comput Med Imaging Graph 2009, 33: 343–352. 10.1016/j.compmedimag.2009.02.006View ArticleGoogle Scholar
- Klinder T, Ostermann J, Ehm M, et al.: Automated model-based vertebra detection, identification, and segmentation in CT images. Med Image Anal 2009, 13: 471–482. 10.1016/j.media.2009.02.004View ArticleGoogle Scholar
- Kaminsky J, Klinge P, et al.: Specially adapted interactive tools for an improved 3D-segmentation of the spine. Comput Med Imag Graph 2004, 28(3):119–127. 10.1016/j.compmedimag.2003.12.001View ArticleGoogle Scholar
- Poay HL, Ulas B, Li B: Introducing Willmore flow into level set segmentation of spinal vertebrae. IEEE Trans on Biomed Eng 2013, 60(1):115–122.View ArticleGoogle Scholar
- Osher S, Sethian JA: Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations. Comput Phys J 1988, 79(1):79.MathSciNetView ArticleGoogle Scholar
- Caselles V, Kimmel R, Sapiro G: Geodesic active contours. Comput Vision J 1997, 22: 61–79. 10.1023/A:1007979827043View ArticleGoogle Scholar
- Chan T, Vese L: Active contours without edges. IEEE Trans Image Process 2001, 10: 266–277. 10.1109/83.902291View ArticleGoogle Scholar
- Morris H, Kallinderis Y: Octree-advancing front method for generation of unstructured surface and volume meshes. AIAA Jounal J 1997, 35(6):976–984. 10.2514/2.206View ArticleGoogle Scholar
- Saxena M, Perucchio R: Parallel FEM algorithms based on recursive spatial decomposition automatic mesh generation. Computers Structures J 1992, 45(5):817–831.View ArticleGoogle Scholar
- Löhner R: A parallel advancing front grid generation scheme. Int J Numer Meth Engng 2001, 51: 663–678. 10.1002/nme.175.absView ArticleGoogle Scholar
- Xu YA, Yang Q, Wu ZZ, et al.: The algorithm of 3D constrained Delaunay triangulation. Chinese J Software J 2001, 12(1):103–110.Google Scholar
- Canvendish JC, Field DA, Frey WH: An approach to automatic three-dimensional finite element mesh generation [J]. Int J Num Methods Engin 1985, 21(2):329–347. 10.1002/nme.1620210210View ArticleGoogle Scholar
- Caselles V, Kimmel R, Sapiro G: Geodesic active contours. Int J Comput Vision 1997, 22: 61–79. 10.1023/A:1007979827043View ArticleGoogle Scholar
- Li C, Xu C, Gui C, et al.: Level set evolution without re-initialization: a new variational formulation. Proceedings of the 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05) 2005, 1: 430–436.View ArticleGoogle Scholar
- Zhang K, Zhang L, Song HH: Active contours with selective local or global segmentation: a new formulation and level set method. Image Vision Comput 2010, 28: 668–676. 10.1016/j.imavis.2009.10.009View ArticleGoogle Scholar
- Li C, Kao CY, Gore JC, Ding Z: Implicit active contours driven by local binary fitting energy. Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR) 2007, 9737883: 1–7.Google Scholar
- Li C, Kao CY, Gore JC, Ding Z: Minimization of region-scalable fitting energy for image segmentation. IEEE Transact Image Process 2008, 17(10):1940–1949.MathSciNetView ArticleGoogle Scholar
- Zhang K, Song H, Zhang L: Active contours driven by local image fitting energy. Pattern Recognition 2010, 43(4):1199–1206. 10.1016/j.patcog.2009.10.010View ArticleGoogle Scholar
- Wang L, Li C, Sun Q, Xia D, Kao CY: Active contours driven by local and global intensity fitting energy with application to brain MR image segmentation. J Comput Med Imaging Graph 2009, 33(7):520–531. 10.1016/j.compmedimag.2009.04.010View ArticleGoogle Scholar
- Liu S, Peng Y: A local region-based Chan–Vese model for image segmentation. Pattern Recognit 2012, 45(1):2769–2779.View ArticleGoogle Scholar
- Asari KV: A fast and accurate segmentation technique for the extraction of gastrointestinal lumen from endoscopic images. Med Engin Phys 2000, 22: 89–96. 10.1016/S1350-4533(00)00015-1View ArticleGoogle Scholar
- Kim DY, Chung SM, Park JW: Automatic navigation path generation based on two-phase adaptive region-growing algorithm for virtual angioscopy. Med Engin Phys 2006, 28: 339–347. 10.1016/j.medengphy.2005.07.011View ArticleGoogle Scholar
- Cristoforetti A, Luca F, Flavia R, et al.: Isolation of the left atrial surface from cardiac multi-detector CT images based on marker controlled watershed segmentation. Med Engin Phys 2008, 30: 48–58. 10.1016/j.medengphy.2007.01.003View ArticleGoogle Scholar
- Otsu N: A threshold selection method from gray-level histograms. IEEE Trans Sys, Man, Cyber 1979, 9(1):62–66.MathSciNetView ArticleGoogle Scholar
- Li BN, Chui CK, Chang S, Ong SH: Integrating spatial fuzzy clustering with level set methods for automated medical image segmentation. Comput Biol Med 2011, 41: 1–10. 10.1016/j.compbiomed.2010.10.007View ArticleGoogle Scholar
- Huh S, Ketter TA, Sohn KH, Lee C: Automated cerebrum segmentation from three-dimensional sagittal brain MR images. Comput Biol Med 2002, 32: 311–328. 10.1016/S0010-4825(02)00023-9View ArticleGoogle Scholar
- Dice L: Measures of the amount of ecologic association between species. Ecology 1945, 26: 297–302. 10.2307/1932409View ArticleGoogle Scholar
- Rockafellar RT, Wets RJB: Variational Analysis. NewYork: Springer-Verlag; 2005.Google Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.