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

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

Guozhi Dong        Michael Hintermüller    Ye Zhang

In this paper, we study damped second-order dynamics, which are quasilinear hyperbolic partial differential equations (PDEs). This is inspired by the recent development of second-order damping systems for accelerating energy decay of gradient flows. We concentrate on two equations: one is a damped second-order total variation flow, which is primarily motivated by the application of image denoising; the other is a damped second-order mean curvature flow for level sets of scalar functions, which is related 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. For the latter, we draw a connection between the equation and some second-order geometric PDEs evolving the hypersurfaces which are described by level sets of scalar functions, and show the existence and uniqueness of the solution for a regularized version of the equation. The latter is used in our algorithmic development. A general algorithm for numerical discretization of the two nonlinear PDEs is proposed and analyzed. Its efficiency is demonstrated by various numerical examples, where simulations on the behavior of solutions of the new equations and comparisons with first-order flows are also documented.

Submitted: 03 May 2019

Keywords: Quasilinear hyperbolic equation, geometric PDEs; total variation flow; mean curvature flow; level set; second-order dynamics; non-smooth and non-convex variational methods; image denoising; displacement error correction.

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 or from a computational viewpoint , 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 , or more specifically a rectangular domain . This is also the particular focus of the current paper, where we consider as an image function, and as given degraded image data. In our practical context, distorted images are (i) subject to some additive noise in which case , where denotes the true image, or (ii) corrupted by displacement errors which gives .

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 , 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 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 to occur on only one direction. In the work of , 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 , the evolution of the level-set MCF produces a regularized solution which remedies the displacement errors in . 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:


where is a general convex functional, and 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 in (1) to be an immersion of a hypersurface representing a level set of a proper function, and denotes the area functional of the hypersurface. For TVF, we can think of as the evolutionary image function, and denotes the total variation of .

More recently, second-order dynamics of the form


have been of great interest in the field of (convex) optimization; see, e.g. . By some of the authors and other colleagues , 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 being a constant is sometimes called a Heavy-Ball-with-Friction system (HBF) in the literature, see, e.g., in . 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 , which moves under the action of the gravity force, and the friction force ( is the associated friction parameter). The introduction of the inertial term 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 . Instead, it decreases the total energy (kinetic+potential). The damping parameter may control the kinetic part. Larger values of in (1) result in more rapid evolution, while smaller values yields (1) more wave-like characteristics. The optimization properties of the HBF system have been studied intensively in  , 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 . 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 is time-variant. In particular in the recent work many associated properties have been carefully analyzed, also a connection to Nestrov’s acceleration algorithm has been revealed.

We note that the standard theory on HBF does typically not apply to PDEs, i.e., when 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 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 Description
with only spatial variable , is time-independent
with only temporal variable , maps to either real values,
or elements in some function spaces that is
is both space- and time-dependent
& first-order and second-order time derivative of , respectively
spatial gradient (including distributional sense) of function
spatial divergence of the vector field
inner product of two elements in Hilbert spaces (mostly space)
inner product of two elements in
Table 1: Notations and abbreviations.

Notation: In Table 1 we summarize some notations and abbreviations which will be frequently used in this paper. The function , 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 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 , who introduced the following nonsmooth variational model for recovering noisy images


Here is a regularization parameter, and is known as the total variation functional. Problem (2.1) is usually referred to as ROF model in the literature. From a practical point of view, 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 here. Let . For a function , the total variation is defined as


where presents the set of infinitely continuously differentiable functions compactly supported in . The space of functions of bounded variation on usually denoted by , is given by


It is well known that is a Banach space, and the Sobolev space is embedded into . We are reminded that for functions the total variation is equally characterized by the norm of the spatial gradient of , that is


In the following, we shall consider a Hilbert space for the function , in particular we assume for the purpose of subsequent studies. Note that for simplicity, in the following we always use the gradient notation instead of for functions also in . Let be the total variation of the function , then

It is not difficult to find that the functional is convex, proper, and lower semi-continuous on the Hilbert space .

Given the following minimization problem


the first-order TVF is nothing but the negative gradient flow for minimizing (2.5), which reads


Note here that has been identified with an element of the subgradient of , which is rather formal. It is important to give a sense to (2.6) as a partial differential “equation”. This was addressed in e.g. . The idea there is to introduce some vector field as an element in the space for all . Then the equation (2.6) is understood in the sense of , where has the form:


which provides a more detailed understanding of (2.6). This also applies later to (2.9), one of our target equations in this paper. For a further mathematical analysis of the first-order TVF, we refer to, e.g., . There, the existence and uniqueness of solutions of the Cauchy problem (2.6) with Neumann/Dirichelet boundary condition on was established. Also, the more general case where is the entire space was studied. These developments are mostly motivated by applications in image denoising. Indeed, setting the initial value , 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 .

A formal connection between the TVF (2.6) and the ROF variational model (2.1) can be drawn as follows. Given the initial value , we consider an implicit time discretization of the TVF (2.6) using the following iterative procedure:


Identifying the time step with the regularization parameter in (2.1), that is , we see that (2.8) is in fact the Euler-Lagrange equation of (2.1). Therefore each iteration in (2.8) can be equivalently approached by solving (2.1), where we take and .

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:


where is the so-called damping parameter, which is assumed to be a constant, and denotes the boundary of the domain which is Lipschitz continuous, is the normal derivative and denotes the outward unit normal vector on .

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

Definition 2.1

A function is called a strong solution of (2.9) provided


given the initial and boundary conditions in (2.9).

Before discussing the existence of solutions for (2.9), we recall the resolvent operator as well as the Yosida approximation operator for the functional. These are standard tools available in many classic textbooks (see e.g. ):

Definition 2.2
  • The resolvent operator is defined by , where is the unique solution of

  • The Yosida approximation operator is defined as

The operators and have the following properties (see also ):

Proposition 2.1
  1. For any fixed , is a Lipschitz continuous mapping, i.e.

  2. is a monotone operator, i.e. for all .

  3. for all .

  4. For all :

  5. For every :

We need the following lemma.

Lemma 2.1

Let , and . Then,

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 besides the functional. Having these results at hand, we proceed to proving the existence of a solution to (2.9).

Theorem 2.1

Given , there exists a solution of (2.9) in the sense of Definition 2.1.

Proof.  We first consider the following approximate problem with fixed :


For simplicity, denote , and introduce the space with scalar product

and the corresponding norm . Note that in the following proof, if there is no specification, then always means the norm. Now, we define and , and then, rewrite (2.12) as a first-order dynamical system in the phase space , i.e.


where .

We show first that is Lipschitz continuous for every fixed . This is true by using the Lipschitz continuity of the Yosida approximation operator , and we have the following inequalities

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

and indicates that .

In the remaining part of the proof, we show that as the function sequence converges to a solution of our original problem (2.9) in the sense of Definition 2.1. We prove this by the following steps.

Step 1. We show that .

According to the definition of and the assumption , we have


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


it is not difficult to show that


by considering (2.12). Integrating both sides in (2.16), we obtain


which yields for all .

Step 2. We prove that both are uniformly bounded.

Since , according to Sobolev’s inequality, see e.g. , a constant independent of exists such that


On the other hand, according to the assertion in Proposition 2.1, a constant exists such that for all ,


Note that it follows from (2.16) that is a non-increasing function. Thus, we obtain together with (2.18) and (2.19) that

which yields for all .

The uniform boundedness of follows from the following inequality:

Step 3. We argue that both are also uniformly bounded.

We have shown in Step 2 that is uniformly bounded for all , and for all . Now we show, by contradiction, that is uniformly bounded . Assume that there exists such that is an unbounded sequence, i.e., . On the other hand, is uniformly bounded for all . Hence, there exists a weakly convergent subsequence, still denoted by with some weak limit in . Note that there must exist a subsequence such that the elements of are not constant functions, otherwise we get which is already a contradiction. We consider in particular this subsequence and use the same notation.

Now, let us define and its smooth approximation , such that for arbitrary we have:

Note that such always exists since and is dense in . Note that and are uniformly bounded in for all . Moreover, is also uniformly bounded as and is compact. Since , we have for every ,

for , which implies

This means that is a bounded sequence, which is a contradiction. Therefore, we have that is uniformly bounded for all .

The uniform boundedness of for follows from the obtained results and equation (2.12).

Step 4. Now we are ready to show that there exists a function which is a solution to (2.9) in the sense of Definition 2.1.

First, we claim that for every sequence with , there exists a uniformly convergent subsequence (here we do not change the notation for the subsequence) and for every of arbitrary , so that


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

are uniformly bounded for all , as well as . Therefore, all elements of both and are Lipschitz continuous thus equicontinuous over . Note that subsequences have to be applied here whenever they are needed.

Furthermore, the uniform boundedness of in implies that there exists a subsequence such that for almost every and arbitrary :


Now we show that for every , and it holds that

We first notice for each that

which means that for arbitrary but fixed , we have

Consequently, for , it holds that

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

It has been shown in Step 3 that is uniformly bounded for all and all . In combination with the uniform convergence of , we have that as uniformly for all . Using the lower semi-continuity of the functional, Fatou’s Lemma, and the convergence (2.20) and weak convergence (2.21), we conclude upon sending that

Thus, when is a Lebesgue point of and , it holds that

for all . Since is bounded, by definition of subgradient we have and

Finally we show that for all , . For every , let and , and . Because of the uniform boundedness of both and , there exists with

For every fixed , we have

Passing to the limit , and due to the continuity of and the lower semi-continuity of , we arrive at

Let , and . Then we have shown that and for all . This concludes the proof.


We remark that for the first-order TVF (2.6), using tools from semi-group theory, the regularity of the initial data can be relaxed to or even to prove the existence and uniqueness of solutions . However, this does not seem to hold true for the second-order TVF (2.9) as it is a nonlinear wave equation, and particularly the semi-group theory does not apply here. Also the initial value 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.2

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

Proof.  Let and both be solutions of (2.9), that satisfy both the initial and boundary conditions. Further and are the function forms in (2.7) corresponding to and , respectively. For every , define for every function

It is not hard to see that , . Let . Compute (2.9) once for and then for , subtract the two PDEs, and then test the resulting equation by to obtain:


Using integration by parts and the initial conditions , equation (2.22) becomes:


Then (2.23) is explicitly written as


with . The second equality holds thanks to the continuity of and the mean value theorem.

In the following, we prove by contradiction that over the time domain . We first notice that because of equation (2.9), if for , and , then for any nonzero constant . As , let be the first occasion such that . If no such exists, then we are done. In case is non-zero immediately after , then we choose a sufficiently small such that .

Then we have

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

for . Note that the inequality is strict as for all and . Recall that . Then, by continuity of , there exists a neighborhood of such that for all , the following relation holds true:


Now we return to the right-hand side of (2.24) with , and find

This implies

which yields a contradiction. Therefore over . Now we repeatedly apply this procedure to the time domains for every . This shows that over . Thus, equation (2.9) admits a unique solution.


Finally, we show a decay rate for the energy when applying the second-order TVF (2.9) as a total variation minimizing flow. For the formulation of the results we use the Landau symbol , i.e., for .

Proposition 2.2

Let be the solution of the second-order TVF (2.9), then

Proof.  We adopt an idea of . Let us first introduce the auxiliary function


By elementary calculations, we derive that