PDE Based Image Diffusion 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 diffusion equation. The paper
shows how a continuous version of the diffusion 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 definitions, the idea here being that
anyone familiar with diffusion 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 defines the number
of channels (e.g. K = 3 in the case of RGB images). It is in this regular domain
where our PDEs are defined. However, in order for us to solve the PDEs, they
need to be discretised first. 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 define 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 defined 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 defined 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.fi
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, first 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
define different 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 different 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 defined by J ∈ N (I). (b) Painted circles denote pixels J that
belong to the neighbourhood defined by J ∈ N
−
(I), while unpainted circles
denote pixels J that belong to the neighbourhood defined by J ∈ N
+
(I). (c)-
(d) Painted circles denote pixels J that belong to the neighbourhood defined by
J ∈ N
−
(I), while unpainted circles denote pixels J that belong to the neigh-
bourhood defined by J ∈ N
+
(I). Processing order is defined 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) defined 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 defined
by this operator depends on the direction of the solver.
3
2 Scalar Valued Diffusion
The basic formula describing a scalar valued diffusion is:
∂I
k
∂t
= DIV
g(x, y, t)∇I
k
(2)
where k refers to the channel in question, g(x, y, t) defines the (scalar) diffusion
weights, and ∇ := [
∂
∂x
,
∂
∂y
] is the spatial gradient. Since g(x, y, t) is a function
of t, the diffusion 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 Difference Discretisation
3.1 Finite Difference Operators
Before going any further, we remind the reader of the ‘standard’ finite difference
operators, which are used for approximating derivatives. Here, the function of
interest, for which we want to find derivatives, is defined as f(x, y) : Ω → R.
1. First order forward difference 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 difference 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 difference 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 difference 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 differentiated. The above formulations can be simplified
further if we assume a uniform grid with ∆x = ∆y = 1. The difference operators
shown here approximate derivatives with respect to x. By switching the roles of
4
x and y, corresponding difference 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 first and second order derivatives, we can
discretise the divergence, DIV , operator. Conceptually, we have two different
cases:
DIV
∇f
DIV
g(x, y, t)∇f
(8)
Here, physical interpretation of the divergence is, in a sense, that of diffusion
[2]. In the case of DIV (∇f), diffusivity is the same in each direction, whereas
in the case of DIV (g(x, y, t) ∇f), diffusivity is defined (or controlled) by the
function g and is not necessarily the same in all the directions. Mathematically,
for a differentiable vector function F = U
i + V
j, divergence operator is defined
as:
DIV
F
=
∂U
∂x
+
∂V
∂y
(9)
In other words, divergence is a sum of partial derivatives of a differentiable
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 finite differences introduced above, one way of
discretising the divergence terms in (10) is using the central difference. First we
apply the central difference and then the forward- and the backward differences
for approximating the corresponding derivatives. The ‘trick’ here is to realise
that (f
x
)(x + 0.5, y) is actually the forward difference D
+
x
f(x), while (f
x
)(x −
0.5, y) is the backward difference 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 difference in the directions given by W, N, E, S.
As it was already mentioned, first we apply first order central difference 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 diffusivity in the directions given by W, N, E, S. As
it can be observed from Equation (12), we need to approximate the diffusivity
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, diffusion: 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 Diffusion 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 different. 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 defined
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 defined 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 ‘diffusion’ 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 different ‘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 diffusing 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 verified. 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 simplified 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 efficient block-wise solvers,
such as TDMA (TriDiagonal Matrix Algorithm).
4.1 TDMA
TDMA. Tridiagonal matrix algorithm (TDMA) is a simplified 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 first (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 diffusion
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 coefficients 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 first 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 first-row coefficients
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 Differential 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
diffusion. 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. Efficient and re-
liable schemes for nonlinear diffusion filtering. IEEE Transactions on Image
Processing, 7:398–410, 1998.
14