1. Introduction
In recent years, active contour model has caught wide attention in the research of medical image segmentation. [1, 2, 3] Active contours have ability in deformation forced by their inherent geodesic character. For segmentation, edge characters are often introduced as an external force field to drive contours[4, 5]. However, various kinds of noises and illumination bias often distort the force fields and weaken the edge characters. Though corresponding assumption for distortion is still avaliable for specific medical images segmentation, the adaptiveness and automation of active contour models are highly impaired. Accordingly, a modified adaptive acitve contour model is proposed in this paper overcoming various nosies and bias. Especially, two novelly designed energy terms guide the contours rapidly to convergence.
One of the typical active contour models describing the force field is C-V model, which together with other multiple-phase C-V models are called Piecewise Constant (PC) model [2, 3, 4]. Chan and Vese modified the Mumford and Shah functional with piecewise constants and achieved better performance. [2, 7] Ideally, Mumford and Shah assumed piecewise smooth functional for N disjoint and homogeneous distributed regions in a whole image domain. By minimizing the functional, curve will converge to the minimized location. But N partial differential equations must be solved for N disjoint regions, which is computationally expensive. With level set method, Chan and Vese simply separated the curves with only two regional constants indicated by the signed distance function (SDF). The numerical implementation was soon efficacious with C-V model.
However, this property makes PC model fails to segment images in the presence of intensity inhomogeneities, even the image with gradually changing intensity levels at edges. In order to avoid the expensive cost of computation and the contours re-initialization in each step, a distance regularized level set evolution(DRLSE) method are independently introduced by Li et al.[18, 12]. Their model consists of an internal energy term and an external term, respectively. The internal energy term is introduced to rectify the deviation of the level set function from a SDF, whereas the external energy term drives the contours to the desirable boundaries in the image.
Effectively, DRLSE modified the inherent property of the active contour evolution. But it is still confusing when segmenting images with bias field and strong noises (e.x. microscope or MR images). Considering the blurry part of the image as an impact of the bias field, Li et al. proposed a level set method by rectifying the piecewise constant with a bias correction and a pre-defined additive noise term. [13] However, the extent of noise on the images is unknown, of which deviation should not be assumed as a constant. Correspondingly, the predefined additive noise obviously shrinks the model’s segmentation ability into only a small group of images with specific noise. Focusing on the similar problem, Wang et al. modelled the level set with an image guided regularization. By pre-featuring image field with discriminative boundaries for a pre-training, contours are guided towards the true boundaries despite the blurry distribution of the intensity level in local regions. Such a guiding strategy is effective to estimate the bias field for medical images. But the hand labeling process is too empirical and tedious, which also biased the automatic segmentation of active contours. [22].
Utilizing the iterative capability of active contours, this paper designs the noise deviation and bias terms directly into the energy functional for various medical images segmentation. The converged variables give a comprehensive description of the image regions meanwhile segmentation adaptively.
For a regional analysis, a localized mutual information clustering method is introduced to re-describe the piecewise constant model at first. Then, according to the information theory, a probabilistic distribution image model is established to rectify the original variantional framework of the PC model. The level set is finally driven by the total mutual information between the piecewise constant and the target image. As an information measurement, the original energy terms are mapped by a logarithm function. Two problems are generated by the logarithm calculation: (1) the increment of the evolution is not adequate enough after logarithm computation; (2) the distance signed function of the level set is distorted heavily by the logarithm calculation which makes the level set need to be re-initialized in each step. Therefore, corresponding regularization is necessary. For problem (1), an expanding energy functional is designed with edge detector using the image information as prior knownledge. For problem (2), an internal energy functional is designed according to the signed distance function defined in our paper. The variance between the zero level set and original image in each step is introduced as the internal functional rectifying the signed distance to an ideal status. Experiments are designed to illustrate the validation and optimization of our model.
2. Related Works
Active contours are mainly driven by the normal curvature of the curve together with the image features as external forces. A balanced status of the contours will then be found by minimizing either the total forces or the energy. Technically, for a scalar computing, the minimization problem is often defined by an energy functional with relative minimizers. A classical active contour model proposed by Caselles et al. is written as follows[1, 2]:
The first two terms indicate the internal energy of the curve while the last term introduces an external energy generated by the objects; ∇u is the gradient of the image u. Obviously, Eq.(1) is aiming to evolve curve C(s) to the place where the image edges reaches the maximum. Consequently, the converged curves will segment the image at the corresponding edges [1, 2].
With the gradient generated from the image as external energy and curvature central of the curve as its internal energy, it is an available method to evolve the curve and complete the segmentation. However, Eq.(1) lacks detailed force field estimation of the image and hard to segment blurry images with intensive noises since the external forces depend excessively on the gradient character. [2, 5, 13]
For more effective image segmentation, Mumford and Shah proposed the following model with region-based clustering analysis to solve the problems above[4]:
where μ, v are positive constants, Ω represents the domain of the image, and curve C indicates the boundaries in the image. u is the image field, and u0 stands for a piecewise smooth function with marked curve C. Mumford-Shah functional EMS(u, C) is targeting to approximate the force field in the image domain with regional edge features. But if there are N piecesmooth regions needed to be recognized, N partial differential equations (PDEs) must be calculated, which is computationally expensive. Therefore, Chan and Vese simplified the Mumford-Shah functional as follow[2, 5, 12]:
where H is a Heaviside function; ϕ indicates the level set. A level set with a marked contour C has mathematical property: . Obviously, the two disjoint regionsΩ1 and Ω2 are discriminated by curve C naturally. The driving energy includes the first two terms, and the regularization energy, the third term, controls evolution of curve C. The minimization problem is then derived to find three minimizers: level set ϕ and the piecewise constants c1 c2. Hence, the method above is also called piecewise constant (PC) model[12].
PC model considering the regional character without edges instead of edge detectors in the Mumford-Shah functionals by introducing the level set method, which makes it a great distribution on the active contour models. Level set ϕ is often initialized as a signed distance function:
Which marks out the distance and direction from each local region to the curve C. However, the signed distance character will be distorted when minimizing functional ECV(c1, c2, ϕ) Also, it is computational expensive to re-initialize the signed distance function in the evolution.
In order to avoid above disadvantage of CV model, C Li et al. proposed a distance regularization level set evolution (DRLSE). A regularization term is introduced into the energy functional of CV model:
Since level set ϕ is with strict step character, should be a constant. Penalty(ϕ) is a mean square error penalty of introduced into the level set methods. It will keep step character of SDF instead of extra re-initialization. Therefore, the energy functional of DRLSE is written as
Where ELES(R, ϕ) denotes the energy functional of level set functionals, minimizer R denotes the regional indicators.
DRLSE preserves the inherent attributes of level set. Then dealing with outer various noises and illumination bias, Li et al proposed an optimized C-V model overcoming the inhomogeneous intensity levels: let J represent an ideal piecewise constant; b indicates the bias field and n indicates the additive noise. For an observed image I:
where b modulates J to simulate the intensity inhomogeneity generated by the image field and n is an additive noise pre-set before the evolution. The CV energy functional is then wrriten as:
Where Kσ is a pre-defined Gassian constant kernel function with standard deviation σ convlved with the driving term in the local region; the bias field b is added as a variable; ∗ denotes the convolution operator. Obviously, the deviation σ of Gaussian kernel must be adjusted beforehand which highly jeopardizes the automatic and adaptiveness of acitve contours. The specified deviation σ also constrained functional ECV-Li(b, c1, c2, ϕ) into a small group of medical image segmentation.
3. Active Contour Model Based on Localized Mutual Information
3.1 Local Probabilistic Model of Medical Images
In order to establish the local probabilistic model between the piecewise constant and the image, some relative characters of medical imaging equipment should be acclaimed at first. One of the most important property is that all the equipment must keep the original information in the imaging object more or less, that means though it’s blurry, each real boundaries of the original image is included in a interfere distribution. Another property is that when different equipments, for example CT and MR, are imaging various objects, there will exist different level of intensity inhomogenities.
In order to simulate the interfere distribution and bias field, the piecewise constant in a sub-domain of the image Ωi can be rectified and modeled as follows[8, 13, 15]:
Additive term n indicates the interfere distribution in the environment and multiplicative term b indicates the bias field. ci is the original piecewise constant and ci' is the re-constructed piecewise constant. Since the environmental interfere is unknown, the additive interfere is assumed to be a Gaussian distribution with deviation σ. Each region of the piecewise constant exists the interfere distribution as the same form, whereas the deviation σand bias field b will change in various status of the piecewise constant when evolution. As the assumed Gaussian noise n, there exists a probability density in each sub-region Ωi.
Where variance vector αi = {b, ci, σi} and I(y) indicates the corresponding image intensity level at position y in the corresponding sub-region Ωi.
In order to uniform the regional computation into the whole image domain, a binary fitting strategy in paper [13] is adopted. But the calculation of Eq.(10) with such strategy will make the piecewise constant ci and bias field b belong to Ωi affect the value at position . According to the impact generated from the probability density, a limitation ρ of the acceptable overlapped location is defined as a joint probability in each corresponding local region to constrain the overlapping effect. Assume ci' = I(x) here x indicates the corresponding positions with y [15]
Window function p(xi, yi) shows a position limitation for sub-region Ωi meanwhile indicates a joint probability density for each corresponding point between the piecewise constant and the image.
3.2 Localized Mutual Information between Piecewise constant and Image Field
Fig. 1 shows a simplified sample of the localized information between piecewise constant and image field with only one contour C in some sub-region. The green line C is the evolving active curve on the green labeled level set [19].
Fig. 1.region label map
The gray levelΩb and the deeper gray level Ωo with a naturally formed edge is the image field and the red circle marked a neighborhood Ox on curve C between the level set ϕ(x) and the image field I(y) where .
The limitation ρ shows the maximum Euclidean distance of the corresponding point between the image field and the piecewise constant. The limitation will play an essential role in keeping the position deviation since their position are not totally overlapped when there exists a bias field.
Introduce a signed distance function (SDF) for the two regions in the image field I(y) Ω1∩Ox and Ω2∩Ox as L(y) : Ω→ {L1, L2} for representing the region Ω1 Ω2 with a mark matrix {L1, L2} :
By introducing the following mutual information between the rectified piecewise constant and the image, the uncertainty of the piecewise constant simulation is measured.
Since the probability exits and the level set has been marked with Eq.(12), first consider the mutual information in a localized region Ω1∩Ox.
where MI(I(x);L(y)) indicates the mutual information between the image field I(y) and the mark matrix L(x):
An entropy for a continuous random variation on the branch set S is defined by [8][13]
Assume each pixel in the image is independent, thus the conditional entropy is the only measurement needed to be considered. Eq.(14) is simplified to
3.3 Active Contour Model and its Level Set Formulation
Consider a two-phase level set:
C is a contour on the level set. Then the mutual information measurement of Eq.(17) in the whole image domain Ωis represented as follows by integrating εx in the whole region Ω in its level set formulation [14]:
M1(ϕ) = H(ϕ), M2(ϕ) = 1-H(ϕ) and H(ϕ) is the Heaviside function. In numerical implementation, the Heaviside function is replaced by a smooth function approximately called the smoothed Heaviside function with step character, which is defined by
3.4 Energy Minimization
As it is proposed in paper [20] and [21], mutual information between the rectified piecewise constant and the image measures the uncertainty of the approximation from piecewise constant to the image. An optimized approximation is got by minimizing the total uncertainty in the image domainΩ which is equivalently considered as minimizing Eq.(19) of the mutual information measurement. In order to minimize the level set formulation of ε(ϕ, α) in Eq.(19) , a gradient flow algorithm is adopted: [7, 23, 22]
After alternative iteration for the variances αi={b, ci, σi} the image will be segmented and the bias field b will be rectified. For fixed ϕ, the optimal ci, bi, σ that minimize the energy is given by
ui = Mi(ϕ(x)) ; Minimize ε(ϕ, α) according to fixed ci, bi, σ, then the evolution partial differential equation can be calculated out as follow[18]：
Accordingly, the Dirac delta function δ(ϕ), the derivation of the Heaviside function H(ϕ), is computed by [7]
4. Regularizations
4.1 Analysis of Localized Mutual Information Active Contour Model
To introduce local statistic character into mutual information measurement only with considering the neighborhood information in a narrowband distribution and the discrepancy of the variance between different textures is an apparent advantage of the adaptive model.
However, the information measurement may cause some negative factors: Firstly, the increment of the evolution is not adequate enough after logarithm computation; Secondly, the computation in the evolution according to the image information will generate discrepancy between the original signed distance function L(x) and zero level set in each step with a nonlinear logarithmic reflection. Accordingly, an expansion energy functional and internal energy functional are designed to impose on the evolution functional E(ϕ, c, K, σ) keeping the smoothness and signed distance function L(x).
4.2 Expansion Energy
Since expanding based on the image’s edge character is an effective method to accelerate the curve evolution, an expansion energy functional is designed as follow [2, 9, 12]:
where
sgn() is a sign function; G is a Gauss filter with a standard deviation d=0.5. Please note that the filtering process of G is introduced to avoid excessive sharp or flat disparity of the intensity distribution in image I, thus the smooth filter G is designed just with a fixed deviation. is the gradient of the image filtered by the Gaussian filter and ΔG∗I is a result of a Laplace operator impacts on the image filtered by the Gaussian filter, α>0 .[2]
Let gbe the edge indicator function [24] defined by
Which is an edge character functional moving the zero level set curve toward the object boundaries.
According to the calculation that the signal of the second derivation is opposite between two sides of an object edge in an image, there is ΔG×I>0, v(I)>0 for the contour outside the object, which motivates the contour evolve to the outside edge of the object and ΔG×I<0, v(I)<0 for the contour inside the object, which induces the contour evolve to the inside edge of the object. Parameter αensures the extents of the impact on the curve evolution.
4.3 Internal Energy
In our level set evolution, the signed distance function is defined as L(x) in Eq.(12). Hence, level set ϕ should satisfy so that the class of Ω1 and Ω2 in Eq.(18) which is marked by L(x) will not be distorted by the iteration computation.
Thus the integral is proposed as follow to indicate the extent of discrepancy between the level set formulation and the sign distance function L(x)[18]:
Parameter μ ensures the extents of its impact on the curve evolution.
4.4 Level Set Formulation with Regularization Energy and its Minimizing
After regularization, the localized mutual information active contour model is represented as level set formulation below:
The PDE which guides the level set formulation to evolve is calculated by minimizing the energy functional with the method of alternative iteration:
where ei indicates the mutual information of the point between the image field and the piecewise constant. When i=1, the information is inside of the curve, i=2 it is outside. indicates the central curvature of the contour and Δϕ indicates the curvature at each point.
5. Algorithms
Step1 Initialize level set formulation ϕ with signed distance function in Eq.(18):
Step2 Initialize c、b、σ according to ϕ(x), where ci is the regional mean value of ϕ in Ωi; bias field b and deviation field σ equals to 1 as the same size as ϕ.
Step3 Calculate v(I) with Eq.(27) and the edge indicator function g(∇I) with Eq. (28);
Step4 update c、b、σ with Eq.(23) according to fixed ϕ :
Step5 evolve ϕ with Eq.(31) according to the fixed c、b、σ with Eq.(23);
Step6 check whether ϕ is convergent, if not, back to Step4.
6. Experimental Results
6.1 Experimental Setups
Experiments are designed to prove the feasibility of the region-based mutual information approach for medical image segmentation in the presence of intensity inhomogeneities at first. Then some comparative experiment will be given for its superiority. Since the experiment is not targeting to compare the effect of different binary fitting functions in the segmentation, we just set ε=1 in the Dirac function and Heaviside function to keep a standard step character of these functions for binary fitting of the different sub-regions Ω={Ω1, … , Ωn). Iterative incremental term Δt of Eq.(23) is set to be 0.1. All the experiments with different active contour models are evolved to convergent status without an iteration time limitation so that the convergence of all the models is comparable. In the comparative experiments, the relative parameter settings of C Li’s level set is according to paper [6] and the relative parameters in piecewise constant model is set according to paper [2].
For the purpose of objectiveness, the practical medical images are obtained from two different institutions: Shanghai Huashan Hospital and the First Affiliated Hospital of Anhui Medical University. The experiments were conducted using Matlab R2011b programs on a Lenovo ThinkPad notebook with Intel(R) Core i5-4200M CPU @2.50GHz 2.50GHz,64bit, RAM4.00GB (3.76GB usable)with Windows7. For the practicality of our segmentation, the experimental results are all with certification from professional physician at the hospital.
6.2 Qualitative Analysis for Adaptive Segmentation
Firstly, by testing the quality of segmentation a group of comparative experiments in Fig.2 among our model, PC model and C. Li model are listed as follows testing the quality of segmentation. The segmentation is targeting to recognize the human brain structures from the CT images pictured from the top of the head in different depth level. The brain images in this experiment are imaged by standard CT equipments from Shanghai HuaShan Hospital.
Fig. 2.Segmentation of denoised CT brain images with 2-phase level set evolution. Row 1: Segmentation results of PC model; Row 2: Segmentation results of C Li model; Row 3: Segmentation results of our model.
As Fig. 2 shows, C Li model a complete description of the brain in the second row. Muscle, skull, cerebral cortex in different levels are basically accurate for diagnosis with distinguishable levels. In contrary, although results of PC model in the second row are also acceptable, less information of cerebral cortex has been recognized than C Li model. Due to the slowly changing field of the image, the piecewise constant is not comprehensive enough to recognize more useful details but to converge earlier. In the third row, our model shows a figuration of the brain with more details inside of the brain levels and cerebral cortex than the others. As cross sections goes deeper, the skull and cerebral cortex are discriminatively recognized in each level. Since different depth levels of the brain structures has been segmented, a reconstruction of an individual human brain in 3-dimension for future diagnosis will be easy to follow.
Secondly, 8 different medical objects imaged by various imaging equipments are utilized to be segmented with our model demonstrating the adaptability of our model in Fig. 3. Most of the necessary boundaries are segmented by our model. The magnetic resonance white and gray matters in the brain are recognized despite the sharp bias field. Both the cytoplasm and cell walls of the two cells are also discriminative in the microscope image and the two atriums are impressively distinguishable with the converged curves of the heart image. These results are obvious with nearly little error of the segmentation though the bias fields in the magnetic resonance images are apparently strong and some of the images are blurry. The detailed recognition of the human skull with distinguishable esophagus, trachea and nerve center and the marrow structure in the second row adaptively illustrate the superiority of our model.
Fig. 3.Segmentation results and segmented regions of 8 medical objects with various imaging equipments. Row 1 and 2 (From left to right): Denoised Computed Tomography human brain structure; Magnetic resonance white and gray matters in the brain; Magnetic resonance germs; Magnetic resonance human vessel; Row 3 and 4 (From left to right): Microscope cells; Magnetic resonance imaged marrow; denoised Computed Tomography human skull; Computed Tomography heart image.
6.3 Quantitative Analysis
6.3.1 Segmentation Precision
In order to measure the segmentation and recognition with our model in terms of efficiency, accuracy and convergence, Fig. 4 shows a quantitative comparison of ours, C Li’s and C-V level set methods in 2-phase with CPU times, mean errors and minimized energy for 14 objects recognition. These objects are imaged by digital camera, monocular vehicle, CT and MR separately. The time consumed to convergence of the evolution is computed as the CPU times and the minimized energy is computed by Eq.(31) with converged parameters ϕ, c, K, σ. Mean error is the average minimized distance from the segmented contour C to the ground truth contour T. This contour-based metric is introduced to evaluate the accuracy of local pixel for a segmentation result given by a contour. For N points {x1, x2, …, xN} on contour C:
Fig. 4.Comparison of our method, C Li method and C-V method in terms of CPU times, mean errors and minimized energy with 8 physical objects imaged by medical imaging equipments. The red solid line with circle marked is results of our method; the blue solid line with asterisk mark is the results of C. Li’s method; the black dotted line with cross mark is the results of C-V method. From left to right: CPU times(s); Mean errors (pixels) (c) Minimized energy (intensity levels)
As Fig. 4 (b) and (c) shows, our method has the lowest mean errors and minimized energy among the three level set methods which correspondingly demonstrates that our method performs better in terms of accuracy and convergence in the segmentation. However, in Fig. 4(a), the CPU times of our method is higher than C-V method and lower than C Li method. Such results occur since the added variable σ in our method cost too much in the evolution computation than the simpler variantional framework in C-V model, but such cost is not as much as the kernel convolution computation in C Li’s method.
6.3.2 Classification Capability
Considering our method as a 2-class classification of the physical object for the medical images, a group of comparative experiments are conducted to test the classification capapbility. Reletive mearsurements are introduced at first.
For a measurment of the taxonomic effect, Tanimoto distance is adopted as one of our criteria:
NFS represents number of foreground points in the set of converged level set; NFG represents number of foreground points in the ground truth and NFSÇFG represents the number of overlapped points in the ground truth and the converged level set. ATdist closer to 0 or 1 indicates a better discrimination.
In addition, in order to evaluate the classification precision, F-score is introduced:
TP, FP, and FN are the true positives, false positives, and false negatives, respectively. A good segmentation algorithm should provide large F1-scores.
Table 1 compares our method with C Li’s latest MICO method[25] for MRI brain white matter segmentation. The MR brain image data is released by standard OASIS database. As Table 1 illustrated, our method has a better recognition than MICO in brain tissue classification for fast recognition, convergence and classification character with lowest information loss. The CPU times is nearly 5 times faster than MICO, F-scores the higher than MICO and Tanimoto distance is more closer to 0 than MICO, respectively. To sum up, our algorithm performs better in segmentation and speed than MICO.
Table 1.Comparison with MICO in terms of classification capability
7. Conclusions
The adaptive active contour clustered by localized mutual information shows a huge superiority in accurately capturing the object boundaries. Comparative experiment shows that the proposed evolution equations can be applied to automatic segmentation more effectively for medical images. In addition, the internal energy have a good command of regularization and the expansion energy evolves the curve directly to the edge in the image, which improves the algorithmic efficaciously to some extents. However, the PDE formulation in our method is too complicated for numerical implementations. Such a design is hugely cost. In addition, the internal enery hasn’t enough preservation ability. Our future work will focus on the PDE simplification and internal energy refining for a more rapid and effective segmentation.
피인용 문헌
- Smoothing splines of apex predator movement: Functional modeling strategies for exploring animal behavior and social interactions vol.11, pp.24, 2015, https://doi.org/10.1002/ece3.8294