Local knowledge of the flora of a region: implications in biodiversity conservation

1.   Introduction

In this paper, we deal with a new variational problem suggested by applications to satellite image segmentation. The satellite images are an important source for extracting landscape boundaries and other vegetation structures, which can provide extremely useful insights for applications in environmental monitoring, agriculture, forestry, and other related fields (see, for instance, ^[1,2,3]). In particular, in agricultural crop field classification, one fundamental problem is to provide a disjunctive decomposition of a fixed domain $ \Omega\subset\mathbb{R}^2 $ onto a finite number of nonempty subsets $ \Omega = \Omega_1\cup\Omega_2\cup \dots\cup \Omega_K $ such that each of these subsets could be associated with a crop that is grown in this area, with forest regions, water zones, and so on, and this correspondence must be established at a rather a high level of accuracy. Up-to-date and accurate crop maps (or crop field classification) are needed to update agricultural statistics, to provide agricultural crop yield prediction, and are often used in environmental modeling. Typically, such an association between a given region and some agricultural crop can be made through the detection and quantitative assessment of green vegetation, which is one of the major applications of remote sensing studies. The information obtained in this way is a source of knowledge used for environmental resources management. One of the ways to get such information is the determination of the so-called vegetation indices (see ^[4] for a review and development). Since over the years many vegetation indices have been proposed for determining the vigor and health of vegetation, the reliability of information about vegetation directly and strictly depends on the fidelity, preciseness, and smoothness of the corresponding vegetation indices within each particular crop field.

The most commonly used vegetation index (Ⅵ) is the so-called slope-based infrared percentage vegetation index (IPVI),

where $ u_{i, d} = u_{i, d}(x_1, x_2) $, $ i = 1, 2 $, with $ (x_1, x_2)\in\Omega $, are functions of two variables representing the intensity of red ($ Red $) and near-infrared ($ NIR $) reflectance of some region $ \Omega $ of $ \mathbb{R}^2 $, respectively.

Thus, each pixel $ x = (x_1, x_2)\in\Omega $ of the original image can be associated with the corresponding IPVI-feature. The problem, which is suggested by application to remote sensing satellite image processing, consists of computing a decomposition

of the domain of the image $ F:\Omega\rightarrow \mathbb{R}^2 $ such that

$ ({\rm{a}}) $ the IPVI-characteristic varies smoothly and/or slowly within each $ \Omega_j $;

$ ({\rm{b}}) $ the IPVI-characteristic varies discontinuously and/or rapidly across most of the boundary $ K $ between different $ \Omega_j $.

The distinguished features of this statement that do not permit to reduce it to the standard settings of the segmentation problem (see, for instance, the Mumford-Shah energy-based model ^[5] or the models proposed by Alvarez ^[6], Guichard ^[7], Lions, Morel ^[8], Caselles ^[9], and others) are the following ones:

● Each region $ \Omega_1 $, $ \Omega_2 $, $ \dots $, $ \Omega_K $ should consist of pixels that can be reasonably grouped according to the IPVI-characteristic. Simultaneously, these regions should be easy to differentiate according to the chosen image feature;

● The respective interiors of image regions should have a more or less simple geometry without gaps. Boundaries of image regions should be smooth enough but also accurate with respect to the chosen image feature;

● The most restrictive obstacle in the construction of such decomposition is the fact that these subdomains should not overlap the borders between fields or contain any fragments of such borders, meaning that they cannot take in even small parts of different fields with arguably different crops.

All of these make the abovementioned segmentation problem rather challenging. It is enough to observe that a precise consideration of this problem demonstrates that the quantitative interpretation of remote sensing information from vegetation is a complex task. Many studies have limited this interpretation by assuming that the extracting vegetation information is uniformly and smoothly distributed within the particular crop fields. However, this assumption is broken when trying to apply these types of vegetation indices on heterogeneous canopies such as plantations with a mixed combination of soil, weeds, and other crops, or plantation where the vegetation of interest has different IPVI-characteristic due to spatial variability. The main idea, we realize in the new setting of the variational problem, can be briefly described as follows. We propose to make use of the so-called $ f $-decomposition instead of the standard Chan-Vese "active contours without edges" model ^[10]. The role of the function $ f:\Omega\rightarrow \mathbb{R} $ in such decomposition of $ \Omega $ has to guarantee that the new objects $ \left\{\Omega_j\right\}_{j = 1}^K $ after the $ f $-decomposition will have homogeneous values of the target function $ f $ within each separate field (a similar point of view can be found in ^[4,11]). In particular, in the case of the agricultural applications, where $ \Omega $ stands for a zone of interest, the IPVI-characteristic can be considered as the main feature of this area, i.e., in this case $ f(x) = \mathrm{IPVI}(x) $ for all $ x\in\Omega $. Thus, the main idea that we push forward is to formulate the segmentation problem as a constrained minimization problem in a special anisotropic functional space, with the "effect of anisotropy" we associate with the structure topology of IPVI-distribution. As a result, the main benefit of such an approach can be briefly described as follows:

$ ({\rm{i}}) $ It prevents the appearance of subdomains containing zones of discontinuity of $ f $ or places where this function tends to change rapidly by utilizing the main characteristic of the given function $ f $ — the unit normal vector field $ \theta:\Omega\rightarrow \mathbb{R}^2 $ to the level sets of $ f $. This characteristic has been used to construct the so-called anisotropic diffusion tensor $ M^f $, which can be defined as a square parametrized matrix function $ M^f(x) = \left[I-\eta^2\theta(x)\otimes\theta(x)\right] $. This matrix plays a central role in the process of $ f $-decomposition of domain $ \Omega $, and we associate with it special anisotropic perimeters of the obtained subdomains (segments).

$ ({\rm{ii}}) $ The second characteristic feature of our approach is the fact that we apply the Jeffreys divergence to replace the standard Euclidean distance in the fidelity term of the objective functional. It is well-known that compared with Euclidean distance, Jeffreys divergence leads to more accurate results in information measurement (for the details of this metric and its advantages, we refer to ^[12,13,14]).

The paper is organized as follows. In Section 2, we give some preliminaries related to the space of functions of bounded variation and other notions. Section 3 is devoted to the description of some specification of the standard $ BV(\Omega) $ space. In particular, we introduce the so-called anisotropic version for the total variation of $ L^1(\Omega) $-functions. At the end of this section, we show that some of the results of Samson et al. ^[15] can be extended to the case of subsets with a finite anisotropic perimeter.

The precise setting of the main constrained minimization problem and its previous analysis are given in Section 4. We show that the proposed minimization problem can be interpreted as a special case of the piecewise-constant Mumford-Shah segmentation problem and the Chan-Vese active contour model without edges. We study this problem in the space of $ L^1(\Omega) $-functions with bounded anisotropic total variation, where the type of anisotropy is closely related to the structure of the image $ f $ which is involved in the segmentation procedure. It is worth emphasizing that the anisotropic perimeter of the segments with uniform distribution of $ IPVI $-characteristic can drastically differ from the standard one because the natural edges of the original image $ f $ can affect it significantly. Despite the "natural" setting of the proposed segmentation problem, the existence of its minimizers seems to be an open issue nowadays. The main reason is that the objective functional is neither coercive nor lower semicontinuous with respect to the weak-$ \ast $ topology of $ BV(\Omega) $ space. This circumstance stimulated us to introduce a special family of unconstrained two-parametric problems to approximate the original one. We show that each of those approximated problems is well-posed and has a nonempty set of minimizers.

Section 5 aims to derive optimality conditions for approximated problems and provide their formal substantiation. In Section 6, we study the asymptotic behavior of the approximated problems and their solutions. The main question is to find out whether the convergence of minima of approximated problems is to minima of the original segmentation problem as small parameters tend to zero. Our main result of this section asserts that: If the original problem has a nonempty set of minimizers, then some of them can be successfully attained by the solution of approximated problems. Otherwise, we can come to a solution of the relaxed version of the original segmentation problem. In Section 7, the implementation of the proposed optimization problem is illustrated, providing numerical experiences with satellite images.

A detailed description of the algorithm (including the method of marching squares for the generation of closed contours in the two-dimensional case) and finite-difference scheme for the proposed approach with the results of numerical simulation using the real-life satellite images will be considered in the forthcoming paper.

2.   Auxiliaries

We denote by $ \mathcal{L}^2 $ the Lebesgue $ 2 $-dimensional measure in $ \mathbb{R}^2 $ and by $ \mathcal{H}^1 $ the $ 1 $-dimensional Hausdorff measure. Let $ \Omega $ be a bounded open subset of $ \mathbb{R}^2 $ with a Lipschitz boundary. For any subset $ E\subset\Omega $, we denote by $ |E| $ its $ 2 $-dimensional Lebesgue measure $ \mathcal{L}^2(E) $. For a subset $ E\subseteq \Omega $, let $ \overline{E} $ denote its closure and $ \partial E $ its boundary. We define the characteristic function $ \chi_E $ of $ E $ by

For a function $ u $, we denote by $ u|_E $ its restriction to the set $ E\subseteq\Omega $, and by $ u^{\partial E} $ its trace on $ \partial E $. Let $ C^\infty_0(\Omega) $ be the infinitely differentiable functions with compact support in $ \Omega $. The $ k $-dimensional Hausdorff measure is denoted by $ \mathcal{H}^k $, and is the restriction of the measure $ \mu $ to the set $ E $. For a Banach space $ X $, its dual is $ X^\ast $ and $ \left < \cdot, \cdot\right > _{X^\ast, X} $ is the duality form on $ X^\ast\times X $. By $ \rightharpoonup $ and $ \stackrel{\ast}{\rightharpoonup} $, we denote the weak and weak$ ^\ast $ convergence in normed spaces.

We remind here of the most common definitions of some functional spaces that we will use later on.

2.1.   Weak compactness criterion in $ L^1(\Omega) $

Throughout the paper, we will often use the concept of weak and strong convergence in $ L^1(\Omega) $. Let $ \left\{f_n\right\}_{n\in \mathbb{N}} $ be a bounded sequence of functions in $ L^1(\Omega) $. We recall that $ \left\{f_n\right\}_{n\in \mathbb{N}} $ is called equi-integrable on $ \Omega $ if for any $ \delta > 0 $ there is a $ \tau = \tau(\delta) $, such that $ \int_S |f_n|\, dx < \delta $ for every measurable subset $ S\subset\Omega $ of Lebesgue measure $ |S| < \tau $. Then, the following assertions are equivalent for $ L^1(\Omega) $-bounded sequences (see, for instance, ^[16,17]):

$ ({\rm{i}}) $ a sequence $ \left\{f_k\right\}_{k\in \mathbb{N}} $ is weakly convergent in $ L^1(\Omega) $;

$ ({\rm{ii}}) $ the sequence $ \left\{f_k\right\}_{k\in \mathbb{N}} $ is equi-integrable.

The following theorem holds.

Theorem 1. ^[16] If a bounded sequence $ \left\{f_k\right\}_{k\in \mathbb{N}}\subset L^1(\Omega) $ is equi-integrable and $ f_k\rightarrow f $ almost everywhere in $ \Omega $, then $ f_k\rightarrow f $ strongly in $ L^1(\Omega) $.

2.2.   Functions of bounded variation

We set $ |E| = \mathcal{L}^2(E) $, the Lebesgue measure of a measurable set $ E\subset\mathbb{R}^2 $. Let $ \mathcal{M}(\Omega; \mathbb{R}^2) $ be the space of all $ \mathbb{R}^2 $-valued Borel measures which is, according to the Riesz theory, the dual of the space $ C_0(\Omega; \mathbb{R}^2) $ of all continuous vector-valued functions $ \varphi(\cdot) $ with compact support in $ \Omega $ and equipped with the uniform norm.

Note that $ \mathcal{M}(\Omega; \mathbb{R}^2) $ is isomorphic to the product space

and that $ \mu = (\mu_1, \mu_2)\in \mathcal{M}(\Omega; \mathbb{R}^2) $ $ \Leftrightarrow $ $ \mu_i\in \left[C_0(\Omega)\right]^\ast $, $ i = 1, 2 $.

Given a vector-valued measure $ \mu:\mathcal{B}(\Omega)\rightarrow \mathbb{R}^2 $, we use the notation $ |\mu| $ for its total variation. We recall that

for all measurable $ E\subseteq \Omega $.

The usual weak-$ \ast $ topology on $ \mathcal{M}(\Omega; \mathbb{R}^2) $ is defined as the weakest topology on $ \mathcal{M}(\Omega; \mathbb{R}^2) $, for which the maps $ \mu\mapsto \sum_{i = 1}^2 \int_{\Omega} \varphi_i\, d\mu_i $ are continuous for every $ \varphi\in C_0(\Omega; \mathbb{R}^2) $.

By $ BV(\Omega) $, we denote the space of all functions $ u\in L^1(\Omega) $, for which their distributional derivatives are representable by finite Borel measures in $ \Omega $, i.e.,

for some $ \mathbb{R}^2 $-valued measure $ Du = (D_1u, D_2u)\in \mathcal{M}^2(\Omega) $. It can be shown that $ BV(\Omega) $, endowed with the norm $ \|u\|_{BV(\Omega)} = \|u\|_{L^1(\Omega)}+|Du|(\Omega) $, is a Banach space, where in view of Eq (2.1), the total variation of $ Du $ in $ \Omega $ can be defined as

We recall that the product topology of the strong topology of $ L^1(\Omega) $ for $ u $ and of the weak-$ \ast $ topology of measures for $ Du $ is called the weak-$ \ast $ topology of $ BV(\Omega) $, and it is denoted $ BV $-$ \ast $. As a result, a sequence $ \{f_k\}_{k = 1}^\infty $ $ \ast $-converges to $ f $ in $ BV(\Omega) $ if, and only if, the two following conditions hold (see [18, p.124]): $ f_k\to f $ strongly in $ L^1(\Omega) $ and $ D f_k\stackrel{\ast}{\rightharpoonup} Df $ weakly-$ ^\ast $ in $ \mathcal{M}(\Omega; \mathbb{R}^2) $, i.e.,

where, in fact, $ Df_k = (D_{x_1}f_k, D_{x_2} f_k)\in \mathcal{M}(\Omega; \mathbb{R}^2) $ and, therefore, the notation $ \int_{\Omega} \phi\, Df_k $ should be interpreted as follows:

Moreover, if $ \{f_k\}_{k = 1}^\infty\subset BV(\Omega) $ converges strongly to some $ f $ in $ L^1(\Omega) $ and $ \sup_{k\in \mathbb{N}}\int_\Omega |Df_k| < +\infty $, then (see, for instance, ^[16] and ^[18])

A simple criterion for the $ BV $-$ \ast $ convergence can be stated as follows (see [18, p.125], [19, Theorem 1.19]):

Proposition 2. A sequence $ \left\{u_k\right\}_{k\in \mathbb{N}}\subset BV(\Omega) $ $ BV $-$ \ast $-converges to $ u $ if, and only if, $ \left\{u_k\right\}_{k\in \mathbb{N}} $ is bounded in $ BV(\Omega) $ and $ \left\{u_k\right\}_{k\in \mathbb{N}} $ converges to $ u $ strongly in $ L^1(\Omega) $.

The following embedding result for the $ BV $-function is very useful with respect to the variational problem that we study in this paper.

Proposition 3. [16, p.378] Let $ \Omega $ be an open bounded Lipschitz subset of $ \mathbb{R}^2 $. Then, the embedding $ BV(\Omega)\hookrightarrow L^2(\Omega) $ is continuous and the embeddings $ BV(\Omega)\hookrightarrow L^p(\Omega) $ are compact for all $ p $ such that $ 1\le p < 2 $. Moreover, there exists a constant $ C_{em} > 0 $, which depends only on $ \Omega $ and $ p $ such that for all $ u $ in $ BV(\Omega) $,

We also make use of the following property concerning to approximation of $ BV $-functions by smooth ones.

Theorem 4. ^[20] Assume $ f\in BV(\Omega) $. Then, there exist a sequence $ \left\{f_k\right\}_{k = 1}^\infty\subset BV(\Omega)\cap C^\infty(\Omega) $ such that

Let $ E $ be an $ \mathcal{L}^2 $-measurable subset of $ \mathbb{R}^2 $ with finite Lebesgue measure. Let $ \chi_E $ be its characteristic function. Following R. Caccioppoli ^[21], we say that $ E $ is a set with a finite perimeter in $ \Omega $ if $ \chi_E\in BV(\Omega) $. This means that the distributional gradient $ D\chi_E $ is a vector-valued measure with finite total variation. The total variation $ |D\chi_E|(\Omega) $ is called the perimeter of $ E $ in $ \Omega $, i.e., $ P(E, \Omega) = |D\chi_E|(\Omega) $ and, therefore,

We also notice that if $ E_t: = \left\{x\in\Omega\ :\ f(x) > t\right\} $ stands for the level set for given $ f\in BV(\Omega) $ and $ t\in\mathbb{R} $, then (see [20, Theorem 5.5.1]) $ E_t $ has a finite perimeter for a.e. $ t\in\mathbb{R} $.

3.   Functions with bounded anisotropic total variation

The main goal of this section is to introduce some specification to the standard space of functions with bounded variation $ BV(\Omega) $. This option is mainly motivated by the natural application in image segmentation problems. In view of this, we introduce the so-called anisotropic version for the total variation of the $ BV $-functions. In principle, the notion of anisotropic total variation is not new in the literature (we refer to ^[22,23] for more details), and our representation for anisotropic total variation can be viewed as some specification of the rule in ^[22]. The main interest is in the application of this concept to the generalization of the well-known results of Samson et al. ^[15]. Namely, we focus on the proof of the following relation:

where $ \left[\chi_E\right]_{{\varepsilon}} $ stands for a smooth approximation of the characteristic function of a given set $ E\subset\Omega $. The main results of this section are presented in the form of Lemma 5 and its corollary.

Let $ \Omega $ be an open bounded and connected subset of $ \mathbb{R}^2 $ with a Lipschitz boundary $ \partial\Omega $. Let $ M:\Omega\rightarrow \mathbb{R}^{2\times 2} $ be a given matrix function such that

for some constant $ \beta > 1 $, i.e., $ M(x) $ is a positive-definite symmetric matrix for each $ x\in\Omega $.

We say that $ u\in L^1(\Omega) $ is a function with bounded anisotropic variation if

It means that there exists a Radon measure $ Du $ such that

Moreover, for the total variation of the measure $ MDu $, we have the following representation:

Then, property (3.1) implies that

Hence, the expression can be viewed as an equivalent norm to standard one $ \|\cdot\|_{BV_M(\Omega)} $ on the space $ BV(\Omega) $. As a result, the main properties of $ BV $-functions (see, for instance, ^[18,19,20]) can be reformulated with respect to the new norm. In particular, let $ M:\Omega\rightarrow \mathbb{R}^{2\times 2} $ be a given matrix function with property (3.1). Then:

$ ({\rm{j}}) $ If $ \left\{u_k\right\}_{k\in\mathbb{N}}\subset BV(\Omega) $ is a bounded sequence, then there exist a subsequence $ \left\{u_{k_i}\right\}_{i\in\mathbb{N}} $ and a function $ u\in BV(\Omega) $ such that

$ ({\rm{jj}}) $ If $ \{u_k\}_{k = 1}^\infty\subset BV(\Omega) $ converges strongly to some $ u $ in $ L^1(\Omega) $ and satisfies $ \sup_{k\in \mathbb{N}}\int_\Omega|MDu_k| < +\infty $, then

$ ({\rm{jjj}}) $ Let $ u\in BV(\Omega) $ be an arbitrary function. Then, there exists a sequence $ \left\{u_k\right\}_{k\in\mathbb{N}}\subset C^\infty(\Omega)\cap BV(\Omega) $ such that $ \|u_k-u\|_{L^1(\Omega)}\rightarrow 0 $ and $ { }\int_{\Omega} |M Du_k|\rightarrow \int_{\Omega} |M Du| $ as $ k\to\infty $.

By analogy with the standard notion, we can also define an anisotropic version of the perimeter, namely, we say that an $ \mathcal{L}^2 $-measurable subset $ U\subset\Omega $ has a finite $ M $-perimeter if $ |M D\chi_U|(\Omega) < +\infty $, where $ \chi_U(\cdot) $ stands for the characteristic function of the set $ U $. In this case, we write down

Moreover, for any $ u\in BV(\Omega) $ and $ M(\cdot)\in C^\infty(\Omega; \mathbb{R}^{2\times 2}) $ with property (3.1), the following anisotropic Coarea formula

holds true (see ^[22]).

It is clear that, for a given level parameter $ l\in\mathbb{R} $, the area of the region $ \left\{x\in\Omega\ :\ \varphi(x) > l\right\} $ can be defined as

Here, $ \chi_{\{\varphi(x)\ge l\}}(x) $ stands for the characteristic function of the set $ \{x\in\Omega:\varphi(x)\ge l\} $. To obtain some approximation of this area, we can substitute $ \chi_E $ by its smooth approximation. With that in mind, for a given small positive parameter $ {{\varepsilon}} $, we fix a positive symmetric mollifier $ \eta\in C^\infty_c(\mathbb{R}^2) $, i.e., $ \eta(x) $ is zero outside a compact set $ B_1 = \left\{x\in\mathbb{R}^2: |x|\le 1\right\} $,

and set

that is,

Then, using the standard properties of mollifiers, it can be shown that

$ ({\rm{i}}) $ $ \left[\chi_E\right]_{{\varepsilon}}\rightarrow \chi_E $ in $ L^1(\Omega) $ for any measurable subset $ E $ of $ \mathbb{R}^2 $ as $ \emptyset\to 0 $;

$ ({\rm{ii}}) $ $ 0\le \left[\chi_E\right]_{{\varepsilon}}(x)\le 1 $ for all $ x\in \mathbb{R}^2 $;

$ ({\rm{iii}}) $ If $ E\subset\mathbb{R}^2 $ is bounded and $ g\in L^1(\mathbb{R}^2) $, then $ { }\int_{\mathbb{R}^2} \left[\chi_E\right]_{{\varepsilon}} g\, dx = \int_{\mathbb{R}^2}\chi_E\left[g\right]_{{\varepsilon}}\, dx $;

$ ({\rm{iv}}) $ If $ {{\varepsilon}}\subseteq \Omega $, then $ \mathrm{supp}\, \left[\chi_E\right]_{{\varepsilon}}\subseteq \Omega_{{\varepsilon}} = \left\{x\in\mathbb{R}^2 : \mathrm{dist}\, (x, \Omega)\le{{\varepsilon}} \right\} $.

The following property will be utilized in our further analysis (see Section 6) and it can be considered as a natural generalization of the results of Samson et al. ^[15] (see also for comparison [19, Proposition 1.15]).

Lemma 5. Let $ E $ be an open set such that $ \overline{E}\subset\Omega $ and $ E $ has a finite $ M $-perimeter $ \mathrm{Per}(E; M;\Omega) $, where $ M:\Omega\rightarrow \mathbb{R}^{2\times 2} $ is a given matrix function with property (3.1). Let $ \left[\chi_E\right]_{{\varepsilon}} $ be the mollified characteristic function described above. Then,

Proof. Taking into account the standard properties of mollifiers, we have $ \left\{\left[\chi_E\right]_{{\varepsilon}}\right\}_{{{\varepsilon}} > 0}\subset BV(\Omega) $ and $ \left[\chi_E\right]_{{\varepsilon}}\rightarrow \chi_E $ in $ L^1(\Omega) $ as $ {{\varepsilon}}\to 0 $. Then, inequality (3.4) implies that

To establish a reverse inequality, we fix an arbitrary function $ \varphi\in C^1_0(\Omega; \mathbb{R}^2) $ with $ |\varphi(x)|\le 1 $. Then, there exists a vector-valued function $ \zeta\in C^1_0(\Omega; \mathbb{R}^2) $ such that $ \zeta = M\varphi $. Therefore, by (iii)-property of mollifiers, we have

Taking then the supremum over all such $ \varphi $, we arrive at the relation

Hence,

Since $ E $ and $ \Omega $ are open sets and $ \overline{E}\subset\Omega $, it follows that

Therefore,

As a result, we obtain

It remains to combine this inequality with Eq (3.9).

Arguing similarly, we can generalize the Eq (3.8) as follows:

Corollary 6. Let $ E\subset\Omega $ and the matrix $ M:\Omega\rightarrow \mathbb{R}^{2\times 2} $ be the same as in Lemma 5. Then,

where $ \left\{M_{{\varepsilon}}\right\}_{{{\varepsilon}} > 0}\subset C^\infty(\Omega; \mathbb{R}^{2\times 2}) $ stands for any smooth approximation of the matrix $ M $ such that

4.   Setting of optimization problem and its previous analysis

Let $ f:\Omega\rightarrow \mathbb{R} $ be a given function. Hereinafter, we associate the function $ f $ with a given gray scale image. We define a smoothed version of the original image using the convolution of $ f $ with a Gaussian kernel

i.e., $ f_{\sigma} = \left(G_\sigma\ast f(\cdot)\right)(x): = { }\int_{\Omega} G_\sigma(x-y) f(y)\, dy $. Here, $ \sigma > 0 $ is a given small positive value.

Since $ f_{\sigma}\in C^\infty(\overline{\Omega}) $, it follows that the boundaries of level sets $ \left\{x\in\Omega\, :\, f_{\sigma}(x)\ge\lambda \right\} $, for all feasible $ \lambda\in[0, C_f] $, can be described by smooth curves with finite length. So, at all points $ x\in\Omega $ of each level sets of $ f_{\sigma} $, we can define a unit normal vector field $ \theta(x) $ following the rule

Then, we associate with the vector field $ \theta:\Omega\rightarrow \mathbb{R}^2 $ the following linear operator $ R:\mathbb{R}^2\rightarrow \mathbb{R}^2 $:

where $ \eta\in (0, 1) $ is a given threshold which should be sufficiently close to $ 1 $. Setting

we see that

i.e., $ M^f $ is a symmetric and positive-definite matrix on $ \Omega $ and it satisfies property (3.1) with $ \beta = (1-\eta^2)^{-1} $.

Remark 7. Let's assume for a moment that $ \xi = \nabla v(x) $, where $ v\in W^{1, 1}(\Omega) $ is a given function and $ x\in\Omega $ is a Lebesgue point of $ f $ in which the original image $ f $ is not expected to change drastically in any direction, i.e., $ x $ is not close to a discontinuity of $ f $ or a zone where $ f $ tends to change rapidly. Then, Eq (4.4) implies that $ M^f $ can be represented at this point as a unit matrix. So, at this point, we have $ M^f\nabla v\approx \nabla v $.

On the other hand, if we consider a point $ x\in\Omega $ that is close to a discontinuity of $ f $, then $ M^f\nabla v $ reduces to $ (1-\eta^2)\nabla v $ if the gradient $ \nabla v(x) $ at this point is colinear to $ \theta(x) $, and to $ \nabla v(x) $ provided $ \nabla v(x) $ is orthogonal to $ \theta(x) $. In view of this, the expression $ M^f\nabla v $ can be interpreted as the directional total variation of $ v $ along the vector field $ \theta $ (see ^[24] for the details).

Definition 8. We say that a gray scale image $ f:\Omega\rightarrow \mathbb{R} $ is feasible for the segmentation procedure using level sets if there exists a value $ \gamma > 0 $ such that

We denote the set of all feasible images by $ \mathcal{F}_\gamma $.

We are now in a position to state the main object of our interest in this paper. Let $ f\in\mathcal{F}_\gamma $ be a given image, and let $ \varphi:\Omega\rightarrow \mathbb{R} $ be a level set function such that $ \varphi\in BV(\Omega) $. We associate with this function a collection of $ m+2 $ distinct level values $ l_0 < l_1 < \dots < l_m < l_{m+1} $ such that $ l_0\le \varphi(x)\le l_{m+1} $ almost everywhere in $ \Omega $.

Then, the constrained optimization problem we are going to consider can be stated as follows:

where $ \alpha > 0 $ is a weight coefficient and the set of feasible solutions is defined as follows:

It is worth noticing that the objective functional $ J $ is well-defined on the set $ \Xi $. Indeed, in this case the assumption $ \varphi\in BV(\Omega) $ implies that the level sets $ E_t = \left\{x\in \Omega: \varphi(x) > t\right\} $ have finite anisotropic perimeter $ \mathrm{Per}(\left\{u > t\right\}; M^f; \Omega) $ for $ \mathcal{L}^1 $ a.a. $ t\in\mathbb{R} $. Since

for each $ j = 1, \dots, m $, it follows that $ J(c, \varphi) < \infty $ for each $ (c, \varphi)\in\Xi $.

Minimization problems (4.6) and (4.7) can be interpreted as a special case of the piecewise-constant Mumford-Shah segmentation problem and the active contour model (see, for instance, ^[8,10,25,26]). We can indicate the following principle features of this statement:

● The problem is investigated in the space of $ L^1(\Omega) $-functions with bounded anisotropic total variation and with additional pointwise constraints, where the matrix of anisotropy is closely related to the structure of the image $ f $, which is involved in the segmentation procedure.

As a result, the anisotropic perimeter of the region $ \left\{x\in\Omega: l_j < \varphi(x) \right\} $, which is given by the term $ \int_{\Omega} |M^f D\chi_{\{\varphi(x) > l_j\}}| $, can drastically differ from the standard one because the natural edges of the original image $ f $ can affect it significantly;

● To find a piecewise-constant approximation of the given image $ f\in\mathcal{F}_\gamma $ in the form

we utilize the Jeffreys divergence between two elements $ f, g\in L^2(\Omega) $ instead of the standard $ L^2 $-norm of their difference $ \|f-g\|_{L^2(\Omega)} $. In spite of the fact that the trick of replacing the squared Euclidean norm with the Jeffreys distance does not alleviate the original problem from the point of view of its solvability and mathematical analysis, it makes the segmentation results more stable to the Poisson noise contaminated images (see ^[12,13,14] for the details);

● We consider the segmentation problems (4.6) and (4.7) as a constrained minimization problem in $ BV(\Omega) $ space with the pointwise constraints $ l_0\le \varphi(x)\le l_{m+1} $ on the set of feasible functions $ \varphi:\Omega\rightarrow \mathbb{R} $.

● In the proposed statement of the segmentation problems (4.6) and (4.7), it is admitted that the set

may have a nonzero $ \mathcal{L}^2 $-measure.

However, the existence of minimizers to the problems (4.6) and (4.7) seems to be an open issue nowadays because the standard application of the direct method of calculus of variation to this problem faces some unsolved challenges. To apply the direct method for proving the existence of minimizers, it is necessary to find a topology for which the functional (4.6) is lower semicontinuous while ensuring compactness of minimizing sequences. In view of the structure of the set of admissible solutions $ \Xi\subset \mathbb{R}^{m+1}\times L^1(\Omega) $ (see Eq (4.7)), the natural topology, in this case, is the product of the norm topology in $ \mathbb{R}^{m+1} $ and the weak-$ \ast $ topology in $ BV(\Omega) $. However, the objective functional $ J:\Xi\rightarrow \mathbb{R} $ is not coercive and lower semicontinuous on $ \mathbb{R}^{m+1}\times L^1(\Omega) $ with respect to the above mentioned topology. Moreover, even if $ \varphi_k\rightarrow \varphi $ strongly in $ L^1(\Omega) $ as $ k\to\infty $, it does not imply the strong convergence in $ L^1(\Omega) $ of $ \chi_{\{\varphi_k(x) > l_j\}} $ to $ \chi_{\{\varphi(x) > l_j\}} $. In particular, it is clear that the implication

may hold true if only $ \mathcal{L}^2\left\{x\in\Omega:\varphi(x) = l_j\right\} = 0 $.

To overcome this difficulty, we make use of the following family of two-parametric approximated problems:

where $ \tau $ and $ {{\varepsilon}} $ are small parameters, which vary within strictly decreasing sequences of positive numbers converging to $ 0 $. The functions $ \left[\chi_{E_j^\tau}\right]_{{\varepsilon}}(\cdot)\in C^\infty(\mathbb{R}) $, $ j = 1, \dots, m $, are defined in Eq (3.7),

and

Thus, a pair $ (c, \varphi) $ sounds as feasible to the problems (4.6) and (4.7) if $ (c, \varphi)\in\Xi_{{\varepsilon}} $, i.e.,

Before proceeding further, we make use of the following property.

Lemma 9. Let $ f\in\mathcal{F}_\gamma $ be a given image, let $ c = (c_0, c_1, \dots, c_m) $ with $ c_j > 0 $, let $ \{\varphi_k\}_{k = 1}^\infty\subset L^1(\Omega) $ be a strongly convergent sequence, and let $ \varphi\in L^1(\Omega) $ be its limit. Then, for each $ j = 1, \dots, m-1 $,

Proof. Since $ L^1(\Omega)\ni\varphi_k\rightarrow \varphi $ as $ k\to\infty $, we may suppose that, up to a subsequence, $ \varphi_k(x)\rightarrow \varphi(x) $ almost everywhere in $ \Omega $. Then, taking into account the continuity of the function $ \left[\chi_{E}\right]_{{\varepsilon}}(\cdot) $ and the fact that the passage in the inequality $ \varphi_k(x)\ge l_j+\tau $ leads to $ \varphi(x)\ge l_j+\tau $ for a.e $ x\in\Omega) $, we see that

Moreover, since $ \left|\left[\chi_{\{\varphi_k(x)\ge l_j+\tau\}}\right]_{{\varepsilon}}(x)\right|\le 1 $ in $ \Omega $, it follows from the Lebesgue dominated theorem that

Since $ f\in\mathcal{F}_\gamma $ and $ c = (c_0, c_1, \dots, c_m) $ with $ c_j > 0 $, it follows that the integrant $ \left(f-c_{j}\right)\log\left(\frac{f}{c_{j}}\right) $ is dominated by some integrable function $ g $ in the sense that

In particular,

As a result, Eq (4.14) is a direct consequence of the Lebesgue dominated theorem. Arguing similarly, we can establish the same assertion for the first and the third terms in Eq (4.9). □

To conclude this section, we give an existence result for the parametrized optimization problems (4.9)–(4.13).

Theorem 10. Let $ f\in\mathcal{F}(\gamma) $ be a given gray scale image, and let

where the vector field $ \theta\in C^\infty(\Omega; \mathbb{R}^{2}) $ is defined by the rule (4.2). Then, for each $ {{\varepsilon}}\in(0, 1) $ and $ \tau > 0 $ small enough, the constrained minimization problems (4.9)–(4.13) admit at least one solution.

Proof. Since the objective functional is bounded from below on $ \Xi_{{\varepsilon}}\subset \mathbb{R}^{m+1}\times BV(\Omega) $, it follows that there exists a minimizing sequence to problems (4.9)–(4.13), i.e.,

where $ C $ stands for a strictly positive constant that can be different from line to line. Without loss of generality, we can suppose that $ C = \zeta+1 $.

Taking into account the fact that

and

we see that

Since

and the function $ K:\mathbb{R}_{\ge}\rightarrow \mathbb{R}_{\ge} $, given by $ K(z) = \left(f-z\right) \log\left(\frac{f}{z}\right) $, is locally continuous and coercive, i.e.,

it follows from Eqs (4.16) and (4.17) that there exists a constant $ C^\ast(\zeta, K, m, |\Omega|, l_i, \delta) > 0 $ such that

Hence, the sequence $ \left\{(c_k, \varphi_k)\right\}_{k\in\mathbb{N}} $ is bounded in $ \mathbb{R}^{m+1}\times BV(\Omega) $. Then, from the compactness property in $ BV $-space and the fact that $ BV(\Omega) $ is compactly embedded in $ L^\delta(\Omega) $, we can deduce the existence of a subsequence of $ \left\{(c_k, \varphi_k)\right\}_{k\in\mathbb{N}} $, that we denote in the same way, and a pair $ (c^0, \varphi^0)\in \mathbb{R}^{m+1}\times BV(\Omega) $ such that

Without loss of generality, we can suppose that each of the functions $ \left\{\varphi_k\right\}_{k\in\mathbb{N}} $ and $ \varphi^0 $ is extended by zero outside of $ \Omega $, and $ \varphi_k(x)\rightarrow \varphi^0(x) $ for a.e. $ x\in \Omega_{2{{\varepsilon}}} $. Then, in view of the standard properties of mollifiers, we have

and

Then, property (3.4) implies that

Besides, in view of the properties (4.18)–(4.20) and the fact that $ \Psi\in C^1_{loc}(\mathbb{R}) $, we have the pointwise convergence

Since

and the sequences $ \left\{|l_0-\varphi_k|^\delta\right\}_{k\in\mathbb{N}} $ and $ \left\{|\varphi_k-l_{m+1}|^\delta\right\}_{k\in\mathbb{N}} $ are bounded in $ L^p(\Omega) $ with $ p = 2/\delta > 1 $ (by the continuous embedding $ BV(\Omega)\hookrightarrow L^2(\Omega) $), it follows from Vitali's lemma that

Hence,

As a result, the lower semicontinuity property of the minimizing sequence

is a direct consequence of relations (4.21)–(4.23) and Lemma 9.

It remains to notice that due to the pointwise convergence (4.18), we have

Hence, the limit pair $ (c^0, \varphi^0) $ is a feasible solution, i.e., $ (c^0, \varphi^0)\in \Xi_{{\varepsilon}} $, and, therefore,

Thus, $ (c^0, \varphi^0) $ is a minimizer to the problems (4.9)–(4.13).

5.   Optimality conditions for approximated problem

This section aims to derive some optimality conditions for the minimization problems (4.9)–(4.13). Let $ {{\varepsilon}} $ and $ \tau $ be given small positive values. With that in mind, we study the differentiability properties of the objective functional $ J_{{{\varepsilon}}, \tau}(c, \varphi) $ in order to specify its local behavior in the immediate vicinity of its minimum point. The corresponding Euler-Lagrange system is presented in Theorem 11.

Let $ (c^0_{{{\varepsilon}}, \tau}, \varphi^0_{{{\varepsilon}}, \tau})\in \Xi_{{\varepsilon}} $ be a local minimizer to problems (4.9)–(4.13). Then,

with $ c^0_{{{\varepsilon}}, \tau, j}\ge 0 $ for each $ j\in\left\{0, \dots, m\right\} $ and $ \varphi^0_{{{\varepsilon}}, \tau}\in \mathcal{G}_{ad} $. In the objective functional $ J_{{{\varepsilon}}, \tau} $, we distinguish three terms

with

From the differentiability of $ \left(f-c_{i}\right)\log\left(\frac{f}{c_i}\right) $ and $ \left[\chi_{A_j^\tau}\right]_{{\varepsilon}} $, it is immediate that the functional $ F_{{{\varepsilon}}, \tau} $ is of the class $ C^1 $. Hence, there exist linear continuous functionals

such that

for any $ h\in BV(\Omega) $ and $ \mu\in\mathbb{R} $, where $ |r(h, \lambda)| = o(|\lambda|) $ and $ |r_j(\mu, \lambda)| = o(|\lambda|) $ as $ \lambda\to 0 $, and

Moreover, making use of the following representations

with $ H(z) = \left\{1, \ z\ge 0\atop 0, \ z < 0\right\} $, we have

where $ H_{{\varepsilon}}(z) = \left[\eta_{{\varepsilon}} \ast \chi_E\right](z) $ stands for the smooth approximation of $ H(z) $ through the mollification.

Then, direct calculations show that

with

Since $ f\in \mathcal{F}_\gamma $, it follows from Eq (5.2) that the unique solution of the system

can be expressed as

with $ c^0_{{{\varepsilon}}, \tau, j}\ge 0 $ for all $ j = 0, \dots, m $, i.e., $ c^0_{{{\varepsilon}}, \tau} \in \mathcal{C}_{ad} $.

Arguing in a similar manner and taking into account that $ \Psi\in C^1_{loc}(\mathbb{R}) $, it can be shown that

for any $ h\in BV(\Omega) $, where $ |r(h, \lambda)| = o(|\lambda|) $ as $ \lambda\to 0 $, and $ D_\varphi \Phi_{{{\varepsilon}}} (\varphi^0_{{{\varepsilon}}, \tau}): BV(\Omega)\rightarrow \mathbb{R} $ is a linear continuous functional with the following representation:

Here,

We are now in a position to establish the main result of this section.

Theorem 11. Given $ {{\varepsilon}} > 0 $ and $ \tau > 0 $ small enough, $ f\in \mathcal{F}_\gamma $, $ \alpha > 0 $, and a collection of $ m+2 $ distinct level values $ l_0 < l_1 < \dots < l_m < l_{m+1} $, let $ (c^0_{{{\varepsilon}}, \tau}, \varphi^0_{{{\varepsilon}}, \tau})\in \Xi_{{\varepsilon}} $ be a local minimizer to problems (4.9)–(4.13). Then, the pair $ (c^0_{{{\varepsilon}}, \tau}, \varphi^0_{{{\varepsilon}}, \tau}) $ satisfies the following Euler-Lagrange system:

where the constants $ c^0_{{{\varepsilon}}, \tau, 0}, c^0_{{{\varepsilon}}, \tau, 1}, \dots, c^0_{{{\varepsilon}}, \tau, m} $ are defined by the rule (5.4).

Proof. Since $ (c^0_{{{\varepsilon}}, \tau}, \varphi^0_{{{\varepsilon}}, \tau})\in \Xi_{{\varepsilon}} $ is a local minimum point of Eqs (4.9)–(4.13), we have that

for all $ \rho > 0 $ small enough and any given $ \mu\ge 0 $ and $ \varphi\in \mathcal{G}_{ad} $.

As a result, inequality (5.7) leads to the Eq (5.3), and, hence, to the representation (5.4), whereas Eq (5.8) together with the convexity of $ j_{{\varepsilon}} $ implies

Now, passing to the limit as $ \rho\to 0 $, we get

Finally, using the expression of the Gateaux derivatives $ D_\varphi F_{{{\varepsilon}}, \tau} (c^0_{{{\varepsilon}}, \tau}, \varphi^0_{{{\varepsilon}}, \tau}) $ and $ D_\varphi \Phi_{{{\varepsilon}}} (\varphi^0_{{{\varepsilon}}, \tau}) $ given by Eqs (5.1) and (5.5), respectively, we immediately arrive at the optimality system (5.6).

6.   Asymptotic analysis of the approximated problems as $ {{\varepsilon}}\to 0 $

The main question we are going to discuss in this section is to find out whether the convergence of minima of Eq (4.9) is to minima of Eq (4.6) as $ {{\varepsilon}} $ and $ \tau $ tend to zero. To this end, we make use of the basic results of the variational convergence of minimization problems and $ \Gamma $-convergence theory (see, for instance, ^[27,28,29]). In particular, in Lemmas 12 and 13, we show that the standard properties of $ \Gamma $-limits hold true for the objective functional $ J_{{{\varepsilon}}, \tau} $ with respect to the weak-$ \ast $ topology of $ BV(\Omega) $ space and the pointwise convergence in $ \mathbb{R}^{m+1} $. Utilizing these characteristic features, we establish the main variational property of the proposed approximation procedure (see Theorem 14). Namely, we prove that any sequence of optimal pairs to the approximated problems (4.9)–(4.13) is compact in the weak-$ \ast $-topology of $ \mathbb{R}^{m+1}\times BV(\Omega) $ and each cluster point is a solution of the problem

where the cost functional is defined in Eq (6.6).

We begin with the following noteworthy result.

Lemma 12. Let $ \tau > 0 $ be a given value such that $ \tau\ll 1 $. Let $ \left\{(c_{{\varepsilon}}, \varphi_{{\varepsilon}})\in\Xi_{{\varepsilon}}\right\}_{{{\varepsilon}}\to 0} $ be a sequence of feasible pairs to the approximated problems (4.9)–(4.13), satisfying the conditions

Then, there exist a subsequence $ \left\{(c_{{{\varepsilon}}_j}, \varphi_{{{\varepsilon}}_j})\right\}_{j = 1}^\infty $ with $ {{\varepsilon}}_j\to 0 $ as $ j\to\infty $, and a pair $ (c, \varphi)\in \mathbb{R}^{m+1}\times BV(\Omega) $ such that

where

Proof. In view of the initial assumptions (see Eq (6.1)$ _2 $), the sequence $ \left\{(c_{{\varepsilon}}, \varphi_{{\varepsilon}})\right\}_{{{\varepsilon}}\to 0} $ is compact with respect to the product of norm topology of $ \mathbb{R}^{m+1} $ and the weak-$ \ast $ convergence in $ BV(\Omega) $. So, there exists a subsequence $ \left\{(c_{{{\varepsilon}}_j}, \varphi_{{{\varepsilon}}_j})\right\}_{j = 1}^\infty $ with $ {{\varepsilon}}_j\to 0 $ as $ j\to\infty $ and a pair $ (c, \varphi)\in \mathbb{R}^{m+1}\times BV(\Omega) $ such that, in addition to Eqs (6.2) and (6.3), we have

Since the matrix $ M^f $ is positive-definite on $ \Omega $, property (6.7) implies the weak-$ \ast $ convergence (6.4). Moreover, from Eqs (6.1) and (6.7), we deduce the existence of a constant $ C > 0 $ independent of $ {{\varepsilon}} $ such that

Hence,

It means that the limit function $ \varphi $ satisfies the pointwise constraints $ l_0\le \varphi(x)\le l_{m+1} $ a.e. in $ \Omega $. Thus, we see that $ (c, \varphi)\in\Xi $.

It remains to establish the inequality (6.5)$ _2 $. With that in mind, we make use of the pointwise convergence Eqs (6.2), (6.7), and properties (ⅰ) and (ⅱ) of the smoothed characteristic functions $ \left[\chi_E\right]_{{\varepsilon}} $. Then, arguing as in the proof of Lemma 9, we deduce that

with some integrable function $ g $. Hence, by the Lebesgue dominated theorem, we have:

To end the proof, we have to show that

With that in mind, we notice that

Taking into account that (see Eq (6.1)$ _2 $), all level sets of the functions $ \left\{\varphi_{{{\varepsilon}}_j}\right\}_{j\in\mathbb{N}} $ have a finite perimeters, and we see that $ \chi_{\{\varphi_{{{\varepsilon}}_j}(x)\ge l_k-\tau\}}\in BV(\Omega) $ for each $ k = 1, \dots, m $. Moreover, in view of Eq (6.7), we have

for each $ k\in\left\{1, \dots, m\right\} $.

Hence,

Together with Eq (6.12), we arrive at the announced inequality (6.11). Thus, the desired property (6.5) is a direct consequence of of relations (6.8)–(6.11). □

Lemma 13. For every feasible pair $ (c, \varphi)\in\, \Xi $, there can be found a sequence $ \left\{(\widehat{c}_{{\varepsilon}}, \widehat{\varphi}_{{\varepsilon}})\right\}_{{{\varepsilon}} > 0} $ satisfying the properties

Proof. Let $ (c, \varphi)\in\, \Xi $ be a given pair. In view of Theorem 4, we can suppose that $ \varphi\in C^\infty(\Omega) $.

Let $ \left\{c_{{\varepsilon}} \right\}_{{{\varepsilon}} > 0} $ be an arbitrary sequence in $ \mathbb{R}^{m+1} $ such that

Then, we define the sequence $ \left\{(\widehat{c}_{{\varepsilon}}, \widehat{\varphi}_{{\varepsilon}})\right\}_{{{\varepsilon}} > 0} $ as follows:

Since $ \varphi\in C^\infty(\Omega) $, it follows that $ \widehat{\varphi}_{{\varepsilon}} $ has a bouded anisotropic total variation and, therefore, $ \left\{(\widehat{c}_{{\varepsilon}}, \widehat{\varphi}_{{\varepsilon}})\in \Xi_{{\varepsilon}}\right)\}_{{{\varepsilon}} > 0} $ are the collection of feasible solutions for the corresponding two-parametric approximated problems (4.9). Then, due to the fact that $ l_0\le \varphi(x)\le l_{m+1} $ in $ \Omega $, we have:

Furthermore, arguing as in the proof of Lemma 6 and taking into account the property (ⅰ) of mollifiers, we see that

Thus, in view of Eqs (6.19)–(6.21), in order to deduce the Eq (6.17), it remains to show that

Observing that (see Lemma 5)

where the last term in Eq (6.23) is equal to zero because $ \{x\in\Omega : \varphi(x)\ge l_j-\tau\} $ is a closed subset of $ \Omega $.

As a result, we obtain

This concludes the proof.

We are now in a position to state the main result of this section.

Theorem 14. Let $ \tau\ll 1 $ be a given positive value. Let $ \left\{(c^0_{{{\varepsilon}}, \tau}, \varphi^0_{{{\varepsilon}}, \tau})\in\Xi_{{\varepsilon}}\right\}_{{{\varepsilon}}\to 0} $ be a sequence of optimal pairs to the approximated minimization problems (4.9)–(4.13). Assume that the sequence $ \left\{\varphi^0_{{{\varepsilon}}, \tau}\right\}_{{{\varepsilon}} > 0} $ is bounded in $ BV(\Omega) $, and the relaxed problem

has a nonempty set of minimizers for the given value $ \tau > 0 $. Then, there exists a pair $ (c^\ast_\tau, \varphi^\ast_\tau)\in\Xi $ such that, up to a subsequence,

where the objective functional $ J_\tau:\Xi\to\mathbb{R} $ is defined in Eq (6.6).

Proof. First, we observe that a given sequence of minimizers for approximating problems (4.9)–(4.13) is compact with respect to the convergences (6.25)–(6.27). Indeed, for an arbitrary test function $ \widehat{\varphi}\in C^\infty_c(\mathbb{R}^2) $ and arbitrary vector $ \widehat{c}\in\mathbb{R}^{m+1} $ with positive components, we have:

Let's assume that, in addition, the function $ \widehat{\varphi} $ satisfies the pointwise constraints $ l_{0}\le\widehat{\varphi}(x)\le l_{m+1} $ in $ \Omega $. Then, $ (\widehat{c}, \widehat{\varphi})\in\Xi_{{{\varepsilon}}} $ for each $ {{\varepsilon}} > 0 $, and, therefore,

Hence,

Thus, for the sequence of minimizers $ \left\{\left(c^0_{{{\varepsilon}}, \tau}, \varphi^0_{{{\varepsilon}}, \tau}\right)\in \Xi_{{{\varepsilon}}}\right\}_{{{\varepsilon}} > 0} $, all preconditions of Lemma 12 are fulfilled. Therefore, there exist a subsequence $ \left\{(c^0_{{{\varepsilon}}_k, \tau}, \varphi^0_{{{\varepsilon}}_k, \tau})\in \Xi_{{{\varepsilon}}_k}\right\}_{k\in\mathbb{N}} $ of the sequence $ \left\{(c^0_{{{\varepsilon}}, \tau}, \varphi^0_{{{\varepsilon}}, \tau})\in\Xi_{{\varepsilon}}\right\}_{{{\varepsilon}}\to 0} $ and a pair $ (c^\ast_\tau, \varphi^\ast_\tau)\in\Xi $, such that $ (c^0_{{{\varepsilon}}_k, \tau}, \varphi^0_{{{\varepsilon}}_k, \tau})\rightarrow (c^\ast_\tau, \varphi^\ast_\tau) $ is the sense of convergences (6.25)–(6.27) and

From this, we deduce that

where $ (c^0_\tau, \varphi^0_\tau) $ is a minimizer for the relaxed problem (6.24).

On the other hand, Lemma 13 implies the existence of a realizing sequence $ \left\{(\widehat{c}_{{\varepsilon}}, \widehat{\varphi}_{{\varepsilon}})\right\}_{{{\varepsilon}} > 0} $ such that $ (\widehat{c}_{{\varepsilon}}, \widehat{\varphi}_{{\varepsilon}})\rightarrow (c^0_\tau, \varphi^0_\tau) $ as $ {{\varepsilon}}\to 0 $ in the sense of relations (6.16), and

Utilizing this fact, we get

From this and Eq (6.30), we deduce that

As a result, we have

Using these relations and the fact that the problem (6.24) is solvable, we may suppose that

Since Eq (6.32) holds for all subsequences of $ \left\{(c^0_{{{\varepsilon}}, \tau}, \varphi^0_{{{\varepsilon}}, \tau})\in\Xi_{{\varepsilon}}\right\}_{{{\varepsilon}}\to 0} $, which are convergent in the sense of relations (6.25)–(6.27), it follows that these limits coincide and, therefore, $ (c^0_\tau, \varphi^0_\tau) $ is the limit of the whole sequence $ \left\{(c^0_{{{\varepsilon}}, \tau}, \varphi^0_{{{\varepsilon}}, \tau})\right\}_{{{\varepsilon}} > 0} $. Then, using the same argument for the entire sequence of minimizers, we finally obtain

and this concludes the proof.

7.   Numerical results

To illustrate the implementation of the proposed optimization problem (4.6) to the domain decomposition that corresponds to the homogeneity zones of a given function $ f:\Omega\rightarrow \mathbb{R} $, we provided numerical experiences with images that have been delivered by satellite Sentinel-2. As input data, we have used an image over the Dnipro area, Ukraine, with a resolution of 10m/pixel (see the left panel in Figure 1). This region represents a typical agricultural area with medium-sized fields of various shapes. As follows from the picture given in Figure 1 (see also the corresponding histogram in Figure 2), the observed data suffer from noise and blurs. So, at the first step, we have realized the denoising and debluring procedure (see the right panel in Figure 1) following the variational approach that has been recently proposed in ^[30]. As Figure 2 indicates, the histogram of the smoothed image has a strongly marked compactly localized spectrum that can be considered as a "good option" for its piecewise constant approximation. To conduct the numerical simulations of the segmentation procedure for the given area, we have set $ f(x) = u_2(x) $ in $ \Omega $, where $ u_2 $ stands for the intensity of the de-blurred image (see Figure 1) in the green spectral channel, and

In accordance with the results of Section 5, we have to solve the system (5.6) and find its solution $ \varphi^0_e $ for the corresponding function $ f $. Since, in practical implementations, it is reasonable to define the solution of the problem (5.6) using a "gradient descent" strategy, we started with some initial level-set function $ \varphi_0\in C(\Omega) $ and passed to the corresponding initial-boundary value problem for quasi-linear parabolic equations with Neumann boundary conditions. For numerical simulations, we set $ {{\varepsilon}} = 0.01 $, $ \tau = 10 $, $ \alpha = 1 $, $ \sigma = 3 $, $ \eta = 0.95 $, and the initial level set function $ \varphi_0\in C(\Omega) $ was defined as follows:

where $ S $ is a circle of radius $ 20 $ with a center at a central point of $ \Omega $, and $ d(x, S) $ denotes the Euclidean distance from the point $ x\in\Omega $ to the circle $ S $. We report the level sets in Figures 3 and 4.

Author contributions

All the authors conceived the idea, designed the methodology and the main proofs. All the authors contributed equally in the writing of the article.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

Peter Kogut thanks the support of University of Salerno.

Conflict of interest

The authors declare there is no conflict of interest.