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
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
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 + gf
(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
i1,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,j1
(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
i1,j
+ τg
t
S
(I
k
)
t+1
i+1,j
+ τg
t
W
(I
k
)
t+1
i,j1
+ τ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
N1
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
i1
a
i
, i = 2, 3, ..., N 1
d
i
=
d
1
b
1
, i = 1
d
i
d
i1
a
i
b
i
c
i1
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