ON OPTIMAL 2 -D DOMAIN SEGMENTATION PROBLEM VIA PIECEWISE SMOOTH APPROXIMATION OF A SELECTIVE TARGET MAPPING

. In this paper we propose a new technique for the solution of the image segmentation problem which is based on the concept of a piecewise smooth approximation of some target functional. We discuss in details the consistency of the new statement of segmentation problem and its solvability. We focus our main intension on the rigor mathematical substantiation of the proposed approach, deriving the corresponding optimality conditions, and show that the new optimization problem is rather ﬂexible and powerful model to the study of variational image segmentation problems. We illustrate the accuracy and eﬃciency of the proposed algorithm by numerical experiences.


Introduction
In this paper we discuss a new coupled variational problem which is suggested by applications to satellite image segmentation. From practical point of view, image segmentation is the process of dividing an image into several areas with features and extracting the target of interest. in particular, in agricultural crop field classification, one of a fundamental problem is to provide a disjunctive decomposition of a fixed domain Ω ⊂ R 2 onto finite number of nonempty subsets Ω = Ω 1 ∪ Ω 2 ∪ · · · ∪ Ω K such that each of these subsets could be associated with a crop that is grown in this area, or with a forest regions, or water zones, and so on, and this correspondence must be established at rather high level of accuracy. Therefore, one important premise of object segmentation is to construct image objects that have homogeneous features [20,25,30]. In the agricultural applications, we consider the IPVI-characteristic as the main feature of such images. At the same time, a precise consideration of this problem demonstrates that the quantitative interpretation of remote sensing information from vegetation is a complex where B i : L 2 (Ω) → L 2 (Ω) is a linear continuous operator.
We Loosely speaking, this problem consists in finding a two-phase decomposition Ω = E 1 ∪ E 2 such that the distribution V (u opt ) ∈ L ∞ (Ω) varies slowly with respect to the L r (E i )-norm within each subset E i , and the distribution V (u opt ) ∈ L ∞ (Ω) varies rapidly (or discontinuously) across most of the boundary between E 1 and E 2 .
It is worth to note that minimization problems like (1.1)-(1.2) with an H −1constraint found growing interest in recent years due to several advantages to the L 2 -constrained problems. In particular, the benefit to involve H −1 -norm in image processing models is well-know (see [6] for a general overview on this topic). From practical point of view, the H −1 (Ω) space is intrinsically more appropriate for the modeling texture or oscillatory pattern and, in fact, it provides the norm which is smaller than the L 2 -norm. As for the nonlinear sequentially continuous operators V i : BV (Ω) → BV (Ω) that have been mentioned before, they can be chosen in a different way. However, for practical implementation of this step, we consider the mappings u → V i (u) as the resolvent operator of Alvarez, Lions, and Morel model [2] that takes the form of some initial-boundary value problem for parabolic equation with a nonlinear anisotropic diffusion operator in its principle part and with the function u in its initial condition.
The paper is organized as follows. In Section 2 we give some preliminaries related with the space of functions of bounded variation. The precise setting of coupled optimization problem and its previous analysis is given in Section 3. In particular, we define the concept of feasible solutions to that problem and show that, under some proper assumptions, this problem has a solution. The aim of Section 4 is to derive optimality conditions and provide their substantiation. In Section 5 we discuss the practical implementation of the considered problem to the satellite image processing related to the monitoring of crop fields. We show that this problem can be considered as a particular case of the proposed couple optimization problem. We also give the results of numerical simulation with the real-life satellite images which illustrate the accuracy and efficiency of the proposed algorithm.

Auxiliaries
Let Ω be a bounded open subset of R 2 with a Lipschitz boundary. For any subset E ⊂ Ω we denote by |E| its 2-dimensional Lebesgue measure L 2 (E). For a subset E ⊆ Ω let E denote its closure and ∂E its boundary. We define the characteristic function χ E of E by χ E (x) := 1, for x ∈ E, 0, otherwise.
For a function u we denote by u| E its restriction to the set E ⊆ Ω, and by u ∂E its trace on ∂E. Let C ∞ 0 (Ω) be the infinitely differentiable functions with compact support in Ω. The k-dimensional Hausdorff measure is denoted by H k , and µ E is the restriction of a measure µ to the set E. For a Banach space X its dual is X * and ·, · X * ,X is the duality form on X * × X. By and * we denote the weak and weak * convergence in normed spaces. Throughout the paper we will often use the concept of weak and strong convergence in L 1 (Ω). Let {f n } n∈N be a bounded sequence of functions in L 1 (Ω). We recall that {f n } n∈N is called equi-integrable on Ω, if for any δ > 0 there is a τ = τ (δ) such that´S |f n | dx < δ for every measurable subset S ⊂ Ω of Lebesgue measure |S| < τ . Then the following assertions are equivalent for L 1 (Ω)-bounded sequences (see, for instance, [5,16]): Theorem 2.1 ( [5]). If a bounded sequence {f k } k∈N ⊂ L 1 (Ω) is equi-integrable and f k → f almost everywhere in Ω, then f k → f strongly in L 1 (Ω).
We remind here the most common definitions of some functional spaces that we will use later on.
Then H 1 0 (Ω) is a Banach space and the norm in H 1 0 (Ω) can be defined by We denote the dual of H 1 0 (Ω) by H −1 (Ω). Then (see [26, p.401]), H −1 (Ω) is isometrically isomorphic to the Hilbert space of all distributions F ∈ D (Ω) satisfying Let's fix an arbitrary element u * of H −1 (Ω). Then there exists a vector- Therefore, it is clear now that On the other hand, due to the Lax-Milgram Lemma, the Dirichlet boundary value problem −∆y = u * in Ω, y = 0 on ∂Ω, has a unique solution y = (−∆) −1 u * ∈ H 1 0 (Ω) for each u * ∈ H −1 (Ω), where stands for the Laplace operator. Moreover, in view of the energy equalitŷ which holds true for the weak solution of Dirichlet problems (2.2), we can deduce the following a priori estimate Combining this result with (2.1), we obtain the following chain of inequalities for the dual norm · H −1 (Ω) in H −1 (Ω): Hence, the standard norm in H −1 (Ω) is equivalent to the following one (see also [27]) Remark 2.1. For the our further analysis, we make use of the following relation. Let y ∈ H 1 0 (Ω) be a weak solution to Dirichlet problem (2.2). Then

Functions of Bounded Variation
Let M(Ω; R 2 ) be the space of all R 2 -valued Borel measures which is, according to the Riesz theory, the dual of the space C 0 (Ω; R 2 ) of all continuous vectorfunctions ϕ with a compact support in Ω, equipped with the uniform norm By BV (Ω) we denote the space of all functions u ∈ L 1 (Ω) for which their distributional derivatives are representable by finite Borel measures in Ω, i.e.
for some R 2 -valued measure Du ∈ M 2 (Ω). It can be shown that BV (Ω), endowed with the norm stands for the total variation of u in Ω. It is clear that |Du|(Ω) =´Ω |∇u| dx if u is continuously differentiable in Ω.
Remark 2.2. In the similar manner, we can also define the space BV (Ω; R M ) as the space of all vector-valued functions u : According to the Radon-Nikodym theorem, if u ∈ BV (Ω) then there exists ∇u ∈ L 1 (Ω; R 2 ) and a measure D s u singular with respect to the 2-dimensional Lebesgue measure L 2 Ω restricted to Ω, such that Du = ∇uL 2 Ω + D s u.
We recall that a sequence {f k } ∞ k=1 converges weakly * to f in BV (Ω) if and only if the two following conditions hold (see [4, p.124]): f k → f strongly in L 1 (Ω) and Df k * Df weakly * in M(Ω; R 2 ), i.e.
The following embedding results for BV -function is very useful in connection with variational problem that we study in this paper (see [5, p.378]). are compact for all p such that 1 p < 2. Moreover, there exists a constant C em > 0 which depends only on Ω and p such that for all u in BV (Ω), We also recall the Poincare-Wirtinger inequality: in two dimensional case, there exists a constant C P W such that, for any u ∈ BV (Ω), we have denotes the mean of u in Ω.

Sets of Finite Perimeter
In this subsection, we remind the main properties of the so-called sets of finite perimeter introduced by R. Caccioppoli in [9]. Definition 2.1. Let E be an L 2 -measurable subset of R 2 with finite Lebesgue measure. Let χ E be its characteristic function. We say that E is a set with finite perimeter in Ω if χ E ∈ BV (Ω). This means that the distributional gradient Dχ E is a vector-valued measure with finite total variation. The total variation Dχ E |(Ω) is called the perimeter of E in Ω, i.e., P (E, Ω) = |Dχ E |(Ω) and, therefore, From the theory of functions of bounded variation it is known that for every set E with finite perimeter in Ω there exists a vector-valued Radon measure Dχ E such that a generalized Gauss-Green formula holds: where ν E is the inner unit normal to E and Dχ E = ν E |Dχ E | is called the polar decomposition of Dχ E .
Since the sets with finite perimeter are not smooth in general, the correct way to represent the measure Dχ E is to introduce the so-called reduced boundary ∂ * E. Definition 2.2. Let E be a set of finite perimeter (in R 2 ). We say that x ∈ ∂ * E if 1. for every r > 0 we have 0 < meas (E ∩ B r (x)) < meas (B r (x)); 2. there exists the limit and |ν E (x)| = 1.
where the set ∂ * E is called the reduced boundary of E.
In this way, for every set E of finite perimeter we have The following properties are well-known: admits a subsequence E k(i) i∈N converging in measure in Ω to some E ⊂ Ω, i.e., |E k(i) ∆E| → 0 as i → ∞; (c) E → P (E, Ω) is lower semicontinuous with respect to the convergence in measure in Ω; (d) For any two sets E 1 and E 2 with finite perimeters in Ω, the relation holds, with equality holding if only the distance between these sets in Euclidean space R 2 is non-zero. domain Ω ⊂ R 2 , it is plausible to suppose that Ω is the set with finite perimeter and such that ∂Ω = ∂ * Ω.
We recall here some auxiliary results concerning the vector fields z ∈ L ∞ (Ω; R 2 ) whose divergence in the sense of distribution is an L 2 (Ω)-function. We denote by L ∞ 2,div (Ω; R 2 ) the space of all such vector-valued fields z. Then the trace of the normal component of the vector field z ∈ L ∞ 2,div (Ω; R 2 ) on ∂Ω can be defined as distribution Tr(z, ∂Ω) in the sense of Anzellotti (see [3]). In particular, having assumed that the original domain Ω ⊂ R 2 is of class C 1 , the trace of the normal component of z on ∂Ω is the distribution defined as follows For example, if z is a piecewise continuous vector field that can be extended continuously in Ω, then (see [13, p. 22 Utilizing the property ∂Ω = ∂ * Ω, we have the following result (the so-called Gauss-Green formula) (see [13,Theorem 5.1] and [21, Proposition 6.12] for the details).
Lemma 2.1. For any u ∈ BV (Ω) and z ∈ L ∞ 2,div (Ω; R 2 ) there exists a Radone measure on Ω denoted by (z, Du) and a function Tr(z, ∂Ω) ∈ L ∞ (∂Ω) such that The measure (z, Du) and | (z, Du) | are absolutely continuous with respect to |Du| and, for any open Ω ⊂ Ω, ∀ ϕ ∈ C ∞ 0 ( Ω), and for all Borel sets Ω ⊂ Ω, we have Moreover, it turns out that in this case the following estimate holds true

Setting of the Coupled Optimization Problem and Its Previous Analysis
We begin with the following assumptions: For the collection of energy functionals . . . . . .
with r 1 and we consider the following coupled variational problem Here, , α i 0, (3.11) and the set of feasible solutions Ξ is given by the rule Remark 3.1. If E opt ⊆ Ω minimizes the cost functional J 0 , then adding to E opt any L 2 -negligible piece, we formally obtain a new optimal solution E ⊆ Ω because such perturbation of E opt inside Ω do not decrease the cost functional J 0 . For this reason, in order to avoid any ambiguity in the sequel, it makes sense to introduce the concept of essential optimal solution [4, p. 319]: it has the property that L 2 (E opt ∩ B(x)) > 0 for any x ∈ E opt and any ball B(x) ⊂ Ω. As a result, it can be shown that any optimal solution E opt induces an essential optimal one E , with E ⊂ E opt .
Let us show that the above mentioned variational problem is consistent, that is, there exists a nonempty Borel subset of E of Ω and a vector-valued function With that in mind we begin with the following result (for comparison, we refer to [1, 11,22,24,28]). Proof. Fixing i ∈ {1, . . . , M − 1}, let us show that the corresponding minimization problem (3.9)-(3.10) is consistent, in particular, J i (u) < +∞ for any u ∈ Ξ. Indeed, since ∂Ω is Lipschitz and we deal with two-dimensional domain Ω, it follows from Poincare inequality (2.9) that BV (Ω) is continuously embedded into L 2 (Ω). In fact, in this particular case we can assert that the injection BV (Ω) → L q (Ω) is compact for all q ∈ [1, 2). On the other hand, by Sobolev embedding theorem we have a continuous injection H 1 0 (Ω) → L 2 (Ω). Hence, by the duality arguments and Riesz representation theorem, we deduce that finalize the remark about consistency, it is enough to observe that B i (u) ∈ L 2 (Ω) for all u ∈ BV (Ω). Hence, J(u) < +∞ provided u ∈ Ξ. Further, we note that the cost functional in (3.9) is non-negative. Hence, the infimum in (3.9)-(3.10) is finite. Let {u i,k } k∈N be a minimizing sequence for (3.14) Now, we prove that there exists a positive constant M such that Then, in view of (2.9), we deduce that On the other hand, by estimate (3.14) we see that Hence, and Then, taking into account the definition of w i,k and assumption (ii), we have Hence, due to condition (3.1), we obtain (3.20) Combining this fact with estimates (3.12)-(3.13), we finally deduce that the minimizing sequence {u i,k } k∈N is bounded in BV (Ω) and H −1 (Ω). Therefore (see [18,Theorem 1.19]), there exists a subsequence, still denoted by the same index, and elements u 0 in Ω by uniqueness of the weak limit in reflexive Banach space. Utilizing these properties together with inequality (2.8) and lower semicontinuity of L 2 -norm with respect to the weak convergence, we get As a result, it follows from the above consideration that Since the set of feasible solutions Ξ is convex and closed in BV (Ω), by Mazur's theorem we have that this set is sequentially closed with respect to the weak * convergence in BV (Ω). Hence, u 0 i 0 almost everywhere in Ω. Thus, u 0 i ∈ Ξ and u 0 i is a minimizer for constrained minimization problem (3.9)-(3.10). It remains to show that u 0 i is a unique minimizer for this problem. Indeed, let us assume the converse. Let u * i ∈ Ξ and v * i ∈ Ξ be two minimizers for the problem (3.9)-(3.10). By the strict convexity of norms · L 2 (Ω) and · H −1 (Ω) (see (2.6)), we obviously have which brings us into conflict with the initial assumptions. Hence, , as a consequence, we obtain Thus, u 0 i is a unique minimizer to the problem (3.9)-(3.10). The proof is complete.
Our next step is to study variational problem (3.7). To begin with, let us show that this problem is consistent. It would be rather provocative to assert that any measurable subset E of Ω can be considered as a feasible solution to the problem (3.7). Therefore, it is reasonable to define the corresponding set of feasible solutions Ξ 0 as follows: E ∈ Ξ 0 if and only if, for a given distribution u opt ∈ BV (Ω; R M ), we have E ⊂ Ω, E is a measurable set, and J 0 (E) < +∞. having set Ω i,n = x ∈ Ω : |u opt i (x)| n , we see that nχ Ω i,n (x) n |u opt i (x)|, ∀ x ∈ Ω i,n . Hence, using the monotonicity of Lebesgue integral, we have and letting n → ∞, it follows that L 2 (Ω i,∞ ) = 0, that is, each of functions u opt i ∈ BV (Ω) is finite L 2 -almost everywhere in Ω. As a result, we have Then, we immediately deduce from (3.11) that V (u opt ) ∈ L ∞ (Ω), and therefore, the last two terms in the energy functional J 0 (see (3.4)) are well defined for any L 2 -measurable subset E ⊂ Ω. It remains to notice that the inclusion χ E ∈ BV (Ω) is equivalent to the condition |Dχ E |(Ω) < +∞. Since P (E, Ω) = |Dχ E |(Ω), it follows that J 0 (E) < +∞ if and only if E is a subset of Ω and it has a finite perimeter in Ω. In view of these observations, the set of feasible solutions Ξ 0 to the problem (3.7) can be redefined as follows To clarify this option, we give the following observation. , let E opt ∈ Ξ 0 be an essential optimal solution of minimization problem (3.7). Then the set E opt is unique in the following sense: there does not exist a subset E 0 ∈ Ξ 0 such that E 0 ∈ Ξ 0 is an essential optimal solution, Proof. Let us assume the converse. Let E 0 ∈ Ξ 0 be a subset of Ω such that E 0 is an essential optimal solution of (3.7) and there exists a subset E * ⊆ X opt ∩ E 0 such that L 2 (E * ) > 0. Let χ E * be the corresponding characteristic function. Then it is clear that Then, taking into account representation (3.4) and the well-know inequality we deduce from (2.7) the following relation which leads us into conflict with optimality of E opt and E 0 .
Theorem 3.1. Let u opt = (u opt 1 , . . . , u opt M ) be a vector-valued minimizer to the constrained minimization problems (3.9)-(3.10). Then there exists a subset E opt ⊂ Ω, with a finite perimeter in Ω, such that Proof. Since 0 J 0 (E) < +∞ for all E ∈ Ξ 0 , it follows that there exists a non-negative value µ 0 such that µ = inf minimizing sequence to the problem (3.7), i.e.
So, we can suppose that J 0 (E k ) µ + 1 for all k ∈ N. Then Therefore, {χ E k } k∈N is a bounded sequence in BV (Ω) and, without loss of generality, we can suppose that there exists a subsequence of {χ E k } k∈N (still denoted by the same index) and a function χ ∈ BV (Ω) ∩ L ∞ (Ω) such that (3.25) Let us show that, in fact, there exists a subset E opt ∈ Ξ 0 such that χ is the characteristic function of E opt , i.e. χ = χ E opt . Indeed, the element χ ∈ BV (Ω) is the characteristic function of a domain if the identity χ(x) (1 − χ(x)) = 0 holds true almost everywhere in Ω. Since the equalities 3.26) are satisfied for all k ∈ N, then passing to the limit in both sides of (3.26) as k → ∞, we deduce from (j)-property that the same identity holds for the limit function χ ∈ L ∞ (Ω). Hence, we can suppose that χ = χ E opt for some measurable subset E opt ⊆ Ω. In order to establish the inclusion E opt ∈ Ξ 0 , it remains to utilize estimate (3.24) together with the lower semi-continuity property (jj). As a result, we obtain: and, therefore, E opt ∈ Ξ 0 .
Thus, in view of compactness properties of minimizing sequence {χ E k } k∈N with respect to the weak * convergence in BV (Ω), we arrive at the following obvious relation Thus, E opt ⊆ Ω is a solution of variational problem (3.7).

Optimality Conditions
In this section, following in general aspects the technique of Temam for the problem of minimal surfaces [14] and duality results from [15], we derive necessary optimality conditions in order to characterize the solution

Optimality Conditions for Constrained Minimization Problem (3.9)-(3.10)
We assume that, for given i ∈ {1, . . . , M − 1}, λ 1 0, λ 2 > 0, u i,d ∈ L 2 (Ω), and v i,d ∈ L 2 (Ω), condition (ii) holds true. For our further analysis, we reformulate problem (3.7)-(3.11) as an equivalent problem on the space X := L 2 (Ω). For this reason we define the following functionals on X: (4.1) Notice that all indicated functional are well-defined on X. In view of this, we extend the energy functional J i to the entire set L 2 (Ω) by the rule Then it is clear that a minimizer u opt i ∈ Ξ of (3.7)-(3.11) is also a minimizer of the modified problem because the inclusion u ∈ Ξ is equivalent to the consistency condition of the problem (4.6), i.e.
A fundamental difficulty that typically appears in deriving of necessary conditions for minimizers of a constrained minimization problem (4.6), is the lack in differentiability of the energy functional (4.5) and the functional G(u) in the side condition. Since the functionals J i : BV (Ω) → R and G : L 2 (Ω) → R are convex, a necessary condition for a minimizer u opt i ∈ Ξ should employ the convex subdifferentials ∂ J i (u opt i ) and ∂G(u opt i ) which are some subsets of the dual spaces [BV (Ω)] * and L 2 (Ω) * = L 2 (Ω), respectively. In spite of the fact that almost nothing is known about the structure of [BV (Ω)] * (see [4] for the details), we have the continuous embedding BV (Ω) → L 2 (Ω). Therefore, the elements of ∂ J i (u opt i ) can be evaluated in L 2 (Ω) because, by the duality arguments, we have We begin with the following technical results.

Proposition 4.2.
Let v i,d ∈ L 2 (Ω) be a given distribution. Then the functional E 3 : L 2 (Ω) → R is convex and continuously differentiable on L 2 (Ω) with where B * i ∈ L(X, X) stands for the adjoint operator.
Proof. Fixing arbitrary elements u ∈ L 2 (Ω) and h ∈ L 2 (Ω), we get which obviously leads us to the desired conclusion.
Before proceeding further, we recall the definition of the subdifferential ∂F (u) of a convex proper functional F : X → R ∪ +∞ at some element u ∈ X. Setting X = L 2 (Ω), we see that X * = X. Then, for given u ∈ X, an element ξ ∈ X belongs to ∂F (u) if and only if, ∀ v ∈ X, Thus, ξ ∈ ∂F (u) if u is a minimizer on X of the following variational problem (4.9) Proposition 4.3. The functional G : L 2 (Ω) → R is convex, lower semicontinuous with respect to the weak convergence in L 2 (Ω), and Gâteaux differentiable on Proof. Since convexity and the lower semicontinuity property of G : L 2 (Ω) → R are the direct consequence of (4.4), it follows that this functional has a nonempty subdifferential at each point u ∈ L 2 (Ω). Let us assume that u is a nonzero element of L 2 (Ω). Let ξ ∈ ∂G(u). Then, by definition of subgradient, we have Setting v = u + λz for any z ∈ L 2 (Ω) and λ > 0, we obtain Hence,ˆΩ and by choosing z = − [4(u − |u|) − ξ], we deduce that It remains to consider the case u ≡ 0. Letting ξ ∈ ∂G(0), we obtain or in other termsˆΩ Since the fulfillment of this relation for all v ∈ L 2 (Ω) is possible if only ξ = 0 in L 2 (Ω), it follows from (4.11) that ∂G(u) is a singleton for all u ∈ L 2 (Ω). Hence, Our next intension is to define the structure of subdifferential for the functional E 1 : L 2 (Ω) → R ∪ +∞ given by the rule (4.1).
Remark 4.1. As immediately follows from (4.13), a vector field z : Ω → [−1, 1] can be formally identified with the quotient Du |Du| provided |Du| is nonzero and well defined at a given point x ∈ Ω. However, due to the Azellotti's theory of pairing, the correct interpretation of the quotient Du |Du| can be made through the equality (z, Du) R 2 = |Du|, where the field z ∈ L ∞ 2,div (Ω; R 2 ) is such that z L ∞ (Ω;R 2 ) 1 (see [13] for the details).
Proof. We set M * is the set of functions ξ ∈ L 2 (Ω) such that ξ = −div z for some z ∈ L ∞ 2,div (Ω; R 2 ), and z L ∞ (Ω;R 2 ) 1 with Tr(z, ∂Ω) = 0. In view of definition of the trace of the normal component of z on ∂Ω, this set is nonempty and well defined. Let us show that M * is closed with respect to the norm convergence in L 2 (Ω). Let {ξ k } k∈N ⊂ M * be a sequence such that ξ k → ξ in L 2 (Ω). Let z k ∈ L ∞ (Ω; R 2 ) k∈N by their prototypes, i.e., Then, for all test functions ϕ ∈ C ∞ 0 (R 2 ) and each k ∈ N, we have the equalitieŝ Since the sequence {z k } k∈N is bounded in L ∞ (Ω; N 2 ), by Banach-Alaoglu theorem, there exists a vector field z ∈ L ∞ (Ω; N 2 ) such that, up to a subsequence, we have the convergence z k * z in L ∞ (Ω; N 2 ) as k → ∞. Taking this fact into account and passing to the limit in (4.15) (for the chosen subsequence), we get i.e., in view of Lemma 2.1, we can deduce that ξ = −div z ∈ L 2 (Ω) and Tr(z, ∂Ω), ϕ = 0 . Making use of the lower semicontinuity of L ∞ -norm with respect to the weak * convergence z L ∞ (Ω;R 2 ) lim inf k→∞ z k L ∞ (Ω;R 2 ) 1, we finally obtain ξ ∈ M * . Thus, M * is a closed subset of L 2 (Ω).
Let I M * be the indicator function of M * . Then the Young-Fenchel transform of I M * is the function on L 2 (Ω) defined by Fixing arbitrary ξ ∈ M * and v ∈ L 2 (Ω), we see that Hence, On the other hand, by definition of the total variation of v in Ω, we get and As a result, we deduce from (4.17) that Taking into account the fact that M * is a closed and convex subset of L 2 (Ω), it follows that its indicator function I M * is convex and lower semicontinuous. Therefore, by the well known results of convex analysis (see [15,Proposition 3.1]), we have Thus, from this we immediately conclude that It remains to notice that the last equality together with estimate (4.16) implies the relation We are now in a position to derive the necessary conditions for a unique minimizer u opt i ∈ Ξ ⊂ BV (Ω) of constrained minimization problem (3.9)-(3.10). Since u opt i is also a minimizer of the modified problem (4.6) and the functionals J i : L 2 (Ω) → R and G : L 2 (Ω) → R are convex, a necessary conditions should employ the convex subdifferentials ∂ J i u opt i and ∂G u opt i . As a result, utilizing Proposition 6.3 in [21], we arrive at the following result.
the derivation of a non-smooth Lagrange multiplier rule for problem (4.6) can be done as in [21,Proposition 6.3]. As a result, using the sum rule for subdifferential of convex functionals saying that we deduce the existence of elements

Optimality Conditions for Segmentation Problem (3.7)
In order to derive the necessary conditions of optimality in optimal segmentation problem (3.7), we make use of the so-called boundary variation method (for the details of this approach, we refer to [8, Section 1.2]). Let E opt ⊆ Ω be an essential optimal solution to variational problem (3.7). Then and χ E opt ∈ SBV (Ω). Without loss of generality, we assume that the set E opt ⊆ Ω has a Lipschitz boundary and, therefore, we can suppose that and where ν E opt (x) stands for the exterior unit normal vector to E opt at x ∈ ∂ * E opt .
Following the boundary variation method, we assume that E opt is an open subset of Ω and its boundary ∂E opt is regular in the following sense: there exists a δ > 0, a finite number of points {ξ k } m k=1 ⊂ ∂E opt , and a collection of C 0,1 - Actually, as it is mentioned in [8], the regularity of ∂E opt does not need to be assumed as a hypothesis but it is a consequence of some suitable conditions on the class of feasible solutions to the problem (3.7). Since this is a quite delicate matter which goes under the name of regularity theory, we do not drill down into this field and we refer the interesting reader to the various books available in the literature (see, for instance, references [18,23]). Thus, in view of the above given assumptions, for H 1 -almost all x 0 ∈ ∂E opt the boundary ∂E opt near a given point x 0 ∈ ∂E opt can be written as the graph of some function φ(x), where x varies in an open subset ω of R. Namely, for a fixed x 0 ∈ ∂E opt , there exists k ∈ {1, . . . , m} such that representations (4.29)-(4.31) remain valid for k = k. So, we can suppose that ω = x 1 ∈ R : |x 1 − ξ k 1 | < δ . Since, without loss of generality, we can suppose that it follows that the corresponding part of the energy functional J 0 (in the region of x 0 ∈ ∂E opt ) can be represented in the Cartesian form as follows where the symbol D stands for the differentiation operator d/dx 1 . We emphasize that because of the Lipschitz property of functions φ k : R → R, they are differentiable for almost all x 1 ∈ ω.
Following the boundary variation method, we perturb of ∂E opt , and hence φ k (x 1 ), by taking a comparison function of the form φ k (x 1 ) + εψ(x 1 ), where ε > 0 is a small parameter, and ψ : ω → R is a smooth function with support in ω.
Since E opt is a minimizer for the problem (3.7), this leads us to the inequality (4.32) For our further analysis, we make use of the following representationŝ In view of the initial assumptions, we see that Hence, by Lebesgue differentiation theorem, almost every point x ∈ Ω is a Lebesgue point of V (u opt ) − c 1 r and V (u opt ) − c 2 r . Therefore, setting we deduce from (4.32) that Integrating by parts and taking into account that ψ ∈ C ∞ 0 (ω), we obtain Since ψ ∈ C ∞ 0 (ω) is an arbitrary function, we finally deduce that φ k ∈ C 0,1 (ω) must satisfy the following differential equation where the term represent the mean curvature of ∂E opt written in Cartesian coordinates. Thus, we have arrived at the following necessary conditions of optimality for a regular solution E opt ⊂ Ω of segmentation problem (3.7): the mean curvature of ∂E opt ∩ Π 2 (ξ k ) is locally equal to for each k = 1, . . . , m.

Application to the Satellite Image Processing Problem
It is well known that the satellite images can be considered as a source of knowledge used for monitoring and evaluating the Earth's vegetative cover in order to estimate the amount of vegetation, its biomass, percent of coverage, and provide the distinctions between bare soils and vegetation zone. One of the standard ways enabling to get such information is determination and timedependent analysis of the so-called vegetation indices. The most commonly used vegetation index (VI) is the so-called slope-based Infrared Percentage Vegetation Index (IPVI), ∈ Ω are functions of two variables representing the intensity of red (Red) and near-infrared (N IR) reflectance of some region Ω of R 2 , respectively.
Let Ω = (a 1 , b 1 ) × (a 2 , b 2 ) be a rectangle domain in R 2 . We associate with this domain the mapping as a two-bands satellite image representing the intensity of red reflectance (Red) u 1,d and near-infrared reflectance (N IR) u 2,d of the region Ω, respectively. The problem, which is suggested by application to remote sensing satellite image processing, consists in computing a decomposition Ω = Ω 1 ∪ Ω 2 ∪ · · · ∪ Ω K such that the IPVI-characteristic (5.1) varies slowly within each Ω j , and this characteristic varies discontinuously and/or rapidly across most of the boundary between different Ω j . Let us show that this problem can be considered as a particular case of the coupled optimization problem (3.2)-(3.7).
With that in mind we notice that the observed imaging data u i,d , i = 1, 2, typically suffer from noise and blurs. So, it is reasonably to consider the deblurring and de-noising problem as the first step of preprocessing for the mapping F : Ω → R 2 in order to recover the original high resolution image which is unknown a priori but it must be, in some sense, as close as possible to the observed data u i,d , i = 1, 2. The core idea of our approach to the de-noising and de-blurring problem of two-band satellite images is to use the total variation as a regularizing term in the following variational problems (Ω) stand for the denoising and de-blurring terms, respectively.
As follows from Proposition 3.1, for given λ 1 0, λ 2 > 0, and u i,d ∈ L 2 (Ω), i = 1, 2, the constrained minimization problems (5.3) admit unique solutions u opt 1 and u opt 2 , respectively. Further, we define the so-called selective smoothing operators {V i : BV (Ω) → BV (Ω)} 2 i=1 as follows: for the de-blurred and de-noised images u opt i ∈ BV (Ω), i = 1, 2, we set up V i (u opt i ) as a steady-state solution of the following initial-boundary value problem with anisotropic diffusion operator Here, G i is the Gauss smoothing kernel, T (u) :=´Ω T (x − y)u(y) dy ∈ W 1,∞ (Ω) is the Steklov smoothing operator with positive compactly supported smooth function T : R 2 → R such that T (x) dx = 1, and T (x) = T (−x), and g i (·) is the edge stopping function that we define as follows g i (s) = 1 1 + k i s 2 . (5.8) It should be noted that the system (5.5)-(5.7) belongs to the class of well-posed problems (existence and uniqueness were proven in [10], for the details we refer to [2, Section 3]). As a result, it can be showed that the above mentioned operators V i are sequentially continuous. Thus, the optimal satellite image segmentation problem via piecewise smooth approximation of Infrared Percentage Vegetation Index (5.1) can be stated as follows: Find a nonempty subset E opt ⊂ Ω such that: E opt ∈ Ξ 0 = {E ⊆ Ω : P (E, Ω) < +∞} , L 2 (E opt ∩ B(x)) > 0 for any x ∈ E opt and any ball B(x) ⊂ Ω, and where It remains to notice that, in view of Remark 3.2 and Theorem 3.1, this segmentation problem admits at least one optimal partition Ω = E opt ∪ Ω \ E opt that can be characterized by necessary optimality conditions (4.33).

Results of Numerical Simulation
We illustrate here the accuracy and efficiency of the proposed algorithm by numerical experiences with images that have been delivered by satellite Sentinel-2. From technical point of view, Sentinel-2 carries a multispectral imager with a swath of 290 km, and delivers high-resolution optical images for land monitoring, emergency response and security services. The imager provides a versatile set of 13 spectral bands spanning from the visible and near infrared to the shortwave infrared, featuring four spectral bands at 10 m, six bands at 20 m and three bands at 60 m spatial resolution.
In what follows, we associate with the original image (see Figure 5.1) the mapping (5.2), where the intensities of red and infrared reflectance u 1,d , u 2,d are presented in .
As follows from the pictures given in Figures 5.1-5.2 (see also the corresponding histogram of IPVI-characteristic for the original image), the observed data u i,d = u i,d (x 1 , x 2 ), i = 1, 2, suffer from noise and blurs.
In accordance with algorithm that we propose in Section 5, we consider the de-blurring and de-noising problem for u i,d = u i,d (x 1 , x 2 ), i = 1, 2, as the first step of the image preprocessing. To do so, we solve the variational problems (5.3) with λ 1 = 20 and λ 2 = 10. Figures 5.5 and 5.6 contain the corresponding solutions u opt 1 and u opt 2 of these problems. The next step is to solve the system (5.5)-(5.7) and find its steady-state solutions V i (u opt i ) for the corresponding deblurred images u opt i . Making use of the standard approach for the numerical simulation of the initial-boundary value problem (see, for instance, [32] for the details), we show that the chosen selective smoothing procedure transforms the original image and its IPVI-characteristic to their 'cartoon' versions as it is shown in Figures 5.7 and 5.8.
Results of optimal segmentation are presented on Figures 5.9-5.14.

Conclusion
In this paper, we have proposed a new setting for the optimal image segmentation problem which is based on the concept of piecewise smooth approximation of some selective target mappings. We have shown that the remote sensing satellite image segmentation problem, based on the analysis of the slope-based vegetation indices, is a particular case of the proposed setting. Focusing mainly on the rigor mathematical substantiation of the proposed approach, we discuss in details the  consistency of the new statement of segmentation problem and its solvability. We derive the corresponding optimality conditions and provide their substantiation. We show that the proposed coupled optimization problem (3.7)-(3.11) is rather flexible and powerful model to the study of variational image segmentation. Mostly motivated by the crop field classification problem and the auto-  Result of optimal segmentation of domain E opt as a solution of piecewise constant approximation problem for the smoothed IPVI-characteristic restricted to E opt . mated computational methodology for the extraction of agricultural crop fields from satellite data, we provide a numerical simulation with a Sentinel-2 remote sensing image. The obtained results confirmed the effectiveness of the proposed algorithm. In particular, the experimental results indicate that IPVI-histograms Result of optimal segmentation of domain Ω \ E opt as a solution of piecewise constant approximation problem for the smoothed IPVI-characteristic restricted to Ω \ E opt .
of the extracted agricultural crop fields, that have been obtained as a result optimal segmentation of the test image, have small and compactly disjointed supports in a given spectrum range. As it has been indicated in many recent publications (see [7,12,29,31]), this fact is crucial for the segmentation of blur and noise corrupted images with heterogeneity of their spectral and shape features. Additionally, the sensitivity analysis show that the evaluation of the smoothed versions of IPVI have a stable relationship with respect to the segmentation results.