#### Abstract

Image segmentation is a fundamental task for many computer vision and image processing applications. There exist many useful and reliable models for two-phase segmentation. However, the multiphase segmentation is a more challenging problem than two phase segmentation, mainly due to strong dependence on initialization of solutions. In this paper we propose a reliable hierarchical algorithm for multiphase texture image segmentation by making full use of two-phase texture models in a fuzzy membership framework. Application of the new algorithm to the synthetic and real medical imaging data demonstrate more satisfactory results than existing algorithms.

#### 1. Introduction

Image segmentation is a fundamental problem in image processing which is a prerequisite to high-level computer vision applications. It aims to divide an image representing a real scene or a synthetic one into classes or categories, corresponding to different objects and the background in the image. In the end, each pixel should belong to one class and only one. In other words, we look for a partition of the image into distinct segments, and each of them shares some features in common such as intensities, color, or texture. In particular, image texture defined by repeated patterns of intensities adds much complication in image processing tasks. A textured image often has several regions with different textures in existence, and the task of segmentation is to locate the texture boundaries, with which this paper is concerned.

Over the past two decades, a variety of different techniques have been developed to solve the problem of image segmentation, ranging from region growing and emerging [1], watershed algorithms [2], minimum description length criteria [3], and active contour models [4, 5] to Mumford-Shah energy minimization model [6].

Historically, image segmentation is known as the process to segment an image into two categories of regions: the foreground and the background. This process nowadays is referred to as a two-phase modeling, whereas multiphase modeling is specifically to deal with segmentation of more than two regions. Earlier work on segmentation attempts to detect the feature boundaries directly by edge detection [7β9]. These methods are susceptible to noise which is often present in real applications, and as such they are not suitable for processing either textured or noisy images, unless one applies it to the transformed image after applying a Gabor-type filter to the original one [10].

For nontexture images, the most influential model, known as MS-model, is proposed by Mumford and Shah [6], where the required boundary set of features as well as the segmented image is determined from minimizing an energy functional. Unfortunately this elegant model is numerically difficult to realize, as such considerable effort has been made in order to alleviate this problem. For instance, Ambrosio and Tortorelli [11] have proposed a more solvable model by approximating the MS-model with elliptic functionals defined on Sobolev spaces. However, the most well-known paper based on [6] is the algorithm proposed by Chan and Vese [12]. Known as the CV-model, this was initially designed for two-phase segmentation of images of approximately piecewise constant intensities, using the framework of level set functions [13].

Segmentation of texture image is intrinsically more challenging than intensity-based ones. For texture images, detecting edges directly is problematic as it would treat texture patterns as image features. One of the strategies is to derive texture features from the image by first applying certain filters to it, and subsequently these features will be treated as a multichannel (i.e., color) segmentation problem. For instance, Sandberg et al. [14] have used Gabor filters for texture image segmentation by extending the original CV-model [12], and a similar strategy has also been adopted by a later work [15]. On the other hand, there have been other powerful models proposed to tackle texture segmentation of texture image directly [16β18]. Especially, both the region competition model of Zhu and Yuille [16] and the nonparametric model of Kim et al. [17] minimized mutual information as fidelity measure plus the regularization term. More recently Ni et al. [18] proposed a more reliable model based on Wasserstein distance as fidelity measure which leads to global optimal solutions to segmentation of texture image, where local histogram information is extensively used in order to divide the image domain into two regions in each of which the difference of the cumulative distribution function from its median is minimal.

Mory and Ardon [19] used the concept of fuzzy region competition to unify and extend it to model texture segmentation problems. Here, each pixel is assigned a probability, instead of a precise membership integer, of belonging to a particular region. The fuzzy region competition has been further generalized to deal with multiphase segmentation problems by Li and coworkers [20, 21], and in particular the multiphase texture segmentation model [20], referred to as LN-model, will be discussed in the following sections.

For the past decade, research in multiphase nontexture segmentation has been active with a focus on grayscale images. The work on multiphase texture segmentation is more recent.

A number of multiphase models have been proposed most of which are generalized from two-phase ones; see [22β30]. Gao and Bui [26] proposed a hierarchical model for multiphase segmentation problem for piecewise smoothing image. Jung et al. [31] formulated a multiphase segmentation model built upon the celebrated phase transition work of Modica and Mortola in material sciences. An efficient algorithm for minimizing the piecewise constant Mumford-Shah functional of image segmentation was proposed [32] based on the threshold dynamics of Merriman et al. [33] for evolving an interface by its mean curvature. More recent work include models based on shape and topological sensitivity [28] and regularization model [34]. Salah et al. [35] recently proposed to use a kernel function to map the original image into data of a higher dimension so that the piecewise constant model becomes applicable. Inclusion of shape constraints into the multiphase segmentation was also explored by Cremers [36]. Several models for multiphase segmentation of image frames have also been proposed [29, 37]. We also remarked that there is a very interesting unsupervised multiphase segmentation model from [30], which can automatically determine the number of regions during the segmentation process.

It has been observed that most generalized models work well only for a small number of problems, and robustness is a major issue. For instance, the above-mentioned CV-model [12] was generalized to multiphases in [22] and then refined in [23]. However as noted by [38, 39], the new model [23] has a strict requirement on its initial guess which contradicts the idea of automatic segmentation. Moreover, the phase number has to be , where is the number of level set functions. To improve this model, the idea of hierarchical implementation of the robust CV-model [12] was first considered in [38] (with time-marching solver a.k.a. additive operator splitting (AOS)) and was further improved by Badshah and Chen [39] with a robust multigrid solver. Along this line, Ni et al. [40] employed the fast time-marching dual algorithm proposed by [41]. These improved hierarchical models have produced convincing results in segmentation of intensity images but not yet applied to segmentation of texture images.

In clear contrast to two-phase segmentation, multiphase texture segmentation is relatively understudied so far. Several feature-based models have been proposed to tackle the problem by working on features derived by different strategies. For instance, Aujol et al. [42] have proposed a model based on features derived from applying wavelet transforms to the image of interest. Similarly, Wei and Xin [43] proposed a supervised model based on contourlet features for aerial image segmentation. These feature-based models have limitation in choosing appropriate descriptors of texture. The most recent Li-Ng (LN) model [20] was based on mutual information which implies that there is no need to detect features first as is required by feature-based models. However, this model is not globally convex, as such the performance is sensitive to initialization; Li and Ng [20] have remarked that this occurs in segmentation problems with three or more phases. As such more advanced segmentation models are needed.

We finally remark that there have been recent advances in efficient solvers for variational models. Different from the past solvers such as gradient descent level set methods, a diversity of fast solvers have been proposed; these include graph cut [44], multigrid [39, 45], dual projection algorithm [41], as well as genetic algorithms [46]. Amongst these, alternating optimization strategies via the elegant dual projection model have been increasingly employed in the literature [20, 40].

The rest of the paper is organized as follows. Section 2 reviews some existing segmentation models, paying special attention to fuzzy region competition models which assign a probability value at each pixel rather than an integer for phases. Section 3 introduces our new hierarchical algorithm based on general two-phase models. Section 4 shows a series of experiments for comparisons and verifications. Some conclusions are drawn in Section 5.

#### 2. Review of Some Segmentation Models

We shall review five models that are reliable for two-phase segmentation. The first two models are designed for non-texture segmentation, and the last three models are for texture segmentation; in particular the fourth model by Ni et al. is the most reliable as well as the most expensive.

Let be a bounded open subset of with its boundary, and let be the given grayscale image. The aim of multiphase segmentation is to partition into regions , , where and .

Mathematically, each phase may be characterized by a distinct constant as defined in [5] and in two phases with , one may simply try to find a piecewise constant solution taking two constant values in (or two other constants for membership identity). However it is often advantageous to look for a solution that takes values in in a fuzzy membership framework [19].

##### 2.1. The Mumford-Shah-Based CV Model

A general form of two-phase segmentation problem can be represented [19] as where is the (unknown) union of all boundaries of , and are two positive and suitable multipliers; often we choose . The two functions , and are respective similarity measures of pixels (or region errors) in domains , and with their representative values.

The well-known CV model [12] chooses in (1). Here, takes values in . The model can be adapted to tackle the multiphase segmentation problem [38β40] when more constants are defined.

The level set formulation (as a means to rewrite the above integrals over as over ) is commonly used to represent this model, and subsequently the resulting partial differential equation (PDE) is solved by a gradient descent time-marching method, as with the celebrated CV model [12]. However, this model (1) is not convex (neither is the CV model), as such the level set methods need a reasonable initialization in order to avoid local minima. This prerequisite and the slow convergence are inherent weakness of this PDE-based method.

##### 2.2. A Global Convex Model

Recently, Bresson et al. [47] have tried to reformulate the above problem into a convex one so that the global minimum becomes easier to compute. The idea is to introduce constraints and to convert the two-phase model into a convex total variational (TV) model as follows: It should be remarked that Bae and Tai [48] have proposed a convexified model for four-phase segmentation.

Further, the Chambolleβs pioneering work [41] of fast dual projection can be used to provide an elegant efficient solver for this equation. More specifically, after an auxiliary variable is added to it, (2) will be formulated as where the convex form here is to force and to be close to each other (while serving the purpose of decoupling the nonlinearity), is a small parameter to penalize the error between and . One hopes that, on convergence (with smaller and smaller ), is obtained. Then the model (3) becomes where the minimization problem (4) can be efficiently solved by a fast dual projection algorithm [41] and (5) is solved explicitly. The derived solution is where can be solved by a fixed point method as follows: where is some suitable time step. The solution of (5) can be derived as

Thus, we have discussed a simple algorithm for and . We note that similar use of an intermediate variable such as here can be seen in other works for example, [49].

##### 2.3. The Fuzzy Region Competition Form of Two-Phase Segmentation

The fuzzy region competition model by Zhu and Yuille [16] or Kim et al. [17] for texture segmentation chooses where denotes the probability density function (usually a Gaussian) for a conditional probability. Although Zhu and Yuille have used parametric representation [16], and Kim et al. used nonparametric representation [17], both of them have used standard PDE solution methods. So these models again are not globally convex; as such solvers may suffer from getting stuck at local minima. Of note, if the probability density is a Gaussian distribution and the variance is known, this model will then reduce to the aforementioned CV-model.

##### 2.4. The Local Histogram-Based Model

Rather than the mutual information-based strategies, Ni et al. [18] proposed a new strategy using local histograms (i.e., Wasserstein distance, instead of image intensities) to define a region in order to segment texture images. First of all, a Wasserstein distance with exponentββ1 (measuring the distance of two histograms , and with the range ) is defined as where for each , is the cumulative histogram for ; usually for intensity images. At each pixel , an intensity-based histogram of a small local ball of radius centered on is defined by for each prescribed . Then in the framework of (1), Ni et al. [18] chooses That is to say, an image is segmented according to how accumulated local histograms of a group of pixels are close to a fixed value. We remark that this approach so far has only been used in two-phase texture image segmentation. Since both and are involved, the model is expensive to apply.

##### 2.5. The LN Multiphase Texture Model

Derived from fuzzy region competition method, the multiphase segmentation problem for texture images can be formulated through the concept of mutual information as follows [20]: where is the nonparametric probability density function (pdf) that is estimated by Parzen window by the intensity values of pixels in region . In the case of , the energy in the second term will be similar to the one by Kim et al. [17], who firstly proposed it for texture segmentation.

This energy can be reformulated in a total variational framework as follows: where is the fuzzy membership vector, and is the corresponding pdf vector and is then approximated by By relaxing , the above energy is further approximated by subject to , for , where .

Then minimization of this energy was solved following the fast dual projection model iteratively and in each step:(i)for fixed , , , ;(ii)for fixed and , , . Here, can be solved by a fixed-point method similar to the two-phase problem, by initializing , (iii)for fixed , and , update as in [20].

In summary, the above two-phase models are reliable (when ), while the multiphase texture models are not (when ). Our idea below is to make full use of the reliable two-phase models for multiphase texture segmentation.

#### 3. A Hierarchical Algorithm for Texture Segmentation

Since most multiphase segmentation models solve a nonconvex minimization problem, the solution methods can be seriously affected by initialization or easily get stuck in a minimizer, when the number of phases is larger than two. It should be remarked that there has been some promising progress recently [48] that can reformulate models into convex ones which have a global minimizer but no local minimizers. To improve on [20] when , we pursue an alternative approach in this paper.

Motivated by [39] for a hierarchical algorithm and [1] for seeded region growing models, for non-texture images, we now propose a hierarchical algorithm for segmenting a multiphase texture image.

Our idea is the following. Assume that each phase contains a seeded point indicating the location of a phase; practically one simply clicks on a few points within the image to define seeds. We recursively apply a two-phase model (3) to split a partitioned region into two further subregions and generates an ordered binary tree to represent the structure of the image, as is illustrated in Figure 1. The entire image domain is initially defined as , representing the initialization, which can be arbitrary, we may conveniently choose the image itself (normalized to range from 0 to 1) as initialization for . Then, our two-phase model (i.e., mutual information or local histogram) is to segment the given image , representing the entire domain, into two phases (a domain and its complement ) using fast dual projection strategy. The partitioned regions and will be stored as the tree nodes. For each of and , there is a decision here to make to determine if they are further segmented. This can be done by checking how many seeds a subregion contains, and further segmentation is only carried out using the two-phase model again if more than 1 seed point is contained. Figure 1(a) illustrated the work flow and at the end the original image was successfully segmented into five regions (i.e., , , , , and ). Of note, region term can either be represented by mutual information or Wasserstein distance, or even an alteration between them may be implemented. Other strategies such as Mory and Ardonβs strategy [50] as cost function might be useful, or mixed form during the iteration might be useful.

**(a) Example 1 (Figure 5(a)) demonstrating segmented phases**

**(b) Example 2 (Figure 3(a)) demonstrating segmented phases**

**(c) Entropy map of example 1**

**(d) Intrinsic texture map of example 1**

**(e) Entropy map of example 2**

**(f) Intrinsic texture map of example 2**

For multiphase segmentation problems, it is always assumed that the number of phases are known and more than two, so at least one 2-phase segmentation will be unconditionally performed. For instance, for a 3-phase problem as illustrated in Figure 1(b), for and , at least one of them has to be further segmented, but the program has to follow a rule to decide which one to continue. On the other hand, for the 5-phase problem Figure 1(a), both and should be further segmented.

In summary, our algorithm can be stated in the following steps.(i) initialization: Initialize and set seed points interactively.(ii) two-phase segmentation by(a) initialization: initialize ;(b) iterations:βupdate by (8), where will be replaced with the strategy adopted, that is, mutual information strategy or local histogram;βupdate by (6);βupdate *region term* defined either as or local histogram distance metrics;(c) termination;(iii) for each detected phase, run the previous step if there are two or more seeds in it.(iv) termination strategy: no further segmentation is possible.

Finally, we show some results on designing alternative and possible indicators that may be used to decide on further segmentation without using seeds. For example, region size, change of uniformity described by standard deviation, or even the difference between the child phases or between them and parent one, are candidate choices for this purpose. None of the indicators we tested are suitable for texture segmentation. Our analysis established that conventional measures that worked well for non-texture segmentation problems cannot be readily applied. In searching for new strategies, a measure that can reflect the uniformity of texture would be particularly attractive. From the many different options we tested, due to the limit of space, here we only discuss two of them: entropy and a new texture feature based on semilocal image information [51]. Maps of both of them were given in Figures 1(c), 1(d), 1(e), and 1(f), respectively. Values of mean and standard deviation (STD) of each region are presented in Tables 1 and 2. For the entropy one, an entropy map was produced by computing entropy values of a small window centred at each pixel in the image. By doing this, we wish we can reuse a concept similar to the work on non-texture images, unfortunately, this does work as stated. For the 3-phase problem, values of mean entropy STD are and for phases and , respectively. If using STD as criteria, should be further split because of its nonuniformity, but this is incorrect. The second measure fails as well for this example. Using these measures for the 5-phase problem are also problematic mainly because that they do not monotonically decrease when the number of different textures becomes smaller. For instance, STD values of are lower than . Combinations of these have also been considered, but this becomes a tuning process for threshold values from image to image. All of them are applicable only to a specific problem and a generic one proved to be nonexistent due to the complexity of texture. This automation is certainly one of the future research directions.

#### 4. Experimental Results

We first tested our algorithms on synthetic texture images and real medical images including optical coherence tomography (OCT) images. Note, here for the testing purpose, we set , , and . For simplicity, we denote by:(i)M1βthe method of CV model [12];(ii)M2βthe method of mutual information by Li and Ng [20];(iii)M3βthe method of local histogram by Ni et al. [40];(iv)M4Iβthe method of our framework using mutual information as fidelity term;(v)M4IIβthe method of our framework using local histogram as fidelity term.

##### 4.1. Numerical Tests and Comparisons

Below, we use 2 sets of experiments to show that direct multiphase segmentation models can be sensitive to initialization, while our new algorithm is robust.

###### 4.1.1. Comparisons of Two-Phase Models

In this first set of tests, we applied the CV-model (M1) and two texture segmentation models (M2 and M3) to synthetic and real texture images, see Figure 2. All three models performed reasonably well on a noisy synthetic grayscale image, but for texture image, the CV-model (M1) was no longer able to detect the object with textures. On the other hand, both texture models (M2 and M3) are good for two phases, while it appears that the mutual information one (M2) performs the best.

**(a) Original image**

**(b) Result with M1**

**(c) Result with M2**

**(d) Result with M3**

**(e) Original image**

**(f) Result with M1**

**(g) Result with M2**

**(h) Result with M3**

**(i) Original image**

**(j) Result with M1**

**(k) Result with M2**

**(l) Result with M3**

**(m) Original Image**

**(n) Result with M1**

**(o) Result with M2**

**(p) Result with M3**

###### 4.1.2. Segmentation of Multiphase Texture Images

In this second set of tests, we show that the M2 is good for three phases but can also fail for these phases if the initial guess is not chosen accurately enough, see Figure 3. Other two models M1 and M3 were not tested as the CV-model (M1) cannot cope with texture as shown in the previous, example while the M3 is not yet developed for multiphase segmentation.

**(a) Original image**

**(b) Random initialisation**

**(c) Results from random initialisation**

**(d) User initialisation**

**(e) Results from user initialisation**

**(f) Original image**

**(g) Random initialisation**

**(h) Results from random initialisation**

**(i) User initialisation**

**(j) Results from user initialisation**

We also applied our new segmentation algorithms M4I and M4II to the same three phase images and a five-phase segmentation problem, respectively, and the results are pleasing as expected; see Figure 4 for the results of three-phase problems and Figure 5 for the five-phase one.

**(a) Original image**

**(b) Seed points**

**(c) Result with M4I**

**(d) Results with M4II**

**(e) Original image**

**(f) Seed points**

**(g) Result with M4I**

**(h) Results with M4II**

**(a) Original image**

**(b) Seed points**

**(c) Result with M4I**

**(d) Result with M4II**

##### 4.2. A Medical Imaging Application Segmentation of OCTs

Optical coherence tomography (OCT) is a noninvasive imaging technology that can reveal cross-sectional information of the biological tissues, and it is a powerful tool assisting diagnosis and management of a wide range of eye diseases. All OCT images are extremely noisy and none of the existing denoising models work well for such OCT images. In this section, we show a useful application of our proposed segmentation strategies to some OCT images, when treated as texture images. A typical OCT image of the retina demonstrates multilayers of the retina tissues, see Figures 6(a) and 6(e). The acceptable results by our M4I and M4II, as one can see from Figure 6, demonstrate the potential capability of our strategies for segmenting multiple retinal layers in the OCT images.

**(a) Original image**

**(b) Seed points**

**(c) Result with M4I**

**(d) Result with M4II**

**(e) Original Image**

**(f) Seed points**

**(g) Result with M4I**

**(h) Result with M4II**

#### 5. Conclusions

Multiphase texture segmentation problem represents a significant challenge in image processing. There is a lack of robust and reliable models in the literature. In this paper, we proposed a reliable hierarchical multiphase texture segmentation strategy based on recursive use of the two phase approaches. We have unified recently proposing two-phase strategies of representing region errors, that is, mutual information-based and Wasserstein distance-based, into the same framework. Our application to the synthetic and real medical imaging data demonstrated satisfactory results. One of our future work is to continue investigating alternative strategies for further segmentation (subdividing). We shall also investigate how to use our new algorithm for OCT layer segmentation for clinical research use.