The aim of image segmentation is to obtain meaningful partitions of an input image into a finite number of disjoint homogeneous objects. Active contour models are popular in the regard. Chan and Vese  proposed an active contour without edges scheme based on the classical work of Mumford and Shah  variational energy minimization model. Since biomedical images typically have multiple regions of interest with different characteristics, deriving a multiphase active contour scheme for efficient segmentation is an important area of research in image processing [19, 48, 31].
In MRI (magnetic resonance image) images, segmentations based on active contours have been used with traditional level set method . Active contours can also be improved using region information , salient features  or mathematical morphology  etc. Traditionally these schemes use a gradient descent formulation to implement the non-convex energy minimization and can stuck in undesired local minima thereby lead to erroneous segmentations. Moreover, traditional level set based implementation is prone to slower convergence due to the well-known re-initialization requirement and discretization errors. More recently quite a lot of interest is being shown in techniques that can obtain a general convex formulation for active contours schemes based on energy minimization which can alleviate the problem of local minima at the same time focussing on the computational complexity [23, 13, 8, 14, 39, 30]. Among other techniques for MRI image segmentation, we mention fuzzy C-means based models [1, 35, 21], fuzzy connectedness , automatic labeling 
, adaptive expectation-maximization (EM), Bayesian EM 
, hidden Markov model EM, kernel clustering , optimum-path clustering , anisotropic diffusion combined with classical snakes model , discriminant analysis 
, and neural networks. We also refer to [42, 43, 36] for reviews about segmentation for medical images in general and  for MR images in particular. The area of MRI image segmentation has seen tremendous research activity and a more detailed review in this particular field can be found in .
In this paper, we consider a globally convex version of the four phase piecewise constant energy functional following the seminal work of Chan et al . By deriving an approximate novel convex functional we change the original formulation into a binary segmentation problem and utilize a dual minimization to solve the relaxed formulation . The proposed global methodology avoids the level set re-initialization constraint and other ad-hoc techniques  used for fixing level set active contour movements throughout the iterations. The proposed approach is used to obtain white matter and gray matter partitions on brain MRI images as can be seen for example in Figure 1. Our scheme does not involve level sets or re-initialization and instead relies on the relaxed globally convex formulation of the Vese and Chan multiphase active contours. Comparison results on the different image sets with varying noise and inhomogeneities show that we can obtain better results than traditional level set multiphase schemes [48, 5, 7, 6, 10] and primal-dual approach of . Moreover, compared to these traditional level set based implementations we achieve faster convergence due to the usage of efficient alternating dual minimization. The proposed approach is general in the sense that we can add domain specific knowledge to improve such active contour schemes further for various tasks [33, 50, 26, 25, 44, 45].
The main contribution of our work is two-fold: 1) a fast four phase active contour model using a relaxed globally convex minimization approximation; 2) using an efficient dual minimization based implementation for performing segmentation on MRI images. The rest of the paper is organized as follows. Section 2 introduces the multiphase variational active contour scheme and provides a globally convex formulation. Section 3 illustrates the segmentation results on various Brain MRI images including comparison of different schemes. Finally, Section 4 concludes the paper.
2 Multiphase active contours model
We first recall the multiphase formulation of Vese and Chan  and restrict ourselves to the piecewise constant four phase model since the general case can be derived similarly. Let be the two level sets. , and , , where is the Heaviside function, representing four regions. Our goal is to solve a minimization problem
where , and the constant mean values can be derived as
Note the the zero level sets , , represent object boundaries and the mean values represent the expected average pixel values in these objects. Vese and Chan  used the corresponding gradient descent equations to implement the active contours . In the numerical implementation of the above PDEs, a non-compactly supported, smooth approximation of the Heaviside function , such that as is utilized. Since the above minimization (1) is non-convex the time discretized gradient descent PDEs usually require large iterations and small time steps to convergence (typically in ’s of iterations). Moreover, the final segmentation result may not correspond to the global minimum of the energy function as the gradient descent scheme can be stuck at a local minima of the corresponding energy functional given in Eqn. (1).
We briefly recall the corresponding gradient descent equations (time dependent Euler-Lagrange equations of Eqn. (1)) for the level sets functions and ,
respectively. Here, the image fitting terms are given by,
Then correspondingly we can derive an energy functional which does not depend on regularized Heaviside functions. Thus, we can solve the following globally convex energy minimization problem,
where Heaviside functions are replaced by which are known as binary partitioning functions. The above modified minimization problem (5) can further be relaxed to the set of functions in order to solve a convex minimization problem. That is, the binary partitioning functions based energy minimization becomes,
The following theorem provides the guarantee of finding a global minimizer for the derived functional (5) in terms of the relaxed version in (6). We follow arguments similar to the work of Chan et al  and  to prove the following result.
We use the standard notation for functions of bounded variation . Since , it follows from the standard total variation based Coarea Formula,
For the image fitting term,
Further similar computations yield,
Defining , it follows that
The final segmentation is obtained by thresholding the functions and with any number in the interval for example at , as shown in Figure 1
(b) and (c). Note that the above modified minimization model does not involve level sets and thus can be solved efficiently. Further, we can prove that the above relaxed minimization problem can be solved in a binary variable minimization formulation to find a global minimum. The existence of minimizers of the modified energygiven in Eqn. (6) is proved using the theory of functions of bounded variation (BV) space .
For a given input gray scale image , there exists a minimizer for the functional in (6) in
Let and be a minimizer sequence for the energy , i.e.,
Since is bounded in , there is a subsequence also denoted by , strongly convergent to an element . Furthermore, . Therefore, it follows that and
Now, considering as a function of , its minimization brings the following two equations,
Since , it follows is uniformly bounded. Hence, there is a subsequence also denoted by
and a constant vectorsuch that
Then, from Fatou’s lemma we get for the suitable sequence :
i.e., is a minimizer of the functional . ∎
Note that the values given in the above theorem are computed in the numerical scheme based on a dual minimization formulation which we describe next.
2.1 Implementation details
The four phase convex minimization problem in (6) is solved in an alternating fashion for the image variables :
First fix , and solve for :
Then fix , and solve for :
where the image region fitting terms are given by,
To solve the above convex optimization problems we use the Chambolle’s dual formulation [16, 12] of the total variation regularization function which occurs as the first term in the energy functional in Eqn. (6). Thus, the new unconstrained minimization problems to consider are (for ):
where and , is chosen to be small and and . We solve the above by further splitting into two sub-problems:
Solve for :
The solution is given by
The vector satisfy the equation
and it is solve by a fixed point method: and
Solve for the auxiliary variable :
for which the solution is given by:
Furthermore, at every few iterations the vector is updated according to the following equations:
The computation of values are similar to the ones in Vese and Chan model  (see Eqn. 1) except that they are now based on the binary partitioning functions and does not involve computing regularized Heaviside functions. We refer to  for more details on this particular form of dual minimization and the motivation for the fixed point method used to derive the solution for the auxiliary variable in the second step.
3 Experimental results
. (e) Color visualization of the segmentation result. Second row: (g) Cumulative distribution function (CDF) of the four regions from (e). (g) Histogram of the four regions showing the separation clearly.
We have used full brain MRI images available at the whole brain atlas 555http://www.med.harvard.edu/aanlib/home.html. The parameters were fixed for the segmentation results reported here. In order to simplify notations we use and we fix in all our experiments as well. Equal weights are used for the four regions to be segmented as we do not want to introduce bias for certain phases. The images presented here are from MRI imaging modality with slice thickness of . Our scheme takes less than seconds (for iterations) on MATLAB2012a on a Mac laptop with Intel Core i7 CPU GHz, GB RAM CPU. Meanwhile, the average computation time for related models compared from the literature are in the region of seconds (for iterations) to converge to the final segmentation.
Figure 2 shows another example segmentation result of our globally convex four phase scheme. The noise (calculated relative to the brightest tissue, and denoted by ””) is set to with intensity non-uniformity (denoted by ””) is of strength . The result of our four phase model is displayed in Figure 2(b) with two contours (Purple, Light-Blue) overlaid on top of the input image. Figure 2(c) and (d) show the two functions computed using our scheme and thresholded at . The function captures the background shape 2(a) (corresponding to level set ) whereas function in Figure 2(b) (corresponding to level set ) contains the white matter. Figure 2(e) we use four different colors (Blue, Green, Yellow, and Maroon) to highlight different phases for better visualization of phase separation and boundary detection of regions. In Figure 2(f) and (g), we show cumulative distribution function (CDF) and the histogram of each of the four regions computed by the proposed method. The histograms highlight separation of different phases/regions indicating the superior performance of our splitting based numerical approach.
Figure 3 shows a comparison result with other multiphase active contour methods from [5, 7, 10, 17] called in short, Ker, Mean, Clust and Primal-Dual respectively, for the same image in Figure 2(a). Note that to make a fair comparison with other models we used the same noise level and intensity non-uniformity for this example image. In Figure 3 bottom two rows we show the cumulative distribution function (CDF) and histograms computed for each of the computed phases respectively. Compared with the histograms shown in Figure 2(f) and (g) for our scheme we see that proposed model provides better separation of regions. The histograms for the other schemes in Figure 3(last row) show nontrivial intersections, highlighting the drawback in using level set based implementations. Moreover, the noise remains as speckles in the segmented regions whereas our model handles it efficiently.
3.1 Error metrics computation
We use the following quantitative error metrics to compare the schemes with gold standard ground truth segmentations. For more details about objective evaluation of image segmentation algorithms and for precise definitions of these metrics we refer to .
The Dice coefficient  is a popular error metric and is used to compare ground truth segmentation with those obtained with automatic multiphase segmentation schemes. By definition, for two binary segmentations and , the Dice coefficient is computed as:
Here the binary segmentation is computed automatically, using the segmentation curves and by thresholding regions obtained by all algorithms. The notation denotes the number of pixels in the set . Note that, a value of indicates perfect agreement. In particular, higher numbers indicate that the results of that particular scheme’s result match the gold standard better than results that produce lower Dice coefficients.
Rand Index: A metric based on a classical nonparametric test and is computed by counting pairs of pixels that have compatible label relationships in the two segmentations to be compared.
Global Consistency Error: A metric which computs the degree of overlap of the cluster associated with each pixel in one segmentation and its closest approximation in the other segmentation. Values to closer to indicate better segmentation results.
Variation of Information: A metric related to the conditional entropies between the class label distribution of the segmentations. This computes a measure of information content in each of the segmentations and how much information one segmentation gives about the other. Values closer to indicate better segmentation results.
Note that all these metrics are for comparing two segmentations, one of which is assumed to be the available ground truth. Table 1 shows the comparison of average Dice values (for images) of different models for different noise and intensity inhomogenieties taken from Brainweb database. As can be seen our scheme performs better in terms of the Dice coefficient compared with other related approaches. Similarly in Table 2 we see that the average RI, GCE and VI for different schemes against our model shows that the proposed globally convex multiphase scheme performs well overall.
Figure 4 shows representative segmentation results for full Brain data-sets (axial slices are shown) with different noise (“”) and non-uniformity (“”) levels for our scheme. Different and are specified in Figure 4 for each row. This illustrates that our scheme preserves the topological changes as we move through the image stack. Moreover, our scheme can handle noise and intensity non-uniformity together effectively. Finally, in Figure 5 we show different segmentation results for a particular image (slice number 79) taken across all noise and inhomogeniety levels for different schemes. The results indicate that Ker and Mean methods can lead to poor separation of different regions whereas noise can affect the result of Clust and Primal-Dual schemes. Meanwhile, our approach performs well and handles higher non-uniformity without degrading the final segmentation results. Further data-sets and extensive comparison results of all the schemes for full brain stacks are available online 666http://dx.doi.org/10.6084/m9.figshare.781297.