Differential Equations in Image Processing and Computer Vision

Assignment T2

Jan Hendrik Dithmar

2031259

Group 2: Thu, 12-14 (Markus Mainberger)

Assignment T2  Dithmar (2031259)

2.1 Isotropic Nonlinear Diffusion: Diffusivities

a.

g(s^{2}) = \frac{1}{1 + \frac{s^{2}}{λ^{2}}}
Φ(s) = s · g(s^{2})
= s · \frac{1}{1 + \frac{s^{2}}{λ^{2}}} = \frac{s}{1 + \frac{s^{2}}{λ^{2}}}
Φ′(s) = \frac{1 · (1 + \frac{s^{2}}{λ^{2}} ) − s · \frac{1}{λ^{2}} · 2s}{(1 + \frac{s^{2}}{λ^{2}} )^{2}}
= \frac{1}{1 + \frac{s^{2}}{λ^{2}}} − \frac{2s^{2}}{(1 + \frac{s^{2}}{λ^{2}} )^{2} · λ^{2}}
= \frac{(1 + \frac{s^{2}}{λ^{2}} ) · λ^{2} − 2s^{2}}{(1 + \frac{s^{2}}{λ^{2}} )^{2} · λ^{2}}
Φ′(s) = \frac{λ^{2} − s^{2}}{(1 + \frac{s^{2}}{λ^{2}} )^{2} · λ^{2}}
⇒ Φ′(s) > 0 ⇔ |s| < λ (forward diffusion)
Φ′(s) < 0 ⇔ |s| > λ (backward diffusion)

b.

g(s^{2}) = \frac{1}{\sqrt{1 + \frac{s^{2}}{λ^{2}}}}

Φ(s) = \frac{s}{\sqrt{1 + \frac{s^{2}}{λ^{2}}}}

Φ′(s) = \frac{(1 + \frac{s^{2}}{λ^{2}} )^{\frac{1}{2}} − s · \frac{1}{2} · (1 + \frac{s^{2}}{λ^{2}} ) · \frac{1}{λ^{2}} · 2s}{1 + \frac{s^{2}}{λ^{2}}}

= \frac{(1 + \frac{s^{2}}{λ^{2}} )^{\frac{1}{2}} − \frac{s^{2}}{λ^{2}} · (1 + \frac{s^{2}}{λ^{2}} )}{1 + \frac{s^{2}}{λ^{2}}}
= \frac{1}{(1 + \frac{s^{2}}{λ^{2}} )^{\frac{1}{2}} − \frac{s^{2}}{λ^{2}}}

Φ′(s) = \frac{λ^{2} − s^{2} · (1 + \frac{s^{2}}{λ^{2}} )^{\frac{1}{2}}}{λ^{2} · (1 + \frac{s^{2}}{λ^{2}} )^{\frac{1}{2}}}

⇒ Φ′(s) > 0 ⇔ |s| < \frac{λ}{(1 + \frac{s^{2}}{λ^{2}} )^{\frac{1}{4}}} (forward diffusion)

Φ′(s) < 0 ⇔ |s| > \frac{λ}{(1 + \frac{s^{2}}{λ^{2}} )^{\frac{1}{4}}} (backward diffusion)

c.

g(s^{2}) = exp (\frac{−s^{2}}{2λ^{2}} )
Φ(s) = s · exp (\frac{−s^{2}}{2λ^{2}} )
Φ′(s) = exp (−\frac{s^{2}}{2λ^{2}} ) + s · exp (−\frac{s^{2}}{2λ^{2}} ) · (−\frac{2s}{2λ^{2}} )
= exp (−\frac{s^{2}}{2λ^{2}} ) − \frac{s^{2}}{λ^{2}} · exp (−\frac{s^{2}}{2λ^{2}} )
Φ′(s) = (1 − \frac{s^{2}}{λ^{2}} ) · \underbrace{exp (−\frac{s^{2}}{2λ^{2}} )}_{>0}

⇒ Φ′(s) > 0 ⇔ |s| < λ (forward diffusion)
Φ′(s) < 0 ⇔ |s| > λ (backward diffusion)

2.2 Isotropic Nonlinear Diffusion: Stencil Notation of Explicit Scheme

∂_{t}u = div (g(|∇u|^{2})∇u)
= ∂_{x} (g(|∇u|^{2})∂_{x}u) + ∂_{y} (g(|∇u|^{2})∂_{y}u)

We know from the lecture (Lecture 5, slide 3) that:

\frac{du_{i,j}}{dt} = \frac{1}{h_{1}} (\frac{g_{i+1,j} + g_{i,j}}{2} \frac{u_{i+1,j} − u_{i,j}}{h_{1}} − \frac{g_{i,j} + g_{i−1,j}}{2} \frac{u_{i,j} − u_{i−1,j}}{h_{1}} )

+\frac{1}{h_{2}} (\frac{g_{i,j+1} + g_{i,j}}{2} \frac{u_{i,j+1} − u_{i,j}}{h_{2}} − \frac{g_{i,j} + g_{i,j−1}}{2} \frac{u_{i,j} − u_{i,j−1}}{h_{2}} )

Furthermore we know that

\frac{du}{dt} = A(u)u
⇔ \frac{u_{i,j}^{k+1} − u_{i,j}^{k}}{τ} = \frac{1}{2h_{1}^{2}} ((g_{i+1,j} + g_{i,j})(u_{i+1,j} − u_{i,j}) − (g_{i,j} + g_{i−1,j})(u_{i,j} − u_{i−1,j}))

+\frac{1}{2h_{2}^{2}} ((g_{i,j+1} + g_{i,j})(u_{i,j+1} − u_{i,j}) − (g_{i,j} + g_{i,j−1})(u_{i,j} − u_{i,j−1}))

⇔ u_{i,j}^{k+1} = \frac{τ}{2h_{1}^{2}} (g_{i+1,j}u_{i+1,j} − g_{i+1,j}u_{i,j} + g_{i,j}u_{i+1,j} − g_{i,j}u_{i,j})

−\frac{τ}{2h_{1}^{2}} (g_{i,j}u_{i,j} − g_{i,j}u_{i−1,j} + g_{i−1,j}u_{i,j} − g_{i−1,j}u_{i−1,j})

+\frac{τ}{2h_{2}^{2}} (g_{i,j+1}u_{i,j+1} − g_{i,j+1}u_{i,j} + g_{i,j}u_{i,j+1} − g_{i,j}u_{i,j})

−\frac{τ}{2h_{2}^{2}} (g_{i,j}u_{i,j} − g_{i,j}u_{i,j−1} + g_{i,j−1}u_{i,j} − g_{i,j−1}u_{i,j−1})

+u_{i,j}

= u_{i,j} (1 − τ \frac{g_{i+1,j} + 2g_{i,j} + g_{i−1,j}}{2h_{1}^{2}} − τ \frac{g_{i,j+1} + 2g_{i,j} + g_{i,j−1}}{2h_{2}^{2}} )

+u_{i+1,j} (τ \frac{g_{i+1,j} + g_{i,j}}{2h_{1}^{2}} ) + u_{i−1,j} (τ \frac{g_{i,j} + g_{i−1,j}}{2h_{1}^{2}} )

+u_{i,j+1} (τ \frac{g_{i,j+1} + g_{i,j}}{2h_{2}^{2}} ) + u_{i,j−1} (τ \frac{g_{i,j} + g_{i,j−1}}{2h_{2}^{2}} )

With this knowledge, we can give the stencil notation:

\begin{matrix}0 & τ \frac{g_{i,j+1} + g_{i,j}}{2h_{2}^{2}} & 0 \\ τ \frac{g_{i,j} + g_{i−1,j}}{2h_{1}^{2}} & 1 − τ ( \frac{g_{i+1,j} + 2g_{i,j} + g_{i−1,j}}{2h_{1}^{2}} + \frac{g_{i,j+1} + 2g_{i,j} + g_{i,j−1}}{2h_{2}^{2}} ) & τ \frac{g_{i+1,j} + g_{i,j}}{2h_{1}^{2}} \\ 0 & τ \frac{g_{i,j} + g_{i,j−1}}{2h_{2}^{2}} & 0 \end{matrix}

2.3 Isotropic Nonlinear Diffusion: Matrix Notation of

Semi-Implicit Scheme

For a 1-D signal of length 4 we get a 4 × 4 matrix A which has the following entries:

a_{11} = − (\frac{g_{1} + g_{0}}{2h^{2}} + \frac{g_{1} + g_{2}}{2h^{2}} ) = −\frac{g_{0} + 2g_{1} + g_{2}}{2h^{2}}
a_{22} = − (\frac{g_{2} + g_{1}}{2h^{2}} + \frac{g_{2} + g_{3}}{2h^{2}} ) = −\frac{g_{1} + 2g_{2} + g_{3}}{2h^{2}}
a_{33} = − (\frac{g_{3} + g_{2}}{2h^{2}} + \frac{g_{3} + g_{4}}{2h^{2}} ) = −\frac{g_{2} + 2g_{3} + g_{4}}{2h^{2}}
a_{44} = − (\frac{g_{4} + g_{3}}{2h^{2}} + \frac{g_{4} + g_{5}}{2h^{2}} ) = −\frac{g_{3} + 2g_{4} + g_{5}}{2h^{2}}
a_{12} = \frac{g_{1} + g_{2}}{2h^{2}}
a_{13} = 0
a_{14} = 0
a_{21} = \frac{g_{1} + g_{2}}{2h^{2}}
a_{23} = \frac{g_{2} + g_{3}}{2h^{2}}
a_{24} = 0
a_{31} = 0
a_{32} = \frac{g_{2} + g_{3}}{2h^{2}}
a_{34} = \frac{g_{3} + g_{4}}{2h^{2}}
a_{41} = 0
a_{42} = 0
a_{43} = \frac{g_{3} + g_{4}}{2h^{2}}

In the case of a 2-D image of size 3 × 4, we get a 12 × 12 matrix A with the following entries (where all the entries not stated here are 0):

a_{1,1} = − (\frac{g_{1} + g_{1}}{2h_{1}^{2}} + \frac{g_{1} + g_{2}}{2h_{1}^{2}} + \frac{g_{1} + g_{1}}{2h_{2}^{2}} + \frac{g_{1} + g_{4}}{2h_{2}^{2}} )
= −\frac{3g_{1} + g_{2}}{2h_{1}^{2}} − \frac{3g_{1} + g_{4}}{2h_{2}^{2}}
a_{2,2} = − (\frac{g_{2} + g_{1}}{2h_{1}^{2}} + \frac{g_{2} + g_{3}}{2h_{1}^{2}} + \frac{g_{2} + g_{2}}{2h_{2}^{2}} + \frac{g_{2} + g_{5}}{2h_{2}^{2}} )
= −\frac{g_{1} + 2g_{2} + g_{3}}{2h_{1}^{2}} − \frac{3g_{2} + g_{5}}{2h_{2}^{2}}
a_{3,3} = −\frac{g_{2} + 3g_{3}}{2h_{1}^{2}} − \frac{3g_{3} + g_{6}}{2h_{2}^{2}}
a_{4,4} = −\frac{3g_{4} + g_{5}}{2h_{1}^{2}} − \frac{g_{1} + 2g_{4} + g_{7}}{2h_{2}^{2}}
a_{5,5} = −\frac{g_{4} + 2g_{5} + g_{6}}{2h_{1}^{2}} − \frac{g_{2} + 2g_{5} + g_{8}}{2h_{2}^{2}}
a_{6,6} = −\frac{g_{5} + 3g_{6}}{2h_{1}^{2}} − \frac{g_{3} + 2g_{6} + g_{9}}{2h_{2}^{2}}
a_{7,7} = −\frac{3g_{7} + g_{8}}{2h_{1}^{2}} − \frac{g_{4} + 2g_{7} + g_{10}}{2h_{2}^{2}}
a_{8,8} = −\frac{g_{7} + 2g_{8} + g_{9}}{2h_{1}^{2}} − \frac{g_{5} + 2g_{8} + g_{11}}{2h_{2}^{2}}
a_{9,9} = −\frac{g_{8} + 3g_{9}}{2h_{1}^{2}} − \frac{g_{6} + 2g_{9} + g_{12}}{2h_{2}^{2}}
a_{10,10} = −\frac{3g_{10} + g_{11}}{2h_{1}^{2}} − \frac{g_{7} + 3g_{10}}{2h_{2}^{2}}
a_{11,11} = −\frac{g_{10} + 2g_{11} + g_{12}}{2h_{1}^{2}} − \frac{g_{8} + 3g_{11}}{2h_{2}^{2}}
a_{12,12} = −\frac{g_{11} + 3g_{12}}{2h_{1}^{2}} − \frac{g_{9} + 3g_{12}}{2h_{2}^{2}}

a_{1,2} = \frac{g_{1} + g_{2}}{2h_{1}^{2}}  a_{1,4} = \frac{g_{1} + g_{4}}{2h_{2}^{2}}  a_{2,1} = \frac{g_{1} + g_{2}}{2h_{1}^{2}}^{a_{2,3} = \frac{g_{2} + g_{3}}{2h_{1}^{2}}  a_{2,5} = \frac{g_{2} + g_{5}}{2h_{2}^{2}}  a_{3,2} = \frac{g_{2} + g_{3}}{2h_{1}^{2}}}^{a_{3,6} = \frac{g_{3} + g_{6}}{2h_{2}^{2}}  a_{4,1} = \frac{g_{1} + g_{4}}{2h_{2}^{2}}  a_{4,5} = \frac{g_{4} + g_{5}}{2h_{1}^{2}}}^{a_{4,7} = \frac{g_{4} + g_{7}}{2h_{2}^{2}}  a_{5,2} = \frac{g_{2} + g_{5}}{2h_{2}^{2}}  a_{5,4} = \frac{g_{4} + g_{5}}{2h_{1}^{2}}}^{a_{5,6} = \frac{g_{5} + g_{6}}{2h_{1}^{2}}  a_{5,8} = \frac{g_{5} + g_{8}}{2h_{2}^{2}}  a_{6,3} = \frac{g_{3} + g_{6}}{2h_{2}^{2}}}^{a_{6,5} = \frac{g_{5} + g_{6}}{2h_{1}^{2}}  a_{6,9} = \frac{g_{6} + g_{9}}{2h_{2}^{2}}  a_{7,4} = \frac{g_{4} + g_{7}}{2h_{2}^{2}}}^{a_{7,8} = \frac{g_{7} + g_{8}}{2h_{1}^{2}}  a_{7,10} = \frac{g_{7} + g_{10}}{2h_{2}^{2}}  a_{8,5} = \frac{g_{5} + g_{8}}{2h_{2}^{2}}}^{a_{8,7} = \frac{g_{7} + g_{8}}{2h_{1}^{2}}  a_{8,9} = \frac{g_{8} + g_{9}}{2h_{1}^{2}}  a_{8,11} = \frac{g_{8} + g_{11}}{2h_{2}^{2}}}^{a_{9,6} = \frac{g_{6} + g_{9}}{2h_{2}^{2}}  a_{9,8} = \frac{g_{8} + g_{9}}{2h_{1}^{2}}  a_{9,12} = \frac{g_{9} + g_{12}}{2h_{2}^{2}}}^{a_{10,7} = \frac{g_{7} + g_{10}}{2h_{2}^{2}}  g_{10,11} = \frac{g_{10} + g_{11}}{2h_{1}^{2}}  g_{11,8} = \frac{g_{8} + g_{11}}{2h_{2}^{2}}}^{g_{11,10} = \frac{g_{10} + g_{11}}{2h_{1}^{2}}  g_{11,12} = \frac{g_{11} + g_{12}}{2h_{1}^{2}}  g_{12,9} = \frac{g_{9} + g_{12}}{2h_{2}^{2}}}^{g_{12,11} = \frac{g_{11} + g_{12}}{2h_{1}^{2}}}

2.4 Linear Systems of Equations and Thomas’ Algorithm

a. We want to decompose the matrix

A = ( \begin{matrix}γ_{1} & δ_{1} & ε_{1} & 0 & 0 \\ β_{1} & γ_{2} & δ_{2} & ε_{2} & 0 \\ α_{1} & β_{2} & γ_{3} & δ_{3} & ε_{3} \\ 0 & α_{2} & β_{3} & γ_{4} & δ_{4} \\ 0 & 0 & α_{3} & β_{4} & γ_{5} \end{matrix} )

into

L := ( \begin{matrix}1 & 0 & 0 & 0 & 0 \\ b_{1} & 1 & 0 & 0 & 0 \\ a_{1} & b_{2} & 1 & 0 & 0 \\ 0 & a_{2} & b_{3} & 1 & 0 \\ 0 & 0 & a_{3} & b_{4} & 1 \end{matrix} )

and

R := ( \begin{matrix}c_{1} & d_{1} & e_{1} & 0 & 0 \\ 0 & c_{2} & d_{2} & e_{2} & 0 \\ 0 & 0 & c_{3} & d_{3} & e_{3} \\ 0 & 0 & 0 & c_{4} & d_{4} \\ 0 & 0 & 0 & 0 & c_{5} \end{matrix} ) .

The matrix A can be written as follows:

A = LR

= ( \begin{matrix}c_{1} & d_{1} & e_{1} & 0 & 0 \\ b_{1}c_{1} & b_{1}d_{1} + c_{2} & b_{1}e_{1} + d_{2} & e_{2} & 0 \\ a_{1}c_{1} & a_{1}d_{1} + b_{2}c_{2} & a_{1}e_{1} + b_{2}d_{2} + c_{3} & b_{2}e_{2} + d_{3} & e_{3} \\ 0 & a_{2}c_{2} & a_{2}d_{2} + b_{3}c_{3} & a_{2}e_{2} + b_{3}d_{3} + c_{4} & b_{3}e_{3} + d_{4} \\ 0 & 0 & a_{3}c_{3} & a_{3}d_{3} + b_{4}c_{4} & a_{3}e_{3} + b_{4}d_{4} + c_{5} \end{matrix} )

LR-Decomposition: We know that e_{i} = ε_{i}, i ∈ {1, 2, 3}. The coefficients c_{i}, b_{i} and a_{i} can be computed by:

c_{1} := γ_{1}
b_{1} := β_{1} / c_{1}
a_{1} := α_{1} / c_{1}
d_{1} := δ_{1}

for i = 2, ..., 5:
    if i ≤ 4:

d_{i} := δ_{i} - b_{i−1}e_{i−1}

c_{i} := γ_{i} - b_{i−1}d_{i−1} - a_{i−2}e_{i−2}

    if i ≤ 4:

b_{i} := β_{i} / c_{i} - a_{i−1}d_{i−1}

    if i ≤ 3:

      a_{i} := α_{i} / c_{i}

b. Forward Elimination: We solve Ly = u for y:

( \begin{matrix}1 & 0 & 0 & 0 & 0 \\ b_{1} & 1 & 0 & 0 & 0 \\ a_{1} & b_{2} & 1 & 0 & 0 \\ 0 & a_{2} & b_{3} & 1 & 0 \\ 0 & 0 & a_{3} & b_{4} & 1 \end{matrix} ) ( \begin{matrix}y_{1} \\ y_{2} \\ y_{3} \\ y_{4} \\ y_{5} \end{matrix} ) = ( \begin{matrix}u_{1} \\ u_{2} \\ u_{3} \\ u_{4} \\ u_{5} \end{matrix} )

y_{1} := u_{1}
y_{2} := u_{2} - b_{1}y_{1}
for i = 3, .., 5:

y_{i} := u_{i} - a_{i−2}y_{i−2} - b_{i−1}y_{i−1}

Backward Substitution: We solve Rx = y for x:

( \begin{matrix}c_{1} & d_{1} & e_{1} & 0 & 0 \\ 0 & c_{2} & d_{2} & e_{2} & 0 \\ 0 & 0 & c_{3} & d_{3} & e_{3} \\ 0 & 0 & 0 & c_{4} & d_{4} \\ 0 & 0 & 0 & 0 & c_{5} \end{matrix} ) ( \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \\ x_{5} \end{matrix} ) = ( \begin{matrix}y_{1} \\ y_{2} \\ y_{3} \\ y_{4} \\ y_{5} \end{matrix} )

x_{5} := y_{5} / c_{5}
x_{4} := (y_{4} - d_{4}x_{5}) / c_{4}
for i = 3, ..., 1:
    x_{i} := (y_{i} - d_{i}x_{i+1} - e_{i}x_{i+2}) / c_{i}