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 [1]. 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) [4] 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) [5] 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 [6] 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 [7] 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 [8] 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 [11] 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 [12] 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 [13] 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. [14] 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 [12] and WFLS method [14]. Finally, some conclusive remarks are included in "Conclusion".

### The related methods

Level set method was developed by Osher and Sethian in 1988 [15], 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 [16] and region-based models [17].

Early level set methods [18–22] mostly belong to edge-based models, which mainly use image gradient to construct an edge detecting function to stop the contour evolution on the object boundary. The popular formulation for level set segmentation is [23]

\frac{\partial \varphi}{\partial t}=g\left|\nabla \varphi \right|\left(\text{div}\left(\frac{\nabla \varphi}{\left|\nabla \varphi \right|}\right)+\nu \right)

(1)

where div(∇*ϕ*/|∇*ϕ*|) approximates mean curvature, *ν* is a balloon force and *φ* is the level set function. The function *g* is image gradient, namely an edge detecting function (*edf*(*I*)), which is defined as

\mathit{edf}\left(I\right)=\frac{\mathit{1}}{\mathit{1}+\left|\nabla {G}_{\sigma}\ast I\right|}

(2)

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.

The regularity of *φ* is very important for stable evolution and accurate computation in level set methods. A common way to reinitialize *φ* is to set |∇*ϕ*| = 1 before the curve deviates from the level set function, so that the curve can evolve stably and accurate segmentation results can be obtained. However, the re-initialization is very complicated and may bring some side effects, e.g., the evolving level set function can deviate remarkably from the signed distance function with a few iterations, especially when the time step chosen is not small enough. In order to overcome the problem, a fast level set formulation was proposed [24]

\frac{\partial \varphi}{\partial t}=\mathit{\mu P}\left(\varphi \right)+\eta \left(g,\varphi \right)

(3)

where *μ*>0 is a parameter controlling the strength with which the deviation of *φ* from a signed distance function is penalized. The first term *P* (*φ*) penalizes the deviation of *φ* from a signed distance function during its evolution and is defined as the following:

P\left(\varphi \right)=\Delta \mathit{\varphi}-\mathrm{div}\left(\frac{\nabla \varphi}{\left|\nabla \varphi \right|}\right)

(4)

The second term *η*(*g*,*φ*) incorporates the image gradient information by

\eta \left(g,\varphi \right)=\mathit{\lambda \delta}\left(\varphi \right)\mathrm{div}\left(g\frac{\nabla \varphi}{\left|\nabla \varphi \right|}\right)+\mathit{\nu g\delta}\left(\varphi \right)

(5)

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. [25] 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.

To overcome the difficulty caused by intensity inhomogeneities, Li et al. proposed the local binary fitting (LBF) model [26, 27], which makes use of the local intensity information. In the LBF model, two spatially varying fitting functions *f*_{1}(*x*) and *f*_{2}(*x*) are introduced to approximate the local intensities on the two sides of the contour, and for a given point *x*∈Ω, the local intensity fitting formulation is:

\frac{\partial \varphi}{\partial t}=\delta \left(\varphi \right)\left(\mu \mathrm{div}\left(\frac{\nabla \varphi}{\left|\nabla \varphi \right|}\right)-{\lambda}_{1}{e}_{1}+{\lambda}_{2}{e}_{2}\right)+\nu \left({\nabla}^{2}\varphi -\mathrm{div}\left(\frac{\nabla \varphi}{\left|\nabla \varphi \right|}\right)\right)

(6)

Where λ_{1} and λ_{2} are positive constant, and *e*_{1} and *e*_{2} are the functions as the following

\left\{\begin{array}{l}\mathit{e}{}_{1}\mathit{=}{\displaystyle {\int}_{\Omega}\mathit{K}{}_{\sigma}\left(\mathit{y}\mathit{-}\mathit{x}\right)\left|\mathit{I}\left(\mathit{x}\right)\mathit{-}{\mathit{f}}_{1}\left(\mathit{y}\right)\right|{}^{2}\mathrm{d}\mathit{y}}\\ \mathit{e}{}_{2}\mathit{=}{\displaystyle {\int}_{\Omega}{\mathit{K}}_{\sigma}\left(\mathit{y}\mathit{-}\mathit{x}\right)\left|\mathit{I}\left(\mathit{x}\right)\mathit{-}{\mathit{f}}_{2}\left(\mathit{y}\right)\right|{}^{2}\mathrm{d}\mathit{y}}\end{array}\right.

(7)

Where *K*_{
σ
}(*y* − *x*)is a Gaussian kernel function, and *f*_{1}(*x*) and *f*_{2}(*x*) are two values that approximate image intensities inside and outside contour C, respectively.

\left\{\begin{array}{l}{f}_{1}\left(x\right)=\frac{{K}_{\sigma}\left(x\right)\left[H\left(\varphi \left(x\right)\right)I\left(x\right)\right]}{{K}_{\sigma}\left(x\right)\left[H\left(\varphi \left(x\right)\right)\right]}\\ {f}_{2}\left(x\right)=\frac{{K}_{\sigma}\left(x\right)\left[\left(1-H\left(\varphi \left(x\right)\right)\right)I\left(x\right)\right]}{{K}_{\sigma}\left(x\right)\left[1-H\left(\varphi \left(x\right)\right)\right]}\end{array}\right.

(8)

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. [28]. In addition, the LBF model is sensitive to initialization to some extent [29], which limits its practical applications. Recently, Liu et al. [30] proposed LRCV model, which have similar capability of handling intensity inhomogeneity as LBF model.