A fast two-stage active contour model for intensity inhomogeneous image segmentation

This paper presents a fast two-stage image segmentation method for intensity inhomogeneous image using an energy function based on a local region-based active contour model with exponential family. In the first stage, we preliminary segment the down-sampled images by the local correntropy-based K-means clustering model with exponential family, which can fast obtain a coarse result with low computational complexity. Subsequently, by taking the up-sampled contour of the first stage as initialization, we precisely segment the original images by the improved local correntropy-based K-means clustering model with exponential family in the second stage. This stage can achieve accurate result rapidly as the result of the proper initialization. Meanwhile, we converge the energy function of two-stage by the Riemannian steepest descent method. Comparing with other statistical numerically methods, which are used to solve the partial differential equations(PDEs), this method can obtain the global minima with less iterations. Moreover, to promote regularity of energy function, we use a popular regular method which is an inner product and applies spatial smoothing to the gradient flow. Extensive experiments on synthetic and real images demonstrate that the proposed method is more efficient than the other state-of-art methods on intensity inhomogeneous images.


Introduction
Image segmentation is still a popular problem in the field of image processing and computer vision [1]. The intensity inhomogeneity exits in images, which is caused by the imperfections of image acquisition, influence of the illumination and other factors of environment, is one of the main issues in subject of image segmentation. The existences of intensity inhomogeneity in images may lead to misjudgments of doctors and researchers. Therefore, the subject for segmenting images with intensity inhomogeneity attracts more and more researchers to study. Recently, active contour models [2] have been a successful branch for image segmentation.
The existing active contour models can be divided into two major types: edge-based models [3][4][5][6][7] and region-based models [1,2,[8][9][10][11][12][13][14][15]. Edge-based models utilize the image gradient a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 clustering (LCK) model [13], which could obtain ideal segmentation results for intensity inhomogeneous images. In addition, in order to further improve the speed of operation and weaken the dependence of LCK model on the initial position, we divide our segmentation algorithm into two stages. The first stage can roughly but fast obtain the coarse contour near the object boundaries with low computational complexity in coarse space. The second stage can easily get the accurate segmentation results with suitable initialization in original space. Meanwhile, the two minimization problems of the energy function with exponential family [29] are iterated by natural gradient method, which is faster than other methods for solving active contour models to find the global minima.
The remainder of this paper is structured as follows. In the next section we introduce the classic region-based active contour modes and the frameworks of region-based active contour models with exponential family. In Section 3 we derive a fast two-stage segmentation model. More precisely, we develop the mathematical implement for local region-based active contour models. In Section 4 we list the experimental results illustrating the performance of the proposed method. Finally, in Section 5 we describe the conclusion and future work.

Background and theoretical foundations
Traditional region-based active contour models Chan-Vese model. Chan and Vese [2] proposed a global region-based active contour model for image segmentation. This model approximates the image intensities outside and inside of contour by the average image intensities outside and inside of contour, respectively.
Where μ � 0, ν � 0, λ 1 � 0 and λ 2 � 0 are fixed constants. c 1 and c 2 are the average image intensities outside and inside contour C, respectively. H � (ϕ) is the smooth approximate of the Heaviside function and the derivative of H � is defined as Chan and Vese minimized the energy function (1) by the standard gradient descent method [31]. The gradient descent flow corresponding to ϕ can be written as As a global active contour model, the curve convolution of the CV model is only related to global characteristic of the image region. Therefore, the CV model can achieve the satisfactory results for images with intensity homogeneity, whereas it can not deal with the images with intensity inhomogeneity. Local binary fitted model. Li et al. [11] proposed the local binary fitting (LBF) model by considering a kernel function into a local region-based model. The LBF model can deal with intensity inhomogeneity images. The energy function of LBF model is defined as where λ 1 , λ 2 � 0, ν � 0 and μ � 0 are fixed constants. K s ðx À yÞ ¼ 1 ð2pÞ ðn=2Þ s n e À jxÀ yj 2 =2s 2 is Gaussian kernel function with kernel width σ, which is introduced to control the spatial distance between the x-th and the y-th; m 1 (x) and m 2 (x) are the local clusters of the x-th pixel. Meanwhile, in order to regularize the level set function ϕ, the third and the forth terms are added as level set regularization term.
The level set function and the local clusters can be updated from Eq (5) by the standard gradient descent method, which can be computed as Where e 1 and e 2 are The LBF model considers the Gaussian kernel function to capture the local region information of images and can segment the images with intensity inhomogeneity. However, the Gaussian kernel function is not sufficient to introduce the image information and the LBF model cannot segment images with severe intensity inhomogeneity. Moreover, this model is very sensitive to the initialization and the iterative approach of it easily falls into local minima.
Local correntropy-based K-means clustering model. Wang et al. [13] presented a more accurate segmentation method (LCK) for images with unknown complex noise. The main difference between the LCK model and LBF model is that the LCK model utilizes the pixel-tocluster distance, which makes this model is more robust to unknown complex noise. The objective function of the LCK model is The weight of the y-th is calculated by Minimizing Eq (7) by using the standard gradient method, the corresponding local clusters and level set formulation is obtained by The LCK model is robust to images with complex noise and intensity inhomogeneity, whereas it is also a little bit sensitive to the initialization. Meanwhile, the iterative method of the LCK model also easily falls into minima and needs much more times to obtain the final segmentation results.

Region-based model with exponential family observation
Due to the fact that the energy functions of region-based active contour model does not have a structure of a vector space, Lecellier et al. [29] proposed the exponential family observation framework for region-based model. This framework is specifically calculated when coping with global region-based information such as statistical image features (histogram, variance and mean).
Exponential family. In this section, we provide the necessary mathematical concepts on the exponential family, which covers most noise models of the image acquisition system. For the given point x 2 < d , the image values can be distributed by Where the m 1 (x) and m 2 (x) are the clusters of the foreground and background respectively and the function f: < p ! < + is the probability density function with exponential family distribution. Definition 0.1 The family of distributions of a random variable, is called a k-parameter canonical exponential family. If there exists canonical parameter vector η = (η 1 , . . ., η k ) T and log-normalizer A(η), and real-valued functions h, T 1 , . . ., T k : < k ! <, the probability density function f(�|m) may be written as here T = (T 1 , . . ., T k ) T is the sufficient statistic. Note the Definition 0.1 gives the general distribution of exponential family that most of them have distribution for signal and image processing. Table 1 shows some common canonical distributions of the exponential family, such as the Gamma, Beta, Poisson, Gaussian, Exponential, and Rayleigh.
Region-based active contour model with exponential family observation. Following the definition of exponential family, the region-based active contour model (this paper takes LCK model as an example)(Eq (7)) with exponential family observation can be rewritten as The functional optimization problem Eq (12) can be solved by alternative minimization method. In this paper, we take the Gaussian distribution as an example and give the iteration with respect to ϕ, m 1 (x) and m 2 (x) as following: Although we give more general distribution for LCK model, it still can not overcome its inherent limitations. The segmentation method of the above LCK models is not convex and the weights of the pixels w x are always decreased faster than other region pixels during the iterated process. Hence, the main difficulty for the region-based active contour model is still the solution of energy function.

Proposed model
Inspired by the work of Wang et al. [13] and Pereyra et al. [28], we proposed a novel two-stage segmentation method, which can fast segment the intensity inhomogeneous images with accuracy results. In the real-world, the accuracy segmentation results for intensity inhomogeneity images are difficult to obtain and also need to spend a large amount of time, because these intensity inhomogeneity images are always big and complex. Due to the above factor, we split the process of the segmentation into two stages which makes the computation with lower computational complexity. In the first stage, we implement the segmentation model in coarse spaces. Following it, in the accurate segmentation stage, we utilize the up-sampled coarse results as initializations and segment the images in fine space. Finally, to get more accurate segmentation results faster, we present the Riemannian steepest descent method for the two-stage segmentation model. Table 1. Some common distribution of the exponential family.

Coarse segmentation
The LCK model is robust to images with severe intensity inhomogeneity and complex noise. However, as a local region-based active contour, this model needs a perfect initialization to obtain the satisfactory segmentation results. In the first stage, we propose a novel LCK model with exponential family observation in coarse space (RDLCK), which can obtain the segmentation results fast. The energy functional of this model can be defined as Where O τ denotes the uniform sub-space of O with the down-sampled factor τ, and function f is exponential family function. w τy is computed by Hð� t ðyÞÞgðkIðyÞ À m t1 ðxÞk 2 Þ þ ð1 À Hð� t ðyÞÞÞgðkIðyÞ À m t2 ðxÞk 2 Þ dy; and g(x) is Gaussian kernel function g(x) = exp(−x 2 /2σ 2 ). m τ1 (x) and m τ2 (x) are the local clusters in the coarse space, which can be calculated by gradient descent method. For f is Gaussian distribution, then m τ1 (x) and m τ2 (x) can be computed by Most traditional region-based active contour models always converged by using the standard gradient descent method, whereas it easily falls into minima and needs a large number of iterations. Therefore, fixing m τ1 (x) and m τ2 (x), we utilize the Riemannian steepest descent method for RDLCK model. The update for level set function ϕ τ is Here γ τt is a positive parameter, r � t E RDLCK ðIðxÞ; � t t Þ is the standard gradient flow, which is similar to Eq (12) Therefore, the key to Eq (16) is that the N t is a positive definite matrix. (refer to the appendix for a detailed computation for N ).
This paper makes the computation of the N t À 1 ð� t t Þr � t E RDLCK ðy; � t t Þ possible. Moreover, the fast image segmentation method based on natural gradient Eq (15) is applied to local regionbased active contour model successfully.
Lastly, to promote the solutions of Eq (13) smooth, we need to regularise ϕ. The traditional region-based active contour model is regularised by adding the penalty term l R O j d½�ðxÞ� j dx, which regularises the energy function by selecting contours with minimal length. In this paper, we promote regularity by using an increasingly popular alternative, which applies spatial smoothing to the natural gradient flow. This method can be written as Where h s ðs; uÞ ¼ 1 2ps 2 exp À s 2 þu 2 2s 2 À � and the low values of σ are applied for images with low noise, which preserves details of images.
We can obtain the coarse contour ϕ τ in the down-sampled space O τ , but the coarse contour ϕ τ is not suitable for the original image due to down-sampling process. Therefore, we need to up-sample the coarse contour ϕ τ and get the up-sampling version ϕ � of the coarse contour ϕ τ . As shown in Fig 1, we give the segmentation process of a heart CT image (the size of the image is 152 × 128).

Accurate segmentation
The coarse segmentation can obtain the coarse contour with low cost of computation, whereas the coarse contour is not accurate due to the loss of region information from down-sampling. Therefore, we further segment the objects by improved LCK model with satisfactory initialization ϕ � . The improved LCK model with exponential family in original space (ROLCK) is defined as where ξ is a small positive parameter and ϕ � is the up-sampling vision of the coarse contour.
The last term promise that the distance between ϕ � and accurate result is small. Using the Riemannian steepest descent method for converging ROLCK model, the iterative scheme for minimizing the Eq (19) is similar with the one for solving Eq (13).
Where γ t is positive constant and w x is the final weight of the x-th pixel, which is calculated by H(ϕ(x)))g(kI(x) − m 2 (y)k 2 ), the computation of N is also similar to the computation of N t ðN ð�ÞÞ ðx;yÞ ¼ K s ðxÞ � The only difference is the iteration of r � E ROLCK � ðIðxÞ; � t Þ, it is computed by At the original space, there is no region information loss. Since the initialization for Eq (19) is very close to the object, the convergence of Eq (19) is fast and the accurate segmentation stage can obtain the final result after a few iterations. An optimization active contour model

Segmentation procedures
The procedures of the proposed method can be summarized in Algorithm. 1. The initial level function of the coarse segmentation is defined as We provide a positive parameter for ρ in this paper. The proposed method can be summarized in Algorithm 1.

Experimental results and comparison
In order to prove and compare the performance with the CV model [2], the GCV model [10], the LBF model [11], the LIC model [12], the LCK model [13], the RCV model [28] and the LGFI model [14], we apply our approach to several images with intensity inhomogeneity in this section. All the experiments are implemented by using MATLAB R2015a and running on a person computer with Intel core i5, 2.5GHz, and 4.00 GB RAM.

Comparisons with state-of-art methods for intensity inhomogeneity images Comparisons with global active contour models (the CV [2] model, GCV [10] model and RCV [28] model.
In this section, we compare our method with some global active contour models. The global active contour models assume that the intensities of images are piecewise constant and thereby these models can not segment the images with intensity inhomogeneity. Fig 2 shows the segmentation results on two synthetic images (The first row is a T-shape image and the fourth row is a synthetic image.) and two medical images (The second row is a blood image and the third row is a CT heart image.) with intensity inhomogeneity by using the CV model [2], the GCV model [10], the RCV model [28] and the proposed method, respectively. The first column shows the segmentation results with CV model. The second column shows the segmentation results with GCV model. The third column shows the segmentation results with RCV model. Finally the last column shows the segmentation results with our method. It can be seen that the results of the CV model [2], the GCV model [10] and the RCV model [28] are unsatisfactory because these models are global active contour models and use the global intensity means of images to fit them. On the other hand, our method can achieve satisfying results due to taking the local image region information into account, which can deal with the intensity inhomogeneous images well.

Comparisons with local active contour models (the LIC [12] model, the LCK [13] model and the LGFI model [14]).
We demonstrate the superior performance of the proposed method to previous local active contour methods on a synthetic image (see the third row of Fig  3) and three medical images (The first and last row give two brain images, whereas the second row gives a X-ray image.) with severe intensity inhomogeneity (see the last column of Fig 3). Fig 3 shows the results of the LIC model [12], the LCK model [13] and the LGFI model [14] and our method, all of which consider the local region information and fit them by local An optimization active contour model intensities means. Therefore, these methods can achieve much better results for images with intensity inhomogeneity than global active contour models. However, The LIC model [12] only uses the local mean information and thereby it fails to segment the objects when the intensity inhomogeneity of images are severe, which can be seen from the first column of Fig  3. Although the LCK model (see the second column of Fig 3) and LGFI model (see the third column of Fig 3) consider the gaussian kernel function as constraint and improve the segmentation results to some extent, these models still fail to discriminate the intensities of the objects. Due to the proper initialization from the coarse segmentation, our method obtain more satisfactory results (see the forth column of Fig 3) than other methods. Table 2 introduces the computation complexity for Fig 3 in terms of total iterations and CPU time. The proposed method achieves the accurate segmentation results with less CPU time and less number of iterations. The LIC model [12], the LCK model [13] and the LGFI model [14] are all solved by standard gradient method, which take a large number of iterations to converge and make the level set function easily fall into local minima. In order to overcome the drawback of the standard gradient method, our model constructs a steepest descent on the model's intrinsic manifold which converges extremely fast. Moreover, we split the process of segment into two stage and the second stage can get a suitable initialization from the first stage. This experiment proves our method can obtain the best segmentation results with the least time.
Results on natural images with intensity inhomogeneity. In Fig 4, we have shown the performance of our method for natural images with different types of region properties. In the    Table 3, from which we can see that our method converges obviously faster than other methods. On the one hand, our method converges with the Riemannian steepest method is much faster. On the other hand, the first stage provides the proper initialization for second stage. Therefore, our method is also applicable to noise images.

Quantitative comparisons with the-state-of-art models
To quantitatively analysis the performance of the proposed method, we adopt the F-score value [13] and Jaccard similarity (JS) [32] metrics as the evaluation framework of segmentation accuracy. The metric JS is defined as where A is the computed object region and B is the ground truth region. Obviously, the closer the value of JS is to 1, the better segmentation result is obtained. Furthermore, to evaluate the results more precise, we also give the value of F-score, which can be calculated by Where TP = A T B (true positive) corresponds the correct segmented regions, FN (false negative) corresponds the false unsegmented regions, FP (false positive) corresponds the undetected segmented regions. Similar to the value of JS, the higher value of F-score means the better segmentation results.
Robustness to the severe intensity inhomogeneity. Fig 6 shows the results of competing method segmentation on five synthetic images with changed intensity inhomogeneity. The first row of Fig 6 shows the five original images with initial contour. As can be seen, the object of each image is gradual difficult to be segmented from top to bottom. The segmentation results obtained by CV model are presented in the second row, GCV model are presented in the third row, RCV model are presented in the forth row, LBF model are presented in fifth row, LCK model are presented in the sixth row and the proposed model are presented in the last row, respectively. We can see that our method provides the best segmentation results, especially when the strength of the intensity inhomogeneity for objects is strong, which proves that the proposed method is more robust to image intensity inhomogeneity. The corresponding Fscore values and JS values are shown in Fig 7, all the models can achieve high F-score values and JS values, when the strength of intensity inhomogeneity is slow (see the first column and second column of Fig 6). Meanwhile, the total iterations and CPU time of all the models, which are shown in Table 4, are close. Whereas, when the strength of intensity inhomogeneity becomes strong, the F-score values and JS values of the proposed method are still close to 1 and higher than other models. Furthermore, our method converges with 2 steps for all images, even the object of image is difficult to be segmented (see the last column of Fig 6). The main difference between our method and other local active contour is the proper initialization and iterative method.  An optimization active contour model the Fig 8, we present the segmentation results of the proposed method for images with different level of speckle noise. As can be seen, the proposed method can achieve the satisfactory results even the noise is very strong. Furthermore, the F-score values and JS values for proposed method, which are shown in Fig 9, are also higher than 0.9, so the accuracy of the proposed method is ideal.

Conclusion
In this paper, we propose a fast two-stage segmentation method for the segmentation problem of synthetic and real-world (including medical) images with complex noise and severe intensity information. The first stage obtains coarse segmentation results in coarse space with low calculation complexity. Following it, the second stage takes the up-sampled contour of the first stage as the initialization and achieves accurate segmentation results in original space quickly. To minimize the energy function with low time complex, we propose a method measuring intensity under the framework of exponential family and converge the model by Riemannian   An optimization active contour model steepest method. Finally, a smooth operator, which is Toeplitz and separable, is used to regularize the level set function to preserve more image details.
The main contribution of this paper is the formulation of two-stage energy function based on LCK model with exponential family for acquiring the intensity information and solving the two-stage method by smooth natural gradient method. Qualitative and quantitative analysis showed that the performance of the proposed method is superior to the LCK model and other state-of-art methods for intensity inhomogeneity images.
There are still some limitations for the proposed method, such as model still depends on the initializations to some extent and model is difficult to segment images with wispy and clutter targets. In the future, we plan to add shape constraints to optimization model and focus on developing the application of the proposed methodology for more types real-life images.   An optimization active contour model ð30Þ if x = y and (G(ϕ)) (x,y) = 0, then N t ð� t Þ ¼ 0. Otherwise, we need to prove the diagonal elements of matrix N t ð� t Þ are strictly positive. In Eq (30), the E(. . .|ϕ τ (x)) is the expectation with respect to the marginal likelihood According to the fact that δ 0 (−x) = −δ 0 (x) and the Eq (30)  Finally, in this paper the function f belongs to exponential family and we can rewritten (33) in terms of Bregman divergences [34] ðN t ð� t ÞÞ ðx;yÞ ¼ K s ðxÞ � fjd 0 � t ðxÞw tx jB f ðm t1 ðxÞjjm t2 ðxÞÞg if � t ðxÞ � 0