Compensated Convexity on Bounded Domains, Mixed Moreau Envelopes and Computational Methods

TL;DR

This paper introduces local compensated convex transforms for functions on bounded domains, relates them to global transforms via mixed Moreau envelopes, and develops numerical methods with convergence analysis for applications in image processing and shape interrogation.

Contribution

It extends compensated convex transforms to bounded domains using specific extensions and relates them to global transforms through mixed Moreau envelopes, along with proposing convergent numerical schemes.

Findings

Local transforms coincide with global transforms on the domain boundary.

Numerical schemes for Moreau envelopes converge with rates depending on function regularity.

Applications demonstrate effectiveness in image processing and shape analysis.

Abstract

Compensated convex transforms have been introduced for extended real-valued functions defined over $\mathbb{R}^n$. In their application to image processing, interpolation, and shape interrogation, where one deals with functions defined over a bounded domain, one was making the implicit assumption that the function coincides with its transform at the boundary of the data domain. In this paper, we introduce local compensated convex transforms for functions defined in bounded open convex subsets $\Omega$ of $\mathbb{R}^n$ by making specific extensions of the function to the whole space, and establish their relations to globally defined compensated convex transforms via the mixed critical Moreau envelopes. We find that the compensated convex transforms of such extensions coincide with the local compensated convex transforms in the closure of $\Omega$. We also propose a numerical scheme for computing Moreau envelopes, establishing convergence of the scheme with the rate of convergence depending on the regularity of the original function. We give an estimate of the number of iterations needed for computing the discrete Moreau envelope. We then apply the local compensated convex transforms to image processing and shape interrogation. Our results are compared with those obtained by using schemes based on computing the convex envelope from the original definition of compensated convex transforms.

Equations

$$C_{\lambda }^l(g)(x)=\mathrm{co}[g+\lambda |\cdot {|}^2](x)-\lambda |x{|}^2,$$

$$C_{\lambda }^u(g)(x)=\lambda |x{|}^2-\mathrm{co}[\lambda |\cdot {|}^2-g](x).$$

$$C_{\lambda }^l(g)(x)=M^{\lambda }(M_{\lambda }(g))(x)\text{and}{C}_{\lambda }^u(g)(x)=M_{\lambda }(M^{\lambda }(g))(x)\text{for all }x\in {\mathbb{R}}^n$$

$$\begin{array}{rl}& M_{\lambda }(g)(x)=\underset{y\in {\mathbb{R}}^n}{\inf }\{g(y)+\lambda |y-x{|}^2\},\\ \text{and}& M^{\lambda }(g)(x)=\underset{y\in {\mathbb{R}}^n}{\sup }\{g(y)-\lambda |y-x{|}^2\},\end{array}$$

$$g^{\infty}:\,x\in\mathbb{R}^{n}\mapsto g(x)=\left\{\begin{array}[]{ll}g(x),&x\in\Omega\\[6.45831pt] +\infty,&x\in\mathbb{R}^{n}\setminus\Omega\,,\end{array}\right.$$

$$g^{-\infty}:\,x\in\mathbb{R}^{n}\mapsto g(x)=\left\{\begin{array}[]{ll}g(x),&x\in\Omega\\[6.45831pt] -\infty,&x\in\mathbb{R}^{n}\setminus\Omega\,,\end{array}\right.$$

$$O_f:=M-m\ge 0.$$

$$\displaystyle f_{\overline{\Omega}}^{-}:\,x\in\overline{\Omega}\mapsto\,f_{\overline{\Omega}}^{-}(x)=\left\{\begin{array}[]{ll}\displaystyle f(x),&\displaystyle x\in\Omega,\\[6.45831pt] \displaystyle\inf_{\Omega}f,&\displaystyle x\in\partial\Omega\,,\end{array}\right.$$

$$\displaystyle f_{\overline{\Omega}}^{+}:\,x\in\overline{\Omega}\mapsto\,f_{\overline{\Omega}}^{+}(x)=\left\{\begin{array}[]{ll}\displaystyle f(x),&\displaystyle x\in\Omega,\\[6.45831pt] \displaystyle\sup_{\Omega}f,&\displaystyle x\in\partial\Omega\,,\end{array}\right.$$

$$\displaystyle f_{\mathbb{R}^{n}}^{-}:\,x\in\mathbb{R}^{n}\mapsto\,f_{\mathbb{R}^{n}}^{-}(x)=\left\{\begin{array}[]{ll}\displaystyle f(x),&\displaystyle x\in\Omega,\\[6.45831pt] \displaystyle\inf_{\Omega}f,&\displaystyle x\in\mathbb{R}^{n}\setminus\Omega\,,\end{array}\right.$$

$$\displaystyle f_{\mathbb{R}^{n}}^{+}:\,x\in\mathbb{R}^{n}\mapsto\,f_{\mathbb{R}^{n}}^{+}(x)=\left\{\begin{array}[]{ll}\displaystyle f(x),&\displaystyle x\in\Omega,\\[6.45831pt] \displaystyle\sup_{\Omega}f,&\displaystyle x\in\mathbb{R}^{n}\setminus\Omega\,.\end{array}\right.$$

$$\begin{array}{rl}{C}_{\lambda }^l(f_{{\mathbb{R}}^n}^{-})(x)& =\mathrm{co}[f_{{\mathbb{R}}^n}^{-}+\lambda |\cdot {|}^2](x)-\lambda |x{|}^2,\\ C_{\lambda }^u(f_{{\mathbb{R}}^n}^{+})(x)& =\lambda |x{|}^2-\mathrm{co}[\lambda |\cdot {|}^2-f_{{\mathbb{R}}^n}^{+}](x),\end{array}$$

$$C_{\lambda }^l(f_{{\mathbb{R}}^n}^{-})(x)=M^{\lambda }(M_{\lambda }(f_{{\mathbb{R}}^n}^{-}))(x),C_{\lambda }^u(f_{{\mathbb{R}}^n}^{+})(x)=M_{\lambda }(M^{\lambda }(f_{{\mathbb{R}}^n}^{+}))(x).$$

$$M_{\lambda ,\Omega }(g)(x)=\underset{y\in \overline{\Omega }}{\inf }\{g(y)+\lambda |y-x{|}^2\},$$

$$M_{\Omega }^{\lambda }(g)(x)=\underset{y\in \overline{\Omega }}{\sup }\{g(y)-\lambda |y-x{|}^2\},$$

$$C_{\lambda ,\Omega }^l(f_{\overline{\Omega }}^{-})(x):=M_{\Omega }^{\lambda }(M_{\lambda ,\Omega }(f_{\overline{\Omega }}^{-}))(x)\text{for }x\in \overline{\Omega },$$

$$C_{\lambda ,\Omega }^u(f_{\overline{\Omega }}^{+})(x):=M_{\lambda ,\Omega }(M_{\Omega }^{\lambda }(f_{\overline{\Omega }}^{+}))(x)\text{for }x\in \overline{\Omega }.$$

$$C_{\lambda ,\Omega }^l(f_{\overline{\Omega }}^{-})(x)=C_{\lambda }^l(f_{{\mathbb{R}}^n}^{-})(x)\text{and}{C}_{\lambda ,\Omega }^u(f_{\overline{\Omega }}^{+})(x)=C_{\lambda }^l(f_{{\mathbb{R}}^n}^{+})(x).$$

$$C_{\lambda ,\Omega }^l(f_{\overline{\Omega }}^{-})(x)\le C_{\lambda ,\Omega }^l(f)(x)\le C_{\lambda }^l(f^{\infty })(x)\le f(x)\le C_{\lambda }^u(f^{-\infty })(x)\le C_{\lambda ,\Omega }^u(f)(x)\le C_{\lambda ,\Omega }^u(f_{\overline{\Omega }}^{+})(x),$$

$$C_{\lambda ,\Omega }^l(f)(x):=M_{\Omega }^{\lambda }(M_{\lambda ,\Omega }(f))(x)\text{and}{C}_{\lambda ,\Omega }^u(f)(x):=M_{\lambda ,\Omega }(M_{\Omega }^{\lambda }(f))(x)$$

$$C_{\lambda ,\Omega }^l(f_{\overline{\Omega }}^{-})(x)=C_{\lambda }^l(f_{{\mathbb{R}}^n}^{-})(x)\text{and}{C}_{\lambda ,\Omega }^u(f_{\overline{\Omega }}^{+})(x)=C_{\lambda }^l(f_{{\mathbb{R}}^n}^{+})(x).$$

$$\begin{array}[]{ll}\displaystyle\operatorname{\mathsf{co}}\left[f\right](x_{0})=\underset{x_{1},\ldots,x_{n+1}}{\inf}\,\Bigg{\{}\sum^{n+1}_{i=1}\lambda_{i}f(x_{i}):&\displaystyle\sum^{n+1}_{i=1}\lambda_{i}=1,\;\sum^{n+1}_{i=1}\lambda_{i}x_{i}=x_{0},\\[6.45831pt] &\displaystyle\lambda_{i}\geq 0,\;x_{i}\in\mathbb{R}^{n}\Bigg{\}}\,.\end{array}$$

$$\mathrm{co}[f](x_0)=\sum _{k=1}^p{\lambda }_k^{\ast }f(y_k^{\ast }).$$

$$\mathrm{co}\left[f\right](x_0)=\sup \left\{\ell (x_0):\ell \text{affine}\text{and}\ell (y)\le f(y)\text{for all }y\in {\mathbb{R}}^n\right\}$$

$$\begin{array}[]{ll}\displaystyle\operatorname{\mathsf{co}}_{\bar{B}(x_{0};r)}\left[f\right](x_{0})=\underset{x_{1},\ldots,x_{n+1}}{\inf}\,\Bigg{\{}\sum^{n+1}_{i=1}\lambda_{i}f(x_{i}):&\displaystyle\sum^{n+1}_{i=1}\lambda_{i}=1,\;\sum^{n+1}_{i=1}\lambda_{i}x_{i}=x_{0},\\[6.45831pt] &\displaystyle\lambda_{i}\geq 0,\;|x_{i}-x_{0}|\leq r,\;x_{i}\in\mathbb{R}^{n}\Bigg{\}}.\end{array}$$

$${\mathrm{co}}_{\bar{B}(x_0;r)}[f](x_0)=\sum _{k=1}^p{\lambda }_{k}f(y_k),$$

$$C_{\lambda }^l(f)(x)\le f(x)\le C_{\lambda }^u(f)(x),x\in {\mathbb{R}}^n,$$

$$C_{\lambda }^l(f)(x)\le C_{\lambda }^l(g)(x)\text{and}{C}_{\lambda }^u(f)(x)\le C_{\lambda }^u(g)(x),x\in {\mathbb{R}}^n.$$

$$\begin{array}[]{l}\displaystyle C^{l}_{\lambda}(f)(x)=\operatorname{\mathsf{co}}\left[\lambda|(\cdot)-x_{0}|^{2}+f\right](x)-\lambda|x-x_{0}|^{2}\,,\\[6.45831pt] \displaystyle C^{u}_{\lambda}(f)(x)=\lambda|x-x_{0}|^{2}-\operatorname{\mathsf{co}}\left[\lambda|(\cdot)-x_{0}|^{2}-f\right](x)\end{array}$$

$$C_{\lambda }^l(f)(x_0)=\mathrm{co}[\lambda |(\cdot )-x_0{|}^2+f](x_0),C_{\lambda }^u(f)(x_0)=-\mathrm{co}[\lambda |(\cdot )-x_0{|}^2-f](x_0).$$

$$\overline{f}(x)=\underset{y\to x}{\mathrm{limsup}}f(y)\text{and}\underline{f}(x)=\underset{y\to x}{\mathrm{liminf}}f(y).$$

$$\begin{array}{rl}{C}_{\lambda }^l(f)(x_0)& ={\mathrm{co}}_{B(x_0;R_{f,\lambda })}[\underline{f}+\lambda |(\cdot )-x_0{|}^2](x_0)\\ C_{\lambda }^u(f)(x_0)& =-{\mathrm{co}}_{B(x_0;R_{f,\lambda })}[\lambda |(\cdot )-x_0{|}^2-\overline{f}](x_0)\end{array}$$

Full text

Keywords:Compensated convex transforms, Moreau envelopes, Computation of Moreau envelopes, Proximity hull, Scattered data approximation, Shape interrogation

2000 Mathematics Subjects Classification number: 90C25, 90C26, 49J52, 52A41, 65D17

Email: kewei.zhang@nottingham.ac.uk, aorlando@herrera.unt.edu.ar, e.c.m.crooks@swansea.ac.uk

1 Introduction

This paper contains two parts. In the first part, we ‘localize’ the global notions of compensated convex transforms defined over $\mathbb{R}^{n}$, which were first introduced in [1, 2], by defining such transforms over bounded convex closed domains in $\mathbb{R}^{n}$ so that their values in the domain agree with the globally defined transforms applied to some special extensions of the function to $\mathbb{R}^{n}$. The motivation for such local definitions is mainly from applications to digital images and data arrays, where we have to consider functions defined on a rectangular box. In the second part, we propose a new scheme for the computation of Moreau envelopes and we prove its convergence. The scheme can be used to compute the Moreau based definitions of the local compensated convex transforms which we introduce in this paper or it can be useful on its own for general applications of Moreau envelopes.

The theory of compensated convexity transforms has been applied to, for example, digital image processing and computational geometry. So far, applications of the theory include the design of multiscale, parametrized, geometric singularity extraction of ridges, valleys and edges from graphs of functions and from characteristic functions of closed sets in $\mathbb{R}^{n}$ [3]. Several robust methods have been developed to date, namely, for the extraction of the set of intersections between two or more smooth compact manifolds [4, 5], for the extraction of the multiscale medial axis from geometric objects [6] and for the interpolation and approximation of sampled functions [7]. By ‘robustness’ here we mean the Hausdorff stability, that is, the error between the outputs obtained for two data samples is controlled by the Hausdorff distance between the two sampled input data sets.

In the applications mentioned above, the data domains are usually represented by boxes in $\mathbb{R}^{n}$. The numerical schemes used in [3, 4, 6, 7, 5] for the evaluation of the compensated convex transforms relied on the availability of schemes for computing the convex envelope of a given function and on the implicit assumption that the transforms coincide with the function at the boundary of the data domain. In those works, we were applying such schemes mainly to demonstrate the numerical feasibility of the transforms rather than (i) designing efficient numerical schemes for their computation and/or (ii) analysing the effect of the boundary assumptions. In this paper we also address the more practical question of accurately and effectively computing the compensated convex transforms for functions defined on a bounded convex domain without using numerical schemes that compute the convex envelope. In order to do so, we will explore the alternative definitions of the compensated convex transforms based on the Moreau envelopes. One of the advantages of this approach is that when it comes to the numerical implementation by the scheme we propose, we obtain an estimate of the number of iterations which provides the exact discrete Moreau envelope. This is different from the application of iterative schemes to compute the convex envelope which can be shown to converge but for which no convergence rate is known to be available [8].

Before we introduce our local versions of compensated convex transforms on a bounded closed convex domain, we recall from [1] the notions of quadratic compensated convex convex transforms in $\mathbb{R}^{n}$. For a function $g:\mathbb{R}^{n}\to\mathbb{R}\cup\{\infty\}$ satisfying the growth condition $g(x)\geq-c_{0}-c_{1}|x|^{2}$, $x\in\mathbb{R}^{n}$, for some constants $c_{0}$, $c_{1}>0$ with $|x|=\sqrt{x\cdot x}$ the Euclidean norm of $x\in\mathbb{R}^{n}$ and $x\cdot y$ the inner product on $\mathbb{R}^{n}$ between $x\in\mathbb{R}^{n}$ and $y\in\mathbb{R}^{n}$, the quadratic lower compensated convex transform (lower transform for short) of $g$ for $\lambda>c_{1}$ is defined for $x\in\mathbb{R}^{n}$ by [1],

[TABLE]

where $\operatorname{\mathsf{co}}[h]$ is the convex envelope of the function $h:\mathbb{R}^{n}\rightarrow\mathbb{R}\cup\{\infty\}$ bounded below. Given $g:\mathbb{R}^{n}\rightarrow\mathbb{R}\cup\{-\infty\}$ such that $g(x)\leq c_{0}+c_{1}|x|^{2}$, $x\in\mathbb{R}^{n}$, the quadratic upper compensated convex transform (upper transform for short) of $g$ of module $\lambda>0$ is defined for $x\in\mathbb{R}^{n}$ by

[TABLE]

The lower and upper transform can be, in turn, characterized in terms of the critical mixed Moreau envelopes as follows (see Appendix A)

[TABLE]

where, in our notation,

[TABLE]

are the lower and upper Moreau envelope of $g$, respectively [9, 10, 11, 12, 13], defined as the inf- and *sup-*convolution of $g$ with quadratic perturbations, respectively.

In mathematical morphology [14, 15, 16], the Moreau lower and upper envelopes can be viewed as ‘greyscale’ erosions and dilations by quadratic structuring elements [17, 18, 19]. Despite the rich structure of these filters, the use of quadratic structuring functions has not however been fully exploited compared, for instance, to the widespread use of the flat structuring elements [15, 16]. This is usually attributed to the fact that the gray–scale dilation and erosion operators based on a flat structuring function are easy to implement (given that they reduce to a local maximum and minimum filter) and result therefore in fast algorithms. By the definition in terms of the compensated convex transforms it is possible however to offer not only an alternative interpretation of the transforms (1.3) as ‘one-step’ morphological openings and closings, but also to provide an easy evaluation of the geometric properties of such filters [3] and to obtain competetively fast algorithms as we illustrate in this paper.

In variational analysis, by contrast, the lower compensated convex transform (1.1) is also known as proximity hull [20, Example 11.4] and it has been redefined as quadratic envelope in [21]. But it is only starting from [1, 4] that the relation of the proximity hull to the convex envelope, and to the envelope of quadratic functions with given curvature established first in [3, Eq. (1.4)] and then in [21], has been studied in a systematic manner and applied as such to image processing, shape interrogation, and function approximation in [6, 22, 23] and to compressive sensing and low rank approximation in [21].

Note that (1.1), (1.2) and their alternative representations (1.3) are given for extended real valued functions. They can thus be applied to functions $g$, defined in proper subsets $\Omega$ of $\mathbb{R}^{n}$ by defining, as a common practice in convex analysis [24], the following extensions of $g$ in $\mathbb{R}^{n}\setminus\Omega$,

[TABLE]

for the definition of $C_{\lambda}^{l}(g^{\infty})$, and

[TABLE]

for the definition of $C_{\lambda}^{u}(g^{-\infty})$. However, such natural and direct definitions of local compensated convex transforms for functions defined in $\Omega$ using (1.3) depend on values of the Moreau envelope at points outside the domain $\Omega$, while convex envelope based methods using (1.1) or (1.2) will rely on an a-priori assumption about the boundary values of the transforms and on the application of convex envelope based schemes which, as far as we know, are neither efficient nor accurate.

These problems lead us to design the following simple local compensated convex transforms based on the mixed Moreau envelope definitions (1.3) of the compensated convex transforms without the need of calculating values of the Moreau envelopes outside the bounded closed domain $\Omega$. Before introducing our local transforms, we introduce some notation and recall some definitions.

Let $\Omega$ be a non-empty bounded open convex subset of $\mathbb{R}^{n}$. We consider bounded functions $f:\Omega\subset\mathbb{R}^{n}\to\mathbb{R}$ satisfying $m\leq f(x)\leq M$ in $\Omega$ for some constants $-\infty<m\leq M<+\infty$. Without loss of generality, and if not otherwise specified, we set $m=\inf_{\Omega}f$ and $M=\sup_{\Omega}f$ and define the oscillation of $f$ in $\Omega$ by

[TABLE]

We consider the auxiliary functions, $f^{-}_{\overline{\Omega}}$ and $f^{+}_{\overline{\Omega}}$ that extend $f$ from $\Omega$ to its closure $\overline{\Omega}$,

[TABLE]

and the auxiliary functions $f_{\mathbb{R}^{n}}^{-}$ and $f_{\mathbb{R}^{n}}^{+}$ that extend $f$ from $\Omega$ to the whole space $\mathbb{R}^{n}$,

[TABLE]

In practice, the extensions $f_{\overline{\Omega}}^{-}$, $f_{\overline{\Omega}}^{+}$ of $f$ to the boundary $\partial\Omega$ correspond to adding a frame of one pixel wide layer on the boundary of the data array and defining $f$ at each point of the frame by the maximum value or the minimum value of the function.

For the function $f_{\mathbb{R}^{n}}^{-}$ and $f_{\mathbb{R}^{n}}^{+}$ we then consider the following transformations which are well defined for $x\in\mathbb{R}^{n}$,

[TABLE]

and their characterization (1.3) in terms of the critical mixed Moreau envelopes:

[TABLE]

For a bounded function $g:\overline{\Omega}\to\mathbb{R}$ defined in the closure of a convex bounded open set $\Omega$, we also introduce for $x\in\overline{\Omega}$ the notation $M_{\lambda,\Omega}(g)(x)$ and $M^{\lambda}_{\Omega}(g)(x)$ to denote the following inf- and *sup-*convolutions of $g$ with quadratic weights,

[TABLE]

which coincide with the restriction to $\overline{\Omega}$ of $M_{\lambda}(g^{\infty})$ and $M^{\lambda}(g^{-\infty})$, respectively.

Motivated by the characterization (1.3) for compensated convex transforms, we now define the local lower compensated convex transform of $f$ in $\overline{\Omega}$, as

[TABLE]

and the local upper compensated convex transform of $f$ in $\overline{\Omega}$, as

[TABLE]

An important feature of the special extensions $f_{\overline{\Omega}}^{-}$, $f_{\overline{\Omega}}^{+}$, $f_{\mathbb{R}^{n}}^{-}$, $f_{\mathbb{R}^{n}}^{+}$ of functions to $\overline{\Omega}$ and to $\mathbb{R}^{n}$ is that we have, for $x\in\Omega$ that

[TABLE]

The equalities (1.14) do not hold, in general, if we had considered the extensions $f^{\infty}_{\overline{\Omega}}$ and $f^{\infty}$ [respect. $f^{-\infty}_{\overline{\Omega}}$ and $f^{-\infty}$] with $f^{\infty}_{\overline{\Omega}}$ [respect. $f^{-\infty}_{\overline{\Omega}}$] that extends $f$ over $\overline{\Omega}$ by $\infty$ [respect. $-\infty$] on $\partial\Omega$. While our extensions can make the local versions equals to the global ones, the infinity versions would not have this property.

After this brief introduction, in the next Section we present background results on convex analysis and the theory of compensated convex transforms and we will also compare the definitions of the compensated convex transforms for the different extensions of $f$ over $\mathbb{R}^{n}$ by showing that for any $x\in\mathbb{R}^{n}$, there holds (Proposition 2.9)

[TABLE]

where $C_{\lambda,\Omega}^{l}(f)$ and $C_{\lambda,\Omega}^{u}(f)$ are defined, respectively, for $x\in\overline{\Omega}$ as

[TABLE]

by means of (1.11a) and (1.11b). The main theoretical results are stated in Section 3. Here, we state that if $\Omega\subset\mathbb{R}^{n}$ is a bounded convex open set and $f:\Omega\to\mathbb{R}$ a bounded function, then for any $x\in\overline{\Omega}$ (see Theorem 3.1)

[TABLE]

Furthermore, we also show that there exists a constant $\kappa(\lambda,\,f)$ which depends on $\lambda$ and $f$ such that at points $x\in\Omega$ with $\mathrm{dist}(x,\partial\Omega)>\kappa(\lambda,\,f)$, the values of $C^{l}_{\lambda,\,\Omega}(f_{\overline{\Omega}}^{-})(x)$ and $C^{u}_{\lambda,\,\Omega}(f_{\overline{\Omega}}^{+})(x)$ depend only on the values of $f$ on $\Omega$. For the special case of the characteristic function $\chi_{K}$ of a compact set $K\subset\Omega$, which represents ‘geometric shapes’ (binary data), we use the natural and simple extension $\chi^{\Omega}_{K}$ defined in (3.9), which is the restriction of the characteristic function $\chi_{K}$ to $\overline{\Omega}$, and $\chi_{K}$ itself, rather than the extensions defined in (1.7) and (1.8), respectively. Under the condition that $\mathrm{dist}^{2}(x,\partial\Omega)>1/\lambda$, we establish $C^{u}_{\lambda,\Omega}(\chi^{\Omega}_{K})(x)=C^{u}_{\lambda}(\chi_{K})(x)$ (see also the remarks about Theorem 3.7 below). We present in Section 4 an algorithm that allows the numerical realization of the Moreau envelope, whereas Section 5 contains numerical experiments which illustrate how to apply our theoretical findings to carry out, for instance, image processing and computational geometry tasks. As examples, we discuss the finding of ridges in the graph of a function (multiscale medial axis map) and intersections of curves in a plane, image inpainting and salt & pepper denoising. Section 6 concludes the paper with some final remarks whereas the Appendix contains the proofs of the main results.

2 Notation and Preliminaries

This section presents a brief overview of some basic results in convex analysis and the theory of compensated convex transforms that will be used in the sequel for the proof of the main results; for a comprehensive account of convex analysis, we refer to Refs. [24, 25], and to Refs. [1, 3] for an account of the theory of compensated convex transforms.

Proposition 2.1**.**

Let $f:\mathbb{R}^{n}\to(-\infty,\,+\infty]$ be coercive in the sense that $f(x)/|x|\to\infty$ as $|x|\to\infty$, and $x_{0}\in\mathbb{R}^{n}$. Denote by $\mathrm{Aff}(\mathbb{R}^{n},\mathbb{R})$ the set of real–valued affine functions on $\mathbb{R}^{n}$. Then

• (i)

The value $\operatorname{\mathsf{co}}\left[f\right](x_{0})$ of the convex envelope of $f$ at $x_{0}\in\mathbb{R}^{n}$ is given by

[TABLE]

If, in addition, $f$ is lower semicontinuous, then the infimum is attained, that is, for some $1\leq p\leq n+1$ there are $\lambda_{k}^{\ast}>0$, $y_{k}^{\ast}\in\mathbb{R}^{n}$ for $k=1,\dots,p$, satisfying $\sum^{p}_{k=1}\lambda^{\ast}_{k}=1$ and $\sum^{p}_{k=1}\lambda^{\ast}_{k}y^{\ast}_{k}=x_{0}$ such that the points $\left(y_{k}^{\ast},f(y_{k}^{\ast})\right)$, $k=1,\ldots,\,p$, lie in the intersection of a supporting plane of the epigraph of $f$, $\mathrm{epi}(f)$, and $\mathrm{epi}(f)$, and

[TABLE]

• (ii)

The value $\operatorname{\mathsf{co}}\left[f\right](x_{0})$, for $f$ taking only finite values, can also be obtained as follows:

[TABLE]

with the $\sup$ attained by an affine function $\ell^{\ast}\in\mathrm{Aff}(\mathbb{R}^{n})$ if $f$ is lower semicontinuous.

We will also introduce the following local version of convex envelope at a point.

Definition 2.2**.**

Let $r>0$ and $x_{0}\in\mathbb{R}^{n}$. Denote by $B(x_{0};r)$ the open ball centered at $x_{0}$ with radius $r$, and by $\overline{B}(x_{0};r)$ the corresponding closed ball. Suppose $f:\overline{B}(x_{0};r)\to\mathbb{R}$ is a bounded function in $\overline{B}(x_{0};r)$. Then the value $\operatorname{\mathsf{co}}_{\bar{B}(x_{0};r)}\left[f\right](x_{0})$ of the local convex envelope of $f$ at $x_{0}$ in $\overline{B}(x_{0};r)$ is defined by

[TABLE]

Remark 2.3**.**

If $f$ is lower semicontinuous, then by using the second part of Proposition 2.1$(i)$, we see that the infimum is attained in $\overline{B}(x_{0};r)$. This means that for some $1\leq p\leq n+1$ there are $\lambda_{k}>0$, $y_{k}\in\mathbb{R}^{n}$ such that $|y_{k}-x_{0}|\leq r$ for $k=1,\dots,p$, satisfying $\sum^{p}_{k=1}\lambda_{k}=1$, $\sum^{p}_{k=1}\lambda_{k}y_{k}=x_{0}$ and

[TABLE]

thus, in this case, $\operatorname{\mathsf{co}}[f](x_{0})$ depends only on the values of $f$ in $\overline{B}(x_{0};\,r)$.

We recall also the following ordering properties for compensated convex transforms which can be found in Ref. [1],

[TABLE]

whereas for $f\leq g$ in $\mathbb{R}^{n}$, we have that

[TABLE]

Proposition 2.4**.**

(Translation invariance property) For any $f:\mathbb{R}^{n}\to(-\infty,\,+\infty]$ bounded below and for any affine function $\ell:\mathbb{R}^{n}\to\mathbb{R}$, $\operatorname{\mathsf{co}}[f+\ell]=\operatorname{\mathsf{co}}[f]+\ell$. Consequently, both $C^{u}_{\lambda}(f)$ and $C^{l}_{\lambda}(f)$ are translation invariant against the weight function, that is

[TABLE]

for all $x\in\mathbb{R}^{n}$ and for every fixed $x_{0}$. In particular, at $x_{0}$, we have

[TABLE]

For both theoretical and numerical developments, the following property on the locality of the compensated convex transforms for Lipschitz functions and bounded functions plays a fundamental role. The result for bounded functions is a slight modification of the locality property stated in Theorem 3.10, Ref. [3]. For a locally bounded function $f:\mathbb{R}^{n}\to\mathbb{R}$, we define the upper and the lower semicontinuous closure $\overline{f}$ and $\underline{f}$[24, 25], respectively, by

[TABLE]

We have the following result.

Proposition 2.5**.**

Suppose $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is bounded. Let $\lambda>0$ and $x_{0}\in\mathbb{R}^{n}$. Then the following locality properties hold:

[TABLE]

with $R_{f,\lambda}=3\sqrt{O_{f}/\lambda}$. If $f$ is Lipschitz continuous with Lipschitz constant $L\geq 0$, then $R_{f,\lambda}=4L/\lambda$.

Remark 2.6**.**

The values of $R_{f,\lambda}$ given here have improved upon those obtained in Theorem 3.10, Ref. [3].

By Remark 2.3, the convex envelope over $B(x_{0};R_{f,\lambda})$ which enters (2.8) is therefore given by

[TABLE]

for some $1\leq p\leq n+1$, $\lambda_{k}>0$, $y_{k}\in\mathbb{R}^{n}$ for $k=1,\dots,p$, satisfying $\sum^{p}_{k=1}\lambda_{k}=1$, $\sum^{p}_{k=1}\lambda_{k}y_{k}=x_{0}$ and

[TABLE]

for all $k=1,\dots,p$. Similar conclusion can be drawn for $\operatorname{\mathsf{co}}_{B(x_{0};R_{f,\lambda})}[\lambda|(\cdot)-x_{0}|^{2}-\overline{f}](x_{0})$.

Next we state the locality properties for the Moreau envelopes.

Proposition 2.7**.**

Let $f:\mathbb{R}^{n}\to\mathbb{R}$.

• $(i)$

If $f$ is bounded, then for any fixed $x\in\mathbb{R}^{n}$, if $y_{x},\,z_{x}\in\mathbb{R}^{n}$ satisfy

[TABLE]

then

[TABLE]

• $(ii)$

If $f$ is a Lipschitz function with Lipschitz constant $L\geq 0$. Then for any fixed $x\in\mathbb{R}^{n}$, if $y_{x},\,z_{x}\in\mathbb{R}^{n}$ satisfy

[TABLE]

then

[TABLE]

Remark 2.8**.**

Proposition 2.7$(i)$ was also established in [13, Lemma 3.5.7] whereas 2.7$(ii)$ is partially contained in [26, Theorem 4.2]. In Appendix A we give a new the proof of these results.

For completeness, we recall the following relationship between the lower and upper Moreau envelope and the lower and upper compensated convex transform given by

[TABLE]

and the ordering properties of the Moreau envelope which we give next for the lower Moreau envelope in the case of $f\leq g$ with $f,\,g$ bounded from below,

[TABLE]

The following result precises the connection among the compensated convex transforms for the different extensions of $f:\Omega\subseteq\mathbb{R}^{n}\to\mathbb{R}$ over $\mathbb{R}^{n}$. We give this relationship only for the lower compensated convex transform given that by (2.15) we can establish similar ones for the upper transform.

Proposition 2.9**.**

Let $\Omega\subset\mathbb{R}^{n}$ be a convex open set and assume $f:\Omega\subset\mathbb{R}^{n}\rightarrow\mathbb{R}$ to be bounded. Denote by $\underline{f}:\overline{\Omega}\subset\mathbb{R}^{n}\rightarrow\mathbb{R}$ the lower semicontinuous envelope of $f$ defined by (2.7) at the interior points $x$ and by

[TABLE]

at the boundary points $x\in\partial\Omega$. The following inequalities hold

[TABLE]

for any $x\in\mathbb{R}^{n}$.

Next we recall from [27] the definition of modulus of continuity of a function along with some of its properties which will be used to analyse the numerical scheme that computes the Moreau envelope.

Definition 2.10**.**

Let $f:\mathbb{R}^{n}\to\mathbb{R}$ be a bounded and uniformly continuous function in $\mathbb{R}^{n}$. Then,

[TABLE]

is called the modulus of continuity of $f$.

Proposition 2.11**.**

Let $f:\mathbb{R}^{n}\to\mathbb{R}$ be a bounded and uniformly continuous function in $\mathbb{R}^{n}$. Then the modulus of continuity $\omega_{f}$ of $f$ satisfies the following properties:

[TABLE]

A function $\omega$ defined on $[0,\,\infty)$ satisfying (2.20) is called a modulus of continuity. A modulus of continuity $\omega$ can be bounded from above by an affine function (see Lemma 6.1 of Ref. [27]), that is, there exist some constants $a>0$ and $b\geq 0$ such that

[TABLE]

We conclude this Section by recalling the following definitions.

Let $C$ be a subset of $\mathbb{R}^{n}$. We define the distance of $x\in\mathbb{R}^{n}$ from $C$ as

[TABLE]

the diameter of the set $C$ as

[TABLE]

the indicator function $i_{C}$ of $C\subset\mathbb{R}^{n}$ as the function defined in $\mathbb{R}^{n}$ such that

[TABLE]

and the characteristic function $\chi_{C}$ of $C\subset\mathbb{R}^{n}$ as the function defined in $\mathbb{R}^{n}$ such that

[TABLE]

It is then not difficult to verify that for any $x\in\mathbb{R}^{n}$

[TABLE]

that is, $\mathrm{dist}^{2}(\cdot,\,C)$ is the inf-convolution of $i_{C}$ with $\|\cdot\|^{2}$, and is proportional to the lower Moreau envelope of $\lambda i_{C}$ with parameter $\lambda$.

3 Main Results

The main results given in this section consist of two parts. In the first part, we establish the relationship between the local Moreau envelopes of (1.7) and the global Moreau envelopes of (1.8), and between the corresponding mixed Moreau envelopes. This relationship is a consequence of the type of auxiliary functions under consideration. In the second part, we give conditions that ensure that the local Moreau envelopes and the corresponding mixed Moreau envelopes depend only on the local values of $f$. The precise meaning of this statement will be specified for each result. Next we consider the case of a bounded function $f$ defined on a bounded domain $\Omega$.

Theorem 3.1**.**

Let $\Omega\subset\mathbb{R}^{n}$ be a bounded open set and $f:\Omega\subset\mathbb{R}^{n}\to\mathbb{R}$ a bounded function. Consider the extensions $f_{\overline{\Omega}}^{-}$ and $f_{\mathbb{R}^{n}}^{-}$ given by (1.7c) and (1.8c), respectively, and the extensions $f_{\overline{\Omega}}^{+}$ and $f_{\mathbb{R}^{n}}^{+}$ given by (1.7f) and (1.8f), respectively. Then, for any $x\in\overline{\Omega}$,

[TABLE]

and

[TABLE]

Consequently, for any $x\in\overline{\Omega}$,

[TABLE]

Furthermore, we have the following locality results:

• $(i)$

If $x\in\Omega$ is such that $\mathrm{dist}^{2}(x,\partial\Omega)>O_{f}/\lambda$ and there is $z_{x}\in\mathbb{R}^{n}$ such that

[TABLE]

then $M_{\lambda,\Omega}(f_{\overline{\Omega}}^{-})(x)$ [resp. $M^{\lambda}_{\Omega}(f_{\overline{\Omega}}^{+})(x)$] is determined by values of $f$ on $\Omega$, in the sense that $z_{x}\in\Omega$.

• $(ii)$

If $x\in\Omega$ is such that $\mathrm{dist}^{2}(x,\partial\Omega)>4O_{f}/\lambda$ and there is a $z_{x}\in\mathbb{R}^{n}$ such that

[TABLE]

then $M^{\lambda}_{\Omega}(M_{\lambda,\Omega}(f_{\overline{\Omega}}^{-}))(x)$ [resp. $M^{\lambda}_{\Omega}(M_{\lambda,\Omega}(f_{\overline{\Omega}}^{+}))(x)$] is determined by values of $f$ on $\Omega$, in the sense that $M_{\lambda,\Omega}(f_{\overline{\Omega}}^{-})(z_{x})=f_{\overline{\Omega}}^{-}(y_{x})+\lambda|y_{x}-x|^{2}$ where $z_{x},\,y_{x}\in\Omega$.

Remark 3.2**.**

The locality properties of Theorem 3.1 state that under the conditions $(i)$ and $(ii)$ on $x\in\Omega$, respectively, the values of $C^{l}_{\lambda,\Omega}(f_{\overline{\Omega}}^{-})(x)$ [resp. $C^{l}_{\lambda,\Omega}(f_{\overline{\Omega}}^{+})(x)$] depend on the values of $f$ on $\Omega$. This means that the values of $C^{l}_{\lambda,\Omega}(f_{\overline{\Omega}}^{-})(x)$ [resp. $C^{l}_{\lambda,\Omega}(f_{\overline{\Omega}}^{+})(x)$] are not influenced by the values of $f$ we defined on $\partial\Omega$ when we define $f_{\overline{\Omega}}^{-}(x)$ [resp. $f_{\overline{\Omega}}^{+}(x)$]. We express this by saying that $C^{l}_{\lambda,\Omega}(f_{\overline{\Omega}}^{-})(x)$ [resp. $C^{l}_{\lambda,\Omega}(f_{\overline{\Omega}}^{+})(x)$] is not affected by boundary values.

As an application of Theorem 3.1 we consider the case where $f$ is the squared Euclidean distance function to a closed set. The following two results are useful, for instance, when we need to compute the multiscale medial axis map [6]. Let $K$ be a nonempty closed set, the quadratic multiscale medial axis map of $K$ with scale $\lambda>0$ is defined in [6, Definition 3.1] for $x\in\mathbb{R}^{n}$ by

[TABLE]

Next we describe how $\mathcal{M}(\lambda;\,K)(x)$ can be expressed in terms of the local lower transform. The first result can be applied to find the multiscale medial axis map of the set $\Omega\setminus K$, where $\Omega$ is an open subset of $\mathbb{R}^{n}$ and $K\subset\Omega$ a compact set.

Corollary 3.3**.**

Let $\Omega\subset\mathbb{R}^{n}$ be an open set and $K\subset\Omega$ a non-empty compact set. Let $f(x):=\mathrm{dist}^{2}(x,K\cup\Omega^{c})$ for $x\in\mathbb{R}^{n}$ and $f_{\overline{\Omega}}^{-}(x)$ be defined by (1.7c). Then for $x\in\overline{\Omega}$,

[TABLE]

Remark 3.4**.**

Equation (3.7) actually gives $\mathcal{M}(\lambda;\,K\cup\Omega^{c})(x)$ for any $x\in\mathbb{R}^{n}$ given that $\mathcal{M}(\lambda;\,K\cup\Omega^{c})(x)=0$ for $x\in\mathbb{R}^{n}\setminus\overline{\Omega}$.

The next result, on the other hand, applies when we need to define the multiscale medial axis map of an open set $A\subset\Omega$.

Corollary 3.5**.**

Let $\Omega\subset\mathbb{R}^{n}$ be an open set and $A\subset\Omega$ a non-empty open set. Let $f(x):=\mathrm{dist}^{2}(x,A^{c})$ and define $f_{\overline{\Omega}}^{-}(x)$ by (1.7c). Then for any $x\in\overline{\Omega}$

[TABLE]

Remark 3.6**.**

In this case we also have that $\mathcal{M}(\lambda;\,A^{c})(x)=0$ for $x\in\mathbb{R}^{n}\setminus\overline{\Omega}$.

Next we consider the behavior of the local upper compensated transform of the characteristic function of a compact set in view of applications that involve the processing of binary images.

Theorem 3.7**.**

Let $\Omega\subset\mathbb{R}^{n}$ be an open set and $K\subset\Omega$ a non-empty compact set. Let $\chi_{K}$ denote the characteristic function of $K$ defined in $\mathbb{R}^{n}$ and $\chi_{K}^{\Omega}$ the restriction of $\chi_{K}$ to $\overline{\Omega}$, that is,

[TABLE]

Then, if $\mathrm{dist}^{2}(K,\,\partial\Omega)>1/\lambda$, for any $x\in\overline{\Omega}$

[TABLE]

and

[TABLE]

Consequently, if we define $C_{\lambda,\Omega}^{u}(\chi_{K}^{\Omega})(x)=M_{\lambda,\Omega}(M^{\lambda}_{\Omega}(\chi_{K}^{\Omega}))(x)$, it follows that

[TABLE]

Remark 3.8**.**

• $(i)$

Compared to Theorem 3.1, Theorem 3.7 states that in the case of binary functions we can establish the equalities (3.10) and (3.11) using $\chi_{K}^{\Omega}$ as defined by (3.9), and $\chi_{K}$, rather than the corresponding auxiliary functions (1.7f) and (1.8f), respectively, which are the type of functions that are used in Theorem 3.1.

• $(ii)$

It is possible to establish a locality result for $M_{\Omega}^{\lambda}(\chi_{K}^{\Omega})$ and $M_{\lambda,\Omega}(M_{\Omega}^{\lambda}(\chi_{K}^{\Omega}))$ in the following sense. If $x\in\overline{\Omega}$ is such that $\mathrm{dist}^{2}(x,\,\partial\Omega)>4/\lambda$ then $\mathrm{dist}^{2}(z_{x},\,\partial\Omega)>1/\lambda$, where $z_{x}\in\overline{\Omega}$ is such that

[TABLE]

Also, if $x\in\overline{\Omega}$, then $\mathrm{dist}^{2}(z_{x},\,\partial\Omega)>1/\lambda$ and $\mathrm{dist}^{2}(y_{x},\,\partial\Omega)>1/\lambda$ where $y_{x},\,z_{x}\in\overline{\Omega}$ are such that

[TABLE]

As a consequence, if $K=\{y\in\overline{\Omega}:\,\mathrm{dist}^{2}(y,\,\partial\Omega)\geq 1/\lambda\}$ and $\mathrm{dist}^{2}(x,\,\partial\Omega)>9/\lambda$, then both $M_{\Omega}^{\lambda}(\chi_{K}^{\Omega})(x)$ and $M_{\lambda,\Omega}(M_{\Omega}^{\lambda}(\chi_{K}^{\Omega}))(x)$ are determined only by $K$, i.e. $y_{x},\,z_{x}\in K$. The proof of these results follows from the locality property of the Moreau envelopes (see Proposition 2.7) and by a similar argument to the proof of Theorem 3.1$(i)$ and $(ii)$.

For applications to scattered data approximation and image inpainting [7, 22], we assume $K\subset\Omega\subset\mathbb{R}^{n}$ to be a compact set, $M>0$ and $f:K\to\mathbb{R}$ to be a bounded function, and introduce the following auxiliary functions

[TABLE]

where $\inf_{K}f$ and $\sup_{K}f$ denote the infimum and supremum of $f$ over $K$, respectively.

Remark 3.9**.**

With the notation given above, for the inpainting problem the set $D=\Omega\setminus K$ will be the non-empty open subset of $\Omega$ representing the damaged region of the image of domain $\Omega$ which must be repaired, whereas for scattered data approximations $K$ will be the set of sample points.

Theorem 3.10**.**

Let $\Omega\subset\mathbb{R}^{n}$ be a bounded open set and $K\subset\Omega$ a compact set. Assume $f:K\to\mathbb{R}$ to be bounded, and, for $M>0$, consider the auxiliary functions defined by (3.13).Suppose $\lambda>0$, then, for any $x\in\overline{\Omega}$,

[TABLE]

and

[TABLE]

Furthermore, let $D=\Omega\setminus K$ and assume that $\Omega^{c}\cap\overline{D}=\varnothing$. Define $O_{f}=sup_{K}f-\inf_{K}f$ and $\mathrm{dist}(\partial D,\,\partial\Omega)=\inf_{x\in\partial D}\{\mathrm{dist}(x,\,\partial\Omega)\}$. If $M>\sup_{\Omega\setminus D}\,f+\lambda\operatorname{\mathsf{diam}}^{2}(D)$ and $\mathrm{dist}^{2}(\partial D,\,\partial\Omega)>O_{f}/\lambda$, the following locality properties hold:

• $(i)$

If $x\in\Omega$ is such that $\mathrm{dist}^{2}(x,\partial\Omega)>O_{f}/\lambda$ then $M_{\lambda,\Omega}(f_{\overline{\Omega},K}^{M})(x)$ $[\text{resp. }M_{\Omega}^{\lambda}(f_{\overline{\Omega},K}^{-M})(x)]$ is determined by values of $f|_{\Omega\setminus D}$. More precisely, if

[TABLE]

for some $z_{x}\in\overline{\Omega}$, then $z_{x}\in\Omega\setminus D$.

• $(ii)$

If $x\in\Omega$ is such that $\mathrm{dist}^{2}(x,\partial\Omega)>4O_{f}/\lambda$ then $M^{\lambda}_{\Omega}(M_{\lambda,\Omega}(f_{\overline{\Omega},K}^{M}))(x)$ [resp. $M_{\lambda,\Omega}(M^{\lambda}_{\Omega}(f_{\overline{\Omega},K}^{-M}))(x)$] is determined by values of $f|_{\Omega\setminus D}$. More precisely, if

[TABLE]

for some $z_{x}\in\overline{\Omega}$, then

[TABLE]

with $z_{x},\,y_{x}\in\Omega\setminus D$.

Remark 3.11**.**

• $(i)$

Given the definition of $f_{\overline{\Omega},K}^{M}$, it is not difficult to show that if $M>\sup_{K}\,f+\lambda\operatorname{\mathsf{diam}}^{2}(\Omega)$ and if for $x\in\overline{\Omega}$, $M_{\lambda,\Omega}(f_{\overline{\Omega},K}^{M})(x)=f_{\overline{\Omega},K}^{M}(z_{x})+\lambda|z_{x}-x|^{2}$ with $z_{x}\in\overline{\Omega}$, then $z_{x}$ must belong to $K\cup\partial\Omega$. The locality property of Theorem 3.10 is de facto making more precise this result by stating that $z_{x}\in K$. The proof that, in general, $z_{x}\in K\cup\partial\Omega$ can be realized by contradiction. Assume $z\in\overline{\Omega}\setminus(K\cup\partial\Omega)$, then $f_{\overline{\Omega},K}^{M}(z)=M$, thus we have

[TABLE]

In particular, the inequality holds also for $y\in K$ such that $|y-x|^{2}=\mathrm{dist}^{2}(x,\,K)$, that is,

[TABLE]

which is a contradiction. A similar argument can be made for $M^{\lambda}_{\Omega}(f_{\overline{\Omega},K}^{-M})(x)$.

• $(ii)$

The first part of Theorem 3.10 on the equality between the local and global Moreau envelopes of our auxiliary functions is, in fact, a consequence of Theorem 3.1 applied to the functions (3.13) which are bounded in $\Omega$ and are of the type (1.7) and (1.8) considered in Theorem 3.1. However, by Theorem 3.1, we could only conclude that, for instance, for $z_{x}$ as in $(i)$, $(ii)$, we have $z_{x}\in\Omega$. Thus, the relevance of the results stated in Theorem 3.10 is in the locality properties, that is, under the conditions on $M$ and $\mathrm{dist}(\partial D,\,\partial\Omega)$ we can conclude that $z_{x}\in K=\Omega\setminus D$.

• $(iii)$

Theorem 3.10$(ii)$ can be used to get an estimate of the value of $\lambda$ so that the globally defined compensated convex transforms via the mixed critical Moreau envelopes coincides with the local compensated convex transforms apart from a boundary layer. For instance, if we want to restric these boundary effects only to a one-pixel wide boundary layer, that is for $x\in\Omega$ such that $\mathrm{dist}(x,\partial\Omega)<1$, of an image with values in the range $[0,\,1]$ so that $O_{f}=1$, we need to take at most $\lambda\geq 4O_{f}/\mathrm{dist}^{2}(x,\partial\Omega)=4$.

4 Numerical Scheme

Given the equalities (3.5), (3.6), the numerical realization of the local compensated convex transforms boils down to computing the Moreau envelope of a function defined in $\mathbb{R}^{n}$. Without loss of generality, given the relation (2.15) between the upper and lower Moreau envelopes, we will refer in the following only to the computation of the lower Moreau envelope. To compute the upper Moreau envelope, it is not difficult to adapt the algorithm proposed, or to use the relation (2.15) between the two envelopes.

The computation of the Moreau envelope is an established task in the field of computational convex analysis [28] and references therein, that has been tackled by various different approaches aimed at reducing the complexity of a direct brute force implementation of the transform. The methods developed in [29, 30, 31], for instance, are based on a dimensional reduction. The authors exploit the property that the Moreau envelope can be factored by $n$ 1d Moreau envelopes and its relationship with the Legendre-Fenchel transform [20, Example 11.26]. The factorization of the Moreau envelope is also used in [32] where the construction of the *inf-*convolution is reduced to the computation, in constant time, of the envelope of parabolas.

Moreau envelopes can also be obtained by mathematical morphology operations which can be particularly useful in the case when $f$ represents an image. Such a class of methods can be obtained by an appropriate modification of the ones that compute the Euclidean distance transform of binary images. Here we develop such a method that generalizes the one used in [33, 16, 34] to compute the discrete Euclidean distance transform.

The fundamental idea is the characterization (2.26) of the Euclidean distance in terms of the Moreau envelope of the characteristic function. The Euclidean distance transform is the erosion of the characteristic function by the quadratic structuring element whereas the Moreau envelope is the erosion of the image $f$. Thus, one can think of generalizing the Euclidean distance transform of binary images, by replacing the binary image by an arbitrary function on a grid. The decomposition of the structuring element which yields the exact Euclidean distance transform [34] into basic ones, yields a simple and fast algorithm where the discrete lower Moreau envelope can be computed by a sequence of local operations, using one-dimensional neighborhoods. We will use the same structuring element as in [33] and show that we recover the exact discrete lower Moreau envelope.

Unless otherwise stated, in the following, $i,\,j,\,k,\,r,\,s,\,p,\,q\in\mathbb{Z}$ denote integers whereas $m,\,n\in\mathbb{N}$ are non-negative integers. Given $n\geq 1$, we introduce grid of points of the space $\mathbb{R}^{n}$ with regular spacing $h>0$ denoted by $x_{k}\in\mathbb{R}^{n}$, $k\in\mathbb{Z}$.

Definition 4.1**.**

Suppose $f:\mathbb{R}^{n}\to\mathbb{R}\cup\{\infty\}$ satisfies $f(x)\geq-c_{0}|x|^{2}-C_{1}$ with $c_{0},\,C_{1}\geq 0$. Let $h>0$, $n\geq 1$ and denote by $x_{k}$ a point of the grid of $\mathbb{R}^{n}$ of size $h$. Then the discrete Moreau lower envelope at $x_{k}\in\mathbb{R}^{n}$ is defined by

[TABLE]

for $\lambda>c_{0}$.

The approximation of Moreau lower envelope by the discrete Moreau lower envelope is quantified by the following estimation result. We consider the case where $f$ is a uniformly continuous function first.

Theorem 4.2**.**

Let $f:\mathbb{R}^{n}\to\mathbb{R}$ be a function with modulus of continuity $\omega_{f}$ such that $\omega_{f}(t)\leq at+b$ with $a>0$ and $b\geq 0$ for any $t\geq 0$. Assume $h>0$ and $n\geq 1$. Then, for any grid point $x_{k}$ of the grid of $\mathbb{R}^{n}$ of size $h$,

[TABLE]

where $d(\lambda)=\sqrt{\omega_{f}\left(a/\lambda+\sqrt{b/\lambda}\right)}$.

Remark 4.3**.**

Since the modulus of continuity $\omega_{f}=\omega_{f}(t)$ tends to zero as $t\to 0+$ (see Proposition 2.20$(i)$), it follows that the error bound in (4.2) tends to zero when $h\to 0+$.

The rate of convergence in (4.2) can be improved for $L-$Lipschitz functions $f$. In this case, we have the following result.

Corollary 4.4**.**

Let $f$ be an $L-$Lipschitz function in $\mathbb{R}^{n}$. Assume $h>0$ and $n\geq 1$. Then, for any grid point $x_{k}$ of the grid of $\mathbb{R}^{n}$ of size $h$,

[TABLE]

Remark 4.5**.**

Under the same regularity conditions, [29] obtains the same converegnce rate for a real function of one real variable. However, in contrast to the scheme [29] we obtain directy the values of the Moreau envelope at the grid points by a scheme which has the same complexity as the distance transform [16] and can be easily implemented and applied to any dimension.

In Definition 4.1, the infimum is taken over infinitely many grid points thus its computation is not practical. Therefore we introduce the $m$-th approximation of the discrete Moreau lower envelope where the order of approximation $m$ is related to the number of nodes that are taken to compute the infimum in the definition of the envelope.

Definition 4.6**.**

Let $n\geq 1$ and denote by $x_{k}$ a point of the grid of $\mathbb{R}^{n}$ of size $h$. The $m$-th approximation of the $n-$th dimensional discrete Moreau lower envelope $M_{\lambda}^{h}(f)$ is given by

[TABLE]

for $\lambda>c_{0}$, where $|r|_{\infty}$ is the infinity norm of $r\in\mathbb{Z}^{n}$.

Given $m\geq 1$, to evaluate $g_{m}(x_{k})$ at any point $x_{k}$ of the grid of $\mathbb{R}^{n}$ of size $h$, we can, in fact, consider the values $f_{m}(x_{k})$ that are obtained by applying Algorithm 1:

The relation between $f_{m}(x_{k})$ and $g_{m}(x_{k})$ is described in the following theorem.

Theorem 4.7**.**

Let $f$ be bounded in $\mathbb{R}^{n}$. Assume $\lambda>0$. Then for all the grid points $x_{k}$ of the grid of $\mathbb{R}^{n}$ of size $h$, and $m=0,1,2,\dots$,

[TABLE]

with $f_{m}(x_{k})$ computed by applying Algorithm 1.

By definition, for every grid point $x_{k}$, $g_{m}(x_{k})$ is decreasing in $m$. Since $g_{m}(x_{k})$ is bounded from below by $M^{h}_{\lambda}(f)(x_{k})$, $g_{m}$ will then converge as $m$ goes to $\infty$. The following result actually shows that it will take only finitely many iterations for $g_{m}(x_{k})$ to reach $M^{h}_{\lambda}(f)(x_{k})$.

Proposition 4.8**.**

Let $f$ be bounded in $\mathbb{R}^{n}$. Assume $h>0$ and denote by $x_{k}$ a point of the grid of $\mathbb{R}^{n}$ of size $h$. Then

[TABLE]

for $m\geq\lfloor\tfrac{1}{h}\sqrt{\tfrac{\mathrm{osc}(f)}{\lambda}}\rfloor+1$, where $\lfloor x\rfloor$ denotes the integer part of $x$.

Remark 4.9**.**

For an $8-$bit image with $h=1$ the pixel size and $\mathrm{osc}(f)=255$, if we take $m\geq\lfloor 16/\sqrt{\lambda}\rfloor+1$, we will have $g_{m}(x_{k})=M_{\lambda}^{h}f(x_{k})$ at any grid point $x_{k}$.

For completeness, we conclude this section by giving the scheme that we use for the implementation of the convex based definition of the compensated transforms. The scheme is a generalization of the one introduced in [8] that is briefly summarized in Algorithm 2 and described below. Given a uniform grid of points $x_{k}\in\mathbb{R}^{n}$, equally spaced with grid size $h$, denote by $S_{x_{k}}$ the $d-$point stencil of $\mathbb{R}^{n}$ with center at $x_{k}$ defined as $S_{x_{k}}=\{x_{k}+hr,|r|_{\infty}\leq 1,r\in\mathbb{Z}^{n}\}$ with $|\cdot|_{\infty}$ the $\ell^{\infty}$-norm of $r\in\mathbb{Z}^{n}$ and $d=\#(S)$, cardinality of the finite set $S$. At each grid point $x_{k}$ we compute an approximation of the convex envelope of $f$ at $x_{k}$ by an iterative scheme where each iteration step $m$ is given by

[TABLE]

with the minimum taken between $f(x_{k})$ and only some convex combinations of $(\operatorname{\mathsf{co}}f)_{m-1}$ at the stencil grid points $x_{i}$ of $S_{x_{k}}$. Likewise [8], the scheme can be shown to converge but there is no estimate of the rate of convergence.

5 Numerical Experiments

In this section we first present a one-dimensional and a two-dimensional prototype example with analytical expression for $C_{\lambda,\Omega}^{l}(f_{\overline{\Omega}}^{-})$ which we use:

• $(i)$

to verify the numerical scheme introduced in Section 4;

• $(ii)$

to compare the compensated convex transforms using the local Moreau envelope based definition and the convex based definition.

• $(iii)$

to analyze different behaviors of the transformations $C_{\lambda,\Omega}^{l}(f_{\overline{\Omega}}^{-})$ and $C_{\lambda}^{l}(f^{\infty})$ at the boundary of the domain.

For computing the compensated convex transforms by using the Moreau envelope based definition, we use the iterative scheme described by Algorithm 1, whereas for the realization of the convex based definition of the compensated convex transforms, we apply Algorithm 2, which has already been employed to carry out numerical examples of [6, 4]. In the application of Algorithm 1, rather than fixing the number $m$ of iterations, we introduce the convergence check on the $\ell^{\infty}$ norm of the error between two succesive iterates such as we do when we apply Algorithm 2. We then describe numerical experiments on applications of the local compensated convex transforms on a bounded closed convex domain for the extraction of the multiscale medial axis map, the extraction of Hausdorff stable intersections of smooth manifolds and finally, for the interpolation and approximation of sampled functions.

Example 5.1**.**

A one-dimensional prototype example. We consider a one-dimensional model problem given by the piecewise affine double well model

[TABLE]

for $x\in\overline{\Omega}:=\{x\in\mathbb{R},\;|x|\leq 2\}$. Let the corresponding $f_{\overline{\Omega}}^{-}$ be defined by (1.7c). For $\lambda\geq 1$, it is not difficult to show that

[TABLE]

where $x_{1}=2-\sqrt{\lambda}/\lambda$, $x_{2}=1/(2\lambda)$ and $A(x)=-\lambda x^{2}+(2\lambda-2\sqrt{\lambda}+1)x-\lambda+2\sqrt{\lambda}-1$. Given the definition of $f_{\overline{\Omega}}^{-}$, at the boundary nodes $\partial\Omega=\{-2,\,2\}$, we have that $M_{\lambda,\Omega}(f_{\overline{\Omega}}^{-})(x)=M_{\Omega}^{\lambda}(M_{\lambda,\Omega}(f_{\overline{\Omega}}^{-}))(x)=\inf f$ for $x\in\partial\Omega$. After choosing a uniform grid of $\Omega$ with grid size $h$ which we denote next as $\Omega_{h}$, we can therefore run Algorithm 1 only at the interior points $x_{k}\in\Omega_{h}$ and assume $f_{i}=\inf f$ when the scheme is applied at the first grid point $x_{k}$ of $\Omega_{h}$ next to $-2$ and $2$. For a given grid size $h$ and for any given $m$, number of iterations used in the application of Algorithm 1, we compute the $\ell^{\infty}$-norm of the error defined as

[TABLE]

Figure 1$(a)$ displays the convergence plot with respect to the number of iterations for different values of the grid size $h$ and for given values of the parameter $\lambda$, $\lambda=1$ and $\lambda=2$, whereas Figure 1$(b)$ shows the convergence plot with respect to the grid size $h$, using for each $h$ the value of $m$ such that the $\ell^{\infty}$ norm of the error between two iterates is not greater than $10^{-7}$ and for different values of $\lambda$. We observe that the number $m$ of iterations to obtain convergence increases as $h$ is reduced, consistently with the theoretical finding of Proposition 4.8, whereas Figure 1$(b)$ exhibits the linear convergence rate of the scheme as predicted by Corollary 4.4. The graph of (5.2) and $M^{\lambda,h}_{\Omega}(M_{\lambda,\Omega}^{h}(f^{-}_{\overline{\Omega}}))$, with the latter corresponding to the grid size $h=0.01$ and $\lambda=2$ are shown in Figure 2.

To compare the Moreau computation of $C_{\lambda,\Omega}^{l}(f^{-}_{\overline{\Omega}})$ given by $M^{\lambda,h}_{\Omega}(M_{\lambda,\Omega}^{h}(f^{-}_{\overline{\Omega}}))$, to the convex based definition which relies on the equality $C_{\lambda,\Omega}^{l}(f^{-}_{\overline{\Omega}})=C_{\lambda}^{l}(f^{-}_{\mathbb{R}^{n}})$, Table 1 reports, for different grid size $h$, the error $\|e\|_{\ell^{\infty}}$ and the number of iterations $m$ that compute $M^{\lambda,h}_{\Omega}(M_{\lambda,\Omega}^{h}(f^{-}_{\overline{\Omega}}))$ and $C_{\lambda}^{l}(f^{-}_{\mathbb{R}^{n}})$ such that the $\ell^{\infty}$ norm of the error between two iterates is not greater than $10^{-7}$. Inspection of Table 1 shows that for a given grid size $h$ the Moreau based computation uses a much lower number of iterations, especially for small $h$, and the discrete Moreau based lower transform is more accurate than the discrete convex based lower transform. Furthermore, for small values of $h$, we also note that the convex based scheme (Algorithm 2) appears not to be numerically stable. The number of iterations $m$ relative to the computation of $M^{\lambda,h}_{\Omega}(M_{\lambda,\Omega}^{h}(f^{-}_{\overline{\Omega}}))$ is the total number of iterations which sums up the iterations for computing the lower Moreau envelope and the upper Moreau envelope.

We conclude this example by looking at the behavior of $C^{l}_{\lambda}(f^{\infty})$ with $f^{\infty}$ defined by (1.5). For $\lambda\geq 1/2$, we have then the following explicit formula for $C^{l}_{\lambda}(f^{\infty})$,

[TABLE]

The graph of $C^{l}_{\lambda}(f^{\infty})(x)$ is displayed in Figure 3$(a)$ along with that of $C^{l}_{\lambda,\Omega}(f^{-}_{\overline{\Omega}})$ given by (5.2). Note that the two transforms differ at the boundary of $\Omega$, as a result of the type of singularity therein introduced by the definition of $f^{\infty}$ and $f^{-}_{\overline{\Omega}}$, respectively.

At first sight, the computation of $C^{l}_{\lambda}(f^{\infty})$ as $M^{\lambda}(M_{\lambda}(f^{\infty}))$ would be faced with the problem of having to evaluate $M_{\lambda}(f^{\infty})$ in $\mathbb{R}^{n}$. In this case, one could therefore think of resorting to the convex based definition by the assumption that $C^{l}_{\lambda}(f^{\infty})=f^{\infty}$ on $\partial\Omega$ and the application of a scheme that computes the convex envelope of a function. Although there are several such schemes in the literature, such as the quickhull algorithm [35] or the one introduced in [8], their application can be quite cumbersome when applied to compute the convex envelope of a function defined in $\mathbb{R}^{n}$ with $n>1$, or it can exhibit a slow and unknown rate of convergence as in [8] which would not allow any prediction for the rate of error reduction. As a result, even in the case of $C^{l}_{\lambda}(f^{\infty})$ one might wonder whether it is possible to use Algorithm 1 to obtain an approximation of $C^{l}_{\lambda}(f^{\infty})$. Given the localization effect of the inf- and sup-convolution with quadratic perturbations, we can also use the Moreau based definition of $C^{l}_{\lambda}(f^{\infty})$ and apply Algorithm 1 provided that one computes $M_{\lambda}(f^{\infty})$ and $M^{\lambda}(M_{\lambda}(f^{\infty}))$ over an extended domain $\Omega_{e}$ that contains $\Omega$. Figure 3(b) displays the graph of the transforms in the case of $\Omega_{e}=\Omega$ and of $\Omega\subset\Omega_{e}$ with

[TABLE]

In the first case, we are actually computing the transform $M_{\Omega}^{\lambda}(M_{\lambda}^{\infty}(f^{\infty}))$ where we have set

[TABLE]

which produces a boundary error, whereas in the second case we obtain an excellent approximation of $C_{\lambda}^{l}(f^{\infty})$ with $\|e\|_{\infty}$=0.0002. Note that how big the domain $\Omega_{e}$ must be to ensure that $C^{l}_{\lambda}(f^{\infty})(x)=M_{\Omega_{e}}^{\lambda}(M_{\lambda,\Omega_{e}}(f^{\infty}))(x)$ for $x\in\Omega$, is an open question.

Example 5.2**.**

A two-dimensional prototype example. Let $\Omega=B(O;\,2)$ be the open ball with center at the origin $O\in\mathbb{R}^{2}$ and radius $r=2$. Consider the squared distance of $x\in\Omega$ to the boundary $\partial B(O;\,1)$ given by

[TABLE]

and the functions $f_{\overline{\Omega}}^{-}$ and $f^{\infty}$ defined by (1.7c) and (1.5), respectively. Given the radial symmetry of $f_{\overline{\Omega}}^{-}$ and $f^{\infty}$, it is not difficult to verify that, for $\lambda\geq 1$,

[TABLE]

and

[TABLE]

where $r=|x|$, $x_{p}=2-1/\sqrt{1+\lambda}$, $x_{s}=1/(1+\lambda)$ and $m=2(1+2\lambda)-2\sqrt{1+\lambda}$. Though for the computation of (5.6) and (5.7) we could exploit the symmetry of $f_{\overline{\Omega}}^{-}$ and $f^{\infty}$ and reduce their evaluation to $1d$ problems, in order to verify our scheme for $2d$ applications, we will not take the symmetry into account and will refer to $f_{\overline{\Omega}}^{-}$ and $f^{\infty}$ as generic functions of $x\in\mathbb{R}^{2}$. Let $D\subset\mathbb{R}^{2}$ be a box that contains $\Omega$, for instance, let us take

[TABLE]

and consider $f_{\overline{\Omega}}^{-}$ and $f^{\infty}$ to be extended over $D$ by setting $f_{\overline{\Omega}}^{-}(x)=0$ for $x\in D\setminus\Omega$, respectively. Since Algorithm 1 is formulated for a square lattice, we will refer to the above extensions over $D$ for the application of the algorithm and still denote them by $f_{\overline{\Omega}}^{-}$ and $f^{\infty}$.

As in the previous example, Figure 4$(a)$ displays the convergence plot of the $\ell^{\infty}$ norm of the error versus the number of iterations for different grid sizes, with a number of iterations that increases by reducing $h$, whereas Figure 4$(b)$ shows the linear convergence of the error with the grid size. Table 2 contains for different grid sizes $h$, the total number of iterations $m$ and the error $\|e\|_{\ell^{\infty}}$ of the approximations of $M^{\lambda,h}_{\Omega}(M^{h}_{\lambda,\Omega}(f_{\overline{\Omega}}^{-}))$ and $C_{\lambda}^{l}(f_{\mathbb{R}^{n}}^{-})$, for the comparison of the Moreau and convex based computation of the local lower transform, respectively. Also observe that the number of iterations needed to compute the local lower transform using the Moreau based definition is much less than the one needed for the computation based on the convex envelope. For small grid size, while the computation of $M^{\lambda,h}_{\Omega}(M^{h}_{\lambda,\Omega}(f_{\overline{\Omega}}^{-}))$ takes just a few seconds, the computing time of $C_{\lambda}^{l}(f_{\mathbb{R}^{n}}^{-})$ using a tolerance for the convergence of the scheme equal to $10^{-7}$ for the $\ell^{\infty}$ norm of the difference between two succesive iterates, is much longer. Furthermore, for a given grid size, the discrete Moreau based lower transform is much more accurate than the convex based lower transform and, as in the previous example, for small values of $h$ the convex based scheme appears to be unstable, showing oscillations in the error.

We conclude the discussion of this example with some observations on the computation of $C_{\lambda}^{l}(f^{\infty})$. By means of numerical simulations, we show that there exists a domain $D_{e}$ that contains $\Omega$ such that $C_{\lambda}^{l}(f^{\infty})(x)=M^{\lambda}_{D_{e}}(M_{\lambda,D_{e}}(f^{M}))(x)$ for $x\in\Omega$ with $f^{M}$ defined below. Assume $a>0$ and consider the extended domain

[TABLE]

Given $M>0$, define the following auxiliary function

[TABLE]

Then we apply Algorithm 1 to the interior grid points of $D_{e}$ to compute the lower and upper Moreau envelope. Table 3 reports the value of the $\ell^{\infty}$ norm of the error between $M^{\lambda}_{D_{e}}(M_{\lambda,D_{e}}(f^{M}))(x)$ and $C_{\lambda}^{l}(f^{\infty})(x)$ respectively in $\Omega$ for different values of the grid size $h$, the domain extension parameter ${}^{\prime}a^{\prime}$ and the parameter $\lambda$. We observe that by choosing $M$ large enough and by a suitable choice of $a$, we get an excellent agreement in $\Omega$. Also in this case, we conjecture that this extension depends on $\lambda$, but obtaining a formula for it is an open issue at present.

Example 5.3**.**

The Multiscale Medial Axis Map. We present an application of Corollary 3.5 to find the multiscale medial axis map of the closed set $K$ represented in Figure 5$(a)$. The open set $\Omega$ is taken, in this case, as the domain of the whole image and $A=\Omega\setminus K$. By Corollary 3.5 the quadratic multiscale medial axis map with scale $\lambda$ of the closed set $K$ can then be computed for $x\in\overline{\Omega}$ as

[TABLE]

where $f_{\overline{\Omega}}^{-}(x)=\mathrm{dist}^{2}(x;\,K^{c})$ for $x\in\overline{\Omega}$. Note that in this case $f_{\overline{\Omega}}^{-}(x)=0$ for $x\in\partial\Omega$.

The support of $\mathcal{M}(\lambda;\,K)(x)$ with all its fine branches is shown in Figure 5$(b)$. Figure 5$(c)$ depicts the stable branches that correspond to the suplevel set of $\mathcal{M}(\lambda;\,K)(x)$, given by $\{x\in\Omega:\mathcal{M}(\lambda;\,K)(x)>t\}$ with $t>0$ measure of the branch height, whereas Figure 5$(d)$ displays the results of the morphological thinning of the suplevel set shown in Figure 5$(c)$ using the algorithm described in [36, page 879, bottom of first column through top of second column] implemented in MATLAB. To produce Figure 5$(c)$ we have taken $t=1$ after normalizing $\mathcal{M}(\lambda;\,K)$ to the range $[0,\,255]$. While the medial representation of the leaf stem is present in the support of $\mathcal{M}(\lambda;\,K)$ (see Figure 5$(b)$), this disappears in the suplevel set relative to $t=1$. The very small values of $\mathcal{M}(\lambda;\,K)$ at such points is the result of the small value of the separation angle, which is, in turn, related to the values of $\mathcal{M}(\lambda;\,K)$ (see the bound (3.11) in [6]). The application of Algorithm 1 and Algorithm 2, when $f$ is an image, simplifies by taking the digitized image as the grid, and $h=1$ equal to the pixel size. To compare the performance of the two algorithms, Figure 6 displays the variation of the $L^{2}-$norm of $C_{\lambda,\Omega}^{l}(f_{\overline{\Omega}}^{-})$ with the number of iterations $m$ used to compute $M^{\lambda}_{\Omega}(M_{\lambda,\Omega}(f_{\overline{\Omega}}^{-}))$ by Algorithm 1 and to compute $C_{\lambda}^{l}(f_{\mathbb{R}^{n}}^{-})$ by Algorithm 2. By referring to the Moreau based definition, convergence of the scheme is achieved after a finite number of iterations ($m=118$) which is much lower than those needed to compute the convex envelope based definition of the lower transform ($m=22783$).

Example 5.4**.**

Shape interrogation. As an application of Theorem 3.7, we consider the computation of the intersection extraction filter $I_{\lambda}(\cdot;\,K)$ introduced in [3] with the digitized set $K$ as the input image. Given a non-empty compact set $K\subset\mathbb{R}^{n}$, and taking $\Omega$ as a reference bounding box such that $K\subset\Omega$ and with $K$ distant enough from the boundary of $\Omega$, by Theorem 3.7, the filter $I_{\lambda}(\cdot;\,K)$ can be expressed in terms of the local transforms as

[TABLE]

For the digitized set $K$ given by the collection of curves shown in Figure 7$(a)$, the local maxima of $I_{\lambda,\Omega}(\cdot;\,K)$ coincide with all the crossing and turning points of the set $K$. The filter defined by (5.9) can also be applied to $3d$ geometries. Due to the Hausdorff stability of (5.9), we can also consider $K$ represented by point clouds. Figure 7$(b)$ displays the intersection between manifolds of different dimensions with each manifold sampled by point clouds.

Example 5.5**.**

Image Inpainting. Let $\Omega$ be the domain of the whole image, $D\subset\Omega$ the set of missing/damaged pixels and $K=\Omega\setminus D$ the set of the true pixels. The image inpainting problem consists in reconstructing the image over $D$ from knowing $f_{K}$, if we denote by $f$ the original image. In [22] it is shown that the restored image can be obtained by the average compensated convex transform, which under the assumptions of Theorem 3.10, can be computed for $x\in\overline{\Omega}$ as

[TABLE]

where $f_{\overline{\Omega},K}^{M}$ and $f_{\overline{\Omega},K}^{-M}$ have been defined in (3.13d) and (3.13l), respectively. We give here an application of (5.10) for the problem of removing scratches overprinted over an image such as the one displayed in Figure 8$(a)$. By using the Moreau based definition of the upper and lower transforms that enter (5.10) and by applying Algorithm 1, we compute with $\lambda=5$ and $M=10^{4}$ the averaged image $A_{\lambda,\Omega}^{M}$ defined by (5.10) and consider the restored image given by

[TABLE]

The restored image $I$ is shown in Figure 8$(b)$. The number of iterations needed to obtain $A_{\lambda,\Omega}^{M}$ is $m=34$ for a tolerance equal to $10^{-7}$ on the $\ell^{\infty}$ norm of the difference between two succesive iterates. Taking the $\mathrm{PSNR}$ (peak-to-signal ratio) as a measure of the quality of the restoration, which is expressed in the units of $\mathrm{dB}$ and, for an $8-$bit image, is defined by

[TABLE]

where $f_{i,j}$ and $r_{i,j}$ denote the pixels values of the original and restored image, respectively, and $m,\,n$ denote the size of the image $f$, we find a value of $\mathrm{PSNR}$ equal to $33.165\,\mathrm{dB}$. Figure 8$c$ displays the restored image based on the convex based transforms, where we use the convex based definition of the transforms and apply Algorithm 2 to compute the convex envelope with a $tol=10^{-7}$ for the $\ell^{\infty}$ norm of the difference between two successive iterates. We assume $C_{\lambda}^{l}(f_{K}^{M})=C_{\lambda}^{u}(f_{K}^{-M})=f_{K}$ on the boundary of the image array. In this case, we needed $m=1063$ iterations and got $\mathrm{PSNR}=33.239\,\mathrm{dB}$, which is slightly higher than the one that uses the Moreau based transforms but at the expenses of an higher number of iterations. Figure 8$d$ displays the TV based restoration obtained by applying the TV inpainting method described in [37] and solved by the split Bregman method described in [38, 39]. In this case $m=3531$ and $\mathrm{PSNR}=33.027\,\mathrm{dB}$. Figure 9 compares the details of the original image and of the restored images near the right eye, respectively, showing that $A_{\lambda,\Omega}^{M}$, either by the Moreau based transforms or by the Convex based transforms, is able to preserve image details and does not introduce unintended effects.

Example 5.6**.**

Salt & pepper noise removal. As a further example of scattered data approximation and application of Theorem 3.10, we consider the restoration of an image corrupted by salt & pepper noise by computing $A_{\lambda,\overline{\Omega}}^{M}(f_{K})$ defined by (5.11). Given the picture displayed in Figure 10$(a)$ with size $512\times 512$ pixels and damaged by $70\%$ salt & pepper noise as shown in Figure 10$(b)$, let us denote by $K$ the set of the true pixels, by $\Omega$ the domain of the whole image and $f_{K}$ the values of free noise pixels. Since for $x\in\Omega$ the values of $C_{\lambda,\Omega}^{l}(f_{\overline{\Omega},K}^{M})(x)$ and $C_{\lambda,\Omega}^{u}(f_{\overline{\Omega},K}^{-M})(x)$ depend on $f_{K}$ if $\mathrm{dist}(x,\partial\Omega)>4O_{f}/\lambda$, to reduce the boundary error due to the redefinition of $f_{K}$ on $\partial\Omega$, we consider an enlarged image and then restrict the restored image to the original domain. The enlarged image is obtained by padding one pixel before the first image element and after the last image element along each dimension, defining $f_{K}$ thereing equal to $\min_{K}f$ or $\max_{K}f$ according to whether we are computing $C_{\lambda,\Omega}^{l}(f_{\overline{\Omega},K}^{M})$ or $C_{\lambda,\Omega}^{u}(f_{\overline{\Omega},K}^{-M})$, respectively. We apply then Algorithm 1 and Algorithm 2 using a value of $tol=10^{-7}$ to compute the Moreau and convex based upper and lower compensated convex transforms that enter (5.10), respectively. For the application of Algorithm 1 the convergence check was on the $\ell^{\infty}$ norm of the error between two successive iterates, whereas in the application of Algorithm 2, to limit the computing time of the convex envelope scheme, we opted to make the convergence check by comparing the value of the $PSNR$ between two successive iterates. For $\lambda=15$ and $M=10^{13}$, the restored images with the Moreau and convex based scheme are displayed in Figure 10$(c)$ and Figure 10$(d)$, respectively. For the Moreau based scheme we got a value of $\mathrm{PSNR}$ equal to $28.917\,\mathrm{dB}$ after $m=30$ iterations for a tolerance on the $\ell{\infty}$ norm of the difference between two succesive iterates equal to $10^{-7}$, whereas for the convex based scheme, we got $\mathrm{PSNR}=29.308\,\mathrm{dB}$ after $m=185$ iterations, giving an $\ell^{\infty}$ norm of the error between two succesive iterations equal to $2.18\cdot 10^{-5}$. We note that the restored image by the convex based scheme presents a slightly higher value of the $\mathrm{PSNR}$ than the one that uses the Moreau based transforms but at the expenses of an higher number of iterations. For both the restored images, there is no visible effect on the boundary of the image. Finally, we consider the case of noise density equal to 99%. Figure 11 displays the restored images by the Moreau and convex based definitions of the compensated transforms for $\lambda=5$ and $M=10^{13}$ obtained after a number of iterations equal to $m=78$ and $m=3081$, respectively. The two restored images present comparable values of $\mathrm{PSNR}$, with the Moreau based scheme obtained after a much lower number of iterations than the one that uses the convex based definition of the lower and upper compensated convex transform. Even in this case there is no visible boundary effects for the two restored images, and as a result of the theoretical findings of [22, Corollary 4.7], we note $A_{\lambda,\overline{\Omega}}^{M}(f_{K})$ to be realized by an almost continuous piecewise affine interpolation.

6 Conclusions

Compensated convex transforms, or also known as proximity hull in the case of the lower transform, or grey scale erosion and dilation with quadratic structuring elements in mathematical morphology, provide a geometric tight–approximation method for general functions that yields novel ways to smooth functions, to identify singularities in functions, and to interpolate and approximate data. Introduced for real valued functions defined in $\mathbb{R}^{n}$, in view of applications to image processing and computer aided geometric design, in this paper we have proposed a definition for functions defined on bounded convex open sets such that the following equality holds

[TABLE]

where $f_{\mathbb{R}^{n}}^{-}$ and $f_{\mathbb{R}^{n}}^{+}$ are suitable extensions of $f$ to $\mathbb{R}^{n}$, $C^{l}_{\lambda}(f_{\mathbb{R}^{n}}^{-})$ and $C^{u}_{\lambda}(f_{\mathbb{R}^{n}}^{+})$ are the ‘global’ compensated convex transforms (1.1) and (1.2), whereas $C^{l}_{\lambda,\,\Omega}(f_{\overline{\Omega}}^{-})$ and $C^{l}_{\lambda,\,\Omega}(f_{\overline{\Omega}}^{+})$ are the local versions defined in terms of the Moreau envelopes as (1.12) and (1.13), respectively. We have thus proposed a new algorithm to compute the Moreau envelope which has linear complexity with respect to the number of grid points and we have given applications to tasks of image processing such as image inpainting and restoration of image with high density of salt & pepper noise, and of shape interrogation such as detection of intersection of sampled geometries and multiscale medial axis map. The performance of the methods and the accuracy of the results show that, when coupled with an efficient numerical scheme, the theory of compensated convex transforms provides a valid and feasible alternative to state-of-art methods especially for processing data without any a prior information or represented by point cloud.

Acknowledgements

The authors are grateful to the anonymous reviewers for their comments and suggestions. KZ wishes to thank The University of Nottingham for its support, AO acknowledges the partial financial support of the Argentinian Research Council (CONICET) through the project PIP 100, the National University of Tucumán through the project PIUNT CX-E625 and the FonCyT through the project PICT 2016 201-0105 Prestamo BID, whereas EC is grateful for the financial support of the College of Science, Swansea University.

Appendix A Appendix: Proofs of main results

Proof of (1.3). Next we prove only the characterization of the lower compensated convex transform. The one for the upper transform is obtained by using the relationship (2.15). From the definition (1.1) of $C_{\lambda}^{(}f)$ and the characterization (2.3) of the convex envelope, we have

[TABLE]

which shows the lower compensated transform to be equal to the envelope of quadratic functions with given curvature $\lambda$. If we represent the elements of $\mathrm{Aff}(\mathbb{R}^{n})$ as

[TABLE]

then

[TABLE]

thus (A.1) reads also as

[TABLE]

where we have used the fact that for any $z\in\mathbb{R}^{n}$, the largest $c\in\mathbb{R}$ such that

[TABLE]

is given by

[TABLE]

which is $M_{\lambda}(f)(z)$. ∎

Proof of Proposition 2.7. Part $(i)$: Since

[TABLE]

then

[TABLE]

which concludes the proof for $M_{\lambda}(f)(x)$. The proof for $M^{\lambda}(f)(x)$ follows similar arguments. ∎

Part $(ii)$: Again we have

[TABLE]

so that $\lambda|y_{x}-x|^{2}\leq f(x)-f(y_{x})\leq L|y_{x}-x|$. Thus $|y_{x}-x|\leq L/\lambda$. The proof for $M^{\lambda}(f)(x)$ follows similar arguments. ∎

Proof of Proposition 2.9: We first observe that by definition, for $z\in\overline{\Omega}$

[TABLE]

as $f^{\infty}(y)+\lambda|y-z|^{2}=+\infty$ if $y\in\mathbb{R}^{n}\setminus{\overline{\Omega}}$, thus for $x\in\overline{\Omega}$, we have

[TABLE]

hence $C^{l}_{\lambda}(f^{\infty})(x)\geq C^{l}_{\lambda,\Omega}(f)(x)$. By the definition of lower transform for $x\in\overline{\Omega}$, we find that

[TABLE]

given that if there exists $x_{i}\notin\overline{\Omega}$ and $\lambda_{j}>0$, then the sum is $+\infty$. Furthermore, since $M_{\lambda,\Omega}(f_{\overline{\Omega}}^{-})=M_{\lambda,\Omega}(\underline{f_{\overline{\Omega}}^{-}})$ and by definition (1.7c), (2.7) and (2.17), for $x\in\overline{\Omega}$ it is $\underline{f^{-}_{\overline{\Omega}}}(x)\leq\underline{f}(x)$, then it is $M_{\lambda,\Omega}(f^{-}_{\overline{\Omega}})(x)\leq M_{\lambda,\Omega}(\underline{f})(x)=M_{\lambda,\Omega}(f)(x)$, hence

[TABLE]

thus $C^{l}_{\lambda,\Omega}(f^{-}_{\overline{\Omega}})(x)\leq C^{l}_{\lambda,\Omega}(f)(x)$. By comparing (A.5), (A.6) and (A.7) we conclude the proof. ∎

Proof of Theorem 3.1. We next present the proofs only for the lower transform, given that the results for the upper transform will follow from (2.15).

Proof of (3.1): Without loss of generality, suppose that $\inf_{\Omega}f=0$. By definition, for $x\in\overline{\Omega}$, we have

[TABLE]

given that the second infimum is taken over a larger set. In order to show that also the opposite inequality holds, first let $\left(y_{k}\right)_{k\in\mathbb{N}}$ be a minimizing sequence of $g(y):=f^{-}_{\mathbb{R}^{n}}(y)+\lambda|y-x|^{2}$ in $\mathbb{R}^{n}$. Clearly $g$ is coercive in $\mathbb{R}^{n}$ in the sense that $g(y)/|y|\to+\infty$ as $|y|\to+\infty$. Thus $\left(y_{k}\right)_{k\in\mathbb{N}}$ is bounded in $\mathbb{R}^{n}$ and, therefore, there exists a $y_{x}\in\mathbb{R}^{n}$ and a convergent subsequence $\left(y_{k_{j}}\right)_{k_{j}\in\mathbb{N}}$ such that $y_{k_{j}}\to y_{x}\in\mathbb{R}^{n}$.

Next we show that $y_{x}\in\overline{\Omega}$ and then we will use $y_{x}$ to prove that there also holds $M_{\lambda}(f^{-}_{\mathbb{R}^{n}})(x)\geq M_{\lambda,\Omega}(f_{\overline{\Omega}}^{-})(x)$.

If $y_{x}\notin\overline{\Omega}$, since $f^{-}_{\mathbb{R}^{n}}(y)=0=\inf_{\Omega}f$ for $y\in\mathbb{R}^{n}\setminus\Omega$, then $f^{-}_{\mathbb{R}^{n}}$ is lower semicontinuous in $\mathbb{R}^{n}\setminus\Omega$. As a result, $g(y):=f^{-}_{\mathbb{R}^{n}}(y)+\lambda|y-x|^{2}$ is lower semicontinuous at $y_{x}$ and $y_{x}$ is a minimizer of $g(y)$ in $\mathbb{R}^{n}$. We would have therefore

[TABLE]

Since $y_{x}\neq x$, along the line segment $[x,\,y_{x}]$, there would then exist a point $\xi_{x}\in\partial\Omega$ such that

[TABLE]

which is a contadiction, and therefore $y_{x}\in\overline{\Omega}$.

Let us now examine separately the case of $y_{x}\in\partial\Omega$ and $y_{x}\in\Omega$. If $y_{x}\in\partial\Omega$, by an argument similar to the above, we conclude that $g(y)$ is lower semicontinuous at $y_{x}\notin\Omega$ and therefore

[TABLE]

We now have to distinguish the two cases of $x\in\partial\Omega$ and $x\in\Omega$. If $x\in\partial\Omega$, from (A.8), (A.10) and the ordering property of the lower Moreau envelope (2.16), we have

[TABLE]

that is, $y_{x}=x$, hence

[TABLE]

If $x\in\Omega$, on the other hand, let us consider the convergent minimizing sequence $(y_{k_{j}})_{j\in\mathbb{N}}\in\mathbb{R}^{n}$ of $f_{\mathbb{R}^{n}}^{-}(y)+\lambda|y-x|^{2}$. We note that for $j\in\mathbb{N}$ such that $y_{k_{j}}\in\overline{\Omega}$, we have

[TABLE]

and for $j\in\mathbb{N}$ such that $y_{k_{j}}\in\mathbb{R}^{n}\setminus\overline{\Omega}$, there exists $z_{j}\in[x,\,y_{k_{j}}]$ and $z_{j}\in\partial\Omega$ so that

[TABLE]

Thus, from (A.13) and (A.14) we get that for any $j\in\mathbb{N}$, it is

[TABLE]

By letting $j\to\infty$, since $(y_{k_{j}})_{j\in\mathbb{N}}$ is a minimizing sequence of $f^{-}_{\mathbb{R}^{n}}(y)+\lambda|y-x|^{2}$, we conclude that

[TABLE]

If $y_{x}\in\Omega$, for $j>0$ sufficiently large, $y_{k_{j}}\in\Omega$, thus

[TABLE]

Since $(y_{k_{j}})_{j\in\mathbb{N}}$ is a minimizing sequence of $g(y)$ in $\mathbb{R}^{n}$, if we let $j\to\infty$, we obtain also in this case

[TABLE]

By comparing (A.15) and (A.16) with (A.8), we conclude that $M_{\lambda}(f^{-}_{\mathbb{R}^{n}})(x)=M_{\lambda,\Omega}(f^{-}_{\overline{\Omega}})(x)$ for $x\in\overline{\Omega}$ and this completes the proof. ∎

Proof of (3.3): We first observe that for $x\notin\overline{\Omega}$, as [math] is the minimum value of $f_{\mathbb{R}^{n}}^{-}$, we have $M_{\lambda}(f_{\mathbb{R}^{n}}^{-})(x)=0$. Thus

[TABLE]

To complete the proof, we need to show that for $x\in\overline{\Omega}$, $M^{\lambda}(M_{\lambda}(f_{\mathbb{R}^{n}}^{-}))(x)=M^{\lambda}_{\Omega}(M_{\lambda,\Omega}(f_{\overline{\Omega}}^{-}(x))$. By definition, for $x\in\overline{\Omega}$, we have that

[TABLE]

Since $f_{\mathbb{R}^{n}}^{-}$ is bounded in $\mathbb{R}^{n}$, then $M_{\lambda}(f_{\mathbb{R}^{n}}^{-})$ is also a bounded function in $\mathbb{R}^{n}$, thus $g(y):=M_{\lambda}(f_{\mathbb{R}^{n}}^{-})(y)-\lambda|y-x|^{2}$ is coercive, in the sense that $\lim_{|y|\to\infty}\,g(y)/|y|=-\infty$. Furthermore, $g=g(y)$ is also continuous because so is $M_{\lambda}(f_{\mathbb{R}^{n}}^{-})(y)$ [12, Proposition 1.1], thus a maximum point, say, $y_{x}\in\mathbb{R}^{n}$ of $g=g(y)$ exists. If $y_{x}\notin\overline{\Omega}$, we have

[TABLE]

as $M_{\lambda}(f^{-}_{\mathbb{R}^{n}})(y_{x})=0$. Thus $y_{x}\neq x$. Again if we take the point $\eta_{x}\in\partial\Omega$ in the line segment $[x,\,y_{x}]$, then

[TABLE]

This is a contradiction. Thus $M^{\lambda}(M_{\lambda}(f^{-}_{\mathbb{R}^{n}}))(x)=M^{\lambda}_{\Omega}(M_{\lambda,\Omega}(f_{\overline{\Omega}}^{-}(x))$ for $x\in\overline{\Omega}$ so that

[TABLE]

for $x\in\overline{\Omega}$ and this concludes the proof. ∎

As for the estimate in $(i)$, it follows as argued above that there exists a $z_{x}\in\mathbb{R}^{n}$ such that

[TABLE]

which yields the estimate

[TABLE]

Thus if $\mathrm{dist}(x,\partial\Omega)>\sqrt{O_{f}/\lambda}$, then by the triangle inequality we conclude that $z_{x}\in\Omega$.

The estimate in $(ii)$ is obtained by a similar argument applied to $M_{\lambda,\Omega}(f_{\Omega}^{-})$. ∎

Proof of Corollary 3.3: Since $f(x)=0$ for $x\not\in\overline{\Omega}$, this follows from the fact that $f=f^{-}_{\mathbb{R}^{n}}$, the definition of $f_{\overline{\Omega}}^{-}$ and (3.5). ∎

Proof of Corollary 3.5: Same arguments as in the previous corollary apply. ∎

Proof of Theorem 3.7: First we show that if $x\in\mathbb{R}^{n}\setminus\overline{\Omega}$, then $M^{\lambda}(\chi_{K})(x)=0$. Suppose that for some $x\in\mathbb{R}^{n}\setminus\overline{\Omega}$, $M^{\lambda}(\chi_{K})(x)\not=0$. By definition (1.4) of $M^{\lambda}(\chi_{K})$ and since $x\not\in K$, we have that

[TABLE]

Since $\chi_{K}$ is upper-semicontinuous and $\chi_{K}(y)-\lambda|y-x|^{2}\to-\infty$ as $|y|\to\infty$, there exists a $y_{0}\in\mathbb{R}^{n}$ such that

[TABLE]

thus $y_{0}\in K$, otherwise $M^{\lambda}(\chi_{K})(x)\leq 0$. So we have that

[TABLE]

thus

[TABLE]

which contradicts that $\mathrm{dist}(K,\partial\Omega)>\tfrac{1}{\sqrt{\lambda}}$.

Next we show that for any $x\in\overline{\Omega}$, $M_{\Omega}^{\lambda}(\chi_{K}^{\Omega})(x)=M^{\lambda}(\chi_{K})(x)$. If, for some $x\in\overline{\Omega}$, there exists $y_{0}\not\in\overline{\Omega}$ such that

[TABLE]

then

[TABLE]

so that

[TABLE]

which contradicts the fact that $M^{\lambda}(\chi_{K})(x)\geq 0$ for any $x\in\mathbb{R}^{n}$. Thus $y_{0}\in\overline{\Omega}$ and so $M^{\lambda}(\chi_{K})(x)=M_{\Omega}^{\lambda}(\chi_{K}^{\Omega})(x)$ for $x\in\overline{\Omega}$.

Now we consider the mixed transform. Let $x\in\overline{\Omega}$. If there exists $y_{0}\not\in\overline{\Omega}$ such that

[TABLE]

then

[TABLE]

given that $y_{0}\not\in\overline{\Omega}$. Then consider the segment

[TABLE]

for $t\in[0,\,1]$. We have $\ell(0)=y_{0}\not\in\overline{\Omega}$ and $\ell(1)=x\in\overline{\Omega}$. So there exists $\xi\in]0,\,1[$ such that

[TABLE]

We then have

[TABLE]

which is a contradiction as $0\leq(1-\xi)^{2}<1$, thus $y_{0}\in\overline{\Omega}$ and

[TABLE]

which concludes the proof. ∎

Proof of Theorem 3.10: The statements (3.14) and (3.15) are a consequence of Theorem 3.1 applied to the auxiliary functions (3.13) which are bounded functions in $\Omega$ of the type (1.7) and (1.8) considered in Theorem 3.1. Therefore, next, we will focus only on the proof of the local properties.

Part $(i)$: Assume $x\in D$ and let $z_{x}$ be such that

[TABLE]

We will first show that $z_{x}\in(\Omega\setminus D)\cup\partial\Omega$, otherwise we have

[TABLE]

Let $y_{x}\in\partial D$ such that $|y_{x}-x|=\mathrm{dist}(x,\Omega\setminus D)$ we have then that

[TABLE]

so that

[TABLE]

which is a contradiction. Now if $z_{x}\in\partial\Omega$, we have

[TABLE]

On the other hand,

[TABLE]

so that,

[TABLE]

Let $u_{x}\in\partial\Omega$ such that $x-u_{x}=\mathrm{dist}(x,\partial\Omega)$ and $x^{\ast}\in\partial D$ such that $x^{\ast}$ lies on the line segment $[x,\,u_{x}]$. Thus,

[TABLE]

hence

[TABLE]

so that

[TABLE]

which contradicts the assumption.

Assume $x\in\Omega\setminus D$ and that $M_{\lambda}(f_{\overline{\Omega},K}^{M})(x)=f_{\overline{\Omega},K}^{M}(z_{x})+\lambda|x-z_{x}|^{2}$ with $z_{x}\in\partial\Omega$. Then, since

[TABLE]

we have that

[TABLE]

thus

[TABLE]

which is a contradiction. Thus $z_{x}\in\Omega\setminus D$.

Part $(ii)$: Assume first $x\in D$ and let $z_{x}$ be such that

[TABLE]

then we need to distinguish the following three cases:

[TABLE]

Case $(a)$: If $z_{x}\in D$, then by Part (i) we conclude that $M_{\lambda}(f_{\overline{\Omega},K}^{M}))(z_{x})$ is determined by $f|_{\Omega\setminus D}$.

Case $(b)$: If $z_{x}\in\partial\Omega$, we have

[TABLE]

and

[TABLE]

for some $u_{x}\in(\Omega\setminus D)\cup\partial\Omega$, thus

[TABLE]

which is a contradiction, thus $z_{x}\not\in\partial\Omega$.

Case $(c)$: Let $z_{x}\in\Omega\setminus D$. We have then that

[TABLE]

and

[TABLE]

for some $u_{x}\in(\Omega\setminus D)\cup\partial\Omega$. Thus

[TABLE]

that is,

[TABLE]

Therefore,

[TABLE]

and from Part $(i)$ we conclude that $M_{\lambda,\Omega}(f_{\overline{\Omega},K}^{M})(z_{x})$ is determined by $f|_{\Omega\setminus D}$.

Now, let us assume that $x\in\Omega\setminus D$, $\mathrm{dist}(x,\partial\Omega)>2\sqrt{\tfrac{O_{f}}{\lambda}}$ and let $z_{x}$ be such that

[TABLE]

then, also in this case, we need to analyze the following three claims:

[TABLE]

Case $(a)$: By Part $(i)$ we conclude that $M_{\lambda,\Omega}(f_{\overline{\Omega},K}^{M})(z_{x})$ is determined by $f|_{\Omega\setminus D}$.

Case $(b)$: As in the similar case seen previusly, we can mke the same arguments and show that $z_{x}\not\in\partial\Omega$.

Case $(c)$: By imitating the arguments made previously for the same case, let $z_{x}\in\Omega\setminus D$ be such that

[TABLE]

and

[TABLE]

for some $u_{x}\in(\Omega\setminus D)\cup\partial\Omega$. Thus,

[TABLE]

Therefore,

[TABLE]

and from Part $(i)$ we conclude that $M_{\lambda,\Omega}(f_{\overline{\Omega},K}^{M})(z_{x})$ is determined by $f|_{\Omega\setminus D}$. ∎

Proof of Theorem 4.2: We have that

[TABLE]

with some $z_{x_{k}}\in\mathbb{R}^{n}$. By Proposition 2.11,

[TABLE]

which yields the following estimate

[TABLE]

Since $\omega_{f}$ is non-decreasing, from (A.23) and (A.24) we obtain

[TABLE]

where we have set $d(\lambda)=\sqrt{\omega_{f}(a/\lambda+\sqrt{b/\lambda})}$. Now, it is not difficult to realize that we can give the following representation for $z_{x_{k}}$ which enters (A.22) in terms of the grid points as follows 111 Let $x_{z_{x_{k}}}$ be for instance the closest grid point to $z_{x_{k}}$ and set $z_{x_{k}}=x_{k}+x_{z_{x_{k}}}-x_{k}+z_{x_{k}}-x_{z_{x_{k}}}$. Noting that $x_{z_{x_{k}}}-x_{k}=hr_{x_{k}}$, $h$ the grid size and $r_{x_{k}}\in\mathbb{Z}^{n}$, and $z_{x_{k}}-x_{z_{x_{k}}}=h\delta$ with $\delta\in\mathbb{R}^{n}$, $0\leq\delta_{i}\leq 1$, $i=1,\ldots,n$, we obtain (A.26).

[TABLE]

$r_{x_{k}}\in\mathbb{Z}^{n}$ and $\delta\in\mathbb{R}^{n}$ such that $|\delta|<\sqrt{n}$. Thus, we have

[TABLE]

where we have used the fact that

[TABLE]

By the Cauchy-Schwartz inequality,

[TABLE]

and since $\omega_{f}$ is non-decreasing, we have that

[TABLE]

As a result, taking into account (A.29) and (A.30) into (A.27) yields

[TABLE]

From (A.26),

[TABLE]

which, combined with (A.25), gives

[TABLE]

which implies that

[TABLE]

Using (A.34) into (A.31), and by taking into account of (A.22), we get

[TABLE]

which concludes the proof. ∎

Proof of Corollary 4.4: The result follows from Theorem 4.2 by taking $\omega_{f}(t)=Lt$. ∎

Proof of Theorem 4.7: For the sake of clarity, in the following, we present the proof only for $n=1$ and $n=2$, given that the proof of the results for $n>2$ build on the arguments required for $n-1$.

Case $n=1$: We prove the statement by induction. Clearly, the statement is true for $m=0,1$ and for all $x_{k}=kh$, $k\in\mathbb{Z}$. Suppose now $f_{m}=g_{m}$ for some $m\geq 1$. Then for $m+1$, we have

[TABLE]

where we have let $A:=f_{m}(x_{k})$, $B:=f_{m}(x_{k}-h)+\lambda\tau_{m+1}h^{2}$ and $C:=f_{m}(x_{k}+h)+\lambda\tau_{m+1}h^{2}$. By the inductive assumption, we have then

[TABLE]

and

[TABLE]

where we have let

[TABLE]

If we change index $j-1=r$ in $B_{1}$, we have

[TABLE]

Now we show that

[TABLE]

as the right hand side of the above inequality is one of the terms to be minimised in $A$ above. This inequality is equivalent to $(r+1)^{2}+\tau_{m+1}\geq r^{2}$ for $-m\leq r\leq m-1$, which is equivalent to $2r+1+2m+1\geq 0$ and thus to $2(m+r)+2\geq 0$. As $-m\leq r\leq m-1$, the last inequality above clearly holds. Thus we have

[TABLE]

Similarly, we can easily see that

[TABLE]

thus,

[TABLE]

It therefore follows that the statement holds for $n=m+1$ and by induction holds for all $n=1,2,\dots$. ∎

Case $n=2$: The proof of the two-dimensional case is very similar to the one-dimensional case but is slightly more involved. Again, we prove the statement $f_{n}=g_{n}$ by induction. For $m=0,1$, the statement is clearly true. Suppose for some $m\geq 1$, we have $f_{m}=g_{m}$. Then for $m+1$ and each $(x_{k},y_{j})$, we have $f_{m+1}(x_{k},y_{j})=\min\{A,\,B\}$, where

[TABLE]

and

[TABLE]

There are eight terms to be minimised in $B$. We consider first two typical cases.

Case $(a)$: $p=0$ and $q=1$. This is similar to the one-dimensional case as we have

[TABLE]

where we have let

[TABLE]

Note that all terms to be minimised in $B_{1}$ are in the definition of $g_{m+1}(x_{k},y_{j})$ but not in $g_{m}(x_{k},y_{j})$. As in the one-dimensional case, every term to be minimised in $C_{1}$ is greater than a corresponding term to be minimised in $A$ as we have seen for the Case $(i)$. Therefore

[TABLE]

Case $(b)$: $p=1$ and $q=-1$. In this case, we have

[TABLE]

In the above, $f_{0}(x_{k+(m+1)},y_{j-(m+1})+2\lambda(m+1)^{2}h^{2}$ is in the definition of $g_{m+1}(x_{k},y_{j})$ which is one of the terms to be minimised that is not in the definition of $g_{m}(x_{k},y_{j})$. We also have

[TABLE]

Again, as in the Case $(i)$, every term involved in the minimisation in $D_{1}$ is greater than a corresponding term in $A$. Therefore, we have

[TABLE]

For the Case $(c)$: $p=0$, $q=-1$; Case $(d)$: $p=1$, $q=0$ and Case $(e)$: $p=-1$, $q=0$, the proofs are similar to that of Case $(a)$ while for the Case $(f)$ $p=1$, $q=1$; Case $(g)$: $p=-1$, $q=-1$ and Case $(h)$: $p=-1$, $q=1$, the proofs are similar to that of Case $(b)$. Therefore, we have

[TABLE]

where

[TABLE]

and

[TABLE]

By comparing the terms involved in $g_{m+1}(x_{k},y_{j})$ and those involved in $A$, $U$ and $V$, we see that

[TABLE]

Thus the statement holds for $n=m+1$. By induction, we conclude then that $f_{n}(x_{k},y_{j})=g_{n}(x_{k},y_{j})$ holds for all $n=0,1,2,\dots$ and for all $x_{k}=kh,\,y_{j}=jh$ with $k,\,j\in\mathbb{Z}$. The proof is finished. ∎

Proof of Proposition 4.8. For any grid point $x_{k}$, we have that there exists a $r_{x_{k}}\in\mathbb{Z}^{n}$ such that

[TABLE]

thus,

[TABLE]

which implies that

[TABLE]

As a result, by taking in (4.4)

[TABLE]

where $\lfloor x\rfloor$ denotes the integer part of $x$, we are guaranteed that $r_{x_{k}}$, which enters (A.57), belongs to the feasible set of (4.4) and is, therefore, attained by Algorithm 1. This concludes the proof. ∎