A total variation based regularizer promoting piecewise-Lipschitz reconstructions

TL;DR

This paper introduces a novel total variation-based regularizer that enforces piecewise-Lipschitz properties in reconstructions, offering an intuitive parameter and spatial adaptivity, with proven regularization properties and competitive performance.

Contribution

A new regularizer in the total variation family promoting spatially adaptive Lipschitz constraints, with theoretical analysis and practical advantages over existing methods.

Findings

Achieves similar performance to total generalized variation.

Provides an intuitive interpretation of the regularizer parameter.

Enables spatially adaptive regularization.

Abstract

We introduce a new regularizer in the total variation family that promotes reconstructions with a given Lipschitz constant (which can also vary spatially). We prove regularizing properties of this functional and investigate its connections to total variation and infimal convolution type regularizers TVLp and, in particular, establish topological equivalence. Our numerical experiments show that the proposed regularizer can achieve similar performance as total generalized variation while having the advantage of a very intuitive interpretation of its free parameter, which is just a local estimate of the norm of the gradient. It also provides a natural approach to spatially adaptive regularization.

Equations

$${{\mathrm{TVL}}^{\mathrm{p}}}_{\alpha ,\beta }(u):=\underset{w\in L^p(\Omega ;{\mathbb{R}}^d)}{\min }\alpha ||Du-w|{|}_{\mathfrak{M}}+\beta ||w|{|}_{L^p(\Omega ;{\mathbb{R}}^d)}.$$

$$\underset{u\in BV(\Omega )}{\min }\frac{1}{2}\Vert u-f{\Vert }_{L^2(\Omega )}^2+\alpha \Vert Du{\Vert }_{\mathfrak{M}},$$

$$\underset{\begin{array}{c}u\in BV(\Omega )\\ g\in \mathfrak{M}(\Omega ,{\mathbb{R}}^2)\end{array}}{\min }\text{leavevmode}\frac{1}{2}\Vert u-f{\Vert }_{L^2(\Omega )}^2+\alpha \Vert g{\Vert }_{\mathfrak{M}}s.t.\text{leavevmode}Du=g.$$

$$\underset{\begin{array}{c}u\in BV(\Omega )\\ g\in \mathfrak{M}(\Omega ,{\mathbb{R}}^2)\end{array}}{\min }\text{leavevmode}\frac{1}{2}\Vert u-f{\Vert }_{L^2(\Omega )}^2+\alpha \Vert g{\Vert }_{\mathfrak{M}}s.t.\text{leavevmode}|Du-g|⩽\gamma$$

$$\underset{\begin{array}{c}u\in BV(\Omega )\\ g\in \mathfrak{M}(\Omega ,{\mathbb{R}}^2)\end{array}}{\min }\text{leavevmode}\frac{1}{2}\Vert u-f{\Vert }_{L^2(\Omega )}^2+\alpha \Vert Du-g{\Vert }_{\mathfrak{M}}s.t.\text{leavevmode}|g|⩽\gamma ,$$

$${\mathrm{TV}}_{pwL}^{\gamma }(u):=\underset{g\in \mathfrak{M}(\Omega ,{\mathbb{R}}^2)}{\inf }\Vert Du-g{\Vert }_{\mathfrak{M}}\text{s.t. ∣g∣⩽γ}$$

$$|g|(E):=\sup \left\{\sum _{i=1}^{\infty }\Vert g(E_i){\Vert }_2\mid E=\bigcup _{i\in \mathbb{N}}{E}_i,\text{Ei pairwise disjoint}\right\}$$

$${\mathrm{TV}}_{pwL}^{\gamma }(u)\text{leavevmode}=\text{leavevmode}\underset{\begin{array}{c}\phi \in {\mathcal{C}}_0^{\infty }(\Omega ;{\mathbb{R}}^2)\\ |\phi {|}_2⩽1\end{array}}{\sup }\text{leavevmode}\left\{{\int }_{\Omega }u\text{leavevmode}\mathrm{div}\text{leavevmode}\phi \text{leavevmode}dx-{\int }_{\Omega }|\phi {|}_{2}d\gamma \right\},$$

$${\mathrm{TV}}_{pwL}^{\gamma }(u)=\underset{\begin{array}{c}g\in \mathfrak{M}(\Omega ,{\mathbb{R}}^2)\\ |g|⩽\gamma \end{array}}{\inf }\Vert Du-g{\Vert }_{\mathfrak{M}}=\underset{\begin{array}{c}g\in \mathfrak{M}(\Omega ,{\mathbb{R}}^2)\\ |g|⩽\gamma \end{array}}{\inf }\underset{\begin{array}{c}\phi \in {\mathcal{C}}_0(\Omega ;{\mathbb{R}}^2)\\ |\phi {|}_2⩽1\end{array}}{\sup }(Du-g,\phi ).$$

$${\mathrm{TV}}_{pwL}^{\gamma }(u)=\underset{\begin{array}{c}\phi \in {\mathcal{C}}_0(\Omega ;{\mathbb{R}}^2)\\ |\phi {|}_2⩽1\end{array}}{\sup }\underset{\begin{array}{c}g\in \mathfrak{M}(\Omega ,{\mathbb{R}}^2)\\ |g|⩽\gamma \end{array}}{\inf }(Du-g,\phi )$$

$$=\underset{\begin{array}{c}\phi \in {\mathcal{C}}_0(\Omega ;{\mathbb{R}}^2)\\ |\phi {|}_2⩽1\end{array}}{\sup }[(Du,\phi )-\underset{\begin{array}{c}g\in \mathfrak{M}(\Omega ,{\mathbb{R}}^2)\\ |g|⩽\gamma \end{array}}{\sup }(g,\phi )]=\underset{\begin{array}{c}\phi \in {\mathcal{C}}_0(\Omega ;{\mathbb{R}}^2)\\ |\phi {|}_2⩽1\end{array}}{\sup }[(Du,\phi )-(\gamma ,|\phi {|}_2)].$$

$${\mathrm{TV}}_{pwL}^{\gamma }(u)=\underset{\begin{array}{c}\phi \in {\mathcal{C}}_0^{\infty }(\Omega ;{\mathbb{R}}^2)\\ |\phi {|}_2⩽1\end{array}}{\sup }[(u,-\mathrm{div}\phi )-(\gamma ,|\phi {|}_2)]$$

$$TV(u)-\gamma (\Omega )⩽{\mathrm{TV}}_{pwL}^{\gamma }(u)⩽TV(u).$$

$$\Vert Du-g{\Vert }_{\mathfrak{M}}⩾\Vert Du{\Vert }_{\mathfrak{M}}-\Vert g{\Vert }_{\mathfrak{M}}⩾\Vert Du{\Vert }_{\mathfrak{M}}-\Vert \gamma {\Vert }_{\mathfrak{M}}=\mathrm{TV}(u)-\gamma (\Omega ),$$

$${\mathrm{TV}}_{pwL}^{\gamma }(u)=\Vert Du-\tilde{g}{\Vert }_{\mathfrak{M}}.$$

$$|Du|(A)={\int }_{A}f\text{leavevmode}d\mathcal{L}$$

$$\Vert u{\Vert }_{\phi }=\sup \left\{\lambda >0:\int \phi (|u|/\lambda )d\mu ⩽1\right\}.$$

$$T{V}_{pwL}^{\gamma }(u)$$

$$\underset{u\in {\mathbb{R}}^N}{\min }\mathcal{J}(u)\text{s.t. ∥u−f∥22⩽σ2⋅N},$$

$$\underset{u\in {\mathbb{R}}^N}{\min }\frac{1}{2}\Vert u-f{\Vert }_2^2+\alpha \mathrm{TV}(u)$$

Full text

1 Introduction

Since it has been introduced in [20], total variation ($\operatorname{TV}$) has been popular in image processing due to its ability to preserve edges while imposing sufficient regularity on the reconstructions. There have been numerous works studying the geometric structure of $\operatorname{TV}$-based reconstructions (e.g., [19, 17, 11, 16, 8]). A typical characteristic of these reconstructions is the so-called staircasing [19, 17], which refers to piecewse-constant reconstructions with jumps that are not present in the ground truth image. To overcome the issue of staircasing, many other $\operatorname{TV}$-type regularizers have been proposed, perhaps the most successful of which being the Total Generalized Variation ($\operatorname{TGV}$) [3, 5]. $\operatorname{TGV}$ uses derivatives of higher order and favours reconstructions that are piecewise-polynomial; in the most common case of $\operatorname{TGV}^{2}$ these are piecewise-affine.

While $\operatorname{TGV}$ greatly improves the reconstruction quality compared to $\operatorname{TV}$, the fact that it uses second order derivatives typically results in slower convergence of iterative optimization algorithms and therefore increases computational costs of the reconstruction. Therefore, there has been an effort to achieve a performance similar to that of $\operatorname{TGV}$ with a first-order method (i.e. a method that only uses derivatives of the first order). In [9, 10], infimal convolution type regularizers $\operatorname{\operatorname{TVL}^{p}}$ have been introduced that use an infimal convolution of the Radon norm and an $L^{p}$ norm applied to the weak gradient of the image. For an $u\in L^{1}(\Omega)$, $\operatorname{\operatorname{TVL}^{p}}$ is defined as follows

[TABLE]

where $D$ is the weak gradient, $\alpha,\beta>0$ are constants and $1<p\leqslant\infty$. It was shown that for $p=2$ the reconstructions are piecewise-smooth, while for $p=\infty$ they somewhat resemble those obtained with $\operatorname{TGV}$.

The regularizer we introduce in the current paper also aims at achieving a similar performance with second-order methods while only relying on first order derivatives. It can be seen either as a relaxiation of $\operatorname{TV}$ obtained by extending its kernel from constants to all functions with a given Lipschitz constant (for this reason, we call this new regularizer $\operatorname{TV}_{pwL}$, with ‘$pwL$’ standing for ‘piecewise-Lipschitz’), or as an infimal convolution type regularizer, where the Radon norm is convolved with the characteristic function of a certain convex set.

We start with the following motivation. Let $\Omega\subset\mathbb{R}^{2}$ be a bounded Lipschitz domain and $f\in L^{2}(\Omega)$ a noisy image. Recall the ROF [20] denoising model

[TABLE]

where $D\colon L^{1}(\Omega)\to\mathfrak{M}(\Omega,\mathbb{R}^{2})$ is the weak gradient, $\mathfrak{M}(\Omega,\mathbb{R}^{2})$ is the space of vector-valued Radon measures and $\alpha>0$ is the regularization parameter. Introducing an auxiliary variable $g\in\mathfrak{M}(\Omega,\mathbb{R}^{2})$, we can rewrite this problem as follows

[TABLE]

Our idea is to relax the constraint on $Du$ as follows

[TABLE]

for some positive constant, function or measure $\gamma$. Here $|Du-g|$ is the variation measure corresponding to $Du-g$ and the symbol $"\leqslant"$ denotes a partial order in the space of signed (scalar valued) measures $\mathcal{M}(\Omega)$. This problem is equivalent to

[TABLE]

which we take as the starting point of our approach.

This paper is organized as follows. In Section 2 we introduce the primal and dual formulations of the $\operatorname{TV}_{pwL}$ functional and prove their equivalence. In Section 3 we prove some basic properties of $\operatorname{TV}_{pwL}$ and study its relationship with other $\operatorname{TV}$-type regularizers. Section 4 contains numerical experiments with the proposed regularizer.

2 Primal and Dual Formulations

Let us first clarify the notation of the inequality for signed measures in (1).

Definition 1

We call a measure $\mu\in\mathcal{M}(\Omega)$ positive if for every subset $E\subseteq\Omega$ one has $\mu(E)\geqslant 0$. For two signed measures $\mu_{1},\mu_{2}\in\mathcal{M}(\Omega)$ we say that $\mu_{1}\leqslant\mu_{2}$ if $\mu_{2}-\mu_{1}$ is a positive measure.

Now let us formally define the new regularizer.

Definition 2

Let $\Omega\subset\mathbb{R}^{2}$ be a bounded Lipschitz domain, $\gamma\in\mathcal{M}(\Omega)$ be a finite positive measure. For any $u\in L^{1}(\Omega)$ we define

[TABLE]

where $||\cdot||_{\mathfrak{M}}$ denotes the Radon norm and $|g|$ is the variation measure [4] corresponding to $g$, i.e. for any subset $E\subset\Omega$

[TABLE]

(see also the polar decomposition of measures [1]).

The $\inf$ in Definition 2 can actually be replaced by a $\min$ since we are dealing with a metric projection onto a closed convex set in the dual of a separable Banach space. We also note that for $\gamma=0$ we immediately recover total variation.

We defined the regularizer $\operatorname{TV}_{pwL}^{\gamma}$ in the general case of $\gamma$ being a measure, but we can also choose $\gamma$ to be a Lebesgue-measurable function or a constant. In this case the inequality is understood in the sense $|g|\leqslant\gamma\leavevmode\nobreak\ d\mathcal{L}$ where $\mathcal{L}$ is the Lebesgue measure, resulting in $|g|$ being absolutely continuous with respect to $\mathcal{L}$. We will not distinguish between these cases in what follows and just write $|g|\leqslant\gamma$.

As with standard $\operatorname{TV}$, there is also an equivalent dual formulation of $\operatorname{TV}_{pwL}$.

Theorem 2.1

Let $\gamma\in\mathcal{M}(\Omega)$ be a positive finite measure and $\Omega$ a bounded Lipschitz domain. Then for any $u\in L^{1}(\Omega)$ the $\operatorname{TV}_{pwL}^{\gamma}$ functional can be equivalently expressed as follows

[TABLE]

where $|\varphi|_{2}$ denotes the pointwise $2$-norm of $\varphi$.

Proof

Since by the Riesz-Markov-Kakutani representation theorem the space of vector valued Radon measures $\mathfrak{M}(\Omega,\mathbb{R}^{2})$ is the dual of the space $\mathcal{C}_{0}(\Omega,\mathbb{R}^{2})$, we can rewrite the expression in Definition 2 as follows

[TABLE]

In order to exchange $\inf$ and $\sup$, we need to apply a minimax theorem. In our setting we can use the Nonsymmetrical Minimax Theorem from [2, Th. 3.6.4]. Since the set $\{g\leavevmode\nobreak\ |\leavevmode\nobreak\ |g|\leqslant\gamma\}\subset\mathfrak{M}(\Omega,\mathbb{R}^{2})=(\mathcal{C}_{0}(\Omega,\mathbb{R}^{2}))^{*}$ is bounded, convex and closed and the set $\{\varphi\leavevmode\nobreak\ |\leavevmode\nobreak\ \|\varphi\|_{2,\infty}\leqslant 1\}\subset\mathcal{C}_{0}(\Omega,\mathbb{R}^{2})$ is convex, we can swap the infimum and the supremum and obtain the following representation

[TABLE]

Noting that the supremum can actually be taken over $\varphi\in\mathcal{C}_{0}^{\infty}(\Omega;\mathbb{R}^{2})$, we obtain

[TABLE]

which yields the assertion upon replacing $\varphi$ with $-\varphi$.

3 Basic Properties and Relationship with other $\operatorname{TV}$-type regularizers

Influence of $\gamma$.

It is evident from Definition 2 that a larger $\gamma$ yields a larger feasible set and a smaller value of $\operatorname{TV}_{pwL}$. Therefore, $\operatorname{TV}_{pwL}^{\gamma}\geqslant\operatorname{TV}_{pwL}^{\bar{\gamma}}$ whenever $0\leqslant\gamma\leqslant\bar{\gamma}$. In particular, we get that $\operatorname{TV}_{pwL}^{\gamma}\leqslant\operatorname{TV}_{pwL}^{0}=\operatorname{TV}$ for any $\gamma\geqslant 0$.

Lower-semicontinuity and convexity.

Lower-semicontinuity is clear from Definition 2 if we recall that the infimum is actually a minimum. Convexity follows from the fact that $\operatorname{TV}_{pwL}$ is an infimal convolution of two convex functions.

Absolute one-homogeneity.

Noting that $\operatorname{TV}_{pwL}$ is the distance from the convex set $\{g\in\mathfrak{M}(\Omega;\mathbb{R}^{2})\leavevmode\nobreak\ |\leavevmode\nobreak\ |g|\leqslant\gamma\}$, we conclude that it is absolute one-homogeneous if and only if this set consists of just zero, i.e. when $\gamma=0$ and $\operatorname{TV}_{pwL}=\operatorname{TV}$.

Coercivity.

We have seen that $\operatorname{TV}_{pwL}^{\gamma}\leqslant\operatorname{TV}$ for any $\gamma\geqslant 0$, i.e. $\operatorname{TV}_{pwL}$ is a lower bound for $\operatorname{TV}$. If $\gamma(\Omega)$ is finite, the converse inequality (up to a constant) also holds and we obtain topological equivalence of $\operatorname{TV}_{pwL}$ and $\operatorname{TV}$.

Theorem 3.1

Let $\Omega\subset\mathbb{R}^{2}$ be a bounded Lipschitz domain and $\gamma\in\mathcal{M}(\Omega)$ a positive finite measure. For every $u\in L^{1}(\Omega)$ we obtain the following relation:

[TABLE]

Proof

We already established the right inequality. For the left one we observe that for any $g\in\mathfrak{M}(\Omega,\mathbb{R}^{2})$ such that $|g|\leqslant\gamma$ the following estimate holds

[TABLE]

which also holds for the infimum over $g$.

The left inequality in Theorem 3.1 ensures that $\operatorname{TV}_{pwL}$ is coercive on $\operatorname{BV}_{0}$, since $\operatorname{TV}(u_{n})\to\infty$ implies $\operatorname{TV}_{pwL}^{\gamma}(u_{n})\geqslant\operatorname{TV}(u_{n})-\gamma(\Omega)\to\infty$. Upon adding the $L^{1}$ norm, we also get coercivity on $\operatorname{BV}$. This ensures that $\operatorname{TV}_{pwL}$ can be used for regularisation of inverse problems in the same scenarios as $\operatorname{TV}$.

Topological equivalence between $\operatorname{TV}_{pwL}$ and $\operatorname{TV}$ is understood in the sense that if one is bounded then the other one is too. Being not absolute one-homogeneous, however, $\operatorname{TV}_{pwL}$ cannot be an equivalent norm on $\operatorname{BV}_{0}$.

Null space.

We will study the null space in the case when $\gamma\in L^{\infty}_{+}(\Omega)$ is a Lebesgue measurable function and the inequality $|g|\leqslant\gamma$ in Definition 2 is understood in the sense that $|g|\leqslant\gamma d\mathcal{L}$ with $\mathcal{L}$ being the Lebesgue measure.

Proposition 1

Let $u\in L^{1}(\Omega)$ and $\gamma\geqslant 0$ be an $L^{\infty}$ function. Then $\operatorname{TV}_{pwL}^{\gamma}(u)=0$ if and only if the weak derivative $Du$ is absolutely continuous with respect to the Lebesgue measure and its $2$-norm is bounded by $\gamma$ a.e.

Proof

We already noted that, since the space $\mathfrak{M}(\Omega;\mathbb{R}^{2})$ is the dual of the separable Banach space $C_{0}(\Omega,\mathbb{R}^{2})$, the infimum in Definition 2 is actually a minimum and there exists a $\tilde{g}\in\mathfrak{M}(\Omega;\mathbb{R}^{2})$ with $|g|\leqslant\gamma d\mathcal{L}$ such that

[TABLE]

If $\operatorname{TV}_{pwL}^{\gamma}(u)=0$, then $Du=\tilde{g}$ and therefore $0\leqslant|Du|\leqslant\gamma\leavevmode\nobreak\ d\mathcal{L}$, which implies that $|Du|$ is absolutely continuous with respect to $\mathcal{L}$ and can be written as

[TABLE]

for any $A\subseteq\Omega$. The function $f\in L^{1}_{+}(\Omega)$ is the Radon-Nikodym derivative $\frac{d|Du|}{d\mathcal{L}}$. From the condition $|Du|\leqslant\gamma\leavevmode\nobreak\ d\mathcal{L}$ we immediately get that $f\leqslant\gamma$ a.e.

Remark 1

We notice that the set $\{u\in L^{1}(\Omega)\colon\operatorname{TV}_{pwL}(u)=0\}$ is not a linear subspace, therefore, we should rather speak of the null set than the null space.

Remark 2

Proposition 1 implies that all functions in the null set are Lipschitz continuous with (perhaps, spatially varying) Lipschitz constant $\gamma$; hence the name of the regularizer.

Luxemburg norms.

For a positive convex nondecreasing function $\varphi\colon\mathbb{R}_{+}\to\mathbb{R}_{+}$ with $\varphi(0)=0$ the Luxemburg norm $\|\cdot\|_{\varphi}$ is defined as follows [18]

[TABLE]

We point out a possible connection between $\operatorname{TV}_{pwL}$ and a Luxemburg norm corresponding to $\varphi(x)=(x-c)_{+}$ with a suitable constant $c>0$. However, we do not investigate this connection in this paper.

Relationship with Infimal Convolution Type regularizers.

We would like to highlight a relationship to infimal convolution type regularizers $\operatorname{\operatorname{TVL}^{p}}$ [9, 10]. Indeed, as already noticed earlier, $\operatorname{TV}_{pwL}$ can be written as an infimal convolution

[TABLE]

where $C_{\gamma}:=\{\eta\in\mathfrak{M}(\Omega;\mathbb{R}^{2})\leavevmode\nobreak\ |\leavevmode\nobreak\ |\eta|_{p}\leqslant\gamma\}$. If $\gamma>0$ is a constant, we obtain a bound on the $2,\infty$-norm of $g$. This highlights a connection to the $\operatorname{TVL}^{\infty}$ regularizer [10]: for any particular weight in front of the $\infty$-norm in $\operatorname{TVL}^{\infty}$ (and for a given $u$), the auxiliary variable $g$ will have some value of the $\infty$-norm and if we use this value as $\gamma$ in $\operatorname{TV}_{pwL}^{\gamma}$, we will obtain the same reconstruction. The important difference is that with $\operatorname{TVL}^{\infty}$ we don’t have direct control over this value and can only influence it in an indirect way through the weights in the regularizer. In $\operatorname{TV}_{pwL}$ this parameter is given explicitly and can be either obtained using additional a priori information about the ground truth or estimated from the noisy image.

Similar arguments can be made in the case of spatially variable $\gamma$ if the weighting in $\operatorname{TVL}^{\infty}$ is also allowed to vary spatially.

4 Numerical Experiments

In this section we want to compare the performance of the proposed first order regularizer with a second order regularizer, $\operatorname{TGV}$, in image denosing. We consider images (or 1D signals) corrupted by Gaussian noise with a known variance and use the residual method [14] to reconstruct the noise-free image, i.e. we solve (in the discrete setting)

[TABLE]

where $f$ is the noisy image, $\mathcal{J}$ is the regularizer ($\operatorname{TV}_{pwL}$ or $\operatorname{TGV}$), $\sigma$ is the standard deviation of the Gaussian noise and $N$ is the number of pixels in the image (or the number of entries in the 1D signal). We solve all problems in MATLAB using CVX [15]. For $\operatorname{TGV}$ we use the parameter $\beta=1.25$, which is in the range $[1,1.5]$ recommended in [12].

A characteristic feature of the proposed regularizer $\operatorname{TV}_{pwL}$ is its ability to efficiently encode the information about the gradient of the ground truth (away from jumps) if such information is available. Our experiments showed that the quality of the reconstruction significantly depends on the quality of the (local) estimate of the norm of the gradient of the ground truth.

The ideal application for $\operatorname{TV}_{pwL}$ would be one where we have a good estimate of the gradient of the ground truth away from jumps, which, however, may occur at unknown locations and be of unknown magnitude. If such an estimate is not available, we can roughly estimate the gradient of the ground truth from the noisy signal, which is the approach we take.

4.1 1D experiments.

We consider the ground truth shown in Figure 1a (green dashed line). The signal is discretized using $N=1000$ points. We add Gaussian noise with variance $\sigma=0.1$ and obtain the noisy signal shown in the same Figure (blue solid line).

To use $\operatorname{TV}_{pwL}$, we need to estimate the derivative of the true signal away from the jumps. Therefore, we need to detect the jumps, but leave the signal intact away from them (up to a constant shift). This is exactly what happens if the image is overregularized with $\operatorname{TV}$. We compute a $\operatorname{TV}$ reconstruction by solving the ROF model

[TABLE]

with a large value of $\alpha$ (in this example we took $\alpha=0.5$). The result is shown in Figure 1a (red solid line). The residual, which we want to use to estimate the derivative of the ground truth, is shown in Figure 1b. Although the jumps have not been removed entirely, this signal can be used to estimate the derivative using filtering.

The filtered residual (we used the build-in MATLAB function ’smooth’ with option ’rlowess’ (robust weighted linear least squares) and $\omega=50$) is shown in Figure 1c. This signal is sufficiently smooth to be differentiated. We use central differences; to suppress the remaining noise in the filtered residual we use a step size for differentiation that is $20$ times the original step size. The result is shown in Figure 1d (reg solid line) along with the true derivative (green dashed line). We use the absolute value of the so computed derivative as the parameter $\gamma$.

The reconstruction obtained using $\operatorname{TV}_{pwL}^{\gamma}$ is shown in Figure 1e. We see that the reconstruction is best in areas where our estimate of the true derivative was most faithful (e.g., between $4$ and $5$). But also in other areas the reconstruction is good and preserves the structure of the ground truth rather well. We notice a small artefact at the value of the argument of around $2$; examining the estimate of the derivative in Figure 1d and the residual in Figure 1b, we notice that $\operatorname{TV}$ was not able the remove the jump at this location and therefore the estimate of the derivative was too large. This allowed the reconstruction to get too close to the data at this point.

We also notice that the jumps are sometimes reduced, with a characteristic linear cut near the jump (e.g., near $x=1;2;3$ and $4$). This can have different reasons. For the jumps near $x=3$ and $4$ we see that the estimate of $\gamma$ is too large (the true derivative is zero), which allows the regulariser to cut the edges. For the jumps near $x=1$ and $2$ the situation is different. At these positions a negative slope in the ground truth is followed by a positive jump. Since $\gamma$ only constraints the absolute value of the gradient, even with a correct estimate of $\gamma$ the regulariser will reduce the jump, going with the maximum slope in the direction of the jump. Functions with a negative slope followed by a positive jump are also problematic for $\operatorname{TV}$, since they do not satisfy the source condition (their subdifferential is empty [7]). In such cases $\operatorname{TV}$ will also always reduce the jump.

Figure 1e shows the reconstruction obtained with $\operatorname{TGV}$. The reconstruction is quite good, although it is often piecewise-affine where the ground truth is not, e.g. between $4$ and $5$ or between $5$ and $7$. As expected, both regularizers tend to push the reconstructions towards their kernels, but, since $\operatorname{TV}_{pwL}^{\gamma}$ with a good choice of $\gamma$ contains the ground truth in its kernel (up to the jumps), it yields reconstructions that are more faithful to the structure of the ground truth.

4.2 2D experiments.

In this Section we study the performance of $\operatorname{TV}_{pwL}$ in denoising of 2D images. We use two images - “cameraman” (Figure 2a) and “owl” (Figure 3a). Both images have the resolution $256\times 256$ pixels and values in the interval $[0,255]$. The images are corrupted with Gaussian noise with standard deviation $\sigma=0.1\cdot 255=25.5$ (Figures 2b and 3b).

The pipeline for the reconstruction using $\operatorname{TV}_{pwL}$ is the same as in 1D. We obtain a piecewise constant image by solving an overregularized ROF problem (4) with $\alpha=500$ (Figures 2c and 3c) and compute the residuals (Figures 2d and 3d). Then we smooth the residuals using a Gauss filter with $\sigma=2$ (Figures 2e and 3e) and compute its derivatives in the $x$- and $y$-directions using the same approach as in 1D (central differences with a different step size; we used step size $3$ in this example). These derivatives are used to estimate $\gamma$, which is set equal to the norm of the gradient. Figures 2f and 3f show $\gamma$ scaled to the interval $[0,255]$ for better visibility. We use the same parameters (for Gaussian filtering and numerical differentiation) to estimate $\gamma$ in both images. Reconstructions obtained using $\operatorname{TV}_{pwL}^{\gamma}$ and $\operatorname{TGV}$ are shown in Figures 2g-h and 3g-h. As in the 1D example, the parameter $\beta$ for $\operatorname{TGV}$ was set to $1.25$.

Comparing the results for both images, we notice that the residual (as well as its filtered version) captures the details in the “cameraman” image much better than in the “owl” image. The filtered residual in the “owl” image seems to miss some of the structure of the original image and this is reflected in the estimated $\gamma$ (which looks much noisier in the “owl” image and, in particular, does not capture the structure of the needles of the pine tree in the upper left corner). This might be due to the segmentation achieved by $\operatorname{TV}$, which seems better in the “cameraman” image (the one in the “owl” image seems to a bit too detailed). This effect might be mitigated by using a better segmentation technique.

This difference is reflected in the reconstructions. While in the “cameraman” image the details are well preserved (e.g., the face of the cameraman or his camera, as well as the details of the background), in the “owl” image part of them are lost and replaced by rather blurry (if not constant) regions; however, in other regions, such as the feathers of the owl, the details are preserved much better, which can be also seen from the estimated $\gamma$ that is much more regular in this area and closer to the structure of the ground truth. We also notice some loss of contrast in the $\operatorname{TV}_{pwL}$ reconstruction. Perhaps, it could be dealt with by adopting the concept of debiasing [13, 6] in the setting of $\operatorname{TV}_{pwL}$, however, it is not clear yet, what is the structure of the model manifold in this case.

The $\operatorname{TGV}$ reconstructions look strikingly similar to those obtained by $\operatorname{TV}_{pwL}$. Structural similarity between these reconstructions (i.e. SSIM computed using on of them as the reference) is $0.98$ for the “cameraman” image and $0.97$ for the “owl” image. Although the $\operatorname{TGV}$ reconstructions depend on the parameter $\beta$ any may differ more from $\operatorname{TV}_{pwL}$ for other values of $\beta$, the one we chose here ($\beta=1.25$) is reasonable and lies within the optimal range reported in [12].

There are two main messages to be taken from these experiments. The first one is that $\operatorname{TV}_{pwL}$ is able to almost reproduce the reconstructions obtained using $\operatorname{TGV}$ with a reasonable choice of the parameter $\beta$, which is a very good performance for a method that does not use higher-order derivatives. The second one is that the performance of $\operatorname{TV}_{pwL}$ greatly depends on the quality of the estimate of $\gamma$. When we are able to well capture the structure of the original image in this estimate, the structure of the reconstructions is rather close to that of the ground truth. To further illustrate this point, we show in Fig. 4 $\operatorname{TV}_{pwL}^{\gamma}$ reconstructions in the ideal scenario when $\gamma$ is estimated from the ground truth as the local magnitude of the gradient. The quality of the reconstructions is very good, suggesting that with a better strategy of estimating the gradient $\operatorname{TV}_{pwL}^{\gamma}$ could achieve even better performance.

5 Conclusions

We proposed a new $\operatorname{TV}$-type regularizer that can be used to decompose the image into a jump part and a part with Lipschitz continuous gradient (with a given Lipschitz constant that is also allowed to vary spatially). Functions whose gradient does not exceed this constant lie in the kernel of the regularizer and are not penalized. By smartly choosing this bound we can hope to put the ground truth into the kernel (up to the jumps) and thus not penalize any structure that is present in the ground truth.

In this paper we presented, in the context of denoising, an approach to estimating this bound from the noisy image. The approach is based on segmenting the image (and compensating for the jumps) using overregularized $\operatorname{TV}$ and estimating the local bound on the gradient using filtering. Our numerical experiments showed that $\operatorname{TV}_{pwL}$ can produce reconstructions that are very similar to $\operatorname{TGV}$, however, the results significantly depend on the quality of the estimation of the local bound on the gradient. Using a more sophisticated estimation technique is expected to further improve the reconstructions. The ideal application for $\operatorname{TV}_{pwL}$ would be one where there is some information about the magnitude of the Lipschitz part of the gradient.

Acknowledgments

This work was supported by the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 777826 (NoMADS). MB acknowledges further support by ERC via Grant EU FP7 – ERC Consolidator Grant 615216 LifeInverse. YK acknowledges support of the Royal Society through a Newton International Fellowship. YK also acknowledges support of the Humbold Foundataion through a Humbold Fellowship he held at the University of Münster when this work was initiated. CBS acknowledges support from the Leverhulme Trust project on Breaking the non-convexity barrier, EPSRC grant Nr. EP/M00483X/1, the EPSRC Centre Nr. EP/N014588/1, the RISE projects CHiPS and NoMADS, the Cantab Capital Institute for the Mathematics of Information and the Alan Turing Institute. We gratefully acknowledge the support of NVIDIA Corporation with the donation of a Quadro P6000 and a Titan Xp GPUs used for this research.