Fuzzy Bit-Plane-Dependence Region Competition

Abstract

This paper presents a novel variational model based on fuzzy region competition and statistical image variation modeling for image segmentation. In the energy functional of the proposed model, each region is characterized by the pixel-level color feature and region-level spatial/frequency information extracted from various image domains, which are modeled by the windowed bit-plane-dependence probability models. To efficiently minimize the energy functional, we apply an alternating minimization procedure with the use of Chambolle’s fast duality projection algorithm, where the closed-form solutions of the energy functional are obtained. Our method gives soft segmentation result via the fuzzy membership function, and moreover, the use of multi-domain statistical region characterization provides additional information that can enhance the segmentation accuracy. Experimental results indicate that the proposed method has a superior performance and outperforms the current state-of-the-art superpixel-based and deep-learning-based approaches.

1. Introduction

Image segmentation, which aims to partition an image into homogeneous regions [1,2,3], is one of the most challenging problems in computer vision and has various applications, such as pattern recognition and medical imaging. Among the existing approaches, variational methods have been extensively investigated and the best-known approach is the classical active contour model [4,5,6], which adopts a boundary function to evolve the curve on the object boundary. Nevertheless, this approach is sensitive to noise, and thus, an unsatisfactory segmentation result may be obtained for a noisy image. In order to obtain a smooth segmentation from a noisy image, Tong et al. [7] proposed a multi-scale approach to suppress various types of noises by minimizing a boundary threshold surface function. Although these methods work well for some natural images, they largely consider edge information without involving other essential image features, and hence, fail to segment images with complex texture structures. As well as boundary-based methods, region-based approaches are also popular methodologies that make use of both the boundary and region information. The most influential approach is the variational model proposed by Mumford and Shah [8], which estimates an image by the piecewise smooth functions with regular boundaries. Chan and Vese [9] studied a particular case of the Mumford-Shah model using the piecewise constant functions and applied their model to a two-phase image segmentation problem. They proposed a curve evolution technique with a level set [10] formulation to minimize the energy functional. The Chan-Vese model was then extended to vector-valued [11] and multiphase [12] cases. However, experimental results indicate that the piecewise constant functions can only be applied to simple images with homogeneous structures and fail to segment images with complex texture patterns. In [13], the authors proposed a non-convex and convex coupling variational segmentation model based on the total generalized p-variation regularizer and Mumford-Shah model to preserve the boundary and detect the structure in the image. Region competition [1] is another widely used variational image segmentation method that penalizes the boundary length and the Bayes error in each image region is characterized by a probability distribution. Based on the classical region competition model, a novel clustering-based region competition approach [14] was proposed for image segmentation. A non-parametric approach using information theory and curve evolution [15] was also studied, with remarkable success.

The aforementioned approaches adopt binary segmentation that gives hard results. Nevertheless, hard-labeling may not be the optimal scheme because natural images usually have a small dynamic range or limited spatial resolution, which blurs the distinction between image regions, leading to the degradation of boundary identification. The fuzzy segmentation method [16,17,18,19,20,21], which applies a membership function valued between zero and one to measure the association degree of each pixel to all regions, is an alternative approach that has been frequently used in data mining, medical images [22,23,24], and so on. The major advantage of a fuzzy approach is that the optimization problem is convex with respect to the membership function. As a result, the model is not sensitive to initialization and the global minimum may be found explicitly. In [2], Li et al. applied a piecewise constant fuzzy region competition model for multiphase image segmentation. The energy functional was solved using the alternating minimization procedure in which the closed-form solutions were obtained. Experimental results for typical gray-scale and color images have shown satisfactory results. A variational multiphase image segmentation model with the fuzzy membership functions and L1-norm fidelity was proposed in [25]. An alternating direction method of multipliers was adopted to solve the energy functional. Experiments demonstrate that this approach is robust to impulse noise and thus provides satisfactory results. As with the Chan-Vese model, however, the segmentation accuracy of the above methods will be lower for complex images because the piecewise constant functions are not suffice in characterizing indispensable image region information. Recently, a robust active contour segmentation by a fractional-order differential method and fuzzy statistical information of boundaries was proposed for vascular image segmentation [26]. Luo et al. [27] developed an unsupervised multi-region method based on fuzzy active contour model for segmenting SAR images. Compared with the level set-based framework, this method is computationally much more efficient and robust to strong noise. A level set model using an optimized fuzzy region clustering technique [28] was also presented for biomedical MRI and CT scan image segmentation. Their proposed algorithm is capable of identifying weak boundaries and can segment the desired components of an image. In [29], the authors proposed an improved fuzzy region competition-based framework via the hierarchical strategy so that the minimization problem is always convex during the iterative calculation. Their method was applied to noisy SAR images, with satisfactory segmentation results.

Although most variational region-based segmentation approaches perform satisfactorily for some natural and textural images, they consider only simple features from a particular image domain as region characterization, and thus, may not perform well for images with sophisticated texture patterns. In such a case, high-level features from various image domains are critical to deal with complex images. In [30], the authors used a geodesic active contour in conjunction with the Gabor features to segment texture images. A vectorial piecewise constant Mumford-Shah model [31] was adopted for the Gabor filtered images, producing a satisfactory performance. Moreover, a level set formulation with a structure tensor and nonlinear diffusion [32] was employed for unsupervised texture segmentation. A new large-scale image segmentation method that combines fuzzy region competition and the Gaussian mixture model [33] was also proposed. Experimental results indicate that this approach has a superior performance on remote sensing images. In addition to the image features mentioned above, the distribution of image variation is one of the most important and useful features which has been widely used in various areas with satisfactory results. Nevertheless, the distribution should not be directly used because the high dimensionality leads to extremely high computational cost for real-time applications. Hence, there are compelling reasons to develop a more powerful representation to replace the distribution while preserving its important properties. Recently, many effective and efficient statistical methods have been proposed to model the distributions. Essentially, the choice of a parametric family of models depends on the applications and the model itself should be mathematically and statistically tractable. For instance, Do and Vetterli [34] proposed the use of a generalized Gaussian density (GGD) to model wavelet subband histograms and successfully applied GGD to texture retrieval. A characteristic GGD based on Kullback-Leibler divergence [35] and the generalized Gamma density [36] were proposed and applied to supervised texture classification with promising results. While these models usually work well for most distributions of image variation, they always assume the distributions have a specific structure (such as symmetry, monotone and periodicity) and cannot model fluctuating distributions. To remedy this shortcoming, the bit-plane-based probability models [19,37,38] were proposed to characterize image variations that do not need to have specific structures. Thus, incorporating these models to characterize image regions into the energy functional would help to enhance the segmentation of texture images and deal with challenging segmentation problem.

Motivation and Contribution

Variational image segmentation methodology is well-established in the literature and fuzzy region competition is by far one of the most popular approaches. Traditional fuzzy region competition models consider simple boundary and region features extracted in a specific image domain, but they are unable to capture essential image information, and thus, hard to segment images with low color contrast and complex textured patterns. In addition, an appropriate and sophisticated probability distribution should be adopted to correctly characterize the image region in order to enhance the segmentation accuracy. Nevertheless, a simple probability distribution (e.g., Gaussian distribution) is typically employed in order to obtain a closed-form solution of the energy functional so that a time-consuming numerical algorithm is avoided to estimate the parameters of the probability distribution. Motivated by the above issues and recent research on statistical image variation modeling, we propose a novel variational approach that can tackle the limitations of existing algorithms. The contributions of this paper are summarized below.

• We propose a novel variational model that integrates the pixel-level color feature and region-level spatial/frequency information into the fuzzy region competition for image segmentation. Specifically, we propose using the windowed bit-plane-dependence probability models to characterize spatial/frequency region information extracted from various image domains to improve the segmentation accuracy;
• A fuzzy bit-plane-dependence region competition algorithm is developed based on the alternating minimization procedure with the Chambolle’s fast duality projection algorithm. In addition, the closed-form solutions for model parameters and fuzzy membership function are obtained and no numerical algorithm is necessary for parameter estimation, thus making the proposed algorithm much faster.

This paper is organized as follows. Section 2 introduces the bit-plane-based probability models for image region characterization. In Section 3, we present the fuzzy bit-plane-dependence region competition model. Section 4 presents the optimization procedure and the overall implementation. Experimental results are shown in Section 5 and Section 6 concludes the paper.

2. Bit-Plane-Dependence Probability Model

The modeling of image variation by a parametric family of statistical distributions plays an important role in many applications. Many studies reveal that using a robust parametric model to represent image variation leads to satisfactory texture classification and retrieval performance. In this section, we present the bit-plane-dependence probability model and use it as a region characterization for our proposed segmentation model as presented in Section 3.

Pi et al. [37] first proposed using the product of Bernoulli distributions (PBD) to model the distributions of image variation. The idea is summarized as follows: given a collection of coefficients of image variation, each coefficient Y is quantized into a nonnegative integer, which is then transformed into the following m binary bit-planes:

$$Y={\displaystyle \sum _{i=0}^{m-1}{2}^i{Y}_i},$$

(1)

where $Y_i\in \{0,1\}$ is a Bernoulli random variable representing the $i\mathrm{th}$ binary bit of the quantized coefficient. Based on (1), the joint probability distribution of the quantized coefficients is

$$P(Y=y)=P(Y_0=y_0,Y_1=y_1,\dots .,Y_{m-1}=y_{m-1}),y=0, 1,\dots , 2^m-1,$$

(2)

where $\left\{y_0,y_1,\dots ,y_{m-1}\right\}$ is a binary representation of y. Let $p_i=P(Y_i=1),$ $i=0,1,\dots ,m-1,$ and assume that $\left\{Y_i:i=0, 1,\dots , m-1\right\}$ are statistically independent, then (2) can be written as PBD, as follows:

$$P^{\mathrm{PBD}}(Y=y)={\displaystyle \prod _{i=0}^{m-1}P(Y_i=y_i)}={\displaystyle \prod _{i=0}^{m-1}{p}_i^{y_i}{(1-p_i)}^{1-y_i}}.$$

(3)

The major advantage of PBD is that it performs excellently in modeling wavelet subbands and provides a satisfactory performance in fitting low-degree fluctuating histograms [38]. Experimental results [37,38] also show that PBD performs better than GGD [34] in supervised texture retrieval. Furthermore, compared with most existing model parameter estimation methods, the bit-plane probabilities $p_i$ can be efficiently extracted via the counting of 1-bit occurrence for each bit-plane.

Note that the PBD model (3) assumes all bit-planes are independent, but such an assumption may not be appropriate because dependencies may exist between bit-planes. This shortcoming can be alleviated by incorporating conditional probabilities between successive bit-planes into (2):

$$P(Y=y)=P(Y_0=y_0){\displaystyle \prod _{i=1}^{m-1}P(Y_i=y_i|Y_{i-1}=y_{i-1})},y=0, 1,\dots , 2^m-1.$$

Let ${\lambda }_i=P(Y_i=1|Y_{i-1}=1),$$i=1,2,\dots ,m-1,$ and assume successive bit-plane dependence. Then, we obtain the following bit-plane-dependence probability model (BDPM):

$$\begin{array}{l}{P}^{\mathrm{BDPM}}(Y=y) =p_0^{y_0}{(1-p_0^{})}^{1-y_0}{\displaystyle \prod _{i=1}^{m-1}{\lambda }_i^{y_{i-1}{y}_i}{(1-{\lambda }_i^{})}^{y_{i-1}(1-y_i)}}\\ \times {\left(\frac{p_i-{\lambda }_i{p}_{i-1}}{1-p_{i-1}} \right)}^{(1-y_{i-1})y_i}{\left(1-\frac{p_i-{\lambda }_i{p}_{i-1}}{1-p_{i-1}} \right)}^{(1-y_{i-1})(1-y_i)}\end{array}$$

(4)

where $\max (0,1-(1-p_i)/p_{i-1})\le {\lambda }_i\le \min (p_i/p_{i-1},1)$ such that the conditional probabilities are between zero and one. Note that (4) reduces to (3) when ${\lambda }_i=p_i$, $i=1,\dots , m-1$. The major advantages of (4) are threefold [19]. Firstly, the maximum likelihood estimators of model parameters are joint sufficient statistics, which capture all possible information about the model parameters that is in the data. Secondly, compared with the current image variation models which assume that the distributions have specific structures such as symmetry, monotone, and periodicity, BDPM provides a universal parametric representation that can be used to model random distributions without enforcing any specific restrictions on the distributions. Thirdly, experiments show that BDPM outperforms PBD in terms of image variation modeling.

3. Fuzzy Bit-Plane-Dependence Region Competition Model

Having introduced the bit-plane-dependence probability model, we present, in this section, the fuzzy bit-plane-dependence region competition model. We start by formulating the segmentation problem as follows: let I be an image with image domain $\mathsf{\Omega }$, the goal of segmentation is to partition $\mathsf{\Omega }$ into N regions $\left\{{\mathsf{\Omega }}_j:j=1,2,\dots ,N\right\}$ by a suitable measure such that ${\mathsf{\Omega }}_i\cap {\mathsf{\Omega }}_j=\varnothing$, $i\ne j$, and $\mathsf{\Omega }={\mathsf{\Omega }}_i\cup \cdots \cup {\mathsf{\Omega }}_N\cup {\Gamma }$ with ${\Gamma }$ representing the boundaries of all regions.

Various successful fuzzy region competition approaches in the literature are generally based on the optimization of an energy functional consisting of both the boundary and region information. The general form of energy functional for segmenting I into N phases can be represented as:

$$E\left(\delta ,U\right)={\displaystyle \sum _{j=1}^N{\displaystyle {\int }_{\mathsf{\Omega }}\left|\nabla u_j\right|dx}-\alpha {\displaystyle \sum _{j=1}^N{\displaystyle {\int }_{\mathsf{\Omega }}{u}_j\log P_j(I|{\delta }_j)dx}}}$$

(5)

with constraints

$$\left(\mathrm{i}\right) {\displaystyle \sum _{j=1}^N{u}_j}=1, \left(\mathrm{ii}\right) 0\le u_j\le 1,j=1,2,\dots ,N,$$

(6)

where $\delta =\left\{{\delta }_j:j=1,2,\dots ,N\right\}$ is a set of parameters for the regions $\left\{{\mathsf{\Omega }}_j:j=1,2,\dots ,N\right\}$, $U=\left\{u_j:u_j\in B{V}_{[0,1]}(\mathsf{\Omega }), j=1,2,\dots ,N\right\}$ with $B{V}_{[0,1]}(\mathsf{\Omega })$ is a space of the functions of the bounded variations taking their values between zero and one, $P_j(I|{\delta }_j)$ is a probability distribution of region ${\mathsf{\Omega }}_j$ characterized by the parameter ${\delta }_j$, and $\alpha$ is a positive constant that balances the two energy terms. In (5), the first term is a total variation that measures the lengths of boundaries for all of the image regions in a fuzzy manner whereas the second term is the sum of the cost for coding the intensity of every pixel inside each region according to a probability distribution $P_j(I|{\delta }_j)$.

Most existing fuzzy region competition models consider only simple image features extracted from a particular image domain, but these features may not capture essential region information, especially for images with low color contrast and complex textured structures. Hence, it is expected that the overall segmentation results would improve if the effective texture features extracted from different image domains could be adopted in conjunction with color for region characterization. As mentioned in Section 2, we propose using the bit-plane-dependence probability model to characterize local region information in various image domains and incorporating (4) into the energy functional (5).

The Proposed Model

Let $I^0$ be a color image and $I^1$ be the associated gray-scale image. We assume that the spatial image components of $I^1$ in a particular region ${\mathsf{\Omega }}_j$ are independent and identically distributed (i.i.d.) samples generated by a family of probability distribution $P_j^1\left(\left.I^1\right|{\varphi }_j^1\right)$, where ${\varphi }_j^1$ is the parameter of the distribution $P_j^1$. To utilize the region information extracted from other image domains, we let $\{{\mathsf{\Psi }}_q:q=2,3,\dots ,Q\}$ be a collection of image transform operators which transform $I^1$ into a series of transformed images $\{I^q:q=2,3,\dots ,Q\}$ in the associated image domains, i.e., $I^q={\mathsf{\Psi }}_q(I^1)$. In the same region ${\mathsf{\Omega }}_j$, we further assume that the region components of the respective image domain are i.i.d. samples generated by a family of probability distribution $P_j^q\left(\left.I^q\right|{\varphi }_j^q\right)$, where ${\varphi }_j^q$ is the parameter of the distribution $P_j^q$. To integrate the model (4) into the energy functional (5), we use BDPM for all probability distributions $\{P_j^q:j=1,2,\dots ,N;q=1,2,\dots ,Q\}$. That is,

$$\begin{array}{l}{P}_j^q(I^q|{\varphi }_j^q)={\left(p_{0j}^q\right)}^{y_0^q}{(1-p_{0j}^q)}^{1-y_0^q}{\displaystyle \prod _{i=1}^{m-1}{\left({\lambda }_{ij}^q\right)}^{y_{i-1}^q{y}_i^q}{(1-{\lambda }_{ij}^q)}^{y_{i-1}^q(1-y_i^q)}}\\ \times {\left(\frac{p_{ij}^q-{\lambda }_{ij}^q{p}_{i-1j}^q}{1-p_{i-1j}^q} \right)}^{(1-y_{i-1}^q)y_i^q}{\left(1-\frac{p_{ij}^q-{\lambda }_{ij}^q{p}_{i-1j}^q}{1-p_{i-1j}^q} \right)}^{(1-y_{i-1}^q)(1-y_i^q)}\end{array}$$

(7)

where ${\varphi }_j^q=\left\{\left(p_{0j}^q,p_{ij}^q,{\lambda }_{ij}^q\right):i=1,2,\dots ,m-1\right\}$. In addition, the CIELAB color feature with the L₂-norm as a dissimilarity measurement is incorporated into (5) to capture pixel-level information. The CIELAB color space is adopted because it is perceptually uniform with respect to human color vision and performs well in various applications. Then, we obtain the following fuzzy bit-plane-dependence region competition energy functional:

$$E\left(\varphi ,C,U\right)={\displaystyle \sum _{j=1}^N{\displaystyle {\int }_{\mathsf{\Omega }}|\nabla u_j|dx}}+{\alpha }_0{\displaystyle \sum _{j=1}^N{\displaystyle {\int }_{\mathsf{\Omega }}{u}_j{\left(I^0-c_j\right)}^{2}dx}}-{\displaystyle \sum _{q=1}^Q{\alpha }_q{\displaystyle \sum _{j=1}^N{\displaystyle {\int }_{\mathsf{\Omega }}{u}_j\log P_j^q\left(\left.I^q\right|{\varphi }_j^q\right)dx}}}$$

(8)

subject to the constraints (6), where $\varphi =\left\{{\varphi }_j^q:j=1,2,\dots ,N;q=1,2,\dots Q\right\}$ is a set of parameters for the probability distributions $\left\{P_j^q:j=1,2,\dots ,N;q=1,2,\dots ,Q\right\}$, $C=\left\{c_j:j=1,2,\dots ,N\right\}$ is a set of CIELAB color parameters, $U=\left\{u_j: j=1,2,\dots ,N\right\}$ is a fuzzy membership function, and $\{{\alpha }_q:q=0,1,2,\dots ,Q\}$ are fixed parameters.

Note that the major drawback of (8) is that the misclassification of an image pixel would occur because of the noise or statistical fluctuation of data in various image domains, as only a single sample is taken from each probability distribution. To overcome this shortcoming, we propose incorporating the neighboring information around each pixel into (8). Let $W_x^M$ be an $M\times M$ window centered at pixel x. For each image domain, the probability distribution $P_j^q\left(\left.I^q\right|{\varphi }_j^q\right)$ is replaced by ${\left({\displaystyle {\prod }_{W_x^M}^{}{P}_j^q\left(\left.I^q\right|{\varphi }_j^q\right)}\right)}^{1/M^2}$, which is the geometric mean of BDPM over $W_x^M$. In fact, $\log P_j^q\left(\left.I^q\right|{\varphi }_j^q\right)$ is replaced by $\left(1/M^2\right){\displaystyle {\sum }_{W_x^M}\log P_j^q\left(\left.I^q\right|{\varphi }_j^q\right)}$, which is the average of log-likelihood functions over $W_x^M$. Hence, we can rewrite the energy functional (8), as follows:

$$\begin{array}{l}E\left(\varphi ,C,U\right)={\displaystyle \sum _{j=1}^N{\displaystyle {\int }_{\mathsf{\Omega }}|\nabla u_j|dx}}+{\alpha }_0{\displaystyle \sum _{j=1}^N{\displaystyle {\int }_{\mathsf{\Omega }}{u}_j{\left(I^0-c_j\right)}^{2}dx}}\\ -{\displaystyle \sum _{q=1}^Q{\alpha }_q{\displaystyle \sum _{j=1}^N{\displaystyle {\int }_{\mathsf{\Omega }}{u}_j\left(\log {\left({\displaystyle {\prod }_{W_x^M}^{}{P}_j^q\left(\left.I^q\right|{\varphi }_j^q\right)}\right)}^{1/M^2}\right)dx}}}\end{array}$$

(9)

We remark that the CIELAB color characterizes pixel-level information whereas the windowed BDPM model captures region-level information. In addition, the energy functional (9) not only provides local region information in the spatial domain, as in the original and existing region competition models, but also provides local region information extracted from other image domains (i.e., when $q=2,3,\dots ,Q$), which is expected to improve the overall segmentation performance.

4. The Optimization Procedure

The minimization of (9) with the BDPM model (7) and the constraints of (6) form a class of constrained nonlinear optimization problems whose solutions are unknown. Various methodologies can be used to minimize the energy functional, such as curve evolution via the level-set formulation or alternating minimization method. In this paper, we shall use the alternating minimization procedure and make use of the fast duality projection algorithm of Chambolle [39]. Note that the data fidelity and regularization terms in the energy functional are coupled, thus we introduce a collection of auxiliary variables $V=\left\{v_j:j=1,2,\dots ,N\right\}$ to decouple these two terms. Specifically, we shall approximate (9) by replacing $u_j$ with $v_j$ in the regularization term and adding convex terms that force $u_j$ and $v_j$ to be sufficiently close:

$$\begin{array}{l}E\left(\varphi ,C,U,V\right)={\displaystyle \sum _{j=1}^N{\displaystyle {\int }_{\mathsf{\Omega }}|\nabla v_j|dx}}+\frac{1}{2\theta }{\displaystyle \sum _{j=1}^N{\displaystyle {\int }_{\mathsf{\Omega }}{\left(v_j-u_j\right)}^{2}dx}}\\ +{\alpha }_0{\displaystyle \sum _{j=1}^N{\displaystyle {\int }_{\mathsf{\Omega }}{u}_j{\left(I^0-c_j\right)}^{2}dx}}-{\displaystyle \sum _{q=1}^Q{\alpha }_q{\displaystyle \sum _{j=1}^N{\displaystyle {\int }_{\mathsf{\Omega }}{u}_j\left(\log {\left({\displaystyle {\prod }_{W_x^M}^{}{P}_j^q\left(\left.I^q\right|{\varphi }_j^q\right)}\right)}^{1/M^2}\right)dx}}}\end{array}$$

(10)

where $\theta$ is a small positive constant such that $u_j$ is sufficiently close to $v_j$ with respect to the L₂-norm. In what follows, the alternating minimization procedure is adopted to perform the optimization of (10).

4.1. CIELAB Color Parameter

Fixing $\varphi ,$ U and $V$, we compute $C=\left\{c_j:j=1,2,\dots ,N\right\}$ by setting

$$\frac{\partial E_{\varphi ,U,V}\left(C\right)}{\partial c_j}={\displaystyle {\int }_{\mathsf{\Omega }}{u}_j\left(I^0-c_j\right)dx}=0$$

It is straight forward to demonstrate that

$${\widehat{c}}_j=\frac{{\displaystyle {\int }_{\mathsf{\Omega }}{u}_j{I}^{0}dx}}{{\displaystyle {\int }_{\mathsf{\Omega }}{u}_{j}dx}},j=1,2,\dots ,N,$$

(11)

which is simply the CIELAB color weighted by the fuzzy membership function.

4.2. Windowed BDPM Parameter

Fixing $C,U,V$, we solve $\varphi =\left\{{\varphi }_j^q:j=1,2,\dots ,N;q=1,2,\dots Q\right\}$ with ${\varphi }_j^q=\left\{\left(p_{0j}^q,p_{ij}^q,{\lambda }_{ij}^q\right):i=1,2,\dots ,m-1\right\}$ by minimizing

$$E_{C,U,V}\left(\varphi \right)={\displaystyle \sum _{q=1}^Q{\alpha }_q{\displaystyle \sum _{j=1}^N{\displaystyle {\int }_{\mathsf{\Omega }}{u}_j\left(\log {\left({\displaystyle {\prod }_{W_x^M}^{}{P}_j^q\left(\left.I^q\right|{\varphi }_j^q\right)}\right)}^{1/M^2}\right)dx}}}.$$

(12)

Taking the partial derivatives of (12) with respect to $p_{ij}^q$ and ${\lambda }_{ij}^q$ followed by setting them to zero leads to

$${\widehat{p}}_{i,j}^q=\frac{{\displaystyle {\int }_{\mathsf{\Omega }}{u}_j\left(\frac{1}{M^2}{\displaystyle {\int }_{W_x^M}{y}_i^{q}dz}\right)dx}}{{\displaystyle {\int }_{\mathsf{\Omega }}{u}_{j}dx}},i=0,1,\dots ,m-1;j=1,2,\dots ,N;q=1,2,\dots ,Q,$$

(13)

$${\widehat{\lambda }}_{i,j}^q=\frac{{\displaystyle {\int }_{\mathsf{\Omega }}{u}_j\left({\displaystyle {\int }_{W_x^M}{y}_{i-1}^q{y}_i^{q}dz}\right)dx}}{{\displaystyle {\int }_{\mathsf{\Omega }}{u}_j\left({\displaystyle {\int }_{W_x^M}{y}_{i-1}^{q}dz}\right)dx}},i=1,2,\dots ,m-1;j=1,2,\dots ,N;q=1,2,\dots ,Q.$$

(14)

In (13) and (14), the estimators are the weighted local sample mean of 1-bit occurrence, and the weighted local average of 1-bit occurrence in successive bit-planes, respectively, within a window for each region in a particular image domain. It is important to remark that the windowed BDPM parameter has a closed-form solution, which implies that it can be obtained efficiently via only the bit-plane extraction. Note that when M = 1, the neighboring information is neglected and only pixel-level bit-plane information is captured, whereas the region-level bit-plane information is extracted through the windowed BDPM when M > 1.

4.3. Total Variation Minimization

Fixing $\varphi ,$ C and U, we compute V by minimizing

$$E_{\varphi ,C,U}\left(V\right)={\displaystyle {\int }_{\mathsf{\Omega }}|\nabla v_j|dx}+\frac{1}{2\theta }{\displaystyle {\int }_{\mathsf{\Omega }}{\left(v_j-u_j\right)}^{2}dx}$$

which can be solved efficiently using the Chambolle’s fast duality projection algorithm [39]. Then, the solution is given by

$${\widehat{v}}_j=u_j-\theta \mathrm{div} k_j,$$

(15)

where div is the divergence operator. The vector $k_j$ can be computed using a fixed-point algorithm:

$$k_j^{n+1}=\frac{k_j^n+\tau \nabla (\mathrm{div} k_j^n-u_j/\theta )}{1+\tau \left|\nabla (\mathrm{div} k_j^n-u_j/\theta )\right|}$$

with initial value $k_j^0=0$ and $0<\tau \le 1/8$ to ensure the convergence of the algorithm [39].

4.4. Fuzzy Membership Function

Fixing $\varphi ,$ C and V, we solve U by minimizing

$$\begin{array}{l}{E}_{\varphi ,C,V}\left(U\right)=\frac{1}{2\theta }{\displaystyle \sum _{j=1}^N{\displaystyle {\int }_{\mathsf{\Omega }}{\left(v_j-u_j\right)}^{2}dx}}+{\alpha }_0{\displaystyle \sum _{j=1}^N{\displaystyle {\int }_{\mathsf{\Omega }}{u}_j{\left(I^0-c_j\right)}^{2}dx}}\\ -{\displaystyle \sum _{q=1}^Q{\alpha }_q{\displaystyle \sum _{j=1}^N{\displaystyle {\int }_{\mathsf{\Omega }}{u}_j\left(\log {\left({\displaystyle {\prod }_{W_x^M}^{}{P}_j^q\left(\left.I^q\right|{\varphi }_j^q\right)}\right)}^{1/M^2}\right)dx}}}\end{array}$$

subject to

$$(\mathrm{i}) {\displaystyle \sum _{j=1}^N{u}_j}=1, \left(\mathrm{ii}\right) 0\le u_j\le 1,j=1,2,\dots ,N.$$

Equivalently, we can minimize the energy

$${\tilde{E}}_{\varphi ,C,V}^{}\left(U\right)=\frac{1}{2\theta }{\displaystyle {\int }_{\mathsf{\Omega }}{\displaystyle \sum _{j=1}^N{\left(H_j-u_j\right)}^2}dx}+K,$$

subject to the same constraints (i) and (ii), where

$$H_j=v_j-{\alpha }_0\theta {\left(I^0-c_j\right)}^2+\theta {\displaystyle \sum _{q=1}^Q{\alpha }_q\left(\log {\left({\displaystyle {\prod }_{W_x^M}^{}{P}_j^q\left(\left.I^q\right|{\varphi }_j^q\right)}\right)}^{1/M^2}\right)}$$

and

$$K=-\frac{1}{2\theta }{\displaystyle \sum _{j=1}^N{\displaystyle {\int }_{\mathsf{\Omega }}{H}_j{}^{2}dx}}$$

Note that both constraints (i) and (ii) apply to membership function independently for each point $x\in \mathsf{\Omega }$. Hence the minimizer $\widehat{U}$ of ${\tilde{E}}_{\varphi ,C,V}^{}$ is also the pointwise minimizer of the function $f\left(U\right)={\displaystyle {\sum }_{j=1}^N{\left(H_j-u_j\right)}^2}$ subject to the same constraints. That is, for each point $x\in \mathsf{\Omega }$, $\widehat{U}(x)$ is the solution to the problem

$$\underset{U}{\min }{\displaystyle \sum _{j=1}^N{\left(H_j(x)-u_j(x)\right)}^2},$$

subject to

$$(\mathrm{i}) {\displaystyle \sum _{j=1}^N{u}_j(x)}=1, \left(\mathrm{ii}\right) 0\le u_j(x)\le 1,j=1,2,\dots ,N.$$

The above minimization problem is exactly the problem of computing the Euclidean projection of the vector $H(x)=[H_1(x),\cdots ,H_N(x)]$ on the probability simplex ${\Delta }^N$. According to [40], the projection z of a vector $y\in R^N$ onto ${\Delta }^N$ can be expressed as $z_i=\max \{y_i+\lambda ,0\},$ where $i=1,2,\dots ,N$, and $\lambda$ is the Lagrange multiplier chosen such that the constraint ${\displaystyle {\sum }_{i=1}^N{z}_i}=1$ holds. An algorithm to compute $\lambda$, with $y=H(x)$, consists of the following steps. Firstly, we reorder the components $\left\{H_j(x):1\le j\le N\right\}$ into $H_{(1)}(x)\ge H_{(2)}(x)\ge \cdots \ge H_{(N)}(x)$, and find the integer $\rho (x)$ such that

$$\rho (x)=\max \left\{j:H_j(x)+\frac{1}{j}\left(1-{\displaystyle \sum _{i=1}^j{H}_i(x)}\right)>0\right\}$$

Then, the Lagrange multiplier $\lambda (x)$ can be computed as

$$\lambda (x)=\frac{1}{\rho (x)}\left(1-{\displaystyle \sum _{i=1}^{\rho (x)}{H}_{(i)}(x)}\right)$$

Finally, the minimizer $\widehat{U}$ of ${\tilde{E}}_{\varphi ,C,V}^{}$, which is also the projection of $H(x)$ on ${\Delta }^N$ for each $x\in \mathsf{\Omega }$, is given by

$${\widehat{u}}_j=\max \left\{H_j+\lambda ,0\right\}, j=1,2,\dots ,N.$$

(16)

4.5. The Overall Implementation

The optimization of (9) is performed by the alternating minimization between $\varphi$, C and U. The fuzzy bit-plane-dependence region competition algorithm (Algorithm 1) is summarized below. The overall implementation of the proposed method is shown in Figure 1.

Algorithm 1 Fuzzy Bit-plane-dependence Region Competition Algorithm
	Input: Input image I, number of regions N, number of bit-planes m, number of image domains Q, regularization parameters ${\alpha }_q$, $q=0,1,\dots ,Q,$ size of window M, parameters of Chambolle’s fast duality projection algorithm $\theta$ and $\tau$, initial fuzzy membership function U 0 and parameter of termination criterion $\epsilon$.
	Output: Optimal fuzzy membership function ${\widehat{U}}^{opt}$
Step 1:	Compute the color and windowed BDPM parameters (C and $\varphi$) using (11), (13) and (14), respectively.
Step 2:	Compute the auxiliary variable and fuzzy membership function (V and U) using (15) and (16), respectively.
Step 3:	Repeat Step 1 and Step 2 till termination. The termination criterion is ${\left|U^{new}-U^{old}\right|}_{\infty }<\epsilon$, where $\epsilon$ is a small positive constant.

In order to visualize the segmentation results and compare our method with the existing approaches, as reported in the following section, the defuzzification process is performed as follows. Given the optimal fuzzy membership function, ${\widehat{U}}^{opt}=\left\{{\widehat{u}}_j^{opt}:j=1,2,\dots ,N\right\}$, obtained from the algorithm above, the final segmentation result is constructed by assigning each pixel a label $j^{*}$ where $j^{*}=\arg \underset{j}{\max }{\widehat{u}}_j^{opt}$.

5. Experiments

5.1. Experimental Setting

The proposed method is applied to real-world natural images. Unless otherwise specified, we fix the parameters, which achieve the best results, in all the following experiments: $m=8$, $\theta =1.0$, $\tau =0.125$, $M=25$, ${\alpha }_0\in [0.0001,0.005]$, ${\alpha }_q=50{\alpha }_0$, $q=1,2,\dots ,Q,$ and $\epsilon ={10}^{-2}$. To characterize local region information in various image domains using the windowed BDPM model, a collection of transformed images $\{I^q:q=2,3,\dots ,Q\}$ should be generated. In this paper, we use 3-level wavelet transform (i.e., Q = 10) followed by the Hilbert transform to construct the transformed images. We simply report results of 3-level wavelet transform since our segmentation results reveal that 3-level decomposition is good enough to capture essential wavelet information and usually performs better than the others. Note that when 4-level decomposition is used, the subband information at the highest level is redundant and does not enhance the segmentation performance. To obtain the transformed images, we first index each wavelet subband $q=\{2,3,\dots ,10\}=\{(l,\theta ):l=1,2,3;\theta =(H,V,D)\}$, where l denotes the level of wavelet decomposition and $\theta$ is the direction with H, V and D representing horizontal, vertical, and diagonal directions, respectively. Evaluation of all possible wavelet transforms, and filter banks is beyond the scope of this paper; thus, we focus on undecimated wavelet transform (UWT) with Symlets filter bank. Dissimilar to the traditional discrete wavelet transform in which downsampling is performed in each decomposition level, UWT is used so that the filtered image is not subsampled, and therefore, each subband is the same size as the original image. It is important to remark that the UWT coefficients may exhibit a wide spectrum of sinusoidal oscillation, which affects the image variation modeling performance, leading to the degradation of segmentation accuracy. Thus, we narrow down the spectrum of sinusoidal oscillation by applying the Hilbert transform [41] to the UWT coefficients and consider their magnitude to obtain the transformed image. Mathematically, let $C_{(l,b)}$ be the UWT coefficients of the corresponding $q\mathrm{th}$ subband and $H(C_{(l,b)})$ be the Hilbert transform of $C_{(l,b)}$. Then the transformed image is given by Figure 2

$$I^q=255\times \frac{|H(C_{(l,b)})|-\min |H(C_{(l,b)})|}{\max |H(C_{(l,b)})|-\min |H(C_{(l,b)})|}.$$

To quantify the effectiveness of the proposed method, we use three measures, namely: segmentation covering [42], F-measure [43] and Jaccard index [44]. The segmentation covering (SC) measures the overlap of regions of the final segmentation and the regions of the ground truth, F-measure (F) is defined as a harmonic mean of precision and recall, while Jaccard index (JI) measures the intersection over the union of the labelled segments.

5.2. Comparative Segmentation Results

We compare the proposed fuzzy bit-plane-dependence region competition (FBRC) algorithm with several state-of-the-art variational-based, superpixel-based and deep-learning-based approaches, namely: fuzzy region competition with Gaussian mixture model (FRCGMM) [33], L1 fuzzy segmentation (L1FS) [25], piecewise constant fuzzy region competition (PCFRC) [2], superpixel-based fast fuzzy c-means clustering (SFFCM) [45], automatic fuzzy clustering framework (AFCF) [46], backpropagation (BP) [47] and differentiable feature clustering (DFC) [48] algorithms using the well-known Microsoft research in Cambridge dataset [49], Berkeley segmentation dataset 500 [50], Weizmann segmentation evaluation (WSE) database [51] and PH2 database [52]. In the implementation of the above algorithms, we have tried different combinations of parameters to obtain the best possible performance.

Figure 2 displays the comparative segmentation results. In this experiment, we shall simply assess their performance based on visual satisfaction. The following are some discussions about the algorithms and main points observed from the figures.

Figure 2. Comparative segmentation results for various algorithms. (a–j) are images from [49,50,51,52]. The 1st to the 8th columns are the results for FRCGMM, L1FS, PCFRC, SFFCM, AFCF, BP, DFC and FBRC, respectively.

Figure 2. Comparative segmentation results for various algorithms. (a–j) are images from [49,50,51,52]. The 1st to the 8th columns are the results for FRCGMM, L1FS, PCFRC, SFFCM, AFCF, BP, DFC and FBRC, respectively.

(1)

The segmentation results of the proposed method for all images are visually satisfying whereas the rest gives unsatisfactory performance for some images. Specifically, the variational-based FRCGMM, L1FS and PCFRC and superpixel-based SFFCM and AFCF perform similarly while the deep-learning-based BP and DFC provide undesirable results. The success of FBRC is mainly due to the fact that it adopts the CIELAB color feature to capture pixel-level information and the windowed BDPM model to characterize region-level information while the rest employs only pixel-level information in a particular image domain. In addition, the proposed method incorporates additional information from various image domains to further enhance the segmentation accuracy;

(2)

The proposed method usually adheres well to object boundaries for all images, but other approaches give approximated object boundaries for some cases, see, for instance, the segmented boundaries of Figure 2d by FRCGMM, Figure 2b by SFFCM, and Figure 2h by DFC. Essentially, the unsatisfactory boundary adherence of the cases reveals the insufficiency of feature representation adopted by these methods. Note that the proposed method uses the windowed BDPM to characterize local region information and the window may straddle multiple objects in the image, and thus, the quality of the neighboring information in the window will degrade. Nevertheless, the use of windowed BDPM to capture region information from multi-domain provides additional information so that our proposed method can still capture the object boundaries well;

(3)

The proposed method performs better than the other algorithms in terms of mis-segmentation of object parts. For instance, L1FS mis-segments the “Sky” in Figure 2d, AFCF mis-segments the “Rake” in Figure 2e, and BP mis-segments the “Cow” in Figure 2c. Similar examples can also be found in other images. The results clearly indicate that the use of windowed BDPM model to characterize texture patterns is shown to be effective in most cases, see, for instance, the “Tree” in Figure 2b.

Figure 3 shows the comparative segmentation results with the provided ground truths, which are obtained by asking human subjects to manually segment the images into various object parts for all of the algorithms. As can be seen, the segmentation results of Figure 3 for FBRC are visually satisfying and the proposed method outperforms other algorithms for all images. As previously mentioned, this is mainly because most algorithms use a single feature in a particular image domain to characterize image regions, but our method adopts both the color and multi-domain windowed BDPM model to capture essential pixel and region levels information. In addition, we also remark that the proposed method adopts BDPM to capture texture information, and thus, our algorithm can segment images with low color contrast and complex textured patterns.

Table 1. Quantitative comparisons of various approaches for WSE database. The best results are highlighted in bold.

Methods	SC	JI	F	Time (s)
FRCGMM	0.799	0.666	0.765	45.282
L1FS	0.758	0.644	0.748	72.296
PCFRC	0.821	0.687	0.773	39.619
SFFCM	0.798	0.671	0.766	0.815
AFCF	0.801	0.657	0.751	0.673
BP	0.726	0.408	0.496	97.358
DFC	0.640	0.420	0.528	79.903
FBRC	0.829	0.701	0.788	31.473

reports the average quantitative measures for all the algorithms using the WSE database. As far as the segmentation accuracy is concerned, we observe that the variational-based and superpixel-based algorithms perform similarly, and FBRC outperforms the recently developed deep-learning-based algorithms. Lastly, it is important to emphasize that the proposed method generally performs better than all other algorithms in the sense that it achieves the best measures in SC, JI and F.

To provide an additional justification of our approach, we shall compare the proposed method with the aforementioned approaches using the PH2 database, which consists of 200 dermoscopic images of melanocytic lesions.

Table 2. Quantitative comparisons of various approaches for PH2 database. The best results are highlighted in bold.

Methods	SC	JI	F	Time (s)
FRCGMM	0.765	0.674	0.778	35.431
L1FS	0.814	0.712	0.800	35.912
PCFRC	0.807	0.665	0.753	26.769
SFFCM	0.723	0.659	0.773	4.561
AFCF	0.690	0.591	0.711	4.609
BP	0.805	0.471	0.543	285.990
DFC	0.626	0.285	0.373	195.900
FBRC	0.821	0.739	0.827	25.783

shows the corresponding quantitative results. As with the results of the WSE database, both variational-based and superpixel-based methods have similar performance, and FBRC performs better than the deep-learning-based algorithms. Compared with the current state-of-the-art methodologies, the results clearly reveal the superior performance of the proposed method.

5.3. Parameter Sensitivity

In this subsection, we analyze the sensitivity of parameters to the segmentation performance. In the Chambolle’s fast duality projection algorithm, we set $\theta =1.0$ and $\tau =0.125$. Our experiments show that changing the values of $\theta$ and $\tau$ will not affect the segmentation results unless the values of $\theta$ and $\tau$ are sufficiently large. In fact, the algorithm of Chambolle may not converge if $\tau >0.125$ [39]. To update the windowed BDPM parameter during the optimization process, the window of size M = 25 is used. Here, we study the effect of window size to the segmentation accuracy by comparing the performance for different values of M. Figure 4 displays some visual examples when M = 1, 25, 50 and 75, and

Table 3. Quantitative comparisons using different window sizes for WSE database. The best results are highlighted in bold.

M	SC	JI	F
1	0.829	0.692	0.779
25	0.829	0.701	0.788
50	0.828	0.689	0.777
75	0.827	0.688	0.776

reports the average quantitative measures using different window sizes for the WSE database. While the segmentation results are generally satisfactory, we notice that the performance is inferior when the window size is sufficiently small/large, see, for instance, Figure 4a when M = 1 and Figure 4b when M = 75. This phenomenon can be understood as follows: the larger the window size, the smaller the segmentation error because larger window provides more statistical information for image variation modeling. Nevertheless, when the window is sufficiently large, the segmentation error increases since the window may straddle multiple image regions, and thus, the information extracted from the windows centered at the pixels near the region boundaries could be inaccurate. On the other hand, when the window size is small, the quality of image variation modeling will degrade which would increase the segmentation error. When M = 1, the neighboring information is neglected, and the promising segmentation results may not be achieved since only a single sample is taken from the probability distribution. As evident in our experiments, setting the window of size M = 25 is suffice for achieving satisfactory segmentation results. It is important to remark that in the case when the neighboring information is neglected (then (9) reduces to (8)), the use of the BDPM model for region characterization from various image domains can provide additional information to discriminate regions with similar spatial statistics, which in turn improves the segmentation performance. Lastly, we shall investigate the effectiveness of the BDPM model in the proposed energy functional by comparing the segmentation performance between ${\alpha }_q=0$ and ${\alpha }_q>0$ for $q\ge 1$. When ${\alpha }_q=0$, the BDPM model in the energy functional is neglected while it is adopted when ${\alpha }_q>0$. Figure 5 displays some visual examples. As can be seen, the results for ${\alpha }_q>0$ performs better than that for ${\alpha }_q=0$, see, for instance, the “neck of goat” in Figure 5a, the “tail of the plane” in Figure 5c and the “hair of lady” in Figure 5d.

Table 4. Quantitative comparisons for the effectiveness of BDPM to the segmentation performance for WSE database. The best results are highlighted in bold.

	SC	JI	F
${\alpha }_q=0$$, q\ge 1$	0.818	0.695	0.780
${\alpha }_q>0$$, q\ge 1$.	0.829	0.701	0.788

reports the average quantitative performance for the effectiveness of BDPM to the segmentation performance for the WSE database. We observe that the three measures, namely, SC, JI and F, for ${\alpha }_q>0$ performs better than the cases when ${\alpha }_q=0$. These results qualitatively and quantitatively justify that the use of BDPM model for region characterization from multi-domain can improve the segmentation performance.

5.4. Computational Cost

All of the experiments have been implemented on the Intel Core i7-6700HQ laptop. Table 1 and Table 2 list the average computation times of all methods for WSE and PH2 databases, respectively. Briefly speaking, the superpixel-based methods are the fastest since they simply perform standard fuzzy c-means clustering on the superpixels while the deep-learning-based algorithms have high computational cost because they require optimizing complex convolutional neural networks. FRCGMM and LIFS are slower than PCFRC. This may be explained by the fact that FRCGMM performs highly precise fitting of data with Gaussian mixture model, and that L1FS introduces two sets of auxiliary variables and requires solving a collection of sub-problems based on the alternating direction methods of multipliers. The computation time for FBRC is slightly shorter than PCFRC. Both methods adopt the alternating minimization procedure and apply the Chambolle’s fast duality projection algorithm to solve the total variation minimization problem to speed up the segmentation process.

6. Conclusions

In this paper, we have presented a novel multiphase image segmentation model based on fuzzy region competition and statistical image variation modeling. In the proposed energy functional, each region is characterized by the color and spatial/frequency information modeled by the windowed bit-plane-dependence probability models in various image domains, and is represented by the fuzzy membership function. We have employed the alternating minimization procedure in conjunction with the Chambolle’s fast duality projection algorithm to minimize the energy functional, where its closed-form solutions are obtained. Comparative experiments have demonstrated the effectiveness of our proposed method.

While the proposed method provides satisfactory segmentation performance, we shall study the following as future work. Firstly, a few parameters in the algorithm need to be manually selected in order to achieve promising results. Nevertheless, our experiments reveal that selecting a value for each parameter in a certain range is good enough to provide similar performance. In fact, some parameters have significant impact to the segmentation accuracy. Specifically, an automatic parameter selection scheme for ${\alpha }_q$, $q\ge 0$ is necessary so that these parameters can be adaptively adjusted during the optimization procedure, which may improve the segmentation results. Secondly, we use the wavelet and Hilbert transforms to construct the transformed images in all our experiments. However, various image domains are available in the literature and can be adopted in our case. Thus, the selection of image domains and their advantages for image segmentation (or in a specific application domains) should be investigated.