A class of secondorder geometric quasilinear hyperbolic PDEs and their application in imaging science
Abstract
In this paper, we study damped secondorder dynamics, which are quasilinear hyperbolic partial differential equations (PDEs). This is inspired by the recent development of secondorder damping systems for accelerating energy decay of gradient flows. We concentrate on two equations: one is a damped secondorder total variation flow, which is primarily motivated by the application of image denoising; the other is a damped secondorder mean curvature flow for level sets of scalar functions, which is related to a nonconvex 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 secondorder 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 firstorder flows are also documented.
Submitted: 03 May 2019
Keywords: Quasilinear hyperbolic equation, geometric PDEs; total variation flow; mean curvature flow; level set; secondorder dynamics; nonsmooth and nonconvex 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 levelset 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 wellknown RudinOsherFatemi (ROF) model , where total variation (TV) is used for removing additive noise from image data. It is associated with a nonsmooth 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 subclasses 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 levelset 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 levelset MCF is in fact a minimizing flow for the perimeter of the level lines of the functions. By setting , the evolution of the levelset MCF produces a regularized solution which remedies the displacement errors in . A proper application of the levelset 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:
(1.1) 
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 levelset 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, secondorder dynamics of the form
(1.2) 
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 secondorder dynamics are supposed to be superior to the firstorder gradient flows. The case of being a constant is sometimes called a HeavyBallwithFriction 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 wellknown zigzag 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 wavelike 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 secondorder damped systems are far more efficient than algorithms for firstorder systems. There are also studies on cases where is timevariant. 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 secondorder dynamics for both TVF and levelset 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 nonlinearity and nonsmoothness in both the secondorder TVF and the secondorder levelset MCF have to be addressed. Second, both the secondorder dynamics are of PDEs of quasilinear hyperbolic type, which are in general subtler than firstorder ones of parabolic type as the maximum principle is out of reach in the former case. Moreover, for the levelset 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 secondorder levelset 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 secondorder 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 secondorder levelset MCF, we find a connection between the equation and another novel secondorder geometric PDE which evolves hypersurfaces. This provides insight into the behavior of solutions of the secondorder levelset MCF if we take the hypersurfaces to be the level sets of our function. The damped secondorder levelset 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 secondorder levelset 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 firstorder levelset MCF is a minimizing flow of the total variation of the initial data. However, it exhibits different behavior than the firstorder TVF. In fact, while the firstorder TVF is known to decrease the contrast, the firstorder levelset MCF shrinks the scale of image features. Their secondorder counterparts as we studied in this work are able to preserve these features. We also note, however, that secondorder equations are numerically superior when compared to their firstorder counterparts.
Notation  Description 

with only spatial variable , is timeindependent  
with only temporal variable , maps to either real values,  
or elements in some function spaces that is  
is both space and timedependent  
&  firstorder and secondorder 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 
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 timeindependent and timedependent 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 firstorder total variation flow
We start by reviewing the firstorder 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
(2.1) 
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
(2.2) 
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
(2.3) 
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
(2.4) 
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 semicontinuous on the Hilbert space .
Given the following minimization problem
(2.5) 
the firstorder TVF is nothing but the negative gradient flow for minimizing (2.5), which reads
(2.6) 
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:
(2.7) 
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 firstorder 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:
(2.8) 
Identifying the time step with the regularization parameter in (2.1), that is , we see that (2.8) is in fact the EulerLagrange 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 secondorder total variation flow
Following the idea of damped secondorder dynamics for gradient flows of convex functionals in Hilbert spaces, we introduce the following secondorder TVF:
(2.9) 
where is the socalled 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
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

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

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

for all .

For all :
(2.11) 
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).
Proof. We first consider the following approximate problem with fixed :
(2.12) 
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 firstorder dynamical system in the phase space , i.e.
(2.13) 
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 CauchyLipschitzPicard (CLP) theorem (see e.g. ) for firstorder 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
(2.14) 
Defining a Lyapunov function of the differential equation (2.12), that is
(2.15) 
it is not difficult to show that
(2.16) 
by considering (2.12). Integrating both sides in (2.16), we obtain
(2.17) 
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
(2.18) 
On the other hand, according to the assertion in Proposition 2.1, a constant exists such that for all ,
(2.19) 
Note that it follows from (2.16) that is a nonincreasing 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
(2.20) 
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 :
(2.21) 
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 semicontinuity 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 semicontinuity of , we arrive at
Let , and . Then we have shown that and for all . This concludes the proof.
We remark that for the firstorder TVF (2.6), using tools from semigroup 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 secondorder TVF (2.9) as it is a nonlinear wave equation, and particularly the semigroup 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:
(2.22) 
Using integration by parts and the initial conditions , equation (2.22) becomes:
(2.23) 
Then (2.23) is explicitly written as
(2.24) 
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 nonzero 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:
(2.25) 
Now we return to the righthand 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 secondorder 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 secondorder TVF (2.9), then
Proof. We adopt an idea of . Let us first introduce the auxiliary function
(2.26) 
By elementary calculations, we derive that