A class of second-order geometric quasilinear hyperbolic PDEs and their application in imaging science

TL;DR

This paper introduces second-order geometric quasilinear hyperbolic PDEs inspired by image processing applications, demonstrating their advantages over first-order methods through theoretical analysis and numerical comparisons.

Contribution

It develops and analyzes a new class of second-order PDEs for image denoising and correction, including existence, uniqueness, and long-term behavior results.

Findings

Second-order PDEs outperform first-order flows in numerical experiments.

Existence and uniqueness of solutions are established for both PDE models.

Analytical solutions are derived for simplified initial data cases.

Abstract

Motivated by important applications in image processing, we study a class of second-order geometric quasilinear hyperbolic partial differential equations (PDEs). This is inspired by the recent development of second-order damping systems associated to gradient flows for energy decaying. In numerical computations, it turns out that the second-order methods are superior to their first-order counter-parts. We concentrate on (i) a damped second-order total variation flow for e.g., image denoising, and (ii) a damped second-order mean curvature flow for level sets of scalar functions. The latter is connected to a non-convex variational model capable of correcting displacement errors in image data (e.g. dejittering). For the former equation, we prove the existence and uniqueness of the solution and its long time behavior, and provide an analytical solution given some simple initial datum. For the latter, we draw a connection between the equation and some second-order geometric PDEs evolving the hypersurfaces, and show the existence and uniqueness of the solution for a regularized version of the equation. Finally, some numerical comparison of the solution behavior for the new equations with first-order flows are presented.

Tables

Table 1. Table 1: Notations and abbreviations.

Notation	Description
$f\left(x\right)$	with only spatial variable $x$, $f:\Omega \to {ℝ}^n$ is time-independent
$f\left(t\right)$	with only temporal variable $t$, $f$ maps to either real values,
	or elements in some function spaces that is $f\left(t\right)=f(\cdot ,t)$
$f(x,t)$	$f$ is both space- and time-dependent
$\dot{f}$ & $\ddot{f}$	first-order and second-order time derivative of $f$, respectively
$\nabla f$	spatial gradient (including distributional sense) of function $f$
$\mathrm{div}\left(𝐯\right)$	spatial divergence of the vector field $𝐯$
$(\cdot ,\cdot )$	inner product of two elements in Hilbert spaces (mostly $L^2$ space)
$⟨\cdot ,\cdot ⟩$	inner product of two elements in ${ℝ}^n$

Equations

$$\dot{w}=-\partial \Phi (w),$$

$$\ddot{w}+\eta (t)\dot{w}=-\partial \Phi (w),$$

$$\underset{u}{\min }\frac{1}{2}{\int }_{\Omega }{\left|u(x)-u^d(x)\right|}^{2}dx+\alpha TV(u).$$

$$\Vert Du\Vert :=\sup \left\{{\int }_{\Omega }u(x)\mathrm{div}(\mathbf{v}(x))dx:\mathbf{v}\in C_0^{\infty }(\Omega ),\left|\mathbf{v}(x)\right|\le 1\text{ for all }x\in \Omega \right\},$$

$$BV(\Omega ):=\left\{u\in L^1(\Omega ):{\Vert u\Vert }_{BV}<+\infty \right\},\text{ where }{\Vert u\Vert }_{BV}:={\Vert u\Vert }_{L^1}+\Vert Du\Vert .$$

$$\Vert Du\Vert ={\int }_{\Omega }\left|\nabla u(x)\right|dx,\text{ for }u\in W_1^1(\Omega ).$$

$$TV(u):=\left\{\begin{array}[]{ll}\left\|Du\right\|,&u\in L^{2}(\Omega)\bigcap BV(\Omega);\\ \infty,&u\in L^{2}(\Omega)\setminus BV(\Omega).\end{array}\right.$$

$$\min \left\{TV(u):u\in L^2(\Omega )\bigcap BV(\Omega )\right\},$$

$$\left\{\begin{array}[]{ll}\dot{u}(t)=\operatorname{div}\left(\frac{\nabla u(t)}{\left|\nabla u(t)\right|}\right),&\text{ in }\Omega\times(0,\infty)\\ u(0)=u_{0},&\text{ in }\Omega\times 0.\end{array}\right.$$

$$p(t)=\left\{\begin{array}[]{ll}\frac{\nabla u(t)}{\left|\nabla u(t)\right|},&\text{ if }\quad\nabla u(t)\neq 0;\\ \gamma(t)\quad\text{ for some }\left|\gamma(t)\right|\leq 1,&\text{ if }\quad\nabla u(t)=0,\end{array}\right.$$

$$u_m-u_{m-1}\in \Delta t\partial TV(u_m),\text{ for }m\in \mathbb{N}.$$

$$\left\{\begin{array}[]{ll}\ddot{u}(t)+\eta\dot{u}(t)-\operatorname{div}\left(\frac{\nabla u(t)}{\left|\nabla u(t)\right|}\right)=0,&\textmd{in}\ \Omega\times(0,\infty),\\ u(0)=u_{0},\quad\dot{u}(0)=0,&\textmd{in}\ \Omega\times 0,\\ \partial_{\nu}u(t)=0,&\textmd{on}\ \partial\Omega\times(0,\infty),\end{array}\right.$$

$$\ddot{u}(t)+\eta \dot{u}(t)+\partial TV(u(t))\ni 0\text{ for all }t>0.$$

$$u+\lambda \partial TV(u)\ni w.$$

$$A_{\lambda }(w):=(w-J_{\lambda }(w))/\lambda .$$

$$\Vert A_{\lambda }(w_1)-A_{\lambda }(w_2)\Vert \le \frac{2}{\lambda }\Vert w_1-w_2\Vert \text{ for all }{w}_1,w_2\in L^2(\Omega ).$$

$$\underset{\lambda >0}{\sup }\Vert A_{\lambda }(w)\Vert \le \left|(\partial TV{)}^0(w)\right|:=\underset{v\in \partial TV(w)}{\min }\Vert v\Vert .$$

$$\underset{\lambda \to 0}{\lim }{J}_{\lambda }(w)=w.$$

$$(\beta ,u)=TV(u).$$

$$\left\{\begin{array}[]{ll}\ddot{u}_{\lambda}(t)+\eta\dot{u}_{\lambda}(t)+A_{\lambda}(u_{\lambda}(t))=0,&\textmd{in}\ \Omega\times(0,\infty),\\ u_{\lambda}(0)=u_{0},\quad\dot{u}_{\lambda}(0)=0,&\textmd{in}\ \Omega\times 0,\\ \partial_{\nu}u_{\lambda}(t)=0,&\textmd{on}\ \partial\Omega\times(0,\infty).\end{array}\right.$$

$$([u_1,v_1{]}^{\top },[u_2,v_2{]}^{\top }{)}_{H\times H}=(u_1,u_2)+(v_1,v_2),$$

$$\left\{\begin{array}[]{ll}\dot{\mathbf{z}}(t)=F(\mathbf{z}(t))&\textmd{in }\ \Omega\times(0,\infty),\\ \mathbf{z}(0)=[u_{0},0]^{\top}&\textmd{in }\ \Omega\times 0,\end{array}\right.$$

$$\Vert F({\mathbf{z}}_1(t),t)-F({\mathbf{z}}_2(t),t){\Vert }_{H\times H}$$

$$=\sqrt{\Vert v_{\lambda ,1}(t)-v_{\lambda ,2}(t){\Vert }^2+\Vert \eta (v_{\lambda ,1}(t)-v_{\lambda ,2}(t))+(A_{\lambda }(u_{\lambda ,1}(t))-A_{\lambda }(u_{\lambda ,2}(t))){\Vert }^2}$$

$$\le \sqrt{(1+2{\eta }^2)\Vert v_{\lambda ,1}(t)-v_{\lambda ,2}(t){\Vert }^2+\frac{8}{{\lambda }^2}\Vert u_{\lambda ,1}(t)-u_{\lambda ,2}(t){\Vert }^2}$$

$$\le \sqrt{1+2{\eta }^2+\frac{8}{{\lambda }^2}}{\Vert [u_{\lambda ,1}(t),v_{\lambda ,1}(t){]}^{\top }-[u_{\lambda ,2}(t),v_{\lambda ,2}(t){]}^{\top }\Vert }_{H\times H}.$$

$$u_{\lambda }\in C^1((0,\infty );H),\text{ and }{\dot{u}}_{\lambda }\in C^1((0,\infty );H),$$

$$TV(J_{\lambda }(u_0))\le TV(u_0)<\infty .$$

$${\mathcal{L}}_{\lambda }(t)=TV(J_{\lambda }(u_{\lambda }(t)))+\frac{1}{2}\Vert {\dot{u}}_{\lambda }(t){\Vert }^2,$$

$${\dot{\mathcal{L}}}_{\lambda }(t)=-\eta \Vert {\dot{u}}_{\lambda }(t){\Vert }^2$$

$${\int }_0^{\infty }\Vert {\dot{u}}_{\lambda }(t){\Vert }^{2}dt\le {\mathcal{L}}_{\lambda }(0)/\eta =TV(J_{\lambda }(u_0))/\eta <\infty ,$$

$$\Vert J_{\lambda }(u_{\lambda })(t)\Vert \le C\cdot TV(J_{\lambda }(u_{\lambda })(t)).$$

Full text

{history}\submitted

03 May 2019

\ccode

AMS Subject Classification: 35L10, 35L70, 35L72, 35L80, 49K20, 49J52, 65M12

1 Introduction

Total variation flow (TVF) and mean curvature flow for level sets of scalar functions (called level-set MCF in what follows) are important nonlinear evolutionary geometric partial differential equations (PDEs) which have been of interest in many fields during the last three decades. In the literature, they have been intensively investigated either analytically [4, 5, 6, 12, 18, 30] or from a computational viewpoint [13, 14, 21, 38, 41], to name just a few. In particular, they both find application in imaging science and geometry processing, and they are of common interest to variational and PDE methods in image processing and analysis. This is due to the fact that an image (or more general data) can be treated as a function defined on a bounded domain in $\mathbb{R}^{n}$, or more specifically a rectangular domain $\Omega\subset\mathbb{R}^{2}$. This is also the particular focus of the current paper, where we consider $u:\mathbb{R}^{2}\to\mathbb{R}$ as an image function, and $u^{d}$ as given degraded image data. In our practical context, distorted images are (i) subject to some additive noise in which case $u^{d}=u+\delta$, where $u$ denotes the true image, or (ii) corrupted by displacement errors $d:\mathbb{R}^{2}\to\mathbb{R}^{2}$ which gives $u^{d}=u(x+d(x))$.

The first case is a fundamental problem in image processing and has been continuously and intensively studied from many perspectives. Mathematical methods are also developed from several different points of view, and many of them are based on the well-known Rudin-Osher-Fatemi (ROF) model [42], where total variation (TV) is used for removing additive noise from image data. It is associated with a non-smooth energy functional, and it has the beneficial property of preserving the discontinuities (edges) of an image, which are often considered important features. Accordingly, TVF, the gradient flow of the TV functional, has been studied in this context and also beyond, see for instance [5, 6, 12, 14] and the references therein.

Problems with displacement errors have, mathematically, been the subject of several recent studies. This kind of error is not linearly separable like additive noise but rather it constitutes a nonlinear phenomenon calling for new ideas for correction. In the literature, studies are mostly focused on specific sub-classes such as, e.g. image dejittering which restricts the error $d:\mathbb{R}\to\mathbb{R}$ to occur on only one direction. In the work of [37, 23, 24], it is found out that the level-set MCF and some of its variants are capable of correcting displacement errors. An intuitive understanding is that the displacement errors interrupt the level lines of image functions, and level-set MCF is in fact a minimizing flow for the perimeter of the level lines of the functions. By setting $u_{0}=u^{d}$, the evolution of the level-set MCF produces a regularized solution which remedies the displacement errors in $u^{d}$. A proper application of the level-set MCF in this context needs, however, an appropriate stopping. Similarly, let the initial data be a noisy image, the TVF is able to decrease the total variation of the noisy image, and thus regularize the image when it is, again, properly stopped.

To summarize, we use the common framework:

[TABLE]

where $\Phi$ is a general convex functional, and $\partial$ denotes the gradient (or subgradient) operator. Throughout this paper, we will use Newton’s notation for the partial derivative with respect to time. In the context for the level-set MCF, we understand $w$ in (1) to be an immersion of a hypersurface representing a level set of a proper function, and $\Phi$ denotes the area functional of the hypersurface. For TVF, we can think of $w$ as the evolutionary image function, and $\Phi(w)$ denotes the total variation of $w$.

More recently, second-order dynamics of the form

[TABLE]

have been of great interest in the field of (convex) optimization; see, e.g. [45]. By some of the authors and other colleagues [10, 47, 15], it has also been applied as regularization methods for solving inverse problems. The damped second-order dynamics are supposed to be superior to the first-order gradient flows. The case of $\eta(t)$ being a constant is sometimes called a Heavy-Ball-with-Friction system (HBF) in the literature, see, e.g., in [11]. This system is an asymptotic approximation of the equation describing the motion of a material point with positive mass, subject to remaining on the graph of $\Phi(w)$, which moves under the action of the gravity force, and the friction force ($\eta>0$ is the associated friction parameter). The introduction of the inertial term $\ddot{w}(t)$ to the dynamic system allows to overcome some of the drawbacks of gradient descent methods, such as the well-known zig-zag phenomenon. However, in contrast to gradient descent methods, the HBF system is not necessarily a descent method for the potential energy $\Phi$. Instead, it decreases the total energy (kinetic+potential). The damping parameter $\eta$ may control the kinetic part. Larger values of $\eta$ in (2) result in more rapid evolution, while smaller values yields (2) more wave-like characteristics. The optimization properties of the HBF system have been studied intensively in  [1, 2, 11], and the references therein. Numerical algorithms based on the HBF system for solving special problems, such as, e.g. large systems of linear equations, eigenvalue problems, nonlinear Schrödinger problems, inverse source problems, etc., can be found in [25, 26, 43, 46]. There we can also see that numerical algorithms for second-order damped systems are far more efficient than algorithms for first-order systems. There are also studies on cases where $\eta(t)$ is time-variant. In particular in the recent work [45, 8, 9] many associated properties have been carefully analyzed, also a connection to Nestrov’s acceleration algorithm [40] has been revealed.

We note that the standard theory on HBF does typically not apply to PDEs, i.e., when $\partial\Phi$ gives rise to a (nonlinear) partial differential operator. In our context, however, we are confronted with quasilinear hyperbolic PDEs. In fact, in this paper, we investigate the damped second-order dynamics for both TVF and level-set MCF. The aim is to understand the new equations and their solutions from a theoretical point of view, on the one hand, and to apply them to two class of imaging problems, on the other hand. In so doing, there are several mathematical challenges to overcome. First of all, difficulties arising due to the non-linearity and non-smoothness in both the second-order TVF and the second-order level-set MCF have to be addressed. Second, both the second-order dynamics are of PDEs of quasilinear hyperbolic type, which are in general subtler than first-order ones of parabolic type as the maximum principle is out of reach in the former case. Moreover, for the level-set MCF, no convex energy has been found so far to be associated to the function $u$ introduced above. Consequently, convex analysis techniques can not be applied here. Compared to TVF, fundamental mathematical questions such as the existence of solutions and also uniqueness of the solution require more efforts, or even need to introduce a new concept for solutions of the second-order level-set MCF. Therefore, the results of this paper will not only provide novel PDE methods for image processing, but also contribute and propose interesting research questions in the fields of PDE and geometric analysis.

The contributions of the paper are twofold. (i) From a mathematical analysis point of view: We prove existence and uniqueness of the solution to the Cauchy problem for the damped second-order TVF. In doing so, we take the advantage of the TV energy functional being convex. We employ Yosida approximation to show the existence of the solution, and develop an iterative scheme, for proving the uniqueness of the solution. For the damped second-order level-set MCF, we find a connection between the equation and another novel second-order geometric PDE which evolves hypersurfaces. This provides insight into the behavior of solutions of the second-order level-set MCF if we take the hypersurfaces to be the level sets of our function. The damped second-order level-set MCF is a fully degenerate quasilinear hyperbolic PDE, for which general theory seems to be elusive at this point. As a first step towards a solution concept, we show the existence and uniqueness of the solution to a regularized version of the damped second-order level-set MCF, which is also used in our numerical development. However, the resolvability of the original equation remains an open problem. (ii) In view of applicability, it is known that the first-order level-set MCF is a minimizing flow of the total variation of the initial data. However, it exhibits different behavior than the first-order TVF. In fact, while the first-order TVF is known to decrease the contrast, the first-order level-set MCF shrinks the scale of image features. Their second-order counterparts as we studied in this work are able to preserve these features. We also note, however, that second-order equations are numerically superior when compared to their first-order counterparts.

Notation: In Table 1 we summarize some notations and abbreviations which will be frequently used in this paper. The function $u$, which is appeared often in the text, is notationally not distinguished between time-independent and time-dependent versions. However, it should be clear in the context which one is using. In most of the paper, we omit to write the spatial variable $x$ for functions which depend on both space and time.

The remainder of the paper is structured as follows: Section 2 provides the mathematical analysis of the total variation flows. Section 3 investigates the level sets mean curvature flows. Section 4 presents an algorithm and the results of numerical comparisons. A convergence analysis of the algorithm is deferred in the Appendix.

2 Total variation flows

2.1 The first-order total variation flow

We start by reviewing the first-order total variation flow (TVF) and its corresponding variational method. Total variation has become a standard tool in mathematical methods for image processing since the final decade of the last century, which is attributed to the seminal work of Rudin, Osher and Fatemi [42], who introduced the following nonsmooth variational model for recovering noisy images

[TABLE]

Here $\alpha>0$ is a regularization parameter, and $TV(\cdot)$ is known as the total variation functional. Problem (3) is usually referred to as ROF model in the literature. From a practical point of view, $TV(\cdot)$ is preferable in image processing to the standard Tikhonov regularization (quadratic smooth regularization) because it is able to keep sharp contrast (edges) in the image.

We recall the definition of $TV(\cdot)$ here. Let $\Omega\subset\mathbb{R}^{2}$. For a function $u:\Omega\rightarrow\mathbb{R}$, the total variation is defined as

[TABLE]

where $C_{0}^{\infty}(\Omega)$ presents the set of infinitely continuously differentiable functions compactly supported in $\Omega$. The space of functions of bounded variation on $\Omega$ usually denoted by $BV(\Omega)$, is given by

[TABLE]

It is well known that $BV(\Omega)$ is a Banach space, and the Sobolev space $W_{1}^{1}(\Omega)$ is embedded into $BV(\Omega)$. We are reminded that for functions $u\in W_{1}^{1}(\Omega)$ the total variation is equally characterized by the $L^{1}$ norm of the spatial gradient of $u$, that is

[TABLE]

In the following, we shall consider a Hilbert space for the function $u$, in particular we assume $u\in L^{2}(\Omega)$ for the purpose of subsequent studies. Note that for simplicity, in the following we always use the gradient notation $\nabla$ instead of $D$ for functions also in $BV(\Omega)$. Let $TV(u)$ be the total variation of the function $u\in L^{2}(\Omega)$, then

[TABLE]

It is not difficult to find that the functional $TV(\cdot)$ is convex, proper, and lower semi-continuous on the Hilbert space $L^{2}(\Omega)$.

Given the following minimization problem

[TABLE]

the first-order TVF is nothing but the negative $L^{2}$ gradient flow for minimizing (7), which reads

[TABLE]

Note here that $-\operatorname{div}\left(\frac{\nabla u}{\left|\nabla u\right|}\right)$ has been identified with an element of the subgradient of $TV(u)$, which is rather formal. It is important to give a sense to (8) as a partial differential “equation”. This was addressed in e.g. [6]. The idea there is to introduce some vector field $p(t)$ as an element in the space $X(\Omega):=\left\{p(t)\in L^{\infty}(\Omega,\mathbb{R}^{2}):\operatorname{div}(p(t))\in L^{2}(\Omega)\right\}$ for all $t\in(0,\infty)$. Then the equation (8) is understood in the sense of $\dot{u}(t)=\operatorname{div}(p(t))$, where $p$ has the form:

[TABLE]

which provides a more detailed understanding of (8). This also applies later to (11), one of our target equations in this paper. For a further mathematical analysis of the first-order TVF, we refer to, e.g., [5, 6, 12]. There, the existence and uniqueness of solutions of the Cauchy problem (8) with Neumann/Dirichelet boundary condition on $\Omega$ was established. Also, the more general case where $\Omega$ is the entire space $\mathbb{R}^{n}$ was studied. These developments are mostly motivated by applications in image denoising. Indeed, setting the initial value $u_{0}=u^{d}$, and running the flow stopped at a proper time, yields a regularized image. Usually the filtering of TVF is less destructive to the edges in images than filtering with a Gaussian, i.e., solving the heat equation with the same initial value $u_{0}$.

A formal connection between the TVF (8) and the ROF variational model (3) can be drawn as follows. Given the initial value $u_{0}$, we consider an implicit time discretization of the TVF (8) using the following iterative procedure:

[TABLE]

Identifying the time step $\Delta t$ with the regularization parameter in (3), that is $\alpha=\Delta t$, we see that (10) is in fact the Euler-Lagrange equation of (3). Therefore each iteration in (10) can be equivalently approached by solving (3), where we take $u=u_{m}$ and $u^{d}=u_{m-1}$.

2.2 The damped second-order total variation flow

Following the idea of damped second-order dynamics for gradient flows of convex functionals in Hilbert spaces, we introduce the following second-order TVF:

[TABLE]

where $\eta>0$ is the so-called damping parameter, which is assumed to be a constant, and $\partial\Omega$ denotes the boundary of the domain $\Omega\subset\mathbb{R}^{2}$ which is Lipschitz continuous, $\partial_{\nu}$ is the normal derivative and $\nu$ denotes the outward unit normal vector on $\partial\Omega$.

In order to study the resolvability of (11) we consider the following concept for its solutions.

Definition 2.1**.**

A function $u\in\mathcal{V}:=L^{\infty}((0,\infty);D(\partial TV))\cap C^{2}((0,\infty);L^{2}(\Omega))$ is called a strong solution of (11) provided

[TABLE]

given the initial and boundary conditions in (11).

Before discussing the existence of solutions for (11), we recall the resolvent operator as well as the Yosida approximation operator for the $TV$ functional. These are standard tools available in many classic textbooks (see e.g. [16, Chapter 7]):

Definition 2.2**.**

• (i)

The resolvent operator $J_{\lambda}:L^{2}(\Omega)\to D(\partial TV)$ is defined by $J_{\lambda}(w):=u$, where $u\in L^{2}(\Omega)\cap BV(\Omega)$ is the unique solution of

[TABLE]

• (ii)

The Yosida approximation operator $A_{\lambda}:L^{2}(\Omega)\to L^{2}(\Omega)$ is defined as

[TABLE]

The operators $J_{\lambda}$ and $A_{\lambda}$ have the following properties (see also [16, Chapter 7]):

Proposition 2.3**.**

1. (i)

For any fixed $\lambda>0$, $A_{\lambda}$ is a Lipschitz continuous mapping, i.e.

[TABLE] 2. (ii)

$A_{\lambda}$* is a monotone operator, i.e. $(A_{\lambda}(w_{1})-A_{\lambda}(w_{2}),w_{1}-w_{2})\geq 0$ for all $w_{1},w_{2}\in L^{2}(\Omega)$.* 3. (iii)

$A_{\lambda}(w)\in\partial TV(J_{\lambda}(w))$* for all $w\in L^{2}(\Omega)$.* 4. (iv)

For all $w\in D(\partial TV)$:

[TABLE] 5. (v)

For every $w\in L^{2}(\Omega)$:

[TABLE]

We need the following lemma.

Lemma 2.4**.**

[6*, Proposition 1.10]**

Let $u\in D(\partial TV)$, and $\beta\in\partial TV(u)$. Then,*

[TABLE]

Note that this is a special case of a general result which states that the above equality still holds true for an arbitrary convex functional homogeneous of degree $1$ besides the $TV$ functional. Having these results at hand, we proceed to proving the existence of a solution to (11).

Theorem 2.5**.**

Given $u_{0}\in D(\partial TV)$, there exists a solution of (11) in the sense of Definition 2.1.

Proof 2.6**.**

We first consider the following approximate problem with fixed $\lambda>0$:

[TABLE]

For simplicity, denote $H=L^{2}(\Omega)$, and introduce the space $H\times H$ with scalar product

[TABLE]

and the corresponding norm $\|[u,v]^{\top}\|_{H\times H}=\sqrt{\|u\|^{2}+\|v\|^{2}}$. Note that in the following proof, if there is no specification, then $\left\|\cdot\right\|$ always means the $H$ norm. Now, we define $v_{\lambda}(t)=\dot{u}_{\lambda}(t)$ and $\mathbf{z}(t)=[u_{\lambda}(t),v_{\lambda}(t)]^{\top}$, and then, rewrite (14) as a first-order dynamical system in the phase space $H\times H$, i.e.

[TABLE]

where $F(\mathbf{z}(t))=[v_{\lambda}(t),-\eta v_{\lambda}(t)-A_{\lambda}(u_{\lambda}(t))]^{\top}$.

We show first that $F$ is Lipschitz continuous for every fixed $\lambda>0$. This is true by using the Lipschitz continuity of the Yosida approximation operator $A_{\lambda}$, and we have the following inequalities

[TABLE]

The existence and uniqueness of the solution of (15) follow from the Cauchy-Lipschitz-Picard (CLP) theorem (see e.g. [16]) for first-order dynamical systems. In particular we can infer that

[TABLE]

and $\dot{u}_{\lambda}\in C^{1}((0,\infty);H)$ indicates that $u_{\lambda}\in C^{2}((0,\infty);H)$.

In the remaining part of the proof, we show that as $\lambda\to 0^{+}$ the function sequence $(u_{\lambda})_{\lambda}$ converges to a solution of our original problem (11) in the sense of Definition 2.1. We prove this by the following steps.

Step 1. We show that $\dot{u}_{\lambda}\in L^{2}((0,\infty);H)$.

According to the definition of $J_{\lambda}$ and the assumption $u_{0}\in D(\partial TV)$, we have

[TABLE]

Defining a Lyapunov function of the differential equation (14), that is

[TABLE]

it is not difficult to show that

[TABLE]

by considering (14). Integrating both sides in (18), we obtain

[TABLE]

which yields $\dot{u}_{\lambda}\in L^{2}((0,\infty),H)$ for all $\lambda\geq 0$.

Step 2. We prove that both $u_{\lambda},\dot{u}_{\lambda}\in L^{\infty}((0,\infty),H)$ are uniformly bounded.

Since $\Omega\subset\mathbb{R}^{2}$, according to Sobolev’s inequality, see e.g. [3, Theorem 3.47], a constant $C$ independent of $\lambda$ exists such that

[TABLE]

On the other hand, according to the assertion $(\text{v})$ in Proposition 2.3, a constant $\lambda_{1}$ exists such that for all $\lambda\in(0,\lambda_{1}]$,

[TABLE]

Note that it follows from (18) that $\mathcal{E}(t)$ is a non-increasing function. Thus, we obtain together with (20) and (21) that

[TABLE]

which yields $u_{\lambda}\in L^{\infty}((0,\infty),H)$ for all $\lambda\in(0,\lambda_{1}]$.

The uniform boundedness of $\dot{u}_{\lambda}$ follows from the following inequality:

[TABLE]

Step 3. We argue that both $A_{\lambda}(u_{\lambda}),\ddot{u}_{\lambda}\in L^{\infty}((0,\infty);H)$ are also uniformly bounded.

We have shown in Step 2 that $J_{\lambda}(u_{\lambda}(t))\in H$ is uniformly bounded for all $\lambda\in(0,\lambda_{1}]$, and for all $t\in(0,\infty)$. Now we show, by contradiction, that $A_{\lambda}(u_{\lambda})\in L^{\infty}((0,\infty);H)$ is uniformly bounded . Assume that there exists $\lambda^{*}\in(0,\lambda_{1}]$ such that $\left(A_{\lambda^{*}}(u_{\lambda^{*}}(t_{k}))\right)_{k\in\mathbb{N}}$ is an unbounded sequence, i.e., $\lim_{k\to\infty}\|A_{\lambda^{*}}(u_{\lambda^{*}}(t_{k}))\|=\infty$. On the other hand, $\left(J_{\lambda}(u_{\lambda}(t_{k}))\right)_{k\in\mathbb{N}}$ is uniformly bounded for all $\lambda\in(0,\lambda_{1}]$. Hence, there exists a weakly convergent subsequence, still denoted by $\left(J_{\lambda^{*}}(u_{\lambda^{*}}(t_{k}))\right)_{k\in\mathbb{N}}$ with some weak limit in $H$. Note that there must exist a subsequence such that the elements of $\left(J_{\lambda^{*}}(u_{\lambda^{*}}(t_{k}))\right)_{k\in\mathbb{N}}$ are not constant functions, otherwise we get $A_{\lambda^{*}}(u_{\lambda^{*}}(t_{k}))\equiv 0$ which is already a contradiction. We consider in particular this subsequence and use the same notation.

Now, let us define $d_{k}:=A_{{\lambda^{*}}}(u_{\lambda^{*}}(t_{k}))/\left\|A_{\lambda^{*}}(u_{\lambda^{*}}(t_{k}))\right\|$ and its smooth approximation $\tilde{d}_{k}\in C_{0}^{\infty}(\Omega)$, such that for arbitrary $\epsilon\in(0,1)$ we have:

[TABLE]

Note that such $\tilde{d}_{k}$ always exists since $A_{\lambda^{*}}(u_{\lambda^{*}}(t_{k}))\in L^{2}(\Omega)$ and $C^{\infty}_{0}(\Omega)$ is dense in $L^{2}(\Omega)$. Note that $d_{k}$ and $\tilde{d}_{k}$ are uniformly bounded in $L^{2}(\Omega)$ for all $k\in\mathbb{N}$. Moreover, $TV(\tilde{d}_{k})$ is also uniformly bounded as $\tilde{d}_{k}\in C_{0}^{\infty}(\Omega)$ and $\Omega\subset\mathbb{R}^{2}$ is compact. Since $A_{\lambda^{*}}(u_{\lambda^{*}}(t_{k}))\in\partial TV(J_{\lambda^{*}}(u_{\lambda^{*}}(t_{k})))$, we have for every $k\in\mathbb{N}$,

[TABLE]

for $\epsilon\in(0,1)$, which implies

[TABLE]

This means that $\left(\left\|A_{\lambda^{*}}(u_{\lambda^{*}}(t_{k}))\right\|\right)_{k\in\mathbb{N}}$ is a bounded sequence, which is a contradiction. Therefore, we have that $A_{\lambda}(u_{\lambda}(\cdot))\in L^{\infty}((0,\infty);H)$ is uniformly bounded for all $\lambda\in(0,\lambda_{1}]$.

The uniform boundedness of $\ddot{u}_{\lambda}\in L^{\infty}((0,\infty);H)$ for $\lambda\in(0,\lambda_{1}]$ follows from the obtained results $\dot{u}_{\lambda},A_{\lambda}(u_{\lambda})\in L^{\infty}((0,\infty);H)$ and equation (14).

Step 4. Now we are ready to show that there exists a function $u\in C^{2}((0,\infty);H)$ which is a solution to (11) in the sense of Definition 2.1.

First, we claim that for every sequence $(\lambda_{k})_{k\in\mathbb{N}}$ with $\lambda_{k}\to 0$, there exists a uniformly convergent subsequence $(u_{\lambda_{k}})_{k\in\mathbb{N}}\in C^{2}((0,T];H)$ (here we do not change the notation for the subsequence) and for every $t\in(0,T]$ of arbitrary $T\in(0,\infty)$, so that

[TABLE]

This follows from the Arzelá-Ascoli theorem by noting that

[TABLE]

are uniformly bounded for all $\lambda\in(0,\lambda_{1})$, as well as $\ddot{u}_{\lambda_{k}}\in L^{\infty}((0,\infty);H)$. Therefore, all elements of both $(u_{\lambda_{k}})$ and $(\dot{u}_{\lambda_{k}})$ are Lipschitz continuous thus equicontinuous over $t\in(0,\infty)$. Note that subsequences have to be applied here whenever they are needed.

Furthermore, the uniform boundedness of $\ddot{u}_{\lambda_{k}}$ in $L^{\infty}((0,\infty);H)$ implies that there exists a subsequence $(\lambda_{k_{j}})_{j\in\mathbb{N}}$ such that for almost every $t\in(0,T]$ and arbitrary $T<\infty$:

[TABLE]

Now we show that $u(t)\in D(\partial TV)$ for every $t\in(0,T]$, and it holds that

[TABLE]

We first notice for each $t>0$ that

[TABLE]

which means that for arbitrary but fixed $w\in H$, we have

[TABLE]

Consequently, for $0<s\leq t\leq\infty$, it holds that

[TABLE]

Using the triangle inequality and Definition 2.2 of the Yosida approximation operator, we get:

[TABLE]

It has been shown in Step 3 that $\left\|A_{\lambda}u_{\lambda}(t)\right\|$ is uniformly bounded for all $\lambda\in(0,\lambda_{1})$ and all $t\in(0,\infty)$. In combination with the uniform convergence of $u_{\lambda}\to u$, we have that $J_{\lambda}(u_{\lambda})\to u$ as $\lambda\to 0$ uniformly for all $t\in(0,T]$. Using the lower semi-continuity of the $TV$ functional, Fatou’s Lemma, and the convergence (22) and weak convergence (23), we conclude upon sending $\lambda=\lambda_{k_{j}}\to 0$ that

[TABLE]

Thus, when $t$ is a Lebesgue point of $\dot{u},\ddot{u}$ and $TV(u)$, it holds that

[TABLE]

for all $w\in H$. Since $-(\ddot{u}(t)+\eta\dot{u}(t))\in H$ is bounded, by definition of subgradient we have $u(t)\in D(\partial TV)$ and

[TABLE]

Finally we show that for all $t>0$, $u(t)\in D(\partial TV)$. For every $t>0$, let $t_{n}\to t$ and $u(t_{n})\in D(\partial TV)$, and $-(\ddot{u}(t_{n})+\eta\dot{u}(t_{n}))\in\partial TV(u(t_{n}))$. Because of the uniform boundedness of both $\dot{u}$ and $\ddot{u}$, there exists $v,\;v^{\prime}\in H$ with

[TABLE]

For every fixed $w\in H$, we have

[TABLE]

Passing to the limit $n\to\infty$, and due to the continuity of $u$ and the lower semi-continuity of $TV$, we arrive at

[TABLE]

Let $\dot{u}(t)=v$, and $\ddot{u}(t)=v^{\prime}$. Then we have shown that $u(t)\in D(\partial TV)$ and $-(\ddot{u}(t)+\eta\dot{u}(t))\in\partial TV(u(t))$ for all $t>0$. This concludes the proof.

We remark that for the first-order TVF (8), using tools from semi-group theory, the regularity of the initial data can be relaxed to $L^{2}(\Omega)$ or even $L^{1}(\Omega)$ to prove the existence and uniqueness of solutions [6]. However, this does not seem to hold true for the second-order TVF (11) as it is a nonlinear wave equation, and particularly the semi-group theory does not apply here. Also the initial value $\dot{u}(0)=0$ is not compulsory for the analytical results here and later, but it is a natural choice from an algorithmic point of view. Now we continue with the uniqueness of the solution.

Theorem 2.7**.**

The problem (11) admits a unique strong solution given the initial and boundary condition there.

Proof 2.8**.**

Let $u$ and $\bar{u}$ both be solutions of (11), that satisfy both the initial and boundary conditions. Further $p$ and $\bar{p}$ are the function forms in (9) corresponding to $u$ and $\bar{u}$, respectively. For every $s\in(0,T]$, define for every function $g\in\mathcal{V}$

[TABLE]

It is not hard to see that $\phi^{s}_{g}(s)=0$, $\dot{\phi}^{s}_{g}(t)=g(t)$. Let $v=u-\bar{u}$. Compute (11) once for $u$ and then for $\bar{u}$, subtract the two PDEs, and then test the resulting equation by $\phi^{s}_{v}(t)$ to obtain:

[TABLE]

Using integration by parts and the initial conditions $v(0)=\dot{v}(0)=0$, equation (24) becomes:

[TABLE]

Then (25) is explicitly written as

[TABLE]

with $t+h_{s}\in(t,s)$. The second equality holds thanks to the continuity of $v(t)$ and the mean value theorem.

In the following, we prove by contradiction that $v\equiv 0$ over the time domain $(0,s)$. We first notice that because of equation (11), if $u(t)\neq\bar{u}(t)$ for $t\in(0,s)$, and $u_{0}\neq 0$, then $u(t)\neq c\bar{u}(t)$ for any nonzero constant $c$. As $v(0)=0$, let $t=\epsilon>0$ be the first occasion such that $v(\epsilon)\neq 0$. If no such $\epsilon$ exists, then we are done. In case $v$ is non-zero immediately after $t=0$, then we choose a sufficiently small $\epsilon>0$ such that $v(\epsilon)\neq 0$.

Then we have

[TABLE]

On the other hand, using the boundary condition and relation (9) we have

[TABLE]

for $v(\epsilon)\neq 0$. Note that the inequality is strict as $u(t)\neq c\bar{u}(t)$ for all $t\in(0,s)$ and $c\neq 0$. Recall that $v(t)\in C^{2}([0,T];H)$. Then, by continuity of $v(t)$, there exists a neighborhood $B(\epsilon,h_{\epsilon}):=(\epsilon-h_{\epsilon},\epsilon+h_{\epsilon})$ of $\epsilon$ such that for all $t\in B(\epsilon,h_{\epsilon})$, the following relation holds true:

[TABLE]

Now we return to the right-hand side of (26) with $s=\epsilon$, and find

[TABLE]

This implies

[TABLE]

which yields a contradiction. Therefore $v(t)\equiv 0$ over $t\in[0,\epsilon]$. Now we repeatedly apply this procedure to the time domains $[n\epsilon,(n+1)\epsilon]$ for every $n\in\mathbb{N}$. This shows that $v(t)\equiv 0$ over $t\in[0,\infty)$. Thus, equation (11) admits a unique solution.

Finally, we show a decay rate for the $TV$ energy when applying the second-order TVF (11) as a total variation minimizing flow. For the formulation of the results we use the Landau symbol $o(\cdot)$, i.e., $\lim_{t\to 0}\frac{o(t)}{t}=0$ for $t>0$.

Proposition 2.9**.**

Let $u$ be the solution of the second-order TVF (11), then

[TABLE]

Proof 2.10**.**

We adopt an idea of [17]. Let us first introduce the auxiliary function

[TABLE]

By elementary calculations, we derive that

[TABLE]

Then we define the entropy functional

[TABLE]

Note that $\dot{\mathcal{E}}(t)=-\eta\|\dot{u}(t)\|^{2}$ which in combination with Lemma 2.4 implies that

[TABLE]

Integrating the above inequality over $[0,T]$ we obtain together with the non-negativity of $\mathcal{E}(t)$,

[TABLE]

On the other hand, by Theorem 2.5, $u(t)$ and $\dot{u}(t)$ are uniformly bounded. Hence, there exists a constant $M$ such that $|h(t)|\leq M$ for all $t$. Letting $T\to\infty$ in (29), we obtain

[TABLE]

Moreover, since $\mathcal{E}(t)$ is non-increasing, we deduce that

[TABLE]

Using (30), the left side of (31) tends to 0 when $t\to\infty$, which implies that

[TABLE]

Hence, we conclude $\lim_{t\to\infty}tTV(u(t))=0$, which yields that $TV(u(t))=o(\frac{1}{t})$.

This concludes our study of the second-order TVF (11). In the next section, we will study another family of nonlinear flows which are also able to decrease the total variation of a function, albeit in a somewhat different manner. It is motivated from the application in imaging science for correcting displacement errors, which is different to TVFs.

3 Mean curvature motion of level sets

Mean curvature flow of level sets of scalar functions has been analyzed first in [18, 30]. The associated equation reads

[TABLE]

In such a form, the flow can overcome the singularity or topological change which may be generated during the evolution of standard hypersurface mean curvature flow. It finds many applications in surface processing and also in image processing. A particular application which has become a research focus recently is concerned with correcting displacement errors in image data [37]. Also, a connection between the level-set MCF and a non-convex energy functional has been identified in [27].

As mentioned in the introduction, the displacement error in image data can be mathematically modeled as follows:

[TABLE]

where $u^{d}:\mathbb{R}^{2}\to\mathbb{R}$ is the measured image, $M$ is some positive real number, and $u$ is the ideal physical acquisition of the image. Assuming that the magnitude of the error bound $M$ is small, following [37, 23, 24] we may consider a first-order Taylor expansion of the function $u$ along the normal direction of the level sets of $u$:

[TABLE]

assuming $\left|\nabla u(x)\right|>0$. Then the magnitude of the displacement error can be approximated as follows:

[TABLE]

In [37], a generalized total variation regularization is employed to recover $u$ given $u^{d}$, which leads to a non-convex and non-smooth energy functional

[TABLE]

where $\alpha>0$ is a regularization parameter. The parameter $q\in(0,2]$ is introduced in order to simultaneously take care of the displacement error $d(x)$ and also the intensity error $\delta(x)$ in the image data. In the case of $q=1$ we observe that the first term of $\mathcal{E}(u;u^{d})$ in (35) is a measure reflecting both the displacement error and the intensity error by using the geometric mean of $\frac{(u(x)-u^{d}(x))^{2}}{\left|\nabla u(x)\right|^{2}}$ and $(u(x)-u^{d}(x))^{2}$.

Using formal calculations, e.g., the semi-group techniques of [27, 37, 44] or the semi-implicit iterative scheme in [24] by identifying $\alpha$ as a discrete time step (the latter is analogous to our explanation involving the ROF model and the first-order TVF in the previous section), we formally derive the following nonlinear flows:

[TABLE]

We can see that (32) emerges for $q=1$ in (36). In [37, 23, 24], it was documented that the nonlinear flows (36) are able to correct small displacement errors and also to denoise image data where the discrete time step and the stopping time play the roles of regularization parameters. This motivates us to consider a second-order damping flow based on the first-order flow, which is the equation below:

[TABLE]

We refer to this new equation as the damped second-order level-set MCF.

3.1 Heuristic observation on the damped second-order level-set MCF

Suppose that each level set of the function $u$ is a hypersurface which is well-defined in $\mathbb{R}^{2}$. In this case, it has been verified in [30] that the level-set MCF (32) is equivalent to the gradient flow of the volume (perimeter) functional for the hypersurface of every single level set, that is the standard hypersurface mean curvature flow. More precisely, let $\Gamma(t)$ be the immersion of the hypersurface into $\mathbb{R}^{2}$. Without loss of generality, we consider it to be the zero level set of $u(t)$ that is $u(\Gamma(t),t)\equiv 0$. Assume $\Gamma(t)$ to be smooth. Then the evolution of $\Gamma(t)$ governed by the first-order equation (32) is in fact characterized by the following hypersurface mean curvature flow:

[TABLE]

where $\nu(t)$ is the unit normal vector associated to the hypersurface of the level set $\Gamma(t)$. Note here that $\operatorname{div}\left(\nu(t)\right)\nu(t)=\partial\mathcal{V}(\Gamma(t))$, where $\mathcal{V}$ is the volume functional, or more precisely the length of the level lines in our case.

The mean curvature flow (38) for hypersurfaces or general manifolds has been a central topic in geometric analysis. In the level set setting, if the spatial gradient $\nabla u(\Gamma(t),t)\neq 0$, the normal field of the hypersurface of every level set can be represented by $\nu(t)=\frac{\nabla u(\Gamma(t),t)}{\left|\nabla u(\Gamma(t),t)\right|}$.

In this context, a relevant question is connected to identifying an evolutionary equation for the hypersurfaces given by the level sets of $u(t)$ associated to the second-order level-set MCF (37). In the following, we give some heuristics based on formal calculations.

Let us again take $\Gamma(t)$ to be the immersion of the zero level set of the function $u(t)$, and consider the following equation:

[TABLE]

where $\operatorname{Id}$ represents the $2\times 2$ identity matrix, and $\otimes$ denotes the tensor product of vectors. While the tensor product term may appear surprising in the context of (39) at first glance, its role will soon become clear. First, we find that

[TABLE]

where $P_{\tau}$ is the projection operator onto the tangent space of $\Gamma(t)$.

Now, we look for the connection between (39) and (37) for the evolution of the level sets of the function $u$. We first notice that $u(\Gamma(t),t)\equiv 0$ (or any other constant) which gives

[TABLE]

Differentiating with respect to time on both sides of (40), we get:

[TABLE]

Note here that we are not calculating the total time derivative of $u$ but rather the partial derivative with respect to $t$. Since $\Gamma(t)$ now follows the trajectory given by (39), we observe that

[TABLE]

leading to

[TABLE]

where we also use (41). Using the fact that $\nu(t)=\frac{\nabla u(\Gamma(t),t)}{\left|\nabla u(\Gamma(t),t)\right|}$ and $\left|\nabla u(\Gamma(t),t)\right|\neq 0$ we verify that

[TABLE]

Assuming for the moment that $u$ has sufficient regularity and interchanging the order of the time and spatial derivatives in the first term of the right hand side of (42), that is $\dot{\nabla u}(\Gamma(t),t)=(\partial_{t}\nabla u(\Gamma(t),t))$, equation (42) turns into (37) restricted to the level set $\Gamma(t)$, i.e.,

[TABLE]

This indicates that every smooth level set $\Gamma(t)$ of the solution of (37) evolves according to equation (39). Basically, (39) is a vectorial form of second-order dynamics for the mean curvature flow of hypersurfaces. However, the damping coefficient has a matrix form and involves the external function $u$. This shows that (39) is not an independent geometric PDE. Rather it needs to be coupled to (37). This is further expanded in the following remark.

Remark 3.1**.**

Consider the following second-order geometric flow for general smooth hypersurfaces of some immersion function $\Gamma(t)$:

[TABLE]

Define an entropy (Lyapunov function) for (43) through

[TABLE]

where $\mathcal{V}(\Gamma(t))$ presents the volume functional of the hypersurface, and $\nu(t)$ is the unit normal field (or also called Gauss map) of the hypersurface. Since

[TABLE]

taking into account (43) and by direct calculations, we deduce that

[TABLE]

This shows that the entropy $\mathcal{M}(t)$ is monotonically decreasing following the trajectory of the flow (43). Now assume $\Gamma(t)$ again to be the hypersurfaces of the level sets of a scalar function $u$. Isolating the term $\ddot{\Gamma}$ in (43) and inserting it into (41), and also taking into account (40) we derive a new equation corresponding to (43) as follows:

[TABLE]

Formally, for $\nabla u\neq 0$, this suggests the following equation for the scalar function $u$

[TABLE]

Equations (45) and (43) appear novel and they seem to be geometrically meaningful to study. Hyperbolic mean curvature flow for hypersurfaces has been studied in [36, 33], even though no damping term has been involved, not to mention the setting of level sets of functions. However, the new equation (45) looks rather more complicated than (37). Since our motivation here is to develop algorithms for image applications, we will skip detailed discussions on the equations (45) and (43), but rather focus on (37) in this paper.

Using the entropy (44) and following the orbit of the equation (39), we infer

[TABLE]

where $P_{\nu}$ is the normal projection operator onto the hypersurface $\Gamma(t)$. The last term makes the monotonicity of $\mathcal{M}$ unclear. It implies that large $\eta$ are preferred for monotonicity, a practical point which we pick up in Section 4 along with the algorithmic development.

3.2 On the solvability of the damped second-order level-set MCF

In order to study the solvability of the equation (37), we rewrite it to obtain the following explicit form (note that $\nabla u(t)=(u_{x_{1}}(t),u_{x_{2}}(t))^{\top}$):

[TABLE]

where $\delta_{ij}$ is the Kronecker delta function, i.e., $\delta_{ij}=1$ for $i=j$, and $\delta_{ij}=0$ for $i\neq j$. For this problem, because of its geometric meaning, it is natural to study the flow in the domain $\mathbb{R}^{2}\times(0,\infty)$ instead of $\Omega\times(0,\infty)$ as in the TVF case with bounded $\Omega$. For the latter, that is $u_{0}$ is compactly supported in $\mathbb{R}^{2}$, the results developed below will still hold by imposing Neumann boundary conditions on the boundary $\partial\Omega$ for $\partial\Omega$ sufficiently regular. As the right hand sides of (37) and (48), respectively, are not related to gradient (or subgradient) of a convex functional, the standard techniques using test functions are not applicable. Thus our previous approach to the second-order TVF is not suitable for this problem. Also, the concept of the viscosity solution, which has been developed for the first-order level-set MCF [30], is not applicable either because of the degenerate hyperbolic structure of the equation (48). Moreover, it can be checked that the nonlinear coefficients in (48), namely for $\left|p\right|>0$,

[TABLE]

satisfy $\sum_{i,j}\bar{a}_{i,j}(p)\zeta_{i}\zeta_{j}\geq 0$ for almost all $\zeta=(\zeta_{1},\zeta_{2})^{\top}\in\mathbb{R}^{2}$. Therefore (48) is a fully degenerate hyperbolic equation (some times also referred to as weakly hyperbolic) in the domain $\mathbb{R}^{2}\times(0,\infty)$.

The singularity of $\bar{a}_{i,j}$ at $p=0$ is another issue which has to be considered. For the purpose of avoiding singularities and also to eliminate the degeneracy in the equation (48), we construct a regularized version which is quasilinear but strictly hyperbolic. This is also motivated by the numerical realization of (48) from a practical algorithmic point of view.

3.3 Solution of a regularized equation

We concentrate on the following quasilinear but strictly hyperbolic equation as an approximation of (48):

[TABLE]

where $0<\epsilon\ll 1$ is fixed.

The approximation (51) can be interpreted as follows. Consider the function $v^{\epsilon}(y,t):=u^{\epsilon}(x,t)-\epsilon x_{3}$, where $y=(x,x_{3})\in\mathbb{R}^{3}$. Since $|\nabla v^{\epsilon}|^{2}=|\nabla u^{\epsilon}|^{2}+\epsilon^{2}$, the equation in (51) becomes

[TABLE]

where $v^{\epsilon}_{0}(y)=u_{0}(x)-\epsilon x_{3}$. A geometric meaning for this approximation of first-order level-set MCF has been given in [30]. There, it is depicted that $v^{\epsilon}$ is a function defined on a higher dimensional domain, whose zero-level set is a graph given by: $\Gamma^{\epsilon}(t)=\left\{y=(x,x_{3})|x_{3}=u^{\epsilon}(x,t)/\epsilon\right\}$. Then it is argued that the complicated and possibly singular evolution of the level sets $\Gamma(t)$ of $u$ is approximated by a family of well behaved smooth evolutions of the level sets $\Gamma^{\epsilon}(t)$ of a function $v^{\epsilon}$ from a higher dimensional space, in the sense that $\Gamma^{\epsilon}(t)\simeq\Gamma(t)\times\mathbb{R}$ for sufficiently small $\epsilon$ at given $t>0$. We adopt the same geometric intuition for the solution of (51) as [30] did for the first-order level-set MCF. This observation justifies the use of such a regularization properly approximating the original solution when $\epsilon$ is small.

To study the existence and uniqueness of solutions to (51), we rely on the results on linear hyperbolic equations as, for instance [35]. In particular, we consider the following equation

[TABLE]

where $w(t)$ is a given function.

To simplify notations, we write $H^{k}:=W^{k,2}(\mathbb{R}^{2})$ for $k\in\mathbb{N}$, and use $D^{k}$ to represent all derivatives with respect to temporal and spatial variables of differential orders between $[0,k]$. In the following we summarize assumptions for existence results on linear hyperbolic PDEs.

Assumption 3.2

• (i)

The functions $a_{i,j}$ are smooth and satisfy $a_{i,j}=a_{j,i}$ for $i,j=1,2$. Moreover there exists a continuous function $a(p)>0$ such that $a(p)\left|q\right|^{2}\leq\sum_{i,j}a_{i,j}(p)q_{i}q_{j}\leq\sigma_{1}\left|q\right|^{2}$ for all $p\in\mathbb{R}^{2}$ and $q\in\mathbb{R}^{2}$, for some constant $\sigma_{1}>0$.

• (ii)

$w\in C^{0}([0,\infty);H^{k})$, and the initial data are properly bounded, i.e. $\left\|D^{k}\zeta(0)\right\|_{L^{2}}\leq M_{0}$, for some $M_{0}>0$.

The first assumption implies that, for all bounded $p\in\mathbb{R}^{2}$ with $\left\|p\right\|\leq P_{0}$, there exists a fixed $\sigma>0$ such that $a(p)\left|q\right|^{2}\geq\sigma\left|q\right|^{2}$ for all $q\in\mathbb{R}^{2}$. Note that $\zeta_{0}\in H^{k}$ and $\zeta_{1}\in H^{k-1}$ have been integrated into $D^{k}\zeta(0)$ because of the second assumption. The existence and uniqueness of the solution for linear hyperbolic equations of type (55) have been established in [22, Chapter 5]. One may also refer to [35, 29] for more results. We summarize the existence, uniqueness and energy estimate here:

Theorem 3.3**.**

Given Assumption 3.2, and supposing that $f\in C^{s}([0,\infty);H^{k-s-1})$ for $k\in\mathbb{N}$ and $s\in[0,k-1]$ ( $k\geq s+1\geq 1$), the linear hyperbolic equation (55) admits a unique solution $\zeta\in C^{0}((0,\infty);H^{k})\cap C^{s}((0,\infty);H^{k-s})$ and the following estimate holds true:

[TABLE]

for all $t\in(0,\infty)$ and every integer $l\in[1,s+1]$, where the case $l=1$ applies if $f(0)\equiv 0$. Here $e(t)$ is an exponential function of $t$. In particular, if $f\equiv 0$, then there exists some positive constant $\bar{t}_{0}>0$ and $c_{0}<1$, such that for $\left\|\zeta(0)\right\|_{H^{k}}\leq c_{0}M$, it holds that

[TABLE]

Connecting to equation (51), where $a_{i,j}(p)=\delta_{ij}-\frac{p_{i}p_{j}}{\left|p\right|^{2}+\epsilon^{2}}$, we have $a(p)=\frac{\epsilon^{2}}{\left|p\right|^{2}+\epsilon^{2}}$. Moreover, it is not hard to verify that all the assumptions on $a_{i,j}$ are fulfilled for this choice. Also it can be checked that $a(p)\left|q\right|^{2}\geq\sigma\left|q\right|^{2}$ for some $\sigma>0$ as soon as $\left|p\right|$ (or in another words the norm $\left\|u^{\epsilon}(t)\right\|_{H^{k}}$) is uniformly bounded from above for $t\in(0,\bar{t}_{0}]$. With this preparation, we establish now a local (short time) existence and uniqueness of solutions of (51). To simplify the presentation, we omit the superscript $\epsilon$ in the following theorem as it is a fixed parameter anyhow.

Theorem 3.4**.**

For every fixed $\epsilon\in(0,1)$, given $u_{0}\in H^{4}$, a constant $0<c_{0}<1$, and $\left\|u_{0}\right\|_{H^{4}}\leq c_{0}M$, there exists $t_{0}>0$ such that equation (51) admits a unique solution $u\in\bigcap_{s=0}^{2}C^{s}((0,t_{0}];H^{4-s})$. Moreover, $\left\|u(t)\right\|_{H^{4}}\leq M$ holds for all $t\in(0,t_{0}]$.

Proof 3.5**.**

The main idea is borrowed from the proof of [39, Theorem 4]. We first define a differential operator of the following form:

[TABLE]

Then , we construct some initial function $u^{0}(t)\in\bigcap_{s=0}^{3}C^{s}((0,\bar{t}_{0}];H^{4-s})$ which satisfies $u^{0}(0)=u_{0}$, $\dot{u}^{0}(0)=0$, and

[TABLE]

Here $D_{x}$ denotes the spatial derivative. Next, for $m\geq 1$, we consider the following equation recursively:

[TABLE]

Using Theorem 3.3 with $f\equiv 0$, we find that $u^{m}\in\bigcap_{s=0}^{3}C^{s}((0,\bar{t}_{0}];H^{4-s})$ for all $m\geq 0$, and $\left\|u^{m}(t)\right\|_{H^{4}}\leq M$ for all $t\in(0,\bar{t}_{0}]$. As $L_{u^{m}}(u^{m+1})=0$, $L_{u^{m}}(u^{m+1}-u^{m})=-L_{u^{m}}(u^{m})$ and

[TABLE]

we arrive at the following equation:

[TABLE]

Let $f(t):=A^{m}(u^{m}(t)-u^{m-1}(t))$. Using Sobolev embedding (see, e.g., [29]) it is not hard to check that $f\in C^{1}((0,\infty);H^{1})$, as both $u^{m},\;u^{m-1}\in C(0,\infty);H^{4})\bigcap C^{1}((0,\infty);H^{3})$, and $H^{4},H^{3},H^{2}\hookrightarrow L^{\infty}(\mathbb{R}^{2})$. Using the estimate (56) from Theorem 3.3 again for the above equation (note $f(0)\equiv 0$), we have in particular the following estimate:

[TABLE]

Note that $u^{m}(t)\in H^{4}$ is uniformly bounded for $t\in(0,\bar{t}_{0}]$. We also have $u^{m}_{x_{i}}(t)\in H^{3}$ and $u^{m}_{x_{i}x_{j}}(t)\in H^{2}$, and the fact that $H^{2},H^{3}\hookrightarrow L^{\infty}(\mathbb{R}^{2})$ . Then using the following relation

[TABLE]

and applying the triangle inequality, we conclude that

[TABLE]

Here $C_{m}$ is a constant depending on semi-norms of $u_{m}$ and $u_{m-1}$, but not on their difference. Since both $\left\|u_{m-1}(t)\right\|_{H^{4}},\;\left\|u_{m}(t)\right\|_{H^{4}}$ are uniformly bounded by $M$ for all $t\in(0,\bar{t}_{0})$, there exists a constant $c$ independent of $m$ such that

[TABLE]

As $c$ is a fixed constant for all $m>1$, and $\left\|u_{m}(t)\right\|_{H^{4}}\leq M$ for all $t\in[0,\bar{t}_{0}]$, there exists a sufficiently small $t_{0}\in(0,\bar{t}_{0}]$ such that the right hand side is always strictly smaller than $1$. A recursive application of this technique then shows that

[TABLE]

Because $u_{m}\in\bigcap_{s=0}^{3}C^{s}((0,t_{0}];H^{4-s})\hookrightarrow L^{\infty}((0,t_{0}];H^{1})$ for all $m$, and the latter is a Banach space, there exists a function $u(t)\in L^{\infty}((0,t_{0}];H^{1})$ such that $\lim_{m\to\infty}\nabla u_{m}(t)=\nabla u(t)$ for $t\in(0,t_{0}]$. That is $u^{m}\to u$ strongly in $L^{\infty}((0,t_{0}];H^{1})$. When we return to (58) with $\nabla w(t)=\nabla u(t)$, we see that it actually constructs a weakly convergent sequence $(u_{m})_{m\in\mathbb{N}}\in L^{\infty}((0,t_{0}];H^{4})$ to $u\in L^{\infty}((0,t_{0}];H^{4})$ which is a solution of the nonlinear equation (51). Now we show that $u$ satisfies the regularity as stated. We use the estimate (56) for every $u^{m}$ of the equation (59), and consider $l=2,3$ there. Then we get the sequence $(\dot{u}_{m})_{m\in\mathbb{N}}$ and $(\ddot{u})_{m\in\mathbb{N}}$ are uniformly bounded over $[0,t_{0}]$ for all $m\geq 0$, respectively, and in particular they are equicontinuous. This allows us to apply the Arzelá-Ascoli theorem to show that there are subsequences of $(u_{m})_{m\in\mathbb{N}}$, still denoted by $(u_{m})_{m\in\mathbb{N}}$, so that $\dot{u}_{m}\to\dot{u}$ in $C^{1}((0,t_{0}];H^{3})$, and $\ddot{u}_{m}\to\ddot{u}$ in $C^{0}((0,t_{0}];H^{2})$, respectively, and the spatial derives $D_{x}^{2}u_{m}\to D_{x}^{2}u$ in $C^{0}((0,t_{0}];H^{2})$, as well. This yields that $u\in\bigcap_{s=0}^{2}C^{s}((0,t_{0}];H^{4-s})$ is a strong solution to (51).

The argument for uniqueness is rather similar. If there exists another solution $\bar{u}$, then let $u^{m-1}=\bar{u}$, $u^{m}=u$ in (59). Using the estimate in (60) we find then

[TABLE]

Taking into account the initial and boundary conditions, we conclude that $\left\|u(t)-\bar{u}(t)\right\|_{L^{2}}=0$ for $t\in(0,t_{0}]$. Therefore the solution is unique over $(0,t_{0}]$.

We point out that the $H^{4}$ regularity of the initial value $u_{0}$ seems necessary for using our current strategy of proof. Particularly, this regularity is required in order to have the estimate (60).

Remark 3.6**.**

Based on the short time solution, we now comment on how to achieve a global solution of (51) for arbitrary $T\in(0,\infty)$ of the domain $(0,T]$ by assuming a sufficiently regular initial value. The idea is to make sure that for arbitrary $T>0$, one has the estimate

[TABLE]

which is an assumption on the the initial data in Theorem 3.4. If this is fulfilled, we see that $u(t_{0})$ satisfies the requirement for the initial data of Theorem 3.4. Let $u(t_{0})$ again be the initial data. Then one can derive the solution for the time domain $(t_{0},2t_{0}]$ using the same technique as in Theorem 3.4 and (61). The procedure can be repeated for the whole time domain $((n-1)t_{0},nt_{0}]$ for every $n\in\mathbb{N}$ and $n\leq T/t_{0}$. This idea has been realized in [39] where more general equations have been considered. It is proven in [39, Lemma 6] that for sufficiently regular initial data, that is $\left\|u_{0}\right\|_{H^{4}}\leq\varepsilon$ for some $\varepsilon>0$ depending on $M$ in (57) sufficiently small, if (51) has the solution $u\in C((0,T],H^{4})$ for arbitrary $T>0$, and if in addition the estimation $\left\|u(t)\right\|_{H^{4}}\leq M$ holds true over the whole temporal domain $[0,T]$, then (61) holds true. Note that $\varepsilon$ does depend on $M$, but not on $T$.

We mention that for problems with higher spatial dimension, i.e. $x\in\mathbb{R}^{d}$ ($d\geq 3$) there are energy estimates of general quasilinear strictly hyperbolic equations available in [19].

In order to study the solution of the original equation (48), there are certain restrictions using the current framework. First we cannot pass $\epsilon\to 0$ as there is no uniform estimate on the approximating solutions. Second, when $\epsilon=0$, it is a degenerated hyperbolic equation in the entire domain $\mathbb{R}^{2}\times(0,\infty)$. Using the current procedure needs some energy estimates for the corresponding degenerated linear hyperbolic PDEs. However, the literature appears very sparse on such and related issues. There is some work in this direction for degenerate linear hyperbolic equations, such as, e.g., [20, 7], but definitely further efforts, or maybe even completely new concepts are required in order to successfully solve the nonlinear problem. We leave this for future work.

4 Algorithms and numerical results

In this section, we consider the numerical aspect of the two proposed damped second-order nonlinear flows. In particular, we focus on applications in image denoising and correcting displacement errors in image data that motivate this study. We first introduce a discretization of the damped second-order flows and provide a numerical algorithm. The convergence analysis of the algorithm is deferred to the appendix. Finally, we present some simulation results on the behavior of solutions. Comparisons of the numerical results obtained by second-order flows and by first-order flows are also provided.

4.1 An algorithm

Considering the evolutionary PDEs as regularization methods, the stopping time is important as it plays the role of the regularization parameter. In principle, the stopping criterion for image problems typically depends on the noise level and the initial data, just as in standard regularization theory [28], the regularization parameters are chosen according to the magnitude of noise. We provide here an automatic stopping rule based on thresholds on the high frequency in Fourier space.

Discrete images consist of pixels and are stored as matrices. For the convenience of theoretical analysis of the algorithm, we represent the matrices by column vectors in the following. However, in practice coding, direct matrix operations can be correspondingly figured out and they are preferable in terms of computational efficiency, particularly in MATLAB${}^{\text{TM}}$.

Definition 4.1**.**

Given a matrix $\mathbf{u}\in\mathbb{R}^{M}\times\mathbb{R}^{N}$, one can obtain a vector $\vec{\mathbf{u}}\in\mathbb{R}^{MN}$ by stacking the columns of $\mathbf{u}$. This defines a linear operator $vec:\mathbb{R}^{M}\times\mathbb{R}^{N}\to\mathbb{R}^{MN}$,

[TABLE]

Note that $vec(\mathbf{u})$ corresponds to a lexicographical column ordering of the components in the matrix $\mathbf{u}$. The symbol $array$ denotes the inverse of the $vec$ operator, i.e.,

[TABLE]

whenever $\mathbf{u}\in\mathbb{R}^{M}\times\mathbb{R}^{N}$ and $\vec{\mathbf{u}}\in\mathbb{R}^{MN}$.

Denote by $\mathbf{u}^{k}$ the reconstructed image at iteration $k$. Then, based on the above definition and the discretization formula (62) in the appendix, the right-hand side of our damped second-order flows, i.e.,

[TABLE]

can be rewritten in an abstract matrix form as $\mathbf{F}^{k}\vec{\mathbf{u}}^{k}$, where the matrix $\mathbf{F}^{k}$ depends on $\vec{\mathbf{u}}^{k}$. The precise form of $\mathbf{F}^{k}$ is given in (63) in the appendix, where its spectral properties and their usage in convergence considerations are discussed.

In order to set up our stopping rule, we adopt a frequency domain threshold method based on the fact that noise is usually represented by high frequencies in the frequency domain [32, Chapter 4]. An associated high frequency energy is defined by

[TABLE]

where $\mathcal{F}(\mathbf{u})$ denotes a 2D discrete Fourier transform of an image $\mathbf{u}$, and $N_{0}$ presents a selected set which contains high frequency indices. For instance, if one uses the Matlab function $\operatorname{fft2}$ as the discrete Fourier transform operator $\mathcal{F}$, then the high frequency coefficients will be the central part of the matrix, and we take $N_{0}:=[\lfloor\rho M\rfloor,M-\lfloor\rho M\rfloor]\times[\lfloor\rho N\rfloor,N-\lfloor\rho N\rfloor]$, where $\rho\in(0,0.5)$ and $\lfloor\cdot\rfloor$ denotes the floor function which is sufficient in our examples. We define a function called relative denoising efficiency (RDE) as follows:

[TABLE]

Then, the value of RDE at every iteration can be used in a stopping criterion. Based on the above preparation, we propose the following algorithm.

Algorithm 4.2

A symplectic type algorithm for discretizing damped second-order dynamics.

• Input: Image data $\mathbf{u}^{\delta}_{0}$. Parameters $\eta>0$ and $\epsilon>0$. Tolerance $tol>0$.
• Initialization: $\vec{\mathbf{u}}_{0}\leftarrow vec(\mathbf{u}_{0})$, $\vec{\mathbf{v}}_{0}\leftarrow 0$, $\Delta t_{0}>0$, $\text{RDE}(0)\leftarrow E_{N_{0}}(\mathbf{u}_{0})$, $k\leftarrow 0$.
• While: $\text{RDE}(k)>tol$

1. i.

$\vec{\mathbf{v}}^{k+1}\leftarrow(1-\eta\Delta t_{k})\vec{\mathbf{v}}^{k}+\Delta t_{k}\mathbf{F}^{k}\vec{\mathbf{u}}^{k}$; 2. ii.

$\vec{\mathbf{u}}^{k+1}\leftarrow\vec{\mathbf{u}}^{k}+\Delta t_{k}\vec{\mathbf{v}}_{k+1}$; 3. iii.

Updating $\Delta t_{k}$ according to (66) 111In practice we can also use uniform time step, then this updating step can be ignored. However the criteria in (66) (see Appendix) of chosen step size is a kind of guideline for stable convergence of the algorithm.; 4. iv.

$k\leftarrow k+1$; 5. v.

$\text{RDE}(k)\leftarrow E_{N_{0}}(\mathbf{u}^{k}))/\max_{i,j}\left|\mathcal{F}(\mathbf{u}^{k})(i,j)\right|$.

• Output: A corrected image $\hat{\mathbf{u}}\leftarrow array(\vec{\mathbf{u}}^{k})$.

The estimate in (66) gives a theoretical bound on the length of $\Delta t_{k}$ for the convergence analysis of the algorithm. If we look deeper into the iterations of the algorithm, then we find

[TABLE]

As $\vec{\mathbf{v}}^{k}=\frac{\vec{\mathbf{u}}^{k}-\vec{\mathbf{u}}^{k-1}}{\Delta t_{k-1}}$, this turns out to be

[TABLE]

Note that if we use a uniform time step, that is $\Delta t_{k}\equiv\Delta t_{0}$ , and choose $\eta=\frac{1}{\Delta t_{0}}$, then the algorithm is equivalent to the steepest descent method with step size $(\Delta t_{0})^{2}$ whenever $\mathbf{F}^{k}\vec{\mathbf{u}}^{k}$ can be interpreted as the negative gradient direction for an associated energy. On the other hand, we can see that the term $\frac{\Delta t_{k}}{\Delta t_{k-1}}(1-\eta\Delta t_{k})(\vec{\mathbf{u}}^{k}-\vec{\mathbf{u}}^{k-1})$ plays a similar role as the correction step in Nesterov’s scheme [40], by which it is supposed to accelerate the steepest descent method. In the following examples, in order to draw comparisons, we will always implement the first-order method using Algorithm 4.2 by setting $\eta\equiv\frac{1}{\Delta t}$. In this sense, we will find that Algorithm 4.2 by the second-order flows also has an acceleration and be favoured over the first-order flows for the applications. One may note that the time step of the first-order method and the second-order method are not equal, as the former is $(\Delta t)^{2}$, but the latter is $\Delta t$. However, we should be aware that for the discretization of evolutionary PDEs, the Courant-–Friedrichs–-Levy (CFL) condition needs to be taken into account for explicit time discretizations. This puts restrictions on the length of the discrete time steps to the numerical implementations, e.g. for $2$D linear equations, second-order flows can have $\Delta t\sim h_{x}$ while first-order flows usually have $\Delta t\sim h_{x}^{2}$. Here $h_{x}$ is the discretization mesh size of spatial variables. In this sense, we argue that discrete time steps of order $(\Delta t)^{2}$ for first-order flows and $\Delta t$ for second-order ones are justified.

4.2 Numerical results

Evolution of a characteristic function

In order to study some fundamental effects of our PDEs on contour and scale of images, in the first example we test with an image resulting from a characteristic function. We start with the image which is a scaled indicator function of a square $Q\subset\Omega$ (left image in Figure 1)

[TABLE]

It is well-known that the first-order TVF will decrease the intensity value of the region $Q$, but intends to preserve its shape, while the first-order level-set MCF will shrink the square $Q$ slowly to a circular shape, thus reducing the perimeter of boundary of the square $Q$, but it will preserve the intensity value. In Figure 1, we present the simulations on the evolution of the second-order flows (11) and (37). We use a square of size $205\times 205$, and fix the domain $\Omega=(0,1)\times(0,1)$; therefore the spatial step size is $\Delta x=\Delta y=1/205$. We use uniform step size for time discretization, and choose $\Delta t_{1}=0.001$, and $\eta=\frac{1}{50\Delta t_{1}}$ for the TVF methods, and choose $\Delta t_{2}=0.0001$, and $\eta=\frac{1}{20\Delta t_{2}}$ for the level-set MCF methods. For both we run $50000$ iterations. From Figure 1, we find that the damped second-order dynamics present exactly the same behavior as their first-order counterparts. In the three images in Figure 1, we take the same pixel at the position $(106,100)$: The intensity value is $255$ in the original square which is the initial value, it is decreased to $242.2$ in the image evolved with respect to the damped second-order TVF, but remains the same in the image evolved according to the damped second-order level-set MCF. On the other hand, the shape of the square is almost not changed by the damped second-order TVF except the sharp corners, but has been shrunk to a circle by the damped second-order level-set MCF.

Note here and also in the other examples that we take $\epsilon=10^{-16}$ which is already much smaller then the temporal and spatial mesh sizes, respectively. However, it seems sufficient for the numerical examples we considered. Numerical diffusion is observed in the second and third images in Figure 1, and the effect grows as time step and iteration numbers get larger. However, in our following image applications, only a small number of iterations is needed, and thus we will not investigate this issue here further since it is out of the scope of our current paper. Rather we refer to existing methods in the setting of total variation minimization and beyond, such as, e.g., [14, 34].

Image de-noising

Our second set of examples are on image de-noising. Due to the previous discussion on the characteristic image, we will test with both the TVFs and the level-set MCFs. For the reason of comparing the results, we will always consider uniform time steps $\Delta t_{k}=\Delta t$ for all the four kinds of discretized flows. The image we tested here are of pixel size $400\times 400$, and we fix it to be in the domain $\Omega=(0,1)\times(0,1)$; therefore the spatial step size is $\Delta x=\Delta y=1/400$. We choose $\Delta t_{1}=0.003$, and $\eta=\frac{1}{50\Delta t_{1}}$ for the TVF methods, and choose $\Delta t_{2}=0.0001$, and $\eta=\frac{1}{10\Delta t_{2}}$ for the level-set MCF methods. Using the $\eta\equiv\frac{1}{\Delta t}$ in Algorithm 4.2 for the first-order flows, the first-order TVF and first-order MCF has a step size $(\Delta t_{i})^{2}$ for $i=1,2$. For setting the stopping criteria of the algorithm we choose $\rho=0.2$, and $tol=1$ in this example. We simulate the noise by a Gaussian distribution of mean [math] and standard deviation $20$, and add it to the pepper image, which results in an image with approximately $15$ percentage of noise.

The results in Figure 2 indicate that all four methods yield competitive results for denoising, while there is a big difference between the first-order flows and the damped second-order flows in terms of CPU time. Note that this has been reflected in the different iteration times to reach the stopping rule. For the first-order TVF and the second-order one, it is 17969 vs. 349 iterations, while for the first-order level-set MCF and the second-order one, it is 4018 vs. 348 iterations. Each scheme is run on a computer with Intel 3687U CPU, 2.10GHz$\times$4, 15.5GB RAM, and using Matlab 2017b. One may notice that the algorithm by level-set MCFs is capable of denoising almost as good as the TVFs.

Displacement error correction

Now, we show the results of correcting displacement errors in Figure 3. We again consider the pepper image of pixel size $400\times 400$ to be in the domain $\Omega=(0,1)\times(0,1)$. We degrade the image by some displacement error yielding a jittered pepper image. We choose the parameters to be $\Delta t_{1}=0.003$, and $\eta=\frac{1}{50\Delta t_{1}}$ for the damped second-order TVF, and the parameters $\Delta t_{2}=0.0001$, and $\eta=\frac{1}{30\Delta t_{2}}$ for the damped second-order level-set MCF. For setting the stopping rule of the algorithm, we set $\rho=0.2,\;tol=0.25$. We notice that the iteration numbers are again different between the algorithms due to the first-order flows and the second-order flows, in order to reach the stopping criteria. The number of iterations of first-order TVF vs. second-order TVF is $63556:1211$, and the number of iterations of first-order level-set MCF vs. second-order level-set MCF is $6625:268$.

We are aware that the results of the TVFs and level-set MCFs turn out to be quite different in this example. This is not surprising, as we have observed from Figure 1 that the second-order TVF has very limited effect in changing the curvature of the level lines in comparison with second-order level-set MCF. The algorithm with TVFs takes larger efforts to correct the jitter error at the price of sacrificing the contrast of the images as we can see clearly from Figure 3.

Simultaneous denosing and correcting displacement errors

In Figure 4, we show the results on dejittering and denoising simultaneously using the algorithm by the damped second-order level-set MCF and comparing it with the results from the damped second-order TVF. We use the jittered pepper image from the last example and then add the same amount of Gaussian noise as we did in the denoising example (appr. $15$ percentage of noise), and also we select the same disretization parameters as before. We do not show the results given by the first-order flows as they are similar to the second-order ones but requiring larger iteration times. From Figure 4, we see that, upto the stopping threshold with ($\rho=0.2,\;tol=0.5$), the second-order level-set MCF performs better in correcting the displacement error than the second-order TVF, while the latter does better for denoising than the former. Overall, the algorithm by second-order level-set MCF outperforms in this example as both the noise and the jitter are significantly reduced simultaneously. The observed iteration times are $581$ for the second-order TVF and $241$ for second-order level-set MCF, respectively.

Remark 4.3**.**

Notice that a smaller value of $\eta$ in Algorithm 4.2 usually results in better efficiency (less number of iterations to achieve the same outcome with respect to the same discrete time step size) in both denoising and dejittering tasks for both the second-order flows. However, there is a trade off as too small $\eta$ causes unstable evolution. This is particularly relevant for the second-order level-set MCF algorithm as explained in Section 3.1. There we have argued that $\eta$ needs to be sufficiently large to provide an energy decay with respect to the level sets evolution.

5 Concluding remarks

This paper has studied two geometric quasilinear hyperbolic partial differential equations, namely the second-order total variation flow and the second-order level set mean curvature flow. For the former equation, we have a relatively complete result on its well-posedness, which is attributed to the convexity of the total variational functional. However, for the latter, we have only obtained a very preliminary result on its well-posedness by considering a regularized version. The main difficulty there comes from the degeneracy of the hyperbolic structure of the equation. Particularly, different to the former problem, there has no associated convex functional been found out for the latter one. Instead, we identified some novel geometric PDEs evolving hypersurfaces to understand the behavior of the solution of the second-order level-set MCF. From an application point of view, we have observed that both second-order flows are able to generate efficient numerical algorithms for the motivating tasks in imaging sciences. The two types of flows have different behaviors, and this has been verified by numerical examples. As a consequence, they have different strengths in our imaging applications. The TVFs are able to remove additive noise efficiently, but cannot properly treat displacement errors, while the MCFs seem to simultaneously deal with these two tasks, at least to some extent. Based on the above observation, several interesting theoretical problems have been identified. For instance, the well-posedness of the original second-order level-set MCF (37), the asymptotic analysis on the solutions of both the second-order flows and their generalizations to higher dimensional spaces in $\mathbb{R}^{N}$ $(N\geq 3)$ are worthwhile to be further pursued. On the other hand, it would also be interesting to conduct a systematic investigation of the newly derived geometric PDE and its corresponding level sets equation pointed out in equations (43) and (45), respectively.

Acknowledgement

GD thanks Mårten Gulliksson for his hospitality when he was visiting Örebro University, and he also thanks Otmar Scherzer for some stimulating discussions. GD and MH acknowledge support of Institut Henri Poincaré (IHP) (UMS 839 CNRS-Sorbonne Université), and LabEx CARMIN (ANR–10–LABX–59–01), while they both got funded for visiting ‘the Mathematics of Imaging’ research trimester held at IHP in 2019, where the manuscript was finalized.

Funding: The work of GD and MH has been supported by a MATHEON Research Center project CH12 funded by the Einstein Center for Mathematics (ECMath) Berlin, and also funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy –-The Berlin Mathematics Research Center MATH+ (EXC–2046/1, project ID: 390685689). The work of YZ is supported by the Alexander von Humboldt foundation through a postdoctoral researcher fellowship.

Appendix A Convergence analysis of Algorithm 4.2

For the sake of simplicity and clarity of statements, let us consider a uniform grid $\Omega_{MN}=\{(x_{i},y_{j})\}^{M,N}_{i,j=1}$, discretizing $\Omega$ with the uniform step size $h=x_{i+1}-x_{i}=y_{j+1}-y_{j}$. Define $\mathbf{u}(t)=[u(x_{i},y_{j},t)]^{M,N}_{i,j=1}$, and denote by $\mathbf{u}^{k}$ the projection of $u(x,y,t)$ onto the spatial grid $\Omega_{MN}$ and time point $t=t_{k}$. Denote by $c^{\epsilon,k}_{i,j}(\mathbf{u}^{k})=1/\left(\epsilon+\frac{1}{h}\sqrt{(\mathbf{u}^{k}_{i+1,j}-\mathbf{u}^{k}_{i,j})^{2}+(\mathbf{u}^{k}_{i,j+1}-\mathbf{u}^{k}_{i,j})^{2}}\right)$. Using forward-backward differences, we obtain

[TABLE]

To put TVFs and MCFs under a common umbrella, we use $b(\mathbf{u}^{k})\textmd{div}\left(c^{\epsilon,k}(\mathbf{u}^{k})\nabla\mathbf{u}^{k}\right)$ to represent the nonlinear part of the equations: $b(\mathbf{u}^{k})\equiv 1$ for TVFs, and

[TABLE]

for MCFs. By applying lexicographical column ordering of $\mathbf{u}^{k}_{i,j}$ and assuming the Neumann boundary condition, we obtain the matrix representation of $b(\mathbf{u}^{k})\textmd{div}\left(c(\mathbf{u}^{\epsilon,k})\nabla\mathbf{u}^{k}\right)$, denoted as $\mathbf{F}^{k}\vec{\mathbf{u}}^{k}$, where

[TABLE]

and $\mathbf{G}^{k}$ is the $MN\times MN$ matrix with $N\times N$ block entries given by

[TABLE]

Here $I^{k}_{j}$ is the $M\times M$ diagonal matrix $I^{k}_{j}=\operatorname{diag}\left(c^{\epsilon,k}_{1,j},\cdots,c^{\epsilon,k}_{M,j}\right)$, $\mathbf{0}$ represents the $M\times M$ zero matrix, and $L^{k}_{j}$ is the $M\times M$ matrix of the form

[TABLE]

where $\tilde{c}^{\epsilon,k}_{i,j}:=2c^{\epsilon,k}_{i,j}+c^{\epsilon,k}_{i-1,j}+c^{\epsilon,k}_{i,j-1}$.

Proposition A.1**.**

All eigenvalues of $\mathbf{G}^{k}$ (for all $k\in\mathbb{N}$) are non-positive.

Proof A.2**.**

By the definition of $\mathbf{G}^{k}$, i.e. (64), it is not difficult to show that $\mathbf{G}^{k}$ is a symmetric and diagonally dominant matrix. Then, all eigenvalues of $\mathbf{G}^{k}$ (for all $k\in\mathbb{N}$) are real and, by Gershgorin’s circle theorem [31], for each eigenvalue $\lambda$ there exists an index $\nu$ such that

[TABLE]

which implies, by definition of the diagonal dominance, $\lambda\leq 0$. Here, $[\mathbf{G}^{k}]_{\nu,\imath}$ denotes the element of the matrix $\mathbf{G}^{k}$ at the position $(\nu,\imath)$.

Denote $\vec{\mathbf{v}}^{k}=\frac{d\vec{\mathbf{u}}^{k}}{dt}$, and recall Algorithm 4.2 where the symplectic Euler scheme is applied to discretize the second-order flow (11) or (37), i.e.,

[TABLE]

where $\vec{\mathbf{u}}^{d}=vec(\mathbf{u}^{d})$ and $\mathbf{u}^{d}$ is the project of $u^{d}(x)$ onto the grid $\Omega_{MN}$.

Now, we are in a position to give a numerical analysis for the scheme (65) in Algorithm 4.2.

Theorem A.3**.**

Let $\eta>0$ be a fixed damping parameter. If the step size is chosen to fulfill

[TABLE]

where $\lambda^{k}_{max}$ is the maximal eigenvalue of $-\mathbf{F}^{k}$, and $b^{k}_{max}:=\max^{M,N}_{i,j=1}b^{k}_{i,j}$, then the scheme (65) is convergent.

Proof A.4**.**

Denote $\mathbf{z}_{k}=(\mathbf{v}^{k};\mathbf{u}^{k})$, and rewrite equation (65) by

[TABLE]

where

[TABLE]

It is well-known that a sufficient condition for the convergence of the iteration scheme (67) is that $A_{k}$ is a contractive operator, i.e., $\|A_{k}\|_{2}<1$.

By the elementary calculation and the decomposition $-\mathbf{G}^{k}=Q\Lambda^{k}Q^{\top}$ with $\Lambda^{k}=\operatorname{diag}(\lambda^{k}_{i}),\lambda^{k}_{i}\geq 0,i=1,...,MN$, we derive that the eigenvalues of $A_{k}$ are

[TABLE]

where $b^{k}_{i}\geq 0$ represents the $i$-th element in the diagonal of matrix $\mathbf{B}^{k}$. Hence, in order to prove $\|A_{k}\|_{2}\leq 1$, it is sufficient to show that for all $i=1,...,MN$: $|\mu^{k}_{i,\pm}|\leq 1$ for the time step size $\Delta t_{k}$, defined in (66).

For each fixed $i$, there are three possible cases: the overdamped case $\left(\Delta t_{k}b^{k}_{i}\lambda^{k}_{i}+\eta\right)^{2}>4b^{k}_{i}\lambda^{k}_{i}$, the underdamped case $\left(\Delta t_{k}b^{k}_{i}\lambda^{k}_{i}+\eta\right)^{2}<4b^{k}_{i}\lambda^{k}_{i}$, and the critical damping case $\left(\Delta t_{k}b^{k}_{i}\lambda^{k}_{i}+\eta\right)^{2}=4b^{k}_{i}\lambda^{k}_{i}$. We consider each case separately.

For the overdamped case, define $a:=\frac{\eta+\Delta t_{k}b^{k}_{i}\lambda^{k}_{i}}{2\sqrt{b^{k}_{i}\lambda^{k}_{i}}}$ ($a>1$). Then,

[TABLE]

Obviously, $\mu^{k}_{i,\pm}<1$ by noting the positivity of the second term on the right-hand side of the equation above. Now, let us show the inequality $\mu^{k}_{i,\pm}>-1$. By the choice of the time step size $\Delta t_{k}$ in (66), we know that $\eta\Delta t_{k}<1$ and $b^{k}_{max}\lambda^{k}_{max}\Delta t^{2}_{k}<1$, which implies that $\eta\Delta t_{k}+b^{k}_{max}\lambda^{k}_{max}\Delta t^{2}_{k}\leq 2$. Therefore, we have

[TABLE]

Since $a=\frac{\eta}{2\sqrt{b^{k}_{i}\lambda^{k}_{i}}}+\frac{\Delta t_{k}}{2}\sqrt{b^{k}_{i}\lambda^{k}_{i}}$, using inequality (69), we deduce that

[TABLE]

which implies that

[TABLE]

Therefore, we conclude that $|\mu^{k}_{i,\pm}|\leq 1$ for the overdamped case.

Now, consider the underdamped case. In this case, since $\left(\Delta t_{k}b^{k}_{i}\lambda^{k}_{i}+\eta\right)^{2}<4b^{k}_{i}\lambda^{k}_{i}$, we have

[TABLE]

which implies $|\mu^{k}_{i,\pm}|<1$ for any fixed pair $(\eta,\Delta t_{k})$ satisfying (66).

Finally, consider the critical damping case. In this case, the eigenvalue for $\mu^{k}_{i,\pm}$ is simply given by

[TABLE]

which yields the desired result according to the argument in the critical damping case. In conclusion, all the eigenvalues of matrix $A_{k}$ are smaller or equal than $1$, therefore it is a contractive operator. Then the scheme (65) is convergent.