Image Denoising Algorithms Based on the Dual Formulation of Total Variation

(1)

Image Denoising Algorithms Based on the Dual Formulation of Total Variation

Master thesis in applied mathematics

Christoffer A. Elo

Department of Mathematics University of Bergen

April 30, 2009

(2)

Forord

I min mastergrad har jeg studert bildebehandling med retning innen støyfjerning.

Jeg vil gjerne takke min hovedveileder Alexander Malyshev, som introduserte meg til dette feltet og har gitt meg verdifull hjelp underveis. Jeg vil ogs˚a takke min biveilder Talal Rahman som alltid har hatt svar p˚a mine spørsm˚al.

En takk g˚ar ogs˚a til Sveinung Fjær og Tor Harald Sandve for all hjelp og faglige diskusjoner.

Avslutningsvis takker jeg min samboer Eldbjørg, som alltid har gitt meg støtte underveis.

Bergen, April 30, 2009 Christoffer A. Elo

(3)

Introduction

Image processing with the aid of computers has been a prominent research subject for the last decades. It started with basic statistic methods early in the seventies with applications for medical imaging, satellite imagery and photo enhancement. Real-time processing and more advanced methods based on partial differential equations (PDEs) started to evolve in the eighties, along with a better data performance and computational efficiency. With todays modern computers and dedicated hardware for image processing, the applications for image processing are still growing and creating more complicated problems that are yet to be solved.

Among the various image processing applications there is a strong need for mathematics. The image processing field covers a lot of different areas in mathematics, and solving image processing applications, using e.g nonlinear PDEs, requires the study of function spaces in order to get meaningful phys- ical solutions. To be able to arrive at the nonlinear PDE equations, different mathematical frameworks may be needed: optimization techniques, functional analysis, convex analysis, variational methods, statistical analysis, parameter estimation, etc. As one can see, image processing based on PDE can be very sophisticated, and when a nonlinear PDE is established, there is a crucial need for fast computational algorithms. This includes important fields such as: numerical linear algebra, numerical analysis, finite difference/elements, multilevel methods, fourier analysis, approximation theory, etc.

The different image processing applications of interest, that are based on PDEs, include image segmentation [32, 30, 18], image registration [7, 22], image in-painting[3, 28, 38] and image denoising [33, 36, 40, 10, 9, 11]. The latter will be studied in this thesis. For further introduction to the different applications for image processing based on PDEs, see the books on the topic [2, 17] and references therein. Other methods are of cource also possible, and statistic methods and wavelet theory have shown to produce promising results, see the book [17] for furher introduction to other ways of solving applications from image processing.

The image denoising problems that are solved in this thesis uses techniques that aim to minimize a given energy functional, that describes the solutions to the problem. These problems are typically ill-posed inverse problems, i.e we observe the input image, and the main challange is to smooth out the noise while preserving the main features of the true image. The observed image is normally

(5)

given by an input image, which is composed by the right-handed terms:

The above illustration shows that finding the true image requires knowing the blurring kernel and the noise data. This is not easily solved, as we only have one equation and three unknown variables. The kernel operator is often known, e.g the blurring kernel, and this can be approximated. But the noise is often only given by some statistical estimates, such that the noise is zero-mean white Gaussian noise of an estimated variance. By the nature of this problem we can only find an approximated solution.

To find any meaningful solutions to the problem, there is a need for addi- tional assumptions. The variational approach is to regularize the solution, and in the early nineties, the total variational was proposed by the authors Rudin, Osher and Fatemi (ROF model). This was a huge breakthrough, since total variation preserves interesting features in the solution, such as sharp edges that are important to recover.

The thesis is arranged in the following way. First we introduce a brief overview of the mathematical framework, this includes the function spaces, differential operators and how to discretize the nonlinear equations that arise from the functional minimization.

The next chapter will give a brief introduction on how to restore a noisy image based on partial differential equations. These PDE methods evolve in time and depend on the the derivatives of different orders. A famous example is the heat equation, that smoothes the image by diffusion. It will be necessary to introduce nonlinearity, in order to keep the edges in the image. Based on the difficulties from the PDEs, the minimization of an edge preserving energy functional is given. The classical ROF model that uses the total variation regularization, with a L² fidelity term, is given with theoretical results such as existence and uniqueness. The corresponding Euler equation, that minimizes the energy functional, is then derived into a set of nonlinear partial differential equations. An overview of the drawbacks that total variation based restoration introduces into the solution, is also precented. Motivated by the numerical dif- ficluties for the nonlinear term, a dual formulation is solved by the Chambolle iteration proposed in [9], that will improve the speed of convergence.

The second-to-last chapter is the main contribution in this thesis. This chapter is devoted to two models that reduce the staircase effect: the Fourth-Order model [10, 12, 26] and the TV-Stokes model [38, 34, 25, 21]. The latter is a two- step approach that decouples the Fourth-Order model into two second-order problems, since higher-order methods tend to have computational difficulties due to very large conditioning. The first step of the TV-Stokes model smoothes the tangential vector field with the condition that the vector field is divergence free. Once the smoothed tangential vector field is obtained, the second step reconstructs the image by fitting it to the normal vector field. Both the models

(6)

are derived in the dual formulation to improve the convergence speed, and solving the dual formulation is the main novelty in this thesis, see [21]. The discrete algorithm for the dual TV-Stokes exploits the fast Fourier transformation to solve the divergence free condition.

Finally, the last chapter is devoted to numerical experiments. This chapter discusses the numerical behaviour of the proposed dual TV-Stokes model and the choice for the smoothening parameter. Restored images are shown for the different denoising algorithms, aswell as difference images and contour plots to illustrate the denoising quality. The fast convergence for the dual TV-Stokes is compared to the very slow primal TV-Stokes [34, 25], that is solved by explicit schemes.

1.1 Mathematical preliminaries

This section will briefly cover some of the necessarily tools that are needed for modern image processing, based on partial differential equations. The following text is a standard topic in functional analysis, and if terms like function spaces are unknown to the reader, (s)he should refer to books on the topic, see for instance [37, 19]. Function spaces are very useful in image processing since they possess different features that represent an image. Lets start by recalling the Euclidean space followed by the inner product and some properties.

Definition 1 (Vector space). Euclidean n−space Rⁿis the set of all the n−tuples of real numbers. If x denotes a vector in Rⁿ then it is written as

x = [x1, . . . , xn]^T (1.1) with the length, called the Euclidean norm, defined as

p(x, x) = kxk = q

x²₁+ . . . + x²_n (1.2) Definition 2 (Inner product space). The quantity (x, y) =Pn

i=1xⁱyⁱ is called the inner product. Let x, x₁, x₂ and y, y₁, y₂ be vectors in Rⁿ then the following properties must hold

1. (x, y) = (y, x) (symmetry)

2. (ax, y) = (x, ay) = a(x, y) (bilinearity) (x1+ x2, y) = (x1, y) + (x2, y)

(x, y1+ y2) = (x, y1) + (x, y2)

3. (x, x) ≥ 0, and (x, x) = 0 iff x = 0 (positive definiteness) Proposition 1 (Properties of the norm). If x, y ∈ Rⁿ and a ∈ R

1. kxk > 0 if x 6= 0 2. |(x, y)| ≤ kxkkyk 3. kx + yk ≤ kxk + kyk 4. kaxk = |a|kxk

(7)

Definition 3 (Hilbert space). A Hilbert space is a complete inner-product space.

The L^pspace is a Hilbert space with the following norm

kf k_p=

Z

Ω

|f |^pdx

^1/p

. (1.3)

Example 1 (L²-space). L² is a Hilbert space with the following inner-product (f, g) =

Z

Ω

f g dx. (1.4)

1.1.1 Linear mappings

Linear mappings are the tools for manipulating images. They may be described as transformations, e.g to rotate an image by 90 degrees. A more useful transformation in this thesis is the differentiation map, which will be the key for minimization later on.

Consider two vector spaces X and Y and a mapping f : X → Y , then this mapping is linear if

f (αx + βy) = αf (x) + βf (y), ∀x, y ∈ X and α, β ∈ R (1.5) If A ⊂ X and B ⊂ Y then the image and inverse image is defined by

f (A) = {f (x) : u ∈ A}, f⁻¹(B) = {x : f (x) ∈ B} (1.6) The mapping is called a functional if Y is the scalar field and an operator if Y is a vector space.

Definition 4 (Differentiation operator). A operator f : Rⁿ → R^m is differentiable at a ∈ Rⁿ if there is a linear transformation λ : Rⁿ → R^m such that

lim

h→0

kf (x + h) − f (x) − λ(h)k

khk = 0 (1.7)

The linear transformation λ = (∇f (x), h) is denoted by Df (x) and is called the derivative of f at x.

By the above definition for f : R → Rⁿ and for any h 6= 0,

f⁰(x; h) = (∇f (x), h), ∀h (1.8) In particular, for j = 1, . . . , n

(∇f (x), e_j) = lim

h→0

f (x + hej) − f (x)

h = ∂f

∂ξj

(x) (1.9)

ej is the vector of the jth row in a n × n-identity matrix, and ξj is the jth component of f (x) so that

∇f (x) = ∂f

∂ξ1

(x), . . . , ∂f

∂ξn

(x)

(1.10)

(8)

which is called the gradient of f . The second-order derivative of f is called the Hessian matrix and is denoted by D². The trace of this matrix gives the Laplacian

∆f (x) = tr(D²f (x)) =

n

X

i=1

∂²f

∂ξ²_i (1.11)

If f : Rⁿ → Rⁿthen the n×n matrix of Df (x) is called the Jacobian matrix.

The trace of the Jacobian matrix is called the divergence div f (x) = tr(Df (x)) =

n

X

i=1

∂fi

∂ξi

. (1.12)

Note that the Laplacian is also given by div ∇f (x) = ∆f (x).

It is important to introduce a space that possesses derivatives of certain order in order to search for meaningful solutions when solving problems based on differential equations.

Definition 5 (The function space H¹(Ω) is called a Sobolev space).

H¹(Ω) =f ∈ L²: Df ∈ L²(Ω)

(1.13) In the above definition the Df is a distributional derivative such that

Z

Ω

Df φ dx = − Z

Ω

f Dφ dx, φ ∈ C₀¹(Ω) (1.14) Theorem 1 (The Sobolev space H¹(Ω) is complete).

If f ∈ L¹(Ω) and Df ∈ L¹(Ω) in a distributional sense, then the finite total variation is given by

Definition 6 (Total variation of f ). If f ∈ L¹(Ω) and Df ∈ L¹(Ω) Z

Ω

|Df | dx = sup

Z

Ω

f div g dx : g = (g1, . . . , gn) ∈ C₀¹(Ω)ⁿ, |g| ≤ 1

(1.15) C₀¹(Ω) is the set of continuously differentiable function with compact support(functions that vanish on the boundary) in Ω.

Thus, total variation can be seen as a measure of the amount of oscillation in the function f .

Functions that have a finite value for theR

Ω|Df | dx norm are said to have bounded variation. The space of all these functions is called the space of bounded variation and has the following definition

Definition 7 (BV (Ω) the space of function of bounded variation).

BV (Ω) =

f ∈ L¹(Ω) : Z

Ω

|Df | dx < +∞

(1.16) If f ∈ C¹(Ω), then with integration by parts in equation (1.15) yields

Z

Ω

f div g dx = − Z

Ω n

X

i=1

∂f

∂xi

g dx. (1.17)

Thus the adjoint of the gradient is the negative divergence

(∇f, g) = (f, −div g). (1.18)

(9)

1.1.2 Optimization and Convexity

Convexity plays an important role in optimization. It gives the problem a sufficient condition for uniqueness. The references [20, 35] give a good introduction on optimization and convex problems.

In an Euclidean space, a set is convex if every point on a straight line lies within the set. This gives the mathematical definition of a convex set

Definition 8 (Convex set). A subset V of a vector space is convex iff for every pair of elements u, v ∈ V their closed segment [u, v] is also in V .

A strict convex space is e.g the k · k₂unit ball, while the unit balls of k · k_∞ and k · k₁ are not strictly convex, since they contain a line segment on the boundary, thus the name borderline convex.

Definition 9 (Convex functional). Let f be a functional from Rⁿ to [−∞, ∞], then f is convex if and only if

f (λu + (1 − λ)v) ≤ λf (u) + (1 − λ)f (v) ∀λ ∈ [0, 1] (1.19) if the inequality is strict, then f is strictly convex.

Note that a sum of convex functions is still convex, and a sum of a strictly convex function and a convex function is strictly convex.

Definition 10 (Epigraph of a function). It is the set of points of V which lies above the graph of f .

{(x, λ : x ∈ V, λ ∈ R, λ ≥ f(x)} (1.20) It can further be shown that f is convex iff its epigraph is convex.

1.1.3 Discretization

Discretization of our problems is really the practical part of the functional analysis, the following books give a good introduction on numerical linear algebra and analysis, [39, 24, 20]. Numerical algorithms and methods are used to find solutions to the mathematical problems that are typically described in the continuous setting. It is therefore natural to introduce a discrete framework to place images. The image is described in a two dimensional space X with a fixed size N × N . We denote X by using the Euclidean space R^{N ×N} and Y will express the vector X × X. The space X will be equipped with the Euclidean norm from (1.2)

kxkX =



 X

1≤i,j≤N

|xi,j|²





1/2

(1.21)

and Y has the following Euclidean norm

kykY =



 X

1≤i,j≤N

y_i,j¹ + y²_i,j

2





1/2

(1.22)

where y = [y¹, y²]^T.

(10)

An easy way to implement differential operators, as discussed in the previous section, is to replace the operators with finite difference, i.e replace h⁺→ 0 to a fixed parameter h > 0, giving the forward difference in the one-dimentional case where f ∈ R^N

D^h₊f (x) = (f + h) − f (x)

h (1.23)

approaching the differential from h⁻→ 0 gives the backward difference D₋^hf (x) = f (x) − f (x − h)

h . (1.24)

Another way to approximate the derivative would be to evaluate both sides, giving the central difference

D₀^hf (x) =D^h₊f (x) + D₋^hf (x)

2 = f (x + h) − f (x − h)

2h (1.25)

To introduce a discrete version of the gradient defined in (1.10), a natural choice is to use the forward difference in both directions, if d ∈ X then the gradient ∇d is a vector in Y = X × X and is given by

(∇d)_i,j=

(∇d)¹_i,j (∇d)²_i,j

(1.26) where

(∇d)¹_i,j= D₊^h^xd(x, y), if i < N (1.27)

(∇d)²_i,j= D^h₊^yd(x, y), if j < N (1.28) It is useful to introduce a differential matrix B that approximates the differential operator. This useful since the notation will not be so cluttered as with finite differences.

(∇d)¹_i,j = dB_x^T (1.29)

and

(∇d)²_i,j= Byd (1.30)

where Bx(By) stands for differentiation in the x (resp. y) direction. The forward differential (N ) × (N − 1) matrix takes the following form

B = 1 h







−1 1

−1 1 . .. . ..

−1 1







(1.31)

and the adjoint of B is defined as the backward difference (N − 1) × (N ) matrix

B^∗= −B^T = 1 h





 1

−1 1 . .. . ..

−1 1







(1.32)

(11)

The discrete version of the total variation given in equation (1.15) is then easily followed by

J (d) = X

1≤i,j≤N

|(∇d)_i,j|. (1.33)

The discrete version of the divergence operator can be defined as an analog in the continuous setting in equation (1.18), and thus for every p ∈ Y and d ∈ X, (−div p, d)X = (p, ∇d)Y. We will use the backward difference to get symmetry for the discrete divergence operator, hence

(div p)_i,j= −p₁B_x− B_y^Tp₂. (1.34) Staggered grid

As one notices from the Taylor expansion of the centred approximation given in (1.25), this yields second order accuracy

D^h₀f (x) = f (x + h) − f (x − h)

2h = Df (x) + O(h²). (1.35) However, this centred approximation evaluates the difference across two cells which can lead to undesirable approximations for short wavelengths (high fre- quency). A better way to approximate first-order derivatives is to consider staggered grid, or half-length approximation x_i+1

2 = xi+^h₂ f (x + h) − f (x)

h = Df

x + h

2

+ O(h²). (1.36) This is just a matter of re-defining the grid as shown in figure 1.1.

Figure 1.1: Staggered grid in one dimension

1.1.4 Definition of an image

A digital image can be defined as a function, d, that is bounded and piecewise smooth on an open subset Ω ⊂ R²where Ω usually is square. d(x, y) represents a pixel at the space coordinates (x, y). Is also useful to see the function d ∈ X as a matrix dij where each component in the matrix is the finite gray-scale value, dij = d(x = i, y = j), varying in the range from 0 (black) to 255 (white). In practice, the gray-values are usually normalized into [0, 1].

To obtain a digital image from the continuous world, there is a need for discretization, also known as sampling and quantization. This is done by su- perimposing a regular grid on a continuous image, and each pixel will have the value of the average value from the continuous image. The regular grid will have a size (resolution) that is important for the quality of the obtained digital image. Low resolution will give blocky images, while the computation with these images will be fast, since the dimension of the matrix will be lower.

Higher resolution of the discretized image will give a better approximation of

(12)

the continuous image. The problem with discretization of a continuous image is the natural superimposed noise due to defects in the sensors, transmission problems, interference, etc. This can be mathematically described as the observed image d0 that consists of the original image d perturbed by an additive unknown random noise variable η,

d0= d + η. (1.37)

The problem is given by the observed image d0 and the assumption that η is Gaussian white noise. The values ηi,jare independent random values, each with Gaussian distribution with zero-mean and an estimated variance kηk²₂ ≈ σ². The problem is to reconstruct d from (1.37) which is an inverse problem, and one can easily see that the problem is ill-posed, meaning that one could only get an approximated solution for the given problem.

Another problem that can occur is e.g a blur in the observed image created by e.g incorrect lens adjustment, this can be mathematically described as

d₀= Rd + η (1.38)

where R is a convolution. In the rest of this thesis, we will consider the convolution as an identical map R = I, i.e pure noise as described in (1.37).

There are many features in an image such as edges (discontinuities), flat areas (zero gradient), smooth areas and textures (highly oscillating patterns).

A good denoising algorithm should try to keep these features in the image, while effectively removing the noise.

(13)

Chapter 2

The ROF model

This chapter will cover some of the pioneering methods that are based on variational approaches. The first section will give the primal formulation of the ROF model [36, 40, 10] and the last section will present the dual formulation [16, 8, 9, 13] which increases the computational efficiency.

We start by motivating for finding a solution to (1.37), a good start is to assume a smoothness of the solution. This is a well known processes called Tikhonov regularization of the inverse problem. The idea is to build a functional with a smoothening term and a data fitting-term and then try to minimize this functional

inf

d f (d) = inf

d

Z

Ω

|∇d|²dx + λ Z

Ω

|d − d0|²dx. (2.1) The minimizer infdf (d) will have a solution in the Sobolev space H¹(Ω), and it is unique since f (d) is strictly convex.

The Euler equation, see for instance [20, Theorem 7.2-4], for this minimization is the following

−∆d + λ(d − d0) = 0 in Ω

∂d

∂n = 0 in ∂Ω. (2.2)

An easy way to solve the above equation is to march forward in time

∂d

∂t = ∆d − λ(d − d0) in Ω, (2.3)

this corresponds to the standard heat equation with a force term. The solution is obtained when the iteration is stationary ^∂d_∂t → 0. However, the resulting image does not preserve interesting features in the image, since the method smoothes out the edges and overblurrs the image, as illustrated in figure 2.1.

This is a huge drawback, as the solution d is C^∞(Ω) for λ = 0, meaning that the heat equation has infinite speed of propagation.

Over the two last decades there has been a lot of research in image processing to restore the image d in (1.37), and one quite effective way is to introduce a nonlinear diffusion to preserve the edges. Perona and Malik introduced the diffusion coefficient g(|∇d|) in their classical paper [33]

∂d

∂t = div (g(|∇d|)∇d) in Ω, (2.4)

(14)

(a) Noisy image (b) λ = 0 (c) λ = 0.1

Figure 2.1: Heat equation with and without a force term

where g will be small when the gradient is large (detecting an edge) and large when g detects smooth areas.

2.1 Primal Formulation

Another and more popular approach is to use variational formulation, which has been successful and is still one of the most active areas in image processing.

Variational formulation normally means minimizing an object functional subject to a fidelity term, and with image denoising this can be formulated as

min

d J (d) subject to kd − d0k²₂= σ², (2.5) where J will be the regularity term and is bounded below (lower semi-continuity).

The last term ensures that the given image d0is close to d and this term should be proportional to the noise-level σ²

kd^∗− d0k²₂≈ σ², (2.6) where d^∗denotes an approximate solution to (2.5). σ is supposed to be a known estimate of the variance to the error data, i.e kηk₂= σ.

Choosing J (d) equal toR

Ω|∇d|²dx (H₂¹ semi-norm) will smooth the image very effectively including the discontinuous lines in the image which results in a blurry restored image. To preserve the edges, Rudin et al. suggested in [36]

to lower the H₂¹to the H₁¹(weak derivatives in L¹), which defines the functions of bounded variation BV (Ω), recall the definition in (1.16). Thus a function in BV (Ω) is a function defined in L¹(Ω), whose distributional derivatives are a finite total variation over Ω. The total variation (1.15) will be restated here

k∇dk1= sup

Z

Ω

d div ξ(x) dx : ξ ∈ C₀¹(Ω, R²), |ξ(x)| ≤ 1∀x ∈ Ω

. (2.7) The above functional does not penalize the discontinuities in d, and therefore recovers the edges from the original image. Integration by parts yields

mind

Z

Ω

|∇d| dx subject to kd − d0k²₂= σ², (2.8) which is the famous formulation that Rudin, Osher and Fatemi published in 1992. This model is of great value to image processing, since images tend to have

(15)

discontinuities. Another important task in image processing, which is based on the similar model, is when an image is damaged by missing information. The process is then to fill the missing parts in the image with information from surrounding areas. This is called image in-painting, and pioneering works can be found in [3, 28].

Chambolle and Lions showed in [10] that (2.8) is naturally linked with the following unconstrained problem

min

d f (u) = min

d

Z

Ω

|∇d| dx +λ

2kd − d0k²₂, (2.9) where f (d) is the objective function and λ is a non-negative Lagrange multiplier.

This formulation is known as the TV-L² model, and is the convexification of (2.8) when the constraint is equal to kd − d₀k²₂≤ σ².

The λ parameter balances between two terms λ → ∞ and λ → 0. If λ → ∞, the solution is given by d = d0, which is obviously not a good result. If λ → 0, then the solution is zero, or any constant, since ∇d = 0 and does not fit the image with respect to d0. The optimal λ should take advantage of σ such that the constraint kd − d0k = σ is fulfilled.

The following results show that the ROF model has a unique solution, which is important when searching for solutions. Uniqueness and existence among other theoretical results concerning bounded variation, can be found in [10, 1].

Proposition 2 (f (d) has a unique solution.). Suppose that f is coercive and that f is lower semi-continuous, then the problem has a solution. The solution is unique if f is strictly convex.

Proof. Let dn be a minimizing sequence of (2.9) such that f (d_n) → inf

d∈L²(Ω)

f (d) = α, n ≥ 1. (2.10)

Suppose kd₀kL² = 0 and that all the constrains are satisfied by d_n. f (d) is coercive: lim f (d) → +∞ for kdk → ∞. Due to Poincar´e inequalities, d_n is bounded in BV (Ω), see e.g [1, Theorem 2.5]. Thus, an extracted sequence from d_n converges weakly to some ¯d ∈ L²(Ω). f is lower semi-continuous on L², hence

f ( ¯d) ≤ lim

n_i→∞f (dn_i) = α, (2.11) where ¯d is a solution of (2.8) and α 6= −∞.

If f is strictly convex, two solutions is impossible due to f d₁+ d2

2

<1

2(f (d1) + f (d2)) = α. (2.12)

To derive a minimization of f (d) one introduces an admissible variation ψ of d, i.e ψ ∈ H¹(Ω), which has no constraints

φ : t ∈ I → φ(t) = f (d + tψ). (2.13) φ will be minimal at t = 0 and φ⁰(0) = f⁰(d)ψ

0 ≥ lim

t→0⁻

φ(t) − φ(0)

t = φ⁰(0) = lim

t→0⁺

φ(t) − φ(0)

t ≥ 0, (2.14)

(16)

which shows that f⁰(d)ψ = 0, and the ψ is arbitrary, therefore f⁰(d) = 0. This is known as the Euler equation, see [20] for further introduction on the topic.

The Euler equation of (2.9) is the following

φ(t) = f (d + tψ) = Z

Ω

|∇d + t∇ψ| +λ

2|d − d0+ tψ|²dx

= Z

Ω

|∇d + t∇ψ| +λ

2|d − d0|²+ λtψ(d − d0) +λ

2t²ψ²dx. (2.15) Using the fact that φ⁰(t) will give minimum at φ⁰(0) = 0

φ⁰(t) = Z

Ω

∇d + t∇ψ

|∇d + t∇ψ|, ∇ψ

+ λψ(d − d0) + λtψ²dx (2.16)

ψ⁰(0) = Z

Ω

∇d

|∇d|, ∇ψ

+ ψλ(d − d0) dx = 0, (2.17) integrating by parts (strong formulation) yields

Z

∂Ω

ψ∂d

∂ndx + Z

Ω

ψ

−div ∇d

|∇d|+ λ(d − d0)

dx = 0. (2.18) Then for any ψ, we reach the nonlinear partial differential equation with homogeneous Neumann boundary condition

div_|∇d|^∇d − λ(d − d0) = 0 in Ω

∂d

∂n = 0 in ∂Ω. (2.19)

Note that for |∇d| = 0, (2.19) is not well-defined, this is a problem since the solutions of the ROF model have large flat areas, i.e where ∇d = 0 in the image.

One way to solve this is to regularize f (d) inf

d f (d)^β= Z

Ω

p|∇d|²+ β + λ

2|d − d0|²dx, (2.20) where β > 0 implies that f (d)_β is differentiable and this yields the regularized Euler equation

( g(d) = div√ ^∇d

|∇d|²+β − λ(d − d0) = 0 in Ω

∂d

∂n = 0 in ∂Ω, (2.21) with homogeneous Neumann boundary condition and initial solution, given by d = d0. g(d) = 0 is a necessary and sufficient condition for d to be a solution of the convex minimization problem (2.20).

2.1.1 Drawbacks

Total variation restoration is very good for recovering images from noise. Unfor- tunately it may fail in some areas to give a satisfactorily result, due to artificial effects that are introduced while performing pure total variation with a L² fidelity term. These drawbacks can be found in the survey paper [11], and they will be restated here.

(17)

Regularization The first disadvantage one notices about the ROF model is the β-regularization. The drawback is that the smaller β responds to more iterations before convergence, and this makes the ROF model rather slow when solving for time marching schemes. Choosing higher β gives faster convergence, but smooths the image needlessly.

Staircase It is well known that the ROF model produces staircase effect for piecewise smooth areas. In figure 2.2 it is easy to see why this happens. The figure consists of a decreasing function from [5, 25], and when noise is added, the function becomes piecewise constant for the decreasing parts. When running the ROF scheme, the nondecreasing parts will be smoothed out as illustrated in the figure below.

(a) Original function (red-line) and noisy function (blue-line).

(b) Denoised function (blue-line) with the staircase effect.

Figure 2.2: Illustration of staircase effect

Loss of texture Another disadvantage when using total variation based methods, is that texture is easily smoothed out during the denoising. In a recently proposed model by Osher et al. in [31], they prevent this by iterative regularization, which means running the ROF model multiple times, and each time adding the difference between the observed noisy image and the previous one, i.e di+1 = di+ (d0− di). Thus, the noise (and the details!) are added back to the image.

Loss of constrast Loss of contrast is also a limitation for the ROF model. If the input image d0(x, y) is a disk, then the solution is given by αd(x, y), where α ∈ [0, 1), i.e the solution is never α = 1. This disadvantage can also be seen in the one-dimensional example illustrated in figure 2.2, where the solution f (x) is in the range (0, 2).

One approach to this problem is given in [14], where they weaken the fidelity norm. They formulated the problem to measure the L¹-norm between the observed and denoised data. The TV-L¹ minimization problem is then given by

infd

Z

Ω

|∇d| +λ

2|d − d0| dx. (2.22)

The above formulation is not strictly convex, which means that there is no unique global minimizer. The reason why this formulation has been studied, is

(18)

that it has some disirable advantages, e.g that the solution is contrast invariant.

This means that if d is a solution for the problem (2.22), then αd is also a solution for αd0, where d0 is the initial input image and α is a scaling factor.

2.1.2 Numerical methods

There are various ways of solving (2.21), and in [36] the authors Rudin et al.

suggested to forward the equation with artificial time to get a gradient descent scheme

dⁿ⁺¹= dⁿ− kg(dⁿ), n = 0, 1, 2, . . . , (2.23) as n → ∞. When steady-state is achieved ^∂d_∂t → 0, the solution dⁿ⁺¹will satisfy the unique solution of the minimizer given in (2.8). The timestep k must be chosen to satisfy the Courant-Friedrichs-Lewy (CFL) condition k ≤ c|∇d|, for some constant c > 0. Due to this condition, the method converges rather slowly, since k must be chosen small. This is due to that |∇d| is approximately zero in flat regions.

The main difficulty that equation (2.21) poses is the linearization of the highly nonlinear term div (∇d/|∇d|β). Thus, one way to overcome the CFL condition is to linearise (2.21) by using the nonlinear term 1/|∇d|β from the previous iteration. This was proposed in [40] and is known as lagged diffusivity fixed point iteration. The term dⁿ⁺¹can then be solved by the following sparse system of linear equations

dⁿ⁺¹= d⁰+1

λdiv ∇dⁿ⁺¹

|∇dⁿ|β

. (2.24)

This gives a faster iteration than (2.23), but reacts more to the β-regularization due to the highly nonlinear term. Vogel et al. also solve the above iteration with multigrid methods, but unfortunately, standard finite differences or finite elements discretization yield disappointing results when the input image d₀ is not sufficiently smooth. Thus, by achieving O(n) complexity, higher β- regularization is needed such that the nonlinear term is well behaved.

2.2 Dual Formulation

As seen in the previous section, a major drawback of the ROF model is that the formulation is non-differentiable in zero, due to the total variation term.

To overcome further regularization, Chambolle in [9], Carter in [8] and Chan et al. in [16], studied the so called dual formulation for the ROF model. We will reduce the above presented primal formulation to the dual formulation. Lets consider the more general definition of the total variation, given in equation (2.7), and use the right-handed definition to obtain the inf − sup problem

infd sup

|p|≤1

h(d, p) = inf

d sup

|p|≤1

Z

Ω

(d, div p) + 1

2λ|d − d0|²dx

. (2.25)

Hence h( ¯d, ¯p) is a saddle point since ¯d is the minimum of the function h( ¯d, ·), and ¯p is the maximum of the function h(·, ¯p), such that

sup

|p|≤1

f ( ¯d, p) = f ( ¯d, ¯p) = inf

d (d, ¯p). (2.26)

(19)

Given any elements ˆd and ˆp gives the following inequality inf

d f (d, ˆp) ≤ f ( ˆd, ˆp) ≤ sup

|p|≤1

f ( ˆd, p), (2.27)

and as infdf (d, ˆp) is a function of the single variable ˆp, and ˆd is the single variable for sup_|p|≤1f ( ˆd, p), this yields that the following inequality exists

sup

|p|≤1

inf

d f (d, p) ≤ f (d, p) ≤ inf

d sup

|p|≤1

f (d, p). (2.28)

The fact that f ( ¯d, ¯p) is a saddle point from equation (2.26), gives that the necessary converse inequality must also hold. This is of great value for solving (2.25), since it does not matter which problem is solved first. Thus interchanging the sup and inf from equation (2.25) yields the inner minimization

sup

|p|≤1

inf

d h(d, p) = sup

|p|≤1

inf

d

Z

Ω

(d, div p) + 1

2λ|d − d0|²dx. (2.29) Hence the equation is piecewise continuous and differentiation w.r.t d is easy

d = d0− λdiv p. (2.30)

The term λdiv p = πλ is the nonlinear projection that is to be solved and inserted into (2.30) instead of solving the primal problem. This projection can be solved by substituting (2.30) into (2.29) and deriving the following maximization problem

sup

|p|≤1

h(p) = sup

|p|≤1

Z

Ω

(d0, div p)−λ

2|div p|²dx = sup

|p|≤1

1 2λ

Z

Ω

|d0|²−|λdiv p − d0|² dx, (2.31) using sup {h} = − inf {−h} yields the following minimization problem for the dual variable

minp kdiv p − λ⁻¹d0k²₂: |p| ≤ 1 . (2.32) In one dimension p is pair of two vectors, and the above minimization is easily solved by e.g least-square methods. In our case, p is a pair of two matri- ces. Chambolle solved the minimization problem by considering the optimality condition, also known as the KKT conditions (due to Karush-Kuhn-Tucker, cf.

[20, Theorem 9.2-3])

−(∇(div p − λ⁻¹d₀)) + αp = 0, (2.33) where α is the Lagrangian multiplier associated with the constraint p that must satisfy

α > 0

|p| = 1 or

α = 0

|p| < 0. (2.34)

He then did an important observation, that in either case

α = |(∇(div p − λ⁻¹d₀))|, (2.35) such that the (2.35) can be substituted into equation (2.33) which yields the equation

(∇(div p − λ⁻¹d₀)) − |(∇(div p − λ⁻¹d₀))|p = 0. (2.36)

(20)

2.2.1 Discrete algorithm

By considering the minimization problem in a discrete framework we have the following analog to (2.32)

minp

nkdiv^hp − λ⁻¹d₀k²_X: |p_i,j| ≤ 1, i, j = 1, . . . No

. (2.37)

Deriving with the same arguments as (2.33)-(2.36) and considering pⁿ⁺¹− pⁿ

k = ∇^h

div^hpⁿ− λ⁻¹d0

− ∇^h

div^hpⁿ− λ⁻¹d0

pⁿ⁺¹ (2.38) yields a semi-implicit scheme known as Chambolle’s iteration that is given in [9]

p⁰= 0, pⁿ⁺¹=

pⁿ+ k∇^h

div^hp − λ⁻¹d₀ 1 + k

∇^h

div^hp − λ⁻¹d₀

. (2.39)

The iteration converges for k ≤ 1/8 which is proven in [9]. In practice, convergence is observed by setting k equal to 1/4. This gives a rapid convergence up to 50-100 iterations, and then the speed of convergence decays.

(21)

Chapter 3

Staircase reducing models

Section 2.1.1 pointed out some of the areas where the ROF model may suffer from an inefficient restoration. This chapter will first give a brief overview of some of the proposals from the literature that may improve these problems.

Then the main focus will be on the staircasing effect. Two models that deal with this effect will be presented: the Fourth-Order model [10, 12, 26] and the TV-Stokes model [38, 34, 25, 21]. The latter reference is the novelty in research literature and it is the main contribution in this thesis. Both models will be given in dual formulation and the next section will give numerical experiments for these models. It should also be noted that these problems are still challenging problems in image processing and remain to be solved.

The ROF model treats flat regions and edges well, but favours piecewise constant functions in the solution, known as staircase effect. This is an adverse effect in images that contain regions with gradual image variations. To prevent this over-sharpening there have been some proposals to minimize different regularizations terms.

A more general version of the ROF model may be given as the following minimization problem

infd

Z

Ω

φ(|∇d|) +λ

2(d − d0) dx, (3.1)

where φ : R⁺ → R⁺ is a smooth function. The above minimization has the following regularized Euler equation

div

φ⁰(p|∇d|²+ β) ∇d

|∇d|

− λ(d − d0) = 0, (3.2) where β > 0 is a small parameter to avoid |∇d| = 0. Pure total variation corresponds to φ(g) = g and φ(g) = g² corresponds to the H¹ norm. The question now is how to combine the two norms to benefit each of them.

In [4] they propose

φ⁰(|∇d|) = |∇d|^p−1 (3.3)

and then choose p close to 1.

(22)

The smooth function φ in the minimization problem (3.1) will be more interesting when endowed with the following conditions

(C1) : φ(0) = φ⁰(0) = 0 (3.4)

(C2) : lim

s→+∞

φ(s)

s = 0. (3.5)

The first condition (C1) implies that the first term in (3.1) is quadratic near the origin, thus smoothes the flat areas (|∇d| ≈ 0) with the H¹ norm. The second condition (C2) implies that the function is sublinear in growth at infinity. This property will enhance the edges since the cost of edges are low. Using the above conditions for the smooth function φ results in an ill-posed problem since the function is obviously non-convex. A simple example would be choosing φ(s) = s²/(1 + s²). However, the non-convexity is irrelevant for the finite- dimensional case as stated in [2] on page 93, and the numerical experiments often give visually good results. This strategy is throughly presented in [23]

where the authors also gave the dual formulation with a weighted total variation term of the non-convex problem.

In [5, 4] they proposed φ(|∇d|) = |∇d|^p(|∇d|)to adapt the behaviour of |∇d|, such that p = 1 near the edges and p = 2 in flat regions. In between these regions they use a fractional norm H_p¹, p ∈ (1, 2). A simple example would be p(|∇d|) = _1+2|∇d|² . However, in [5, 4], examples are given in just the one dimensional case.

3.1 Fourth-Order denoising in dual formulation

One of the early approaches to remove staircasing in image denoising, was to introduce higher-order derivatives or norms into the energy. The first approach to combine the TV norm and H¹ norm (square norm of the gradient) is to consider the inf-convolution as the authors Chambolle and Lions proposed in [10]. The resulting minimization is equivalent to

inf

d

Z

|∇d|≥

|∇d| dx + 2

Z

|∇|<

|∇d|²dx +λ 2 Z

Ω

|d − d0|²dx, (3.6) where is a given threshold to separate the smooth and discontinuous parts of the image.

The same method can be used to minimize the total variation of the gradient, i.e a second-order functional. Some of the earliest variants of higher-order methods for image denoising can be found in the same paper [10], where they decompose the image into two parts, d = d1+ d2, where d1consists of discontinuous parts and d2consists of the smooth areas of the image, and then minimize the following problem

inf

d₁,d₂

Z

Ω

|∇d₁| + α|∇(∇d₂)| +λ

2|d₁+ d₂− d₀|²dx, (3.7) such that each component is assigned to the appropriate norm, i.e d₁∈ BV (Ω) and d2∈ H₁¹(Ω) with ∇d2∈ BV (Ω). The results given in [10] are very good for the one-dimensional case, but the model suffers from over-smoothening of the edges in two dimensions.

(23)

This is evident, since the problem is how to combine the two different norms, which are dependent on the gradient, and this operation is unstable.

One of the first approaches to minimize only the higher-order derivaives in the energy was proposed by the the authors Lundervold, Lysaker and Tai (LLT) in [26], where they gave the following minimization problem

Z

Ω

(|dxx| + |dyy|) dx dy subject to kd − d0k2≤ σ (3.8) and another functional that is rotational invariant

Z

Ω

q

|dxx|²+ |dxy|²+ |dyy|²+ |dyx|²dx dy subject to kd − d0k2≤ σ. (3.9) The Euler equations for the two above functionals will result into a set of fourth- ordered nonlinear PDEs.

We will further focus on the Fourth-Order model, and derive a dual formulation of the following minimization problem

inf

d

Z

Ω

(|∇²d|) + 1

2λ(d − d0)²dx, (3.10) where the regularization term is defined by |∇²d| =p|dxx|²+ |dyy|². This term is usually replaced by |∇²d|β=p|∇²d|²+ β², which is a necessary regularization to avoid the non-differentiability. The model shares the same numerical difficulities as the ROF model, and a dual formulation will now be given to overcome the β-regularization and the slow iteration. Similar works can be found in [15].

We can define the regularization term as an analog to the total variation defined in (2.7), hence

k∇²dk₁= sup

Z

Ω

d ∆ξ(x) dx : ξ ∈ C₀²(Ω, R²), |ξ(x)| ≤ 1∀x ∈ Ω

. (3.11) Thus, we use the right-handed definition of the above equation and deduce

inf

d sup

|p|≤1

Z

Ω

(d, ∆p) + 1

2λ(d − d₀)²dx, (3.12) where p is the dual variable. By the convexity of d and the concave dual function p, we can use the same argument as in (2.25)-(2.29) to switch the inf and sup in (3.12)

sup

|p|≤1

inf

d

Z

Ω

(d, ∆p) + 1

2λ(d − d0)²dx. (3.13) Solving for the inner minimization with respect to d yields

d = d₀− λ∆p. (3.14)

Thus, by the simple equation above we find our restored image by looking for the projection λ∆p. We deduce a minimization problem for the dual variable by substituting equation (3.14) into (3.13) and setting sup{·} = − inf{−·}, and end up with the following dual problem

infp k∆p − λ⁻¹d0k²₂: |p| ≤ 1 . (3.15)

(24)

The dual minimization problem is solved by considering the optimality condition for equation (3.15)

∇² ∆p − λ⁻¹d0 + αp = 0, (3.16) where either α = 0 and |p| ≤ 1 or α > 0 and |p| = 1, both cases yield

α = |∇² ∆p − λ⁻¹d0 |. (3.17) Thus, inserting equation (3.17) into (3.16) gives

∇² ∆p − λ⁻¹d0 + |∆p − λ⁻¹d0|p = 0, (3.18) which may be solved in the same manner as the second-order Chambolle iteration

pⁿ⁺¹= pⁿ+ τ ∇² ∆pⁿ− λ⁻¹d0

1 + τ |∇²(∆pⁿ− λ⁻¹d0)| . (3.19)

3.1.1 Discrete algorithm

The continuous iteration from equation (3.19) is discretized by using the matrix operators from section 1.1.3. Thus, the second-order operator ∇²may have the following discrete definition

∇²_hd =

−dB_x^TB_x

−B_y^TB_yd

, (3.20)

and the discrete laplacian operator in the same manner

∆^hd = −dB^T_xBx− B_y^TByd. (3.21) The full discrete version of (3.19) is then given by using the above discrete operators

pⁿ⁺¹=pⁿ+ τ ∇²_h ∆^hpⁿ− λ⁻¹d0

1 + τ |∇²_h(∆^hpⁿ− λ⁻¹d0)| . (3.22) The following result gives a condition on τ for fast convergence.

Theorem 2. Let τ ≤ 1/64. Then λ∆^hpⁿ converges to solution λ∆^hp of the minimization problem given in (3.15) as n → ∞.

A proof of the theorem can be found in [15] where they propose a dual algorithm for the minimization problem given in (3.7).

As an alternative argument for convergence and a better condition for τ , we can consider the following

pⁿ⁺¹= pⁿ+ τ ∇²_h ∆^hpⁿ− λ⁻¹d0 = (1 + τ ∇²_h∆^h)pⁿ− τ ∇²_h∆^hλ⁻¹. (3.23) By denoting the above equation as pⁿ⁺¹ = φ(pⁿ) and requiring φ⁰(pⁿ) ≤ 1 (1-Lipschitz) we deduce

kτ ∇²_h∆^hk − k1k ≤ k1 + τ ∇²_h∆^hk ≤ 1. (3.24) Hence, we get the following condition for τ

τ k∇²_h∆^hk ≤ 2. (3.25)

(25)

To estimate the norm of the operator ∇²_h∆^h, we note that kBk_∞= 2 and the SVD of the operator ∇²_h∆^h have the following form in x (resp. y) direction

∇²_h∆^h= −C^T

0 Σ⁴_x

+

0 Σ⁴_xy

C. (3.26)

Thus, the following bound in x (resp. y) direction can be estabilished

λmax ∇²_h∆^h = λmax Σ⁴_x +λmax Σ⁴_xy ≤ kBk⁴₂+ kBk⁴₂≤ kBk⁴_∞+ kBk⁴_∞= 32, (3.27) giving the optimal condition for τ to be equal to 1/16. Note that the operator ∇²_h∆^h needs to be calculated in each iteration, and this operation is ill-conditioned, i.e the condition number κ(∇²_h∆^h) is close to 10¹⁸.

3.2 TV-Stokes denoising in dual formulation

3.2.1 Motivation

To tackle the numerical problems in higher-order methods, one can split the problem into two steps, and solve a second-order problem in each step. There has been some research activity to formulate a two step method, and Lysaker, Osher and Tai (LOT model) proposed in [27] to smooth the normal vector field in the first step with the following minimization

inf

|n|=1

Z

Ω

|∇n| +δ

2|n − n₀|² dx. (3.28) This will find the smoothed flow field, i.e directions of the gradient: n =

∇d/|∇d| and n0 = ∇d₀/|∇d₀|. Once the normal vector field is obtained, the image can be reconstructed by a second step, which matches the normal vector field by the following minimization problem

inf

d

Z

Ω

|∇d| − (n⁰, ∇d) +µ

2|d − d0|²dx, (3.29) where n⁰= [n1, n2]. The LOT method gives a better result than the ROF, since it is more edge preserving, cf. [27]. However, it still suffers from the staircase effect, since it can not recover images containing affine regions of noise.

Another way to smooth the normal vector fields in step one, which is the crucial part of the two-step algorithm, is to notice that the isophote lines in the image must be constant, see Figure 3.1. Hence, the tangential vector field is divergence free. The first step should then minimize with respect to the tangential vector field given by

τ = ∇^⊥d = [−dy, dx]^T = [v, u]^T, (3.30) subject to the divergence of the tangential vector field is equal to zero, known as incompressibility in fluid mechanics

div τ = 0. (3.31)

This is formulated in [34] by the authors Rahman, Tai and Osher. The model will be called the primal TV-Stokes model to indicate that the energy

(26)

Figure 3.1: Normal and tangent vectors to the isophote lines

functionals are solved with respect to the primal variables τ and d. The first step of the primal TV-Stokes model has the following minimization problem

infτ

Z

Ω

|∇τ | +δ

2|τ − τ₀|²dx subject to div τ = 0. (3.32) where τ0 is given by τ0 = [v0, u0]^T. Once the smoothed tangential vector field is obtained, the corresponding normal vector field is then calculated by n = [u, −v]. The normal vector field can be used to reconstruct the noisy image by fitting it with the following second step minimization

inf

d

Z

Ω

|∇d| −

∇d, n

|n|

dx subject to Z

Ω

(d − d0)²dx = σ². (3.33) The authors then solve the problem by finding the Euler-Lagrange equations for (3.32-3.33) and then iterating the nonlinear PDEs with an explicit gradient descent method until steady-state is reached. This, however, is very slow due to the nature of explicit forward schemes. A modified TV-Stokes model is proposed by Litvinov, Rahman and Tai in [25], where they establish existence and uniqueness to the minimization problems for both steps of the model.

The next sections will reformulate the primal equations given by (3.32) and (3.33) into a dual formulation, which improves the speed drastically.

3.2.2 First step: tangent field smoothing

To derive a dual formulation of the primal problem (3.32) we introduce the definition of the total variation of the tangential vector field

Z

Ω

|∇τ | dx = sup

p∈K

Z

Ω

(τ, div pi) dx : i = 1, 2

(3.34) and K is closure of the convex set

div ξi : ξ_i ∈ C₀¹(Ω, R²), |ξ_i(x)| ≤ 1, ∀x ∈ Ω, i = 1, 2 . (3.35) The dual variable p is now a vector composed by the two dual variables for both directions in the tangential vector field. The divergence of this field then

(27)

given by

div p = (div p1, div p2)^T, div p1= ∂p¹₁

∂x +∂p²₁

∂y, div p2= ∂p¹₂

∂x +∂p²₂

∂y. (3.36) This definition of the divergence is similar to the vectorial dual norm from [6]

for vectorial images, e.g. color images.

Using the dual formulation of the total variation norm, the primal problem (3.32) can be written as minimization of the following dual problem

inf

div τ =0 sup

|p|≤1

Z

Ω

(τ, div p) + 1

2δ|τ − τo|²dx, (3.37) interchange sup and inf with the same argument as given in the section 2.2

sup

|p|≤1

inf

div τ =0

Z

Ω

(τ, div p) + 1

2δ|τ − τo|²dx. (3.38) Let us introduce the orthogonal projection ΠPonto the constrained subspace P = {τ : div τ = 0}. This projection may be given as

ΠP

π1

π2

=

π1

π2

− ∇∆⁺div

π1

π2

. (3.39)

Thus, finding the minimum for the inner problem, without constraint div τ = 0, is obtained by using the property from orthogonal projection, (τ, div p) = (ΠPτ, div p) = (τ, ΠPdiv p), and the Euler equation. The equation (3.38) with the orthogonal projection

sup

|p|≤1

inf

div τ =0

Z

Ω

(τ, Π_Pdiv p) + 1

2δ|τ − τ_o|²dx (3.40) is solved by obtaining the following derivation for the Euler equation

φ(t) = E(τ + tψ, λ) = Z

Ω

(τ + tψ, ΠPdiv p) + 1

2δ|τ − τ0+ tψ|²dx