Automatic Image Segmentation Using Active Contours with Univariate Marginal Distribution

This paper presents a novel automatic image segmentation method based on the theory of active contour models and estimation of distribution algorithms. The proposed method uses the univariate marginal distribution model to infer statistical dependencies between the control points on different active contours. These contours have been generated through an alignment process of reference shape priors, in order to increase the exploration and exploitation capabilities regarding different interactive segmentation techniques.This proposedmethod is applied in the segmentation of the hollow core inmicroscopic images of photonic crystal fibers and it is also used to segment the humanheart and ventricular areas fromdatasets of computed tomography andmagnetic resonance images, respectively. Moreover, to evaluate the performance of the medical image segmentations compared to regions outlined by experts, a set of similarity measures has been adopted. The experimental results suggest that the proposed image segmentation method outperforms the traditional active contour model and the interactive Tseng method in terms of segmentation accuracy and stability.


Introduction
Automatic image segmentation is an important and challenging problem in computer vision and medical image analysis.The objective of image segmentation is to separate objects of interest from a given image based on different attributes such as shape, color, intensity, or texture.In recent years, several techniques have been reported for this purpose including graph cut [1,2], improved watershed transform [3], suppressed fuzzy c-means [4], supervised fuzzy clustering [5], multithreshold based on differential evolution [6], and active contour model (ACM), which has been applied in different areas such as intravascular ultrasound images [7], automatic urban buildings [8], and natural images [9], to name a few.
The Active Contour Model is an energy-minimizing spline curve composed of discrete control points called snaxels.The curve is attracted towards features as edges of a target object through the evaluation of internal and external forces.The classical implementation of ACM is prone to be trapped into local minima problem and it is also highly sensitive to initialization of the control points because they require being close to the target object; otherwise failure of convergence will occur.
Since ACM was introduced by [10], many researchers have suggested adapting different techniques to work together with ACM in order to overcome its shortcomings.The suggested improvements of the classical ACM including the introduction of prior knowledge such as active shape models [11], shape prior applied on human cerebellum [12], ACM based on level set method [13], populationbased methods such as genetic algorithms [14], differential evolution [15], and particle swarm optimization [16].The performance of these population-based methods working together with ACM is robust in local minima problem and according to the tests, these methods present a more stable and efficient image segmentation within an appropriate computational time.

Mathematical Problems in Engineering
The population-based methods are an effective way to solve discrete optimization problems.Recently, a new approach known as Estimation of Distribution Algorithms (EDAs) from the family of Evolutionary Algorithms has begun to attract more attention for solving global optimization problems with a fast convergence.EDAs are stochastic methods that incorporate statistical knowledge to solve optimization problems [17].These algorithms consist of a set of potential solutions called population, where each potential solution is referred to as individual, and in general, EDAs work with truncation selection and binary encoding to explore the search space.In the EDAs strategy a subset of individuals is selected, and a probabilistic model of these individuals is constructed.The new individuals will be generated from this model, and the algorithm evolves until a stopping criterion is satisfied.Since EDAs are highly efficient in solving optimization problems, they have been successfully applied in a wide range of applications such as the side chain placement problem [18], dynamic optimization [19], cancer chemotherapy optimization [20], and multiobjective knapsack problem [21].
In this paper, we introduce a novel automatic image segmentation method based on the theory of Active Contour Model and Estimation of Distribution Algorithms.The proposed method uses the Univariate Marginal Distribution Algorithm (UMDA) to infer statistical dependencies between snaxels belonging to a population, in order to increase the exploration and exploitation capabilities regarding the classical Active Contour Model.To establish the initial positioning of the proposed method, a shape prior and the alignment process proposed in [22] to construct a target object template are used.This template is discretized and it is also scaled to different size in order to generate the initial populations of individuals and assuming that the target object is confined within them.This proposed method addresses the problem of segmenting the hollow core in microscopic images of photonic crystal fibers and the human heart and ventricular areas from Computed Tomography and Magnetic Resonance images, respectively.
The remainder of this paper is organized as follows.In Section 2, the fundamentals of Active Contour Model and Univariate Marginal Distribution Algorithm are introduced.In Section 3, the proposed image segmentation method is presented, along with a set of similarity metrics to evaluate its performance.The experimental results are discussed in Section 4, and conclusions are given in Section 5.

Background
In this section, the fundamentals of the Active Contour Model and the Univariate Marginal Distribution Algorithm are described in detail.

Active Contour Models. The traditional Active Contour
Model, also known as snake, is a parametric curve that can move within a spatial image domain where it was defined [10].The curve is defined by (, ) = ((, ), (, )),  ∈ [0, 1], where it evolves through time  to minimize the total energy function given by the following: This energy function is composed of two different energies: the external energy  ext defined by the particular gradient features of the image, and the internal energy  int , which is used to control the shape modification of the curve and to maintain the search within the spatial image domain.On the other hand, the discrete computational implementation of the classical ACM consists of a set of  discrete points {  |  = 1, 2, . . ., }, and the local energy function is given by (2), which is iteratively evaluated in order to minimize, for the actual discrete point, the   index in the   searching window using (3): Since the classical ACM has the drawbacks of initialization and local minima problem, Chan and Vese [13] incorporated a shape prior constraint within the traditional ACM.This method is defined by using the following: represents the total energy function composed of the energies  1 ,  2 , and  3 with their weighting factors  1 ,  2 , and  3 , where  1 represents the active contour or snake.The energy  2 represents the shape energy defined by the difference between the active contour and the shape template and it is expressed as follows: where Ω is the image domain, (⋅) is the Heaviside function,  is a signed distance function,   the deformed template, and   the transformation matrix consisting of translation [  ,   ]  , scaling [], and rotation [] parameters, as follows: where   and   are the translation parameters in the horizontal and vertical axes, respectively.() is the scaling factor, and  represents the rotation angle parameter.Finally, the third energy  3 is the image-based force with an image intensity  and the gradient operator ∇ computed as follows: These three energies are iteratively evaluated, until the difference between the previous and current segmented Mathematical Problems in Engineering 3 area becomes stable.Although the Chan and Vese method is suitable to solve the initialization disadvantage of the classical ACM, this method is prone to be trapped into local minima.An alternative to overcome this drawback is to use population-based methods such as Estimation of Distribution Algorithms, which are described in the following Section 2.2.[23][24][25] are population-based stochastic algorithms that incorporate statistical information to solve optimization problems.Although EDAs are based on the principles of evolutionary computation since they use a population of individuals, binary encoding, and selection operators, EDAs replace the application of the crossover and mutation operators by building probabilistic models of promising solutions based on global statistical information.The probabilistic model that EDAs build in each generation is designed to infer statistical dependencies between the variables in order to generate new solutions.In our work, we will focus on the Univariate Marginal Distribution Algorithm, which is one of EDAs that works perfectly for linear problems and for applications without many significant dependencies [26,27].UMDA works on binary strings and uses a probability vector p = ( 1 ,  2 , . . .,   )  to construct the probabilistic model and then create new solutions for each variable independently, where   represents the probability of obtaining a 1 (binary encoding) in position .The main idea behind UMDA is that it approximates the actual probability distribution of the individuals in P  as the product of the univariate frequencies calculated from the selected population and assuming that all variables are independent [28].In general, UMDA iterates the steps of selection, estimation of probability distribution and the creation of new individuals.In the first step, after the individuals in the search space Ω have been sorted according to fitness, the selection probability  is calculated by using the following proportional selection:

Overview of Estimation of Distribution Algorithms (EDAs). Estimation of Distribution Algorithms
In the second step a joint probability P is computed through the following: where  = ( 1 ,  2 , . . .,   )  is the binary value of the th bit in the binary string (chromosome) and   represents the th component of the random vector .The last step of UMDA generates new individuals from the estimated distribution which will be evaluated according to the fitness function in the next generation.These three steps are iteratively performed until the termination criteria are satisfied.According to the above description, the UMDA algorithm can be implemented through the following.
(2) Select a subpopulation  of  ≤  individuals according to a selection method. (

Proposed Image Segmentation Method
The proposed image segmentation method based on Active Contour Model and the Univariate Marginal Distribution Algorithm is described in Section 3.1.Additionally, to assess the performance of the proposed method, the Jaccard and Dice indices are explained in Section 3.2.

Active Contour Model with Univariate Marginal Distribution.
Due to the two main shortcomings of the traditional ACM discussed above, the Univariate Marginal Distribution Algorithm has been adopted to solve the local minima drawback by building probabilistic models, and the initialization disadvantage is addressed by using scaled templates obtained from an alignment process.Since the methodology of the proposed method makes it possible to apply the UMDA strategy directly in the segmentation task, the advantages of low computational time, efficiency, and robustness are inherently acquired.The procedure of the proposed method is illustrated in Figure 1, and it is described below.This procedure is similar to [15] in the final segmentation step, while the main difference is the introduction of shape templates avoiding the user interaction via seed point and also the use of an Euclidean distance parameter between individuals instead of constrained polar sections.The first step of the proposed method consists of the construction of a shape template through the alignment of a training set of selected reference images, which leads to differences in position, direction, and scale.The alignment process is performed using the technique developed in [22], by estimating the parameters [, , , ]  as follows: where (, ) is the  (horizontal) and  (vertical) translation matrix, () is the scale matrix, and () is the rotation matrix.The product of the three matrices maps the coordinates (, ) ∈ R 2 to coordinates ( x, ỹ) ∈ R 2 , which is used to apply the gradient descent method iteratively in order to minimize the following energy function: where Ω is the image domain and Ĩ is the transformed image.
The final process of this step involves obtaining the final aligned template by superimposing all transformed images and then acquires it through the maximum shape boundary.
Step  In the second step of the method, a preprocessing stage is required, where we first remove the noise from the test image by using a 2D median filter (3 × 3 window size), followed by an edge detection between the background and the target object through applying the Canny edge detector, which is experimentally tuned to  = 1.3,   = 10.0, and  ℎ = 30.0, in order to preserve the real edges in the test image since these can affect the segmentation result.The final step in this preprocessing, is to compute the Euclidean distance map (EDM) according to [29].The EDM is used as potential surface to perform the optimization process, because it assigns high potential values to the image pixels located far from the target object, and low potential values (ideally zero) to pixels located close to the object.The automatic initialization procedure on the resulting EDM is performed by using the maximum mutual information between the template and the current test image.Subsequently, the  initial active contours are created by scaling the final shape template acquired from the previous alignment process.The number of scaled templates has to be considered assuming that the target object is confined within them.After the  contours are defined, these contours must be discretized by a number  of equidistant control points, this parameter has to be considered to adapt and smooth the segmentation result to the shape of the target object.The control points are assigned as individuals to conform a number  of populations , where each population is composed of individuals of different contours with the same position label.The third step of the proposed method involves the numerical optimization and the image segmentation result.The numerical optimization is performed by using 8-bit representation (binary encoding) instead of the intensity of the EDM (real encoding in the range [0, 255]) which is used as fitness function in the optimization process.The UMDA strategy is applied for each population   separately in order to be placed on the nearest edge solution.All the individuals are iteratively evaluated according to the fitness function and the best individual of each population is updated only if a best solution is found considering a maximum distance  max between best individuals.Finally, when the optimization process for each population   is finished, the resulting segmented object is acquired by connecting the best individual of each population to each other.
The procedure of the proposed image segmentation method is described as follows.
(1) Align reference shapes according to [22] and obtain final template after alignment.
(2) Perform maximum mutual information to positioning the template.
(3) Initialize number of active contours  and control points .
(4) Initialize the UMDA parameters: number of generations, number of binary bits, and maximum distance  max .
(5) Generate  populations assigning the control points as individuals.

Evaluation Metrics.
To evaluate the performance of the proposed method on medical images, Jaccard and Dice indices have been adopted to analyze the segmentation results between the regions obtained by computational methods and the regions outlined by experts.The Jaccard index (, ) and Dice index (, ) are similarity measures used for binary variables and defined in the range [0, 1], which are computed using ( 12) and ( 13), respectively.In our tests,  represents the reference segmented object outlined by experts and  represents the automatic segmented object by computational methods [12]: In both indices, when the regions  and  are completely different the obtained result is 0 and is 1 when these two regions are completely superimposed.
In Section 4, the segmentation results obtained from the proposed method on microscopic images of photonic crystal fibers and medical images are presented and analyzed by the evaluation metrics.

Experimental Results
In this section, the proposed method is applied firstly, on the segmentation of the hollow core in microscopic images of photonic crystal fibers, and secondly, to segment the human heart and ventricular areas from computed tomography and magnetic resonance images.The computational implementations presented in this section are performed using the GNU Compiler Collection (C++) version 4.4.5 running on Debian GNU/Linux 6.0, Intel Core i3 with 2.13 GHz and 4 GB of memory.
Figure 2 introduces an image of size 512 × 512 pixels consisting of a microscopic image of a hollow core photonic crystal fiber, where the aim is to separate the hollow core from the image.In Figure 2(a) the original test image is presented and its resulting Euclidean distance map is illustrated in Figure 2(b).On the other hand, in Figure 2(c) the segmentation result by using the classical ACM implementation is presented, in which the noise sensitivity and fitting problem are shown.The ACM parameters were statistically determined as 35 control points,  = 0.01,  = 0.9, and  = 0.05 with an executing time of 0.131 s.In Figure 2(d) the segmentation result obtained with the interactive Tseng method is presented.The parameters of this method were set as 35 control points, window size as 30 × 30 pixels according to [16], and 9 particles for each swarm, obtaining an executing time of 0.160 s.Finally, in Figure 2(e) the segmented image by using the proposed method shows an appropriate hollow core segmentation avoiding the local minima problem and locates the hollow core boundary accurately.In this simulation, the parameters were set as generations = 15,  max = 25, number of contours = 10, and number of control points = 35, obtaining an executing time of 0.177 s. Figure 3 presents a microscopic image of size 512 × 512 pixels consisting of another kind of hollow core photonic crystal fiber.In Figures 3(a   From the dataset of computed tomography images of the human heart described above, in Table 1 the average of the segmentation results obtained from the classical ACM, interactive Tseng method, and the proposed method is compared with the manual delineations by experts.The comparative analysis reveals that the proposed method is promising for the human heart segmentation on CT images. In Figure 5 the segmentation results on a subset of magnetic resonance images containing the human ventricular area are presented.The whole set of MR images consists of 19 images with size 256 × 256 pixels.In Figure 5(a) human ventricular areas outlined by experts are illustrated.According to the human ventricular area dataset of MR images previously described, in Table 2 the average of the segmentation results obtained by the traditional ACM, interactive Tseng method, and our proposed method is compared with the regions outlined by experts.This similarity analysis suggests that the proposed method is more robust in ventricular area segmentation with respect to the comparative computational methods based on the regions outlined by experts.
The use of the Univariate Marginal Distribution Algorithm in the proposed method provides robustness, accuracy and stability in the segmentation problem.Although the computational time of the optimization process is appropriate regarding the comparative computational methods, the proposed image segmentation method improve the segmentation results avoiding the local minima and sensitivity of initialization disadvantages of the classical active contour model.

Conclusions
In this paper, a novel image segmentation method based on the theory of active contour models and estimation of distribution algorithms has been proposed.This segmentation method has introduced two important advantages with respect to different interactive segmentation techniques: firstly, the automatic initialization by using scaled shape templates obtained from an alignment process, in order to overcome the sensitivity to initial contour position and secondly, the incorporation of statistical information of the control points for addressing the local minima problem.Moreover, this proposed method was applied to segment the hollow core in microscopic images of photonic crystal fibers and it was also used to segment the human heart and ventricular areas from CT and MR images.The experimental results demonstrated that the proposed method can lead to more accuracy and efficiency than the traditional active contour model and the interactive Tseng method.In addition, the experimental results have also shown that the exploitation and exploration capabilities of the proposed method, are highly efficient for different applications according to the evidence showed by the set of similarity measures within a competitive computational time.

Figure 1 :
Figure 1: Process of the proposed image segmentation method.

( 6 )
For each population   , (a) apply restriction of the search space to ignore improper individuals; (b) evaluate each individual in fitness function derived from the Euclidean distance map; (c) select a subpopulation of individuals according to selection method; (d) compute the probabilistic model (univariate marginal distribution); (e) generate  new individuals based on the probabilistic model; (f) stop if the convergence criterion is satisfied (e.g., stability or number of generations); otherwise go to step (a).

Figure 2 :
Figure 2: Hollow core photonic crystal fiber: (a) microscopic test image, (b) Euclidean distance map of test image, (c) segmentation result of classical ACM, (d) segmentation result of interactive Tseng method, and (e) segmentation result of proposed method.
) and 3(b) the original test image and its resulting Euclidean distance map are illustrated.In Figure3(c) the segmentation result obtained through the traditional ACM is shown, where due to the noise it cannot adjust to the hollow core boundary accurately.The ACM parameters were experimentally determined as 42 control points,  = 0.013,  = 0.845, and  = 0.195, with an executing time of 0.149 s.In Figure3(d) the segmentation result by using the interactive Tseng method is presented.The parameters of this simulation were set as 42 control points, window size as 30 × 30 pixels, and 9 particles for each swarm, obtaining an executing time of 0.192 s.In addition, Figure3(e)shows the segmentation result obtained from the proposed method, where it can adjust to the hollow core boundary accurately.The parameters of this experiment were set as generations = 15,  max = 20, number of contours = 10, and number of control points = 42, with an executing time of 0.215 s.

Figure 4
Figure 4 presents the segmentation results on a subset of CT images containing the human heart, where the whole dataset consists of 144 CT images of size 512 × 512 pixels obtained from different patients.In Figure 4(a) the human heart outlined by cardiologists is presented.Figure 4(b) illustrates the segmentation results obtained through the classical ACM, where the noise sensitivity and fitting problem are clearly shown.The ACM parameters were determined according to [16] as 45 control points,  = 0.017,  = 0.86, and  = 0.45 with an executing time of 0.161 s.Figure 4(c) illustrates the segmentation results obtained with the interactive Tseng method.The parameters of this simulation were set as 45 control points, window size as 30 × 30 pixels, and 9 particles for each swarm, given an average executing time of 0.205 s per image.Finally, in Figure 4(d) the segmentation results obtained by using the proposed method show an appropriate human heart segmentation avoiding the local minima problem.The parameters of this experiment are set as generations = 15,  max = 15, number of contours = 9, and number of control points = 45, with an executing time of 0.209 s.From the dataset of computed tomography images of the human heart described above, in Table1the average of the segmentation results obtained from the classical ACM, interactive Tseng method, and the proposed method is compared with the manual delineations by experts.The comparative

Figure 4 (
Figure 4 presents the segmentation results on a subset of CT images containing the human heart, where the whole dataset consists of 144 CT images of size 512 × 512 pixels obtained from different patients.In Figure 4(a) the human heart outlined by cardiologists is presented.Figure 4(b) illustrates the segmentation results obtained through the classical ACM, where the noise sensitivity and fitting problem are clearly shown.The ACM parameters were determined according to [16] as 45 control points,  = 0.017,  = 0.86, and  = 0.45 with an executing time of 0.161 s.Figure 4(c) illustrates the segmentation results obtained with the interactive Tseng method.The parameters of this simulation were set as 45 control points, window size as 30 × 30 pixels, and 9 particles for each swarm, given an average executing time of 0.205 s per image.Finally, in Figure 4(d) the segmentation results obtained by using the proposed method show an appropriate human heart segmentation avoiding the local minima problem.The parameters of this experiment are set as generations = 15,  max = 15, number of contours = 9, and number of control points = 45, with an executing time of 0.209 s.From the dataset of computed tomography images of the human heart described above, in Table1the average of the segmentation results obtained from the classical ACM, interactive Tseng method, and the proposed method is compared with the manual delineations by experts.The comparative

Figure 3 :Figure 4 :
Figure 3: Hollow core photonic crystal fiber: (a) microscopic test image, (b) Euclidean distance map of test image, (c) segmentation result of classical ACM, (d) segmentation result of interactive Tseng method, and (e) segmentation result of proposed method.

Figure 5 (Figure 5 :
Figure 5: MR images (ventricular area segmentation): (a) regions outlined by experts, (b) results of classical ACM, (c) results of interactive Tseng method, and (d) segmentation results of proposed method.

Table 1 :
Average similarity measure with the Jaccard and Dice indices among the regions segmented by the traditional ACM, interactive Tseng method, our proposed method, and the regions outlined by experts of the CT dataset.

Table 2 :
Average similarity measure with the Jaccard and Dice indices among the ventricular areas segmented by the traditional ACM, interactive Tseng method, our proposed method, and the regions outlined by experts of the MR dataset.The parameters of this simulation were set as generations = 15,  max = 12, number of contours = 8, and number of control points = 25, with an average executing time of 0.118 s.