Spatially Adaptive Regularization in Image Segmentation

Abstract

We present a total-variation-regularized image segmentation model that uses local regularization parameters to take into account spatial image information. We propose some techniques for defining those parameters, based on the cartoon-texture decomposition of the given image, on the mean and median filters, and on a thresholding technique, with the aim of preventing excessive regularization in piecewise-constant or smooth regions and preserving spatial features in nonsmooth regions. Our model is obtained by modifying a well-known image segmentation model that was developed by T. Chan, S. Esedoḡlu, and M. Nikolova. We solve the modified model by an alternating minimization method using split Bregman iterations. Numerical experiments show the effectiveness of our approach.

1. Introduction

Image segmentation is a fundamental problem in image processing and computer vision. Its goal is to divide the given image into regions that represent different objects in the scene, which can be identified taking into account different features. The research on image segmentation has made several advances in the last decades and various approaches have been developed to satisfy the requirements coming from different applications. In this work, we focus on variational models for image segmentation, which have been widely investigated, proving to be very effective on different images—see, e.g., [1,2,3,4]. Other recent approaches to image segmentation include learning-based methods, which often exploit deep-learning techniques [5,6,7,8]. However, in this case, a large amount of data must be available to train learning networks, thus making those approaches impractical in some  applications.

Roughly speaking, when a variational approach is used, the segmentation may be obtained by minimizing a cost functional that linearly combines regularization and data fidelity terms:

$$\underset{u\in \Omega }{min}E\left(u\right)\equiv F\left(u\right)+\lambda G\left(u\right),$$

(1)

where $u:\Omega \subset {\mathbb{R}}^2\to \mathbb{R}$ provides a representation of the segmentation and $\lambda >0$ is a parameter that controls the weight of the fidelity term G versus the regularization term F. A widely-used segmentation model is the two-phase partitioning model that was introduced by Chan, Esedoḡlu, and Nikolova [3], which we refer to as the CEN model:

$$\begin{array}{cc}{\displaystyle \underset{u,c_1,c_2}{minimize}}\hfill & {\displaystyle {\int }_{\Omega }|\nabla u|dx+\lambda {\int }_{\Omega }\left({(c_1-\overline{u}\left(x\right))}^{2}u\left(x\right)+{(c_2-\overline{u}\left(x\right))}^2(1-u\left(x\right))\right)dx,}\hfill \\ \mathrm{s}.\mathrm{t}.\hfill & 0\le u\le 1,\hfill \\ \hfill & c_1,c_2>0.\hfill \end{array}$$

(2)

Here, $\overline{u}$ denotes the image to be segmented, which is assumed to take its values in $[0,1]$. The CEN model allows us to obtain one of the two domains defining the segmentation, $\Sigma$ and $\Omega \backslash \Sigma$, by setting

$$\Sigma =\left\{x:u\left(x\right)>\alpha \right\}\mathrm{for} \mathrm{a}.\mathrm{e}.\alpha \in (0,1),$$

where u is the solution of problem (2). Note that (2) is the result of a suitable relaxation of the Chan–Vese model [2] leading to a convex formulation of that problem for any given $(c_1,c_2)$.

In this work, we start from a discrete version of the CEN model. Let

$$\Theta =\left\{(i,j):0\le i\le m-1,0\le j\le n-1\right\}$$

be a discretization of $\Omega$ consisting of an $m\times n$ grid of pixels and let

$$|{\nabla }_x{u|}_{i,j}=|{\delta }_x^{+}{u|}_{i,j},|{\nabla }_y{u|}_{i,j}={\left|{\delta }_y^{+}u\right|}_{i,j}$$

where ${\delta }_x^{+}$ and ${\delta }_y^{+}$ are the forward finite-difference operators in the x- and y-directions, with unit spacing, and the values $u_{i,j}$ with indices outside $\Theta$ are defined by replication. We consider the following discrete version of (2) with anisotropic discrete total variation (TV):

$$\begin{array}{cc}{\displaystyle \underset{u,c_1,c_2}{minimize}}\hfill & {\displaystyle \sum _{i,j}(|{\nabla }_x{u|}_{i,j}+{\left|{\nabla }_{y}u\right|}_{i,j})+\lambda \sum _{i,j}\left({\left(c_1-{\overline{u}}_{i,j}\right)}^2{u}_{i,j}+{\left(c_2-{\overline{u}}_{i,j}\right)}^2(1-u_{i,j})\right),}\hfill \\ \mathrm{s}.\mathrm{t}.\hfill & 0\le u_{i,j}\le 1,\hfill \\ \hfill & c_1,c_2>0.\hfill \end{array}$$

(3)

After the minimization problem (3) is solved, the segmentation is obtained by taking

$$\Sigma =\left\{(i,j)\in \Theta :u_{i,j}>\alpha \right\},$$

(4)

for some $\alpha \in (0,1)$. Problem (3) is usually solved by alternating the minimization with respect to u and the minimization with respect to $c_1$ and $c_2$. In the sequel, we denote $E(u,c_1,c_2)$ the objective function in (3), and $F\left(u\right)$ and $G(u,c_1,c_2)$ its discrete regularization and fidelity terms, respectively.

The selection of a parameter $\lambda$ able to balance $F\left(u\right)$ and $G(u,c_1,c_2)$ and to produce a meaningful solution is a critical issue. Too large values of $\lambda$ may produce oversegmentation, while too small values of $\lambda$ may produce undersegmentation [9]. Furthermore, a constant value of $\lambda$ may not be suitable for the whole image, i.e., different regions of the image may need different values. In recent years, spatially adaptive techniques have been proposed, focusing on local information, such as gradient, curvature, texture, and noise estimation cues—see, e.g., [10]. Space-variant regularization has been also widely investigated in the context of image restoration, using TV and its generalizations—see, e.g., [11,12,13,14,15].

In this work, we propose some techniques for setting the parameter $\lambda$ in an adaptive way based on spatial information, in order to prevent excessive regularization of smooth regions while preserving spatial features in nonsmooth areas. Our techniques are based on the so-called cartoon-texture decomposition of the given image, on the mean and median filters, and on a thresholding technique. The resulting locally adaptive segmentation model can be solved either by smoothing the discrete TV term—see, e.g., [16,17]—and applying optimization solvers for differentiable problems, such as spectral gradient methods [18,19,20,21], or by using directly optimization solvers for nondifferentiable problems, such as Bregman, proximal and ADMM methods [22,23,24,25,26,27,28]. In this work, we use an alternating minimization procedure exploiting the split Bregman (SB) method proposed in [24]. The results of numerical experiments on images with different characteristics show the effectiveness of our approach and the advantages coming from using local regularization.

The rest of this paper is organized, as follows. In Section 2 we propose our spatially adaptive techniques. In Section 3, we describe the solution of the segmentation model using those techniques by the aforementioned SB-based alternating minimization method. In Section 4 we discuss the results obtained with our approaches on several test images, also performing a comparison with the CEN model and a segmentation model developed for images with texture. The results show the effectiveness of our approach and the advantages coming from the use of local regularization. Some conclusions are provided in Section 5.

2. Defining Local Regularization by Exploiting Spatial Information

The regularization parameter $\lambda$ plays an important role in the segmentation process, because it controls the tradeoff between data fidelity and regularization. In general, the smaller the parameter in (3) the smoother the image content, i.e., image details are flattened or blurred. Conversely, the larger the parameter the more the enhancement of image details and, hence, noise may be retained or amplified. Therefore, $\lambda$ should be selected according to local spatial information. A small value of $\lambda$ should be used in the smooth regions of the image to suppress noise, while a large value of $\lambda$ should be considered to preserve spatial features in the nonsmooth regions. In other words, a matrix $\Lambda =\left({\lambda }_{i,j}\right)$ should be associated with the image, where ${\lambda }_{i,j}$ weighs pixel $(i,j)$, as follows:

$$\begin{array}{cc}{\displaystyle \underset{u,c_1,c_2}{minimize}}\hfill & {\displaystyle \sum _{i,j}(|{\nabla }_x{u|}_{i,j}+{\left|{\nabla }_{y}u\right|}_{i,j})+\sum _{i,j}{\lambda }_{i,j}\left({\left(c_1-{\overline{u}}_{i,j}\right)}^2{u}_{i,j}+{\left(c_2-{\overline{u}}_{i,j}\right)}^2(1-u_{i,j})\right),}\hfill \\ \mathrm{s}.\mathrm{t}.\hfill & 0\le u_{i,j}\le 1,\hfill \\ \hfill & c_1,c_2>0.\hfill \end{array}$$

(5)

Furthermore, in order to avoid oversegmentation or undersegmentation, it is convenient to fix a minimum and a maximum value for the entries of $\Lambda$, so as to drive the level of regularization in a reasonable range, depending on the image to be segmented.

We define $\Lambda$ as a function of the image $\overline{u}$ to be segmented:

$$f:{\overline{u}}_{i,j}\to {\lambda }_{i,j}\in [{\lambda }_{min},{\lambda }_{max}],$$

(6)

where $0<{\lambda }_{min}<{\lambda }_{max}<\infty$. We propose three choices of f, detailed in the next subsections.

We note that problem (5) still has a unique solution for any fixed $(c_1,c_2)$, as a consequence of the next proposition.

Proposition 1.

For any fixed $(c_1,c_2)$, problem (5) is a convex problem.

Proof. 

Because the CEN model is convex for any fixed $(c_1,c_2)$, the thesis immediately follows from the fact that the parameters ${\lambda }_{i,j}$ are constant with respect to u.  □

2.1. Regularization Based on the Cartoon-Texture Decomposition

We define $f\left({\overline{u}}_{i,j}\right)$ by using the Cartoon-Texture Decomposition (CTD) of the image discussed in [29]. CTD splits any image u into the sum of two images, w and v, such that w represents the cartoon or geometric (piecewise-smooth) component of u, while v represents the oscillatory or textured component, i.e., v contains essentially textures, oscillating patterns, fine details, and noise. The algorithm for computing CTD acts as described next. For each image pixel, we decide whether it belongs to the cartoon or the textural part by computing a local indicator associated with an image window around the pixel. The main feature of a textured region is its high TV, which decreases very fast under low-pass filtering. This leads to the following definition of local total variation (LTV) at a pixel $(i,j)$:

$$LT{V}_{\sigma }{\left(u\right)}_{i,j}=(L_{\sigma }\ast {|\nabla u|}_{i,j}),$$

(7)

where $L_{\sigma }$ is a low-pass filter, $\sigma$ is a scale parameter, ${|\nabla u|}_{i,j}=\sqrt{{\left({\nabla }_{x}u\right)}_{i,j}^2+{\left({\nabla }_{y}u\right)}_{i,j}^2}$, and ∗ denotes the convolution product. The relative reduction rate of LTV,

$${\left({\rho }_{\sigma }\right)}_{i,j}=\frac{LT{V}_{\sigma }{\left(u\right)}_{i,j}-LT{V}_{\sigma }{(L_{\sigma }\ast u)}_{i,j}}{LT{V}_{\sigma }{\left(u\right)}_{i,j}},$$

(8)

gives the local oscillatory behavior of the image. A value of ${\left({\rho }_{\sigma }\right)}_{i,j}$ close to 0 means that there is little relative reduction of LTV by the low-pass filter, thus the pixel $(i,j)$ belongs to the cartoon. Conversely, ${\left({\rho }_{\sigma }\right)}_{i,j}$ close to 1 indicates that the relative reduction is large and, hence, the pixel belongs to a textured  region.

We use (8) for defining the weights ${\lambda }_{i,j}$. The basic idea is that a large regularization parameter is needed if the pixel $(i,j)$ belongs to the cartoon, while the parameter must be reduced in texture regions. Therefore, we set the function f in (6) as

$$f\left({\overline{u}}_{i,j}\right)\equiv f^{\mathrm{CTD}}\left({\overline{u}}_{i,j}\right)=max\left\{\frac{{\lambda }_{min}}{{\lambda }_{max}},1-{\left({\rho }_{\sigma }\right)}_{i,j}\right\}{\lambda }_{max}$$

(9)

and ${\left({\rho }_{\sigma }\right)}_{i,j}$ is defined by using the given image $\overline{u}$. Following [30], we set $L_{\sigma }$ equal to the Gaussian filter.

2.2. Regularization Based on the Mean and Median Filters

We define a technique that is based on spatial filters that are commonly used to enhance low-frequency details or to preserve edges [31,32]. More precisely, we combine the mean and median filters; the former aims at identifying smooth regions, where the regularization parameter can take small values, the latter aims at identifying edges, where the parameter must remain large. Mean filtering is a simple and easy-to-implement method for smoothing images, i.e., for reducing the intensity variation between a pixel and its neighbors. It also removes high-frequencies components due to the noise and the edges in the image, so the mean filter is a low-pass filter. The median filter preserves edges and useful details in the image.

Based on these considerations, we define a weight function, as follows:

$${\omega }_{i,j}=\left\{\begin{array}{cc}|{\overline{u}}_{i,j}-{(L_{h1}\ast \overline{u})}_{i,j}|& \mathrm{if} |{(L_{h1}\ast \overline{u})}_{i,j}-{(M_{h2}\ast \overline{u})}_{i,j}|<t,\hfill \\ 1& \mathrm{otherwise}.\hfill \end{array}\right.$$

(10)

where $h_1$ is the window size of the mean filter $L_{h1}$, $h_2$ is the window size of the median filter $M_{h2}$, and t is a threshold value acting as a cutoff between the two filters. Note that $0\le {\omega }_{i,j}\le 1$ and the pixels in homogeneous regions have ${\omega }_{i,j}$ close to 1. The function f in (6) is set as

$$f\left({\overline{u}}_{i,j}\right)\equiv f^{\mathrm{MM}}\left({\overline{u}}_{i,j}\right)=max\left\{\frac{{\lambda }_{min}}{{\lambda }_{max}},1-{\omega }_{i,j}\right\}{\lambda }_{max},$$

(11)

where MM stands for “Mean and Median”.

2.3. Regularization Based on Thresholding

This approach implicitly exploits the definition of $\Sigma$ in (4) in order to set ${\lambda }_{i,j}$. The idea is to use large values of ${\lambda }_{i,j}$ when $u_{i,j}$ is close to 1 and small values when $u_{i,j}$ is close to 0. Therefore, the parameters ${\lambda }_{i,j}$ are not defined in terms of the given image $\overline{u}$ only. If the function u identifying the segmentation were known a priori, we could define ${\lambda }_{i,j}=f\left(u_{i,j}\right)$ as follows:

$$f\left(u_{i,j}\right)\equiv f^{\mathrm{THR}}\left(u_{i,j}\right)={10}^{{\eta }_{i,j}}\in [{\lambda }_{min},{\lambda }_{max}],$$

(12)

where

$${\eta }_{i,j}=emax-(1-u_{i,j})(emax-emin),$$

$emax={log}_{10}{\lambda }_{max}$ and $emin={log}_{10}{\lambda }_{min}$.

Because the function u must be computed by minimizing (3) and this is done by using an iterative procedure, we decided to update ${\lambda }_{i,j}$ at each iteration, while using the current value of $u_{i,j}$. On the other hand, in this case evaluating f is computationally cheaper than in the previous approaches, which apply two-dimensional convolution operators; thus, the cost of the iterative update of ${\lambda }_{i,j}$ is practically negligible.

3. Solution by Split Bregman Iterations

As mentioned in Section 1, we use an alternating minimization method to solve problem (5). Given $u_{i,j}^{k-1}$, by imposing the first-order optimality conditions with respect to $c_1$ and $c_2$, we get

$$c_1^k=\frac{{\sum }_{i,j}{\lambda }_{i,j}{\overline{u}}_{i,j}{u}_{i,j}^{k-1}}{{\sum }_{i,j}{\lambda }_{i,j}{u}_{i,j}^{k-1}},c_2^k=\frac{{\sum }_{i,j}{\lambda }_{i,j}{\overline{u}}_{i,j}(1-u_{i,j}^{k-1})}{{\sum }_{i,j}{\lambda }_{i,j}(1-u_{i,j}^{k-1})}.$$

(13)

For the solution of (5) with respect to u, we use the split Bregman (SB) method [24]. Let

$$G_{loc}(u,c_1,c_2)=\sum _{i,j}{\lambda }_{i,j}\left({\left(c_1-{\overline{u}}_{i,j}\right)}^2-{\left(c_2-{\overline{u}}_{i,j}\right)}^2\right)u_{i,j}=\sum _{i,j}{r}_{i,j}{u}_{i,j},$$

(14)

where

$$r_{i,j}={\lambda }_{i,j}\left({\left(c_1-{\overline{u}}_{i,j}\right)}^2-{\left(c_2-{\overline{u}}_{i,j}\right)}^2\right).$$

Following [33], we reformulate the minimization problem, as follows:

$$\begin{array}{cc}{\displaystyle \underset{u,d_x,d_y}{minimize}}& {\displaystyle \parallel d_x{\parallel }_1+{\parallel d_y\parallel }_1+G_{loc}(u,c_1,c_2),}\hfill \\ \mathrm{s}.\mathrm{t}.& 0\le u_{i,j}\le 1,\hfill \\ & d_x={\nabla }_{x}u,\hfill \\ & d_y={\nabla }_{y}u.\hfill \end{array}$$

(15)

Given $c_1^k$ and $c_2^k$, the SB method applied to (15) reads:

$$\begin{array}{cc}\hfill u^k& =\underset{0\le u_{i,j}\le 1}{argmin}{G}_{loc}(u,c_1^k,c_2^k)+\frac{\mu }{2}\parallel d_x^{k-1}-{\nabla }_{x}u-b_x^{k-1}{\parallel }_2^2+\frac{\mu }{2}{\parallel d_y^{k-1}-{\nabla }_{y}u-b_y^{k-1}\parallel }_2^2,\hfill \\ \hfill d_x^k& =\underset{d_x}{argmin}\parallel d_x{\parallel }_1+\frac{\mu }{2}{\parallel d_x-{\nabla }_x{u}^k-b_x^{k-1}\parallel }_2^2,\hfill \\ \hfill d_x^k& =\underset{d_y}{argmin}\parallel d_y{\parallel }_1+\frac{\mu }{2}{\parallel d_y-{\nabla }_y{u}^k-b_y^{k-1}\parallel }_2^2,\hfill \\ \hfill b_x^k& =b_x^{k-1}+\mu ({\nabla }_x{u}^k-d_x^k),\hfill \\ \hfill b_y^k& =b_y^{k-1}+\mu ({\nabla }_x{u}^k-d_y^k),\hfill \end{array}$$

(16)

where $\mu >0$.

Closed-form solutions of the minimization problems with respect to $d_x$ and $d_y$ can be computed using the soft-thresholding operator:

$$d_x^k=\mathcal{S}\left({\nabla }_x{u}^k+b_x^k,\frac{1}{\mu }\right),d_y^k=\mathcal{S}\left({\nabla }_y{u}^k+b_y^k,\frac{1}{\mu }\right),$$

where, for any $v=\left(v_{i,j}\right)$ and any scalar $\gamma >0$,

$$\mathcal{S}(v,\gamma )=z=\left(z_{i,j}\right),z_{i,j}=\frac{z_{i,j}}{|z_{i,j}|}max\left\{|z_{i,j}|-\gamma ,0\right\}.$$

Finally, an approximate solution to the minimization problem with respect to u can be obtained by applying Gauss–Seidel (GS) iterations to the following system, as explained in [33]:

$$-\Delta u_{i,j}=\frac{r_{i,j}}{\mu }+\left({\nabla }_x{(b_x^k-d_x^k)}_{i,j}\right)+\left({\nabla }_y{(b_y^k-d_y^k)}_{i,j}\right),$$

(17)

where $\Delta$ is the classical finite-difference discretization of the Laplacian. If the solution to (17) lies outside ${[0,1]}^{m\times n}$, then it is projected onto that set. We denote ${\mathcal{P}}_{[0,1]}$ the corresponding projection  operator.

The overall solution method is outlined in Algorithm 1. Note that, when the approach in Section 2.3 is used, the values ${\lambda }_{i,j}$ must be updated at each iteration k, using (12) with $u=u^k$.

Algorithm 1 Split Bregman (SB)-based method for spatially adaptive segmentation
Input:  $\overline{u},{\lambda }_{min},{\lambda }_{max},f,\mu ,\alpha$ (with f defined in (9) or (11) or (12))
Output:  $u,c_1,c_2$
Set $u^0=\overline{u},d_x^0=0,d_y^0=0,b_x^0=0,b_y^0=0$
Compute $\Lambda =f\left(u^0\right)$
for $k=1,2,\dots$ do
	Compute $c_1^k$ and $c_2^k$ by (13)
	Compute $u^k$ by applying GS iterations to system (17)
	$u^k={\mathcal{P}}_{[0,1]}\left(u^k\right)$
	$d_x^k=\mathcal{S}({\nabla }_x{u}^k+b_x^k,1/\mu )$
	$d_y^k=\mathcal{S}({\nabla }_y{u}^k+b_y^k,1/\mu )$
	Update $b_x^k$ and $b_y^k$ according to (16)
end
Set $\Sigma =\left\{(i,j)\in \Theta :u_{i,j}^k>\alpha \right\}$

4. Results and Comparisons

The three spatially adaptive regularization techniques were implemented in MATLAB, while using the Image Processing Toolbox. Algorithm 1 was implemented by modifying the Fast Global Minimization Algorithm for Active Contour Models by X. Bresson, available from http://htmlpreview.github.io/?https://github.com/xbresson/old_codes/blob/master/codes.html. This is a C code with a MATLAB MEX  interface.

$L_{\sigma }$ in (7) is defined as a rotationally symmetric Gaussian filter with size 3 and standard deviation $\sigma =2$. The mean and median filters use windows of size 3 and 7, respectively, and the parameter t in (10) is set as $t=0.5$. The parameter $\alpha =0.5$ is used to identify the domain $\Sigma$ according to (4).

In the original and modified codes, the SB iterations are stopped, as follows:

$$|{\mathtt{diff}}^k-{\mathtt{diff}}^{k-1}|\le \mathtt{tol}\mathrm{and}k>\mathtt{maxit},$$

(18)

where

$${\mathtt{diff}}^l=\frac{\mathtt{sd}\left(u^l\right)}{\mathtt{sd}\left(u^l\right)·\mathtt{sd}\left(u^{l-1}\right)},\mathtt{sd}\left(u^l\right)=\sum _{i,j}{(u_{i,j}^l-u_{i,j}^l)}^2,$$

$\mathtt{tol}$ is a given tolerance and $\mathtt{maxit}$ denote the maximum number of outer iterations. The stopping criterion for the GS iterations is

$$E\left(u^k\right)=1-\frac{|{\mathtt{msd}}^k-{\mathtt{msd}}^1|}{{\mathtt{msd}}^1}\le {\mathtt{tol}}_{\mathtt{GS}}\mathrm{and}k>{\mathtt{maxit}}_{\mathtt{GS}}$$

where

$${\mathtt{msd}}^l=\frac{1}{mn}\sum _{i,j}{(u_{i,j}^l-u_{i,j}^{l-1})}^2,$$

and ${\mathtt{tol}}_{\mathtt{GS}}$ and ${\mathtt{maxit}}_{\mathtt{GS}}$ are the tolerance and the maximum number of iterations for the GS method, respectively. In our experiments we set $\mathtt{tol}={10}^{-6}$, $\mathtt{maxit}=30$, ${\mathtt{tol}}_{\mathtt{GS}}={10}^{-2}$ and ${\mathtt{maxit}}_{\mathtt{GS}}=50$.

The adaptive models are compared with the original CEN model on different images widely used in image processing tests, as listed in

Table 1. Test images and their sizes.

Image	Size (pixels)
bacteria	$233\times 256$
bacteria2	$380\times 391$
brain	$210\times 210$
cameraman	$204\times 204$
flowerbed	$321\times 481$
ninetyeight	$300\times 225$
squirrel	$167\times 230$
tiger	$321\times 481$

and shown in Figure 1, Figure 2, Figure 3 and Figure 4. In particular, the images bacteria, bacteria2, brain, cameraman, squirrel, and tiger are included in Bresson’s code distribution, flowerbed has been downloaded from the Berkeley segmentation dataset [34] available from https://www2.eecs.berkeley.edu/Research/Projects/CS/vision/grouping/resources.html (image #86016), and ninetyeight is available from https://tineye.com/query/2817cf0d186fbfe263a188952829a3b5e699d839. We note that tiger is included in Bresson’s code as a test problem for a segmentation model specifically developed for images with texture [35], which is also implemented in the code. This model uses the well-known Kullback–Leibler divergence function regularized with a TV term. The model is solved with the SB method. We perform the segmentation of tiger with the CEN model, the textural segmentation model, and our approaches, in order to investigate whether our locally adaptive model can be also suitable for textural images. The cameraman image is perturbed with additive Gaussian noise, with zero mean and two values of standard deviation, $\sigma =15,25$, with the aim of evaluating the behavior of our adaptive approaches on noisy images. The noisy images are called cameraman15 and cameraman25.

For the images provided with Bresson’s code, the values of $\lambda$ and $\mu$ associated with the original CEN model are set as in that code. For the remaining images, the values of $\lambda$ and $\mu$ are set by trial and error, following the empirical rule reported in Bresson’s code. The values of ${\lambda }_{min}$ and ${\lambda }_{max}$ used in the spatially adaptive approaches are chosen, so that the corresponding $\lambda$ in the original CEN model is in $[{\lambda }_{min},{\lambda }_{max}]$, with few exceptions. The associated values of $\mu$ are set as in the non-adaptive case. The values of $\lambda$, ${\lambda }_{min}$, ${\lambda }_{max}$ and $\mu$ used for each image are specified in

Table 2. Segmentation of six test images by the CEN and spatially adaptive models: parameters and iterations. The value of $\mu$ is the same for all the models, thus it is reported only once per image.

	$\mathit{\lambda }$	$\mathit{\mu }$	$\mathtt{it}$	${\mathtt{it}}_{\mathtt{GS}}$	${\mathit{\lambda }}_{\mathit{min}}$	${\mathit{\lambda }}_{\mathit{max}}$	$\mathbf{it}$	${\mathbf{it}}_{\mathbf{GS}}$
Image	CEN				CTD
brain	$0.7\times {10}^3$	$0.1\times {10}^4$	5	19.2	$0.9\times {10}^3$	$0.4\times {10}^5$	5	14.8
squirrel	$0.118\times {10}^3$	$0.1\times {10}^4$	6	50	$0.118\times {10}^3$	$0.1\times {10}^4$	15	50
					$0.2\times {10}^3$	$0.1\times {10}^4$	5	45
	CEN				MM
bacteria2	$0.1\times {10}^3$	$0.1\times {10}^4$	2	50	$0.1\times {10}^3$	$0.1\times {10}^4$	10	50
ninetyeight	$0.2\times {10}^2$	$0.1\times {10}^3$	4	38.5	$0.1\times {10}^2$	$0.5\times {10}^2$	3	50
	CEN				THR
bacteria	$0.5\times {10}^4$	$0.1\times {10}^5$	6	27.5	$0.4\times {10}^4$	$0.8\times {10}^4$	7	30.7
flowerbed	$0.1\times {10}^2$	$0.1\times {10}^4$	2	27	$0.3\times {10}^2$	$0.1\times {10}^3$	2	16.5

,

Table 3. Segmentation of tiger by the textural [35], CEN, CTD, and MM models: parameters and  iterations.

Image	$\mathit{\lambda }$	$\mathit{\mu }$	$\mathtt{it}$	${\mathtt{it}}_{\mathtt{GS}}$	${\mathit{\lambda }}_{\mathit{min}}$	${\mathit{\lambda }}_{\mathit{max}}$	$\mathit{\mu }$	$\mathbf{it}$	${\mathbf{it}}_{\mathbf{GS}}$
	Textural Model				CTD
tiger	$0.3\times {10}^1$	$0.1\times {10}^3$	5	12	$0.1\times {10}^1$	$0.1\times {10}^2$	$0.1\times {10}^2$	8	50
	CEN				MM
tiger	$0.3\times {10}^2$	$0.1\times {10}^3$	7	49.6	$0.1\times {10}^1$	$0.1\times {10}^2$	$0.1\times {10}^2$	8	50

and

Table 4. Segmentations of cameraman and its noisy versions by using CEN and the spatially adaptive models: parameters and iterations. The value of $\mu$ is the same for all of the models, thus it is reported only once per image.

	CEN				CTD
Image	$\lambda$	$\mu$	$\mathtt{it}$	${\mathtt{it}}_{\mathtt{GS}}$	${\lambda }_{min}$	${\lambda }_{max}$	$\mathtt{it}$	${\mathtt{it}}_{\mathtt{GS}}$
cameraman	$0.8\times {10}^3$	$0.1\times {10}^3$	4	2.8	$0.5\times {10}^3$	$0.1\times {10}^4$	4	3
cameraman15	$0.3\times {10}^3$	$0.1\times {10}^3$	5	4.2	$0.3\times {10}^3$	$0.1\times {10}^4$	5	4.2
cameraman25	$0.17\times {10}^3$	$0.1\times {10}^3$	6	7	$0.17\times {10}^3$	$0.8\times {10}^3$	6	7
	MM				THR
Image	${\lambda }_{min}$	${\lambda }_{max}$	$\mathtt{it}$	${\mathtt{it}}_{\mathtt{GS}}$	${\lambda }_{min}$	${\lambda }_{max}$	$\mathtt{it}$	${\mathtt{it}}_{\mathtt{GS}}$
cameraman	$0.5\times {10}^3$	$0.1\times {10}^4$	4	3	$0.5\times {10}^3$	$0.1\times {10}^4$	4	2.8
cameraman15	$0.3\times {10}^3$	$0.1\times {10}^4$	4	4.2	$0.3\times {10}^3$	$0.1\times {10}^4$	5	3.6
cameraman25	$0.17\times {10}^3$	$0.5\times {10}^3$	6	7	$0.17\times {10}^3$	$0.8\times {10}^3$	7	5.7

. The values of ${\lambda }_{i,j}$ in the adaptive models are initialized, as specified in (9), (11), and (12). As described in Section 2.3, in the strategy based on thresholding those values change at each iteration k of Algorithm 1. It is worth noting that our adaptive approaches also simplify the choice of the regularization parameter, which may be a time-consuming task.

The initial approximation $u^0$ is set equal to the given image ${\overline{u}}_{i,j}$, which takes its values in $[0,1]$, as specified in Section 1. This is used to compute the starting values of $c_1$ and $c_2$, and $s_{i,j}={\left(c_1-{\overline{u}}_{i,j}\right)}^2-{\left(c_2-{\overline{u}}_{i,j}\right)}^2$ for all $(i,j)$. In the original non-adaptive code, the value of $\lambda$ is scaled using the difference between the maximum and the minimum value of $s_{i,j}$; the same scaling is applied to ${\lambda }_{min}$ and ${\lambda }_{max}$ in the implementations of the adaptive approaches.

We run the tests on an Intel Core i7 processor with clock frequency of 2.6 GHz, 8 GB of RAM, and a 64-bit Linux system.

For each of the six images bacteria, bacteria2, brain, flowerbed, ninetyeight, and squirrel, we show the results that were obtained with the spatially adaptive strategy yielding the best segmentation for that image. The corresponding (unscaled) values of $\lambda$, ${\lambda }_{min}$, ${\lambda }_{max}$, the value of $\mu$, as well as the number of outer iterations and the mean number of GS iterations per outer iteration are reported in Table 2. Note that, for squirrel, we use two values of ${\lambda }_{min}$, one equal to the value of $\lambda$ in the CEN model and the other greater than that, obtained by trial and error. Both values of ${\lambda }_{min}$ produce the same segmentation, but the larger value of ${\lambda }_{max}$ reduces the number of outer and GS iterations. The segmentations corresponding to the data in Table 2 are shown in Figure 1, Figure 2 and Figure 3. For squirrel, we display the CTD segmentation computed by using the larger value of ${\lambda }_{min}$.

We see that, on selected problems, CTD reduces the number of outer and GS iterations with respect to the CEN model; on the other hand, the setup of the regularization parameters is computationally more expensive. In terms of iterations, there is no clear winner between CEN and MM and between CEN and THR. The models based on the spatially adaptive techniques are slightly more expensive than the CEN model in this case too. A significant result is that the segmentations obtained with the adaptive techniques appear better than those obtained with the non-adaptive model, i.e., the spatially adaptive models can better outline boundaries between objects and foreground. This is clearly visible by looking at the segmentations of brain, bacteria2 (see the upper right corner), ninetyeight, bacteria (see the shape of the bacteria). It is also worth noting that the adaptive model based on THR removes textural details that do not belong to the flowerbed in the homonymous image.

The latter observation is confirmed by the experiments performed on tiger. Table 3 reports the corresponding model and algorithmic details, while the segmentations are shown in Figure 4 along with (visual) information on quantities used to define ${\lambda }_{i,j}$ (see (8) and (10)). We see that the CTD and MM strategies produce satisfactory results, although they have been obtained by generalizing a non-textural model.

Finally, we show the results obtained on cameraman and its noisy versions by using CEN, CTD, MM, and THR. The methods perform comparable numbers of inner and GS iterations (see Table 4), but the spatially adaptive model THR yields some improvement over the CEN model in the segmentation of the images that are affected by Gaussian noise (Figure 5).

5. Conclusions

We introduced spatially adaptive regularization in a well-established variational segmentation model with the aim of improving the segmentation of images by suitably taking into account their smooth and nonsmooth regions. To this aim, we introduced three techniques, based on the application of suitable spatial filters or thresholding. The locally adaptive models, solved via an alternating minimization method using split Bregman iterations, showed the effectiveness of our approaches on several images, also including textural and noisy images. We also believe that the proposed models may simplify the setup of the regularization parameter.

Future work can include the extension of our spatially adaptive strategies to other segmentation models and the development of further adaptive techniques.