Definition of an image

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

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 ob-served 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 convo-lution 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.

Chapter 2

The ROF model

This chapter will cover some of the pioneering methods that are based on vari-ational 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 minimiza-tion 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)

(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

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 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

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ψ) = integrating by parts (strong formulation) yields

Z Then for any ψ, we reach the nonlinear partial differential equation with homo-geneous 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 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).