In image segmentation, the Mumford-Shah model proposed by Mumford D and Shah J  is regarded as the most significant region-based model and has been applied to many applications. In 2001, the two-phase Chan-Vese (CV) model  was proposed to detect objects in a given image. With the increasing complexity of images, the multiphase segmentation models [3, 4, 5] were proposed and these models mainly represent different regions by using the level set functions [6, 13, 15, 17]. In order to reduce the number of level set functions, Chan et al. proposed a multiphase segmentation model , which is a generalization of CV model.
were also established subsequently according to different noise distribution. They obtained the characteristic information contained in images by estimating corresponding parameters. When dealing noisy images, it is known that many segmentation problems need a suitable noise model, e.g., synthetic aperture radar, positron emission tomography, electron micrograph or medical ultrasound imaging, etc. Especially when the data were collected with poor statistics, it is necessary to consider the influence of the noise probability distribution in segmentation implementation.
Recently, authors in [11, 12, 13] made some progress in achieving illusory contour recovery while doing segmentation, which can identify absent boundaries or missing shapes successfully without necessary region features. In detail,  employed the fitting terms of two-phase CV model and the Euler s elastica term 
as the regularization. Its major contribution was that the missing boundaries were interpolated automatically without specifying the regions. improved the segmentation with depth problem [14, 15] and achieved acceleration via the strategies of model simplification and constraint projection. Its significant performance enhancements included shape reconstruction of occluded objects and determination of their ordering relation in a specific scene based on only one single image. Many other works illustrated that the curvature-related terms have played crucial roles in the boundary reconstruction [16, 17] and image restoration [18, 19] with the capacity of producing excellent edge and corner preservation results. All of these researches show the significant potential for curvature-based methods.
However, the current segmentation models mentioned above cannot be directly applied for noisy images when the type of noise is unknown or there are more than one type of noise in the image. The reason is that in these models there exists a one-to-one mapping relationship between the parameters to be evaluated and the noisy images with some probability density distribution. Besides, the curvature-related terms will bring extra computational complexity due to the existence of nonlinear higher-order derivatives. This issue also appears in other variational models such as the non-texture image inpainting and image denoising 
with features (edge, corner, smoothness, contrast, etc.) preservation. Hence it is essential to take some mathematical optimization techniques, for example with global solution, stability guarantee and calculation acceleration, into consideration in the process of algorithm design.
For the difficulty of dealing with stochastic noises which are inherently generated by the acquisition procedure of imaging due to various issues, [22, 23] have given us a great source of inspiration. The authors extended the range of applications for progressive hedging algorithm (PHA) in multistage stochastic variational inequality problems and explored stochastic complementarity problems in a two-stage formulation. Convergence of the algorithm was proved in detail and how the PHA performs was validated through numerical experiments as well as its stability and practicability. One of our major motivations is to embed stochastic property into segmentation energy functional in images with unknown noises or arbitrary damages and then implement PHA to solve it. Furthermore, Euler s elastica term will be also considered as the regularization in our variational formulations design since its better properties in dealing with image feature information. Last but not least, in order to improve the computational efficiency and solve problems caused by the non-convex, non-differentiable, nonlinear and higher order terms involved in Euler s elastica related functional, fast algorithms of alternating direction method of multipliers (ADMM) [13, 17, 19, 24] and curvature weighted approach [20, 25, 26] will be systematically designed in the PHA algorithm framework as a fusion for energy minimization problems. The creative introduction of Euler s elastica term along with the nontrivial task of its analytical study will be another main goal of our research.
Our contributions in this paper can be summarized in the following aspects:
We intend to propose novel formulations for image segmentation with stochastic noises for various applications by transforming the original minimization problems into the optimization framework of stochastic programming.
Euler s elastica term is employed to realize completion of meaningful missing boundaries and reconstruction of occluded structures of objects, which further enhances the segmentation performance.
Our novel variational formulations will be applied in the problems of two-phase Euler s elastica based segmentation and segmentation with depth in gray and color spaces respectively.
A general numerical algorithm based on PHA is proposed to calculate the novel formulations. Fusion of ADMM and curvature weighted approach (ADMM-C) is designed for the minimization of the Euler s elastica energy related sub variational problems. The minimization problems derived from ADMM-C will be then efficiently solved by Fast Fourier transform (FFT) [26, 34] and analytical soft threshold formulas [17, 19].
The rest of this paper is structured as follows. Section 2 briefly reviews the related approaches in this field. Our proposed approach and algorithm framework are presented in Section 3. The experiments conducted with performance evaluation and comparison are described in Section 4 followed by the conclusion in Section 5.
2 Research background
For the purpose of clarifying the motivations in this paper, some related works will be briefly reviewed before presenting our contributions clearly in Sections 3.
2.1 Euler s elastica based segmentation
Illusory contour capture and shape reconstruction [27, 28] is a challenging problem which aims to complete the missing boundaries or fuzzy areas for existing objects in an image. It is a very common phenomenon in human vision. Part of illusory contours consists of objects actual boundaries and the other part is made of missing perceptual edges. However, current computer technique can only deal with closed boundaries. It is extremely difficult for computers to identify illusory contours automatically. Euler s elastica based segmentation techniques have one place in this problem due to their crucial roles in boundary reconstruction and image restoration.
Two-phase Euler s elastica based segmentation: Zhu, Tai, and Chan  proposed the Chan-Vese-Euler (CVE) model designed for the foreground shape recovery problem by combining the CV model  and Euler elastica regularizer . This model could recover the illusory contours and form a complete meaningful object, even without requiring initialization of fixed points. According to their work, the energy functional is defined as
where are positive penalty parameters, is a binary level set representation supposed to take on either 0 or 1. The last term of this functional is the classic Euler s elastica term. denotes the curvature represented as .
Segmentation with depth information: In , Nitzberg, Mumford and Shiota defined the problem of segmentation with depth information as a problem of recovering occluded shapes and their ordering relations based on a 2D image. The variables defined in this problem are in three folds: 1) the shapes of the regions to which different objects belong; 2) the ordering relations among objects; 3) the pixel intensities of objects. Without loss of generality, one can assume that the objects in an image are in ascending order, i. e., is the nearest object to the observer while is the farthest one (i.e. background). Let be defined as the visible part of , i.e., , , (). In addition, is defined as the visible background. Based on the above assumptions and definitions, the level set based energy functional is formulated as in 
where are two positive penalty parameters, is the pixel intensity of the -th object, and is the image to be processed. denotes the curvature of boundary for region . Here is chosen to substitute the square in Euler s elastica term with the reason that the the object corners can be preserved when becomes large. The level set function is represented by a continuous signed distance function. and are Heaviside function and Dirac delta function described in detail in [2, 7].
2.2 Segmentation models incorporating noise distributions
investigated noisy image segmentation problems by using specific parameter estimation based on different noise distributions. All the related parameters are calculated via the maximum a posteriori probability (MAP) estimation from the viewpoint of Bayesian probability models. For example, estimation of variance information is used for images degraded with Gaussian noise; The square of image intensity value with capacity of enhancing weak properties is incorporated in the Rayleigh model; Models with great segmentation performance of dealing with Poisson and Gamma noises are built on the standard deviation and average. The general variational model is written as follows
where are positive penalty parameters, is a binary level set as defined in functional (1). Specific representations of function derived from the maximum likelihood method and the computation of their related parameters are given in Table 1. refers to the corresponding parameters of Function need to be estimated.
|TABLE 1 Potential functions of different noise distributions|
|Functions||Gaussian noise||Rayleigh noise|
|Functions||Poisson noise||Gamma noise|
2.3 Progressive hedging algorithm in stochastic programming
As stated in , is a finite set with scenarios . Each scenario has a known probability , and the sum of these probabilities is 1. Let be the mappings that designate responses to in the following form
Here is used to represent the linear space consisting of all such mappings from to . denote iteration steps. Elements and are derived from elements and . denotes for each a variable in . can be determined by solving a separate problem for each scenario to obtain as follows:
and are the projection mappings onto the sub-spaces and of
. The authors indicated that the vectorexists and can be uniquely determined with the reason that the proximal term guarantees the functional being minimized to be strongly convex. More details in terms of theorems and proofs can be found in .
3 Novel formulations for different segmentation purposes incorporating influence of stochastic noises via progressive hedging
Motivated by the research using PHA to solve the minimization problem of stochastic programming, we aim to propose novel formulations tackling different segmentation issues in consideration of the advantages of Euler s elastica term and the influence of unknown noises. In this way, we not only can fulfill general segmentation tasks as well as the classic model, but also can deal with worse situations such as low quality images with large noises, absent boundaries, missing shapes or occlusion. Then we show how to implement PHA with developed ADMM-C algorithm to obtain the optimal solutions efficiently.
3.1 Two-phase segmentation based application
Based on two-phase Euler s elastica based segmentation formulated in (1), we propose a novel segmentation model incorporating influence of different noises expressed as the following stochastic programming (SP) form. The reconstructed contour can be obtained by minimizing the following energy functional with respect to .
represent different noise distributions and contain the stochastic information need to be estimated. is the optimal solution of (3.1) under distribution . Here we continue the definition of in  using binary representation which can also be explained as a substitution . As they described, this binary representation was originally used for finding the global minimizer. And it can also reduce the computational complexity to some extent such as avoiding the necessary calculation associated with level sets. The last two terms added to sub-problem of in functional (3.1) can guarantee the mathematical convergence strictly.
According to the segmentation model for vector-valued images proposed in , the averages of the data terms over all channels are used for coupling. Let be a original color image defined on a domain . Then function should be also in the multichannel form . In fact, our proposed model used to solve color image segmentation can be stated as follows:
where denote the number of layers of a vector-valued image. In this way, we obtain novel Euler s elastica based formulations embedding stochastic noises for two-phase segmentation. In the following section we shall implement the calculation under PHA, which is one useful and effective tool for solving above multistage stochastic programming problem.
3.2 PHA with developed ADMM-C algorithm for two-phase segmentation application
In order to demonstrate the precise numerical procedure of PHA in solving (3.1) and (3.1), we focus on the general model integrating the gray space and color space cases. Detailed proofs for convergence of this algorithm are provided in . The original minimization problems (3.1) and (3.1) are based on separate sub optimization problems,
where refer to noise distributions. For each distribution has a known probability which can be set empirically through experiments, and the sum has to be guaranteed. According to PHA approach, can be obtained by the following steps.
I. First the minimization problems need to be solved separately, which are given in the right side of (8).
where . For the gray space issue, there is only one layer of the image information to be calculated, which means . For the color space issue, should be substituted with the coupling terms . is the potential function for specific noise distribution in each channel of the image. Table 1 shows the representations of and estimations of parameters . is the optimal solution of (3.2) in distribution .
II. Next all of the obtained optimal solutions are utilized to gain the final optimum .
III. At last, the sub problems solutions and the Lagrangian multipliers need to be updated at the end of each iteration.
Then updated , and the parameters derived from (3.2) are passed to the next iteration from step I.
To solve the minimization problems (3.2) separately, both simplification and effectiveness of the algorithm should be considered. In fact, there are three main computational difficulties in the functional as listed below, followed by the corresponding algorithm design in response.
One main challenge of the Euler s Elastica based functional is due to the non-smoothness and non-convexity of . As described in [20, 25, 26], it is more efficient when the curvature term is computed separately from the functional (3.2). Inspired by the concept of curvature weighted approach, we can rewrite functional (3.2) as the following simplified version
The proposed approach essentially reduces the minimization problem (3.2) to a total variation type . Here the case of division by zero in should be avoided. In practice, the denominator is often replaced by ( is a arbitrarily small positive parameter). Then is represented as
Note that the binary constraint for will also cause non-convexity in sub problems (3.2).  demonstrated that certain non-convex minimization problems can be equivalent to the following convex minimization problems
This convex minimization scheme could find global minimizers for (3.2) by thresholding the solution of (3.2), which was classed as a continuous min-cut algorithm. Together with its equivalent form, the continuous max-flow algorithm, the min-cut algorithm has been proved to be an exact convex relaxation of the original problem as stated in [32, 33].
Another critical issue for solving (3.2) is the inevitable high order derivatives in numerical implementation, which is proven to be tedious and prone to errors. Then a developed ADMM-C algorithm is designed to each sub problem by introducing auxiliary variables, Lagrangian multipliers and an alternating directional optimization strategy.
Here the detailed implementation on solving each sub-problem by the ADMM-C algorithm will be presented. Firstly, some auxiliary variables are introduced, i.e., with property and the Lagrangian multipliers . Based on this observation, we can transform (3.2) into the following augmented Lagrangian functional
where is a positive penalty parameter. It is worth noting that this kind of simple structure of (3.2) requires fewer variables and parameters compared with other works [11,18,19] using ADMM to deal with the curvature term directly. After the initialization of , and , a minimization problem is carried out in each step with respect to one variable while keeping other variables fixed temporarily. When alternative optimization for all the variables is completed, the Lagrangian multipliers will be updated subsequently. This gives
To obtain : In the step of the proposed ADMM-C, the average image intensity values as well as variances in the foreground and background can be obtained by using the standard variational method for (15). Table 1 gives all the solutions for distributions .
To obtain : Optimal value of for a certain distribution is obtained by solving the minimization of (3.2) with respect to . We can get the update rule based on the corresponding Euler-Lagrange equations
where . Like in , equation (19) is a screened Poisson equation for which Fast Fourier transform (FFT) [26, 34] is a well-known solver with very low computational cost for imaging problems. Here FFT is applied for further improving the calculation efficiency. Equation (19) can be rewritten as
where and is the discrete inverse Fourier transform. Then we can obtain as follows
For clarity, we present the overall algorithm for the two-phase Euler s elastica based segmentation in stochastic programming in a pseudo code format as follows.
|Algorithm 1. Computing framework for (3.1) and (3.1) via PHA|
|for , do the following steps recurrently|
|1: Obtain via Algorithm 2|
|2: Update via Equation (10)|
|3: Update , via Equation (11)|
|4: if some stopping criteria (given in Section 3.5) are satisfied break|
|Return optimal value after thresholding|
|Algorithm 2 Detailed implementation for step 1 in Algorithm 1 via ADMM-C|
|else solve the following problems alternatively|
|1: Update according to distribution laws|
|2: Update via minimization problem (3.2)|
|3: Update via minimization problem (3.2)|
|4: Update via (18) using gradient ascent method|
The idea of the ADMM-C method is intentionally applied for the main challenge of the non-convex, non-smooth and non-linear problems in Euler s elastics and it has attracted extensive attention. ADMM algorithm has been given analytical properties in [35, 36] and many other applications [37, 38] are provided for which many similar algorithms are developed and successfully used to achieve excellent performances via solving a variety of non-convex problems. In this paper, we adopt a similar idea as in [20, 26] to design a new algorithm to deal with the sub-problems derived from PHA framework.
3.3 Segmentation with depth based application
In this part, we intend to use similar stochastic programming skills as the ones applied in the two-phase issue. Based on original segmentation with depth model for gray space reviewed in Section 2.1, we can establish the energy functional by introducing random noise set as follows
where denote the number of objects in the image, and is set only for consistency of description. We can still extend above segmentation with depth incorporating stochastic noises model to multichannel case. Analog to the coupling approach used in (3.1), the formulation is written as
In order to describe the calculation procedure explicitly, we plan to conduct on the general model integrating the gray space and color space cases. According to the curvature-weighted approach used in (3.2), the simplified version can be directly written as
where and the definition of is given in (). is for gray space issue and
for color space issue. And the characteristic function for the-th region reads . Particularly, representing the -th region: background. Next section we shall calculate above general formulation under PHA with detailed implementation.
3.4 PHA with developed ADMM-C algorithm for segmentation with depth application
For the numerical part, we shall try to use similar ideas as were used in the sophisticated two-phase imaging tasks. The main program loop for PHA is shown below
I. First we need to solve the sub minimization problems of (25) separately, which gives