PDE Based Image Diﬀusion and AOS

Jarno Ralli

October 20, 2012

Abstract

This technical paper shows how AOS (Additive Operator Splitting)

scheme, together with TDMA (TriDiagonal Matrix Algorithm), can be

used for solving a non-linear, scalar valued diﬀusion equation. The paper

shows how a continuous version of the diﬀusion equation is transferred

into a ‘discrete’ version in order to be solved numerically. Much of the

text is based on [3], and is presented here only for the sake of readabil-

ity. Most part of the text deals with deﬁnitions, the idea here being that

anyone familiar with diﬀusion should be able to understand the used no-

tation and, thus, the formulations that follow. This text is supposed to

be accompanied by a Matlab implementation of the code.

1 Used Notation

We consider the image to be a continuous function with I(x, y, k) : R × R →

R

K+

, where the domain of the image is Ω = R × R and K deﬁnes the number

of channels (e.g. K = 3 in the case of RGB images). It is in this regular domain

where our PDEs are deﬁned. However, in order for us to solve the PDEs, they

need to be discretised ﬁrst. Also, the kind of images that we are dealing with

are, in fact, discretised versions that we receive from the imaging devices, such

as digital- or thermal cameras. We deﬁne a discretisation grid as given in (1):

G

h

:= {(x, y) | x = x

i

= ih

x

, y = y

j

= jh

y

; i, j ∈ Z} (1)

where h = (h

x

, h

y

) is a discretisation parameter. With the discretisation grid,

the domain of the discretised images can be deﬁned as Ω

h

= Ω ∩ G

h

. Instead

of I(x, y) = I(ih

x

, jh

y

), we typically use I

i,j

when pointing to the pixels.

Sometimes positions on a grid are given using both cardinal- and inter-

cardinal directions, as deﬁned in Figure 1. The idea is to simplify the notation

of the discretised versions of the equations which can be rather messy.

NW N NE

W C E

SW S SE

Figure 1: Directions on a grid. To simplify the notation, both cardinal- and

inter-cardinal directions are used. Here W, N, E, S, and C refer to west, north,

east, south, and centre, respectively.

www.jarnoralli.ﬁ

www.jarnoralli.com

1

1.1 Pixel Numbering Schemes

As it was already mentioned, we can refer to any pixel in the image by using I

i,j

,

where 0 ≤ i ≤ m, 0 ≤ j ≤ n. Here the discretisation parameter h = (h

x

, h

y

)

has been chosen so that the discretised domain is Ω

h

: [1, m] × [1, n]. Another

particularly useful way is to think of the discretised image as a vector I ∈ R

N

.

Now, the components of the vector are I

I

where I ∈ {1, . . . , N } and N is

the number of pixels in image. This second numbering scheme is useful for

algorithmic descriptions, and in matrix notation, as will be shown later on.

Figure 2 depicts both column- and row-wise traversing order of pixels. De-

pending on the solver, either column- or row-wise traversing, or a combination of

both (for example, ﬁrst column- and then row-wise), can be used. An interested

reader is pointed to [4] for more information.

1

2

3

4

5

6

7

8

9

(a) Column wise.

1

4

7

2

5

8

3

6

9

(b) Row wise.

Figure 2: Column- and row-wise pixel orderings. Here, I ∈

{1, 2, 3, 4, 5, 6, 7, 8, 9}.

.

1.2 Pixel Neighbourhoods

In order to simplify the notation, for example in algorithmic descriptions, we

deﬁne diﬀerent kinds of pixel neighbourhoods. With pixel neighbourhood we

mean image pixels I

J

that are neighbours of a pixel of interest I

I

. The neigh-

bourhoods are slightly diﬀerent depending on if we are talking about a element-

or a block-wise solver. By element wise solver we mean a Jacobi or a Gauss-

Seidel type iterative solvers, that search for the solution for a single element at

a time. On the other hand, block type solvers search for a solution for a group

of elements (or a block). However, since the neighbourhoods have the same

function in both the above mentioned cases, we use the same neighbourhood

operator to denote the neighbours. It should be clear from the structure of the

solver which kind of a neighbourhood is in question. J ∈ N(I) denotes the set

of neighbours J of I, as seen in Figure 3 (a).

2

❥❡① ①

①

①

❥❡①

①

❥❡①

①

①

(a) 4-neighbourhood with eliminated

boundary conditions.

❥❡① ❤

①

❤

❥❡①

①

❥❡①

①

❤

(b) Element wise neighbourhood with

eliminated boundary conditions.

❥❡③ ❥

❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄

(c) Block wise, column ordering.

❥❡

③

❥

✲

✲

✲

✲

✲

✲

✲

✲

✲

(d) Block wise, row ordering.

Figure 3: Pixel neighbourhoods, where central pixel, I, is denoted with a dou-

ble circle. (a) Painted circles denote neighbouring pixels J that belong to the

neighbourhood deﬁned by J ∈ N (I). (b) Painted circles denote pixels J that

belong to the neighbourhood deﬁned by J ∈ N

−

(I), while unpainted circles

denote pixels J that belong to the neighbourhood deﬁned by J ∈ N

+

(I). (c)-

(d) Painted circles denote pixels J that belong to the neighbourhood deﬁned by

J ∈ N

−

(I), while unpainted circles denote pixels J that belong to the neigh-

bourhood deﬁned by J ∈ N

+

(I). Processing order is deﬁned by the arrows.

Due to the eliminated boundary conditions ‘scheme’, the pixel neighbourhood

operators only point to valid neighbours, as shown in (a) and (b).

Pixel wise. J ∈ N

−

(I) denotes the set of neighbours J of I with J < I

(painted circles in Figure 3 (b)), and J ∈ N

+

(I) denotes the neighbours J of

I with J > I (unpainted circles in Figure 3 (b)).

Block wise. J ∈ N

−

(I) denotes the neighbours J of I with J < I (painted

circles in Figure 3 (c)-(d)). Since we run the block wise solver both column- and

row-wise, depending on the direction, the neighbour(s) deﬁned by this neigh-

bourhood operator changes. J ∈ N

+

(I) denotes the neighbours J of I with

J > I (painted circles in Figure 3 (c)-(d)). Again, the actual neighbours deﬁned

by this operator depends on the direction of the solver.

3

2 Scalar Valued Diﬀusion

The basic formula describing a scalar valued diﬀusion is:

∂I

k

∂t

= DIV

g(x, y, t)∇I

k

(2)

where k refers to the channel in question, g(x, y, t) deﬁnes the (scalar) diﬀusion

weights, and ∇ := [

∂

∂x

,

∂

∂y

] is the spatial gradient. Since g(x, y, t) is a function

of t, the diﬀusion is non-linear.

2.1 Implicit Scheme

For discretisation, we use Euler backward, semi-implicit method, and obtain the

following:

(I

k

)

t+1

− I

t

k

τ

= DIV

g

t

∇(I

k

)

t+1

(3)

where g

t

= g(x, y, t). In other words, we use values of g, available at time t,

when resolving for a new value at time t + 1.

3 Finite Diﬀerence Discretisation

3.1 Finite Diﬀerence Operators

Before going any further, we remind the reader of the ‘standard’ ﬁnite diﬀerence

operators, which are used for approximating derivatives. Here, the function of

interest, for which we want to ﬁnd derivatives, is deﬁned as f(x, y) : Ω → R.

1. First order forward diﬀerence is given by:

D

+

x

f(x, y) = f

+

x

(x, y) =

f(x + ∆x, y) − f(x, y)

∆x

(4)

2. First order backward diﬀerence is given by:

D

−

x

f(x, y) = f

−

x

(x, y) =

f(x, y) − f(x − ∆x, y)

∆x

(5)

3. First order central diﬀerence is given by:

D

0

x

f(x, y) = f

0

x

(x, y) =

f(x +

1

2

∆x, y) − f (x −

1

2

∆x, y)

∆x

(6)

4. Second order central diﬀerence is given by:

DD

0

x

f(x, y) = f

0

xx

(x, y) =

f(x + ∆x, y) − 2f(x, y) + f(x − ∆x, y)

∆x

2

(7)

where in f

x

the sub-index (here x) indicates with respect to which variable

the function has been diﬀerentiated. The above formulations can be simpliﬁed

further if we assume a uniform grid with ∆x = ∆y = 1. The diﬀerence operators

shown here approximate derivatives with respect to x. By switching the roles of

4

x and y, corresponding diﬀerence operators for approximating derivatives with

respect to y can be obtained easily. In the case of images containing several

channels (e.g. RGB-images), derivatives are approximated separately for each

channel.

3.2 Discretisation of DIV Operator

Now that we know how to approximate ﬁrst and second order derivatives, we can

discretise the divergence, DIV , operator. Conceptually, we have two diﬀerent

cases:

DIV

∇f

DIV

g(x, y, t)∇f

(8)

Here, physical interpretation of the divergence is, in a sense, that of diﬀusion

[2]. In the case of DIV (∇f), diﬀusivity is the same in each direction, whereas

in the case of DIV (g(x, y, t) ∇f), diﬀusivity is deﬁned (or controlled) by the

function g and is not necessarily the same in all the directions. Mathematically,

for a diﬀerentiable vector function F = U

i + V

j, divergence operator is deﬁned

as:

DIV

F

=

∂U

∂x

+

∂V

∂y

(9)

In other words, divergence is a sum of partial derivatives of a diﬀerentiable

vector function. Therefore, in our case, we have the following.

DIV

∇f

=

∂

∂x

f

x

+

∂

∂y

f

y

=

∂

2

f

∂x

2

+

∂

2

f

∂y

2

= ∆f

DIV

g(x, y, t)∇f

=

∂

∂x

g(x, y, t)f

x

+

∂

∂y

g(x, y, t)f

y

= ∇g ∇f + g∆f

(10)

Now, by simply using the ﬁnite diﬀerences introduced above, one way of

discretising the divergence terms in (10) is using the central diﬀerence. First we

apply the central diﬀerence and then the forward- and the backward diﬀerences

for approximating the corresponding derivatives. The ‘trick’ here is to realise

that (f

x

)(x + 0.5, y) is actually the forward diﬀerence D

+

x

f(x), while (f

x

)(x −

0.5, y) is the backward diﬀerence D

−

x

f(x). Equations (11), and (12) show the

discretisations for DIV (∇f), and DIV (g(x, y, t)∇f), respectively. This is the

same discretisation as in the famous paper by Perona and Malik [2].

5

∂

∂x

f

x

(x, y) +

∂

∂y

f

y

(x, y) =

f

x

(x + 0.5, y) −

f

x

(x − 0.5, y)

+

f

y

(x, y + 0.5) −

f

y

(x, y − 0.5)

=f(x + 1, y) − f(x, y) + f(x − 1, y) − f(x, y)

+ f(x, y + 1) − f(x, y) + f(x, y − 1) − f (x, y)

=∇

E

f + ∇

W

f + ∇

S

f + ∇

N

f

(11)

where ∇

{W,N,E,S}

f denotes the diﬀerence in the directions given by W, N, E, S.

As it was already mentioned, ﬁrst we apply ﬁrst order central diﬀerence on

f

x

(x, y), and thus obtain D

0

x

f

x

(x, y) =

f

x

(x + 0.5, y) −

f

y

(x − 0.5, y).

∂

∂x

gf

x

(x, y) +

∂

∂y

gf

y

(x, y) =

gf

x

(x + 0.5, y) −

gf

x

(x − 0.5, y)

+

gf

y

(x, y + 0.5) −

gf

y

(x, y − 0.5)

=g(x + 0.5, y)

f(x + 1, y) − f(x, y)

+ g(x − 0.5, y)

f(x − 1, y) − f(x, y)

+ g(x, y + 0.5)

f(x, y + 1) − f(x, y)

+ g(x, y − 0.5)

f(x, y − 1) − f(x, y)

=g

E

∇

E

f + g

W

∇

W

f + g

S

∇

S

f + g

N

∇

N

f

(12)

where g

{W,N,E,S}

denotes diﬀusivity in the directions given by W, N, E, S. As

it can be observed from Equation (12), we need to approximate the diﬀusivity

between the pixels. A simple ‘2-point’ approximation would be the average be-

tween neighbouring pixels, for example g(x + 0.5, y) = [g(x + 1, y) + g(x, y)]/2.

A more precise approximation, leading to better results, is a ‘6-point’ approxi-

mation of Brox [1].

As it was already mentioned above, physical interpretation of the divergence

is, in a sense, diﬀusion: the divergence operator introduces a ‘connectivity’

between the pixels. This simply means, as will be shown later on, that a solu-

tion at any position (i, j) will depend on the solution at neighbouring positions.

Because of this kind of a dependency of the solution between the adjacent posi-

tions, variational correspondence methods are said to be ‘global’. This kind of

connectivity is problematic at image borders, where we do not have neighbours

anymore. In order to deal with this problem, we use a scheme called eliminated

boundary conditions, shown in Figure 4.

6

❡❥❤ ❤

❤

❤

(a)

❡❥

❤

❤

(b)

Figure 4: Double circle denotes the position of interest while single circles are

the neighbouring positions W, N, E, S; (b) shows the eliminated boundary

conditions.

3.3 Discretised Diﬀusion Equation

In order to solve Equation 3, we need to discretise the divergence operator. We

start by marking the positions of the pixels of interest with (i, j):

(I

k

)

t+1

i,j

− (I

k

)

t

i,j

τ

= DIV

g

t

∇I

t+1

k

(13)

After this we discretise the DIV operator, as given by Equation 12, and

obtain the following:

(I

k

)

t+1

i,j

− (I

k

)

t

i,j

=τg

t

N

(I

k

)

t+1

i−1,j

− (I

k

)

t+1

i,j

+ τg

t

S

(I

k

)

t+1

i+1,j

− (I

k

)

t+1

i,j

+ τg

t

W

(I

k

)

t+1

i,j−1

− (I

k

)

t+1

i,j

+ τg

t

E

(I

k

)

t+1

i,j+1

− (I

k

)

t+1

i,j

(14)

The only thing left to do, is arrange the terms:

(I

k

)

t+1

i,j

1 + τ (g

t

N

+ g

t

S

+ g

t

W

+ g

t

E

)

=(I

k

)

t

i,j

+ τg

t

N

(I

k

)

t+1

i−1,j

+ τg

t

S

(I

k

)

t+1

i+1,j

+ τg

t

W

(I

k

)

t+1

i,j−1

+ τg

t

E

(I

k

)

t+1

i,j+1

(15)

3.4 Matrix Format

While Equation (15) shows equation to be solved for a pixel position (i, j),

we can write the system equations in matrix/vector format, covering the whole

7

image, as given in (17). Now, the components of the vector are I

I

where I ∈

{1, . . . , N} and N is the number of pixels in image. We use I, J also to mark

the positions in the system matrix A given in (17). This is done in order to

convey clearly the idea that the domains of the discretised images and the system

matrices are diﬀerent. If the domain of the discretised image is Ω

h

: [1, m]×[1, n]

(discrete image with m columns and n rows), the system matrix A is deﬁned

on [m · n] × [m · n] (here · denotes multiplication). Now, we can write the Euler

forward, semi-implicit formulation in a vector/matrix format as follows:

(I

k

)

t+1

− (I

k

)

t

τ

= A

(I

k

)

t

I

t+1

k

(16)

where I

k

:= (I

k

)

I

with I = [1 . . . N ], and k refers to the channel in question

(e.g. R, G or B). The system matrix A

(I

k

)

t

is deﬁned as follows:

A

(I

k

)

t

= [a

t

I,J

]

a

t

I,J

:=

−τg

t

J ∼I

[J ∈ N(I)],

1 +

J ∈N

−

(I)

J ∈N

+

(I)

τg

t

J ∼I

(J = I),

0 (otherwise)

(17)

where g

t

J ∼I

refers to the ‘diﬀusion’ weight between pixels J and I at time t.

In other words, these refer to the g

{W, N, E, S}

seen previously. Equation (18)

gives an example of how the system matrix A would look like for a 3 × 3 size

image. Here C and N are block matrices that refer to the ‘central’ and the

‘neighbouring’ matrices, correspondingly.

A =

C N 0

N C N

0 N C

C =

1 +

J ∈N

−

(I)

J ∈N

+

(I)

τg

t

J ∼I

−τg

t

J ∼I

0

−τg

t

J ∼I

1 +

J ∈N

−

(I)

J ∈N

+

(I)

τg

t

J ∼I

−τg

t

J ∼I

0 −τg

t

J ∼I

1 +

J ∈N

−

(I)

J ∈N

+

(I)

τg

t

J ∼I

N =

−τg

t

J ∼I

0 0

0 −τg

t

J ∼I

0

0 0 −τg

t

J ∼I

(18)

From (17) we can see how the matrix A looks like for a 3 × 3 size image: it is

a block tridiagonal square matrix, of size 9 × 9, that has non-zero components

only on the main diagonal and on the diagonals adjacent to this. Therefore,

unless the image is very small, it is infeasible to solve the system by inverting

A directly. Instead, we search for a solution using iterative methods.

8

In the following we can see how A would look like using diﬀerent ‘notation’.

It is interesting to see how changing the traversing order changes the neighbours

on the diagonals next to the main diagonal. Sub matrices with grey background

refer to the sub matrix C given in Equation 18

C S 0 E 0 0 0 0 0

N C S 0 E 0 0 0 0

0 N C 0 0 E 0 0 0

W 0 0 C S 0 E 0 0

0 W 0 N C S 0 E 0

0 0 W 0 X C 0 0 E

0 0 0 W 0 0 C S 0

0 0 0 0 W 0 N C S

0 0 0 0 0 W 0 N C

(a) Column-wise traversing.

C E 0 S 0 0 0 0 0

W C E 0 S 0 0 0 0

0 W C 0 0 S 0 0 0

N 0 0 C E 0 S 0 0

0 N 0 W C E 0 S 0

0 0 N 0 X C 0 0 S

0 0 0 N 0 0 C E 0

0 0 0 0 N 0 W C E

0 0 0 0 0 N 0 W C

(b) Row-wise traversing.

Figure 5: Table like representation of the system matrix A for both column-

and row-wise ordering. Here W, N, E, S and C refer to the neighbours given

by the cardinal directions (West, North, East and South), while C refers to the

‘Centre’.

4 AOS

With the vector/matrix format in place, we can now formulate the ‘additive

operator splitting’ scheme by Weickert et al. [5]. In order to simplify the

notation, we write A instead of A

(I

k

)

t

. Id refers to the identity matrix.

Therefore, we have:

(I

k

)

t+1

= (I

k

)

t

+ τAI

t+1

k

(19)

From which (I

k

)

t+1

can be solved as follows:

(I

k

)

t+1

=

Id − τ A

−1

I

t

k

(20)

Now, we ‘decompose’ A so that A =

m

l=1

A

l

, which allows as to write the above

equation as:

9

(I

k

)

t+1

=

m

l=1

1

m

Id − τ

m

l=1

A

l

−1

I

t

k

(21)

where m is the number of dimensions (in our case m = 2). Previous equation

can be written, using only a single summation operator, as:

(I

k

)

t+1

=

m

l=1

1

m

Id − τ mA

l

−1

I

t

k

(22)

The idea of writing A as A =

m

l=1

A

l

is based on the traversing order of A as

shown in Fig. 5. When constructing each A

l

(here l can be though of referring

to the traversing order based on the dimension), we only take into account the

neighbours on the diagonals next to the main diagonal and discard the rest.

Physically this can be interpreted as diﬀusing separately along the vertical and

horizontal directions. When constructing the matrices A

l

, one must be careful

with the ‘central’ elements, Equation (18), so that only the neighbours on the

diagonals next to the main diagonal are taken into account. This is apparent

from Equation (30).

Equation (22) has interesting ‘form’ in the sense that the ‘system matrix’

is decomposed. The problem is that the decomposed system matrix is inside

the

−1

operator. Instead, we would like to construct the solution in parts as

follows:

(I

k

)

t+1

=

m

l=1

1

m

Id − τ mA

l

−1

I

t

k

(23)

The problem is that the right hand sides of Equations (22) and (23) are not

equal, as can be easily veriﬁed. Therefore, we pose the question if there exists

a simple variable x, when used to multiply the right hand side of (23), would

make these equal:

m

l=1

1

m

Id − τ mA

l

B

−1

I

t

k

= x

m

l=1

1

m

Id − τ mA

l

B

−1

I

t

k

(24)

The above can be simpliﬁed into:

B

−1

= xm

2

B

−1

(25)

And, thus we have:

x =

1

m

2

(26)

Based on this, in order to use the ‘additive operator splitting’ scheme given

by Equation (23), we multiply the right hand side with

1

m

2

, and obtain the

following:

10

(I

k

)

t+1

=

1

m

2

m

l=1

1

m

Id − τ mA

l

−1

I

t

k

(27)

which is the same as:

(I

k

)

t+1

=

m

l=1

mId − τ m

2

A

l

−1

I

t

k

(28)

As an example, if l = 2 (2D), then we would have:

(I

k

)

t+1

=

2Id − 4τ A

1

−1

I

t

k

+

2Id − 4τ A

2

−1

I

t

k

(29)

Now, as it can be understood, the whole idea of this scheme is to bring the

equations to a ‘simpler’ form, allowing us to use eﬃcient block-wise solvers,

such as TDMA (TriDiagonal Matrix Algorithm).

4.1 TDMA

TDMA. Tridiagonal matrix algorithm (TDMA) is a simpliﬁed form of Gaussian

elimination, that can be used for solving tridiagonal systems of equations. In

matrix/vector format this kind of a system can be written as in (30)

b

1

c

1

0 0 0

a

2

b

2

c

2

0 0

0 a

3

b

3

.

.

.

0

0 0

.

.

.

.

.

.

c

N−1

0 0 0 a

N

b

N

A

x

1

x

2

x

3

.

.

.

x

N

x

=

d

1

d

2

d

3

.

.

.

d

N

d

(30)

The algorithm consists of two steps: the ﬁrst (forward) sweep eliminates the

a

i

, while the second (backward) sweep calculates the solution. Equation (31)

introduces the forward sweep, while Equation (32) shows the backward sweep.

c

i

=

c

1

b

1

, i = 1

c

i

b

i

− c

i−1

a

i

, i = 2, 3, ..., N − 1

d

i

=

d

1

b

1

, i = 1

d

i

− d

i−1

a

i

b

i

− c

i−1

a

i

, i = 2, 3, ..., N

(31)

x

N

= d

N

x

i

= d

i

− c

i

x

i+1

, i = N − 1, N − 2, ..., 1

(32)

Physical interpretation of the terms a

i

and b

i

is that they are diﬀusion

weights, i.e. how much the neighbouring solutions are taken into account.

11

Figures 6 and 7 show both the horizontal- and vertical traversing orders, and

how the coeﬃcients a

i

and c

i

are chosen.

1 2 3 4 5 6 7 8 9

❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄

(a)

1

2

3

4

5

6

7

8

9

✲

✲

✲

✲

✲

✲

✲

✲

✲

(b)

Figure 6: Column- and row-wise pixel ordering. In AOS we apply, for example,

TDMA ﬁrst along all the columns and then along all the rows. (a) Column wise

ordering; (b) row wise ordering.

❥❡

❤

❤

a

i

c

i

❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄

(a)

❥❡ ❤❤

a

i

c

i

✲

✲

✲

✲

✲

✲

✲

✲

✲

(b)

Figure 7: Pixel positions a

i

and c

i

of b

i

, depending on the pixel ordering: (a)

column wise ordering; (b) row wise ordering.

4.2 Matlab Implementation

Following is a Matlab implementation, based on the Wikipedia Matlab code

1

,

of the TDMA algorithm. Since each of the ‘columns’ (column-wise traversing)

and ‘rows’ (row-wise traversing) are independent of each other, these can be

processed is parallel, if needed.

1

http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm

12

Algorithm 1 Tridiagonal matrix algorithm. a, b, c are the column vectors for

the compressed tridiagonal matrix, d is the right vector.

function x = TDMA(a, b, c, d)

[rows cols dim] = size(a);

%Modify the ﬁrst-row coeﬃcients

c(1, :) = c(1, :) ./ b(1, :);

d(1, :) = d(1, :) ./ b(1, :);

%Forward pass

for i = 2:rows-1

temp = b(i, :) - a(i, :) .* c(i-1, :);

c(i, :) = c(i, :) ./ temp;

d(i, :) = (d(i, :) - a(i, :) .* d(i-1, :)) ./ temp;

end

%Backward pass

x(rows, :) = d(rows, :);

for i = rows-1:-1:1

x(i, :) = d(i, :) - c(i, :) .* x(i + 1, :);

end

end function

5 Acknowledgements

I would like to thank everyone who helped me reviewing my PhD thesis (on

which this text is based on). Especially, I would like to thank Dr. Florian

Pokorny and Dr. Karl Pauwels for their help with the mathematical notation

(which can become quite clumsy) and Pablo Guzm´an for reviewing this partic-

ular version of the text.

13

References

[1] T. Brox. From Pixels to Regions: Partial Diﬀerential Equations in Image

Analysis. PhD thesis, Saarland University, Saarbr¨ucken, Germany, 2005.

[2] P. Perona and J. Malik. Scale-space and edge detection using anisotropic

diﬀusion. IEEE Transactions on Pattern Analysis and Machine Intelligence,

12:629–639, 1990.

[3] J. Ralli. Fusion and Regularisation of Image Information in Variational

Correspondence Methods. PhD thesis, University of Granada, Spain, 2011.

[4] U. Trottenberg, C. Oosterlee, and A. Sch¨uller. Multigrid. Academic Press,

2001.

[5] J. Weickert, B. M. Ter Haar Romeny, and M. A. Viergever. Eﬃcient and re-

liable schemes for nonlinear diﬀusion ﬁltering. IEEE Transactions on Image

Processing, 7:398–410, 1998.

14