\pagestyle{fancy}
\cleardoublepage\chapter{Diffusion\index{diffusion}}\label{diffusion}
\Quote{All science is either physics or stamp collecting.}{Lord Kelvin}

\noindent Diffusion\index{diffusion} processes are motivated from physics. Diffusion\index{diffusion} equilibrates concentration differences by redistributing mass while preserving the total mass. This property is expressed by \emph{Fick's law}\index{Fick's law}
\begin{equation}\label{Ficks_Law}
	j = -D \cdot \nabla u.
\end{equation}
Equation \eqref{Ficks_Law} means that the concentration gradient $\nabla u$ causes a flux $j$ which aims to compensate the gradient. The \emph{diffusion tensor}\index{diffusion!tensor} $D$ -- usually a positive definite symmetric matrix -- describes the relation between $\nabla u$ and $j$. We differentiate between two cases: \emph{isotropic}\index{diffusion!isotropic} and \emph{anisotropic}\index{diffusion!anisotropic} \emph{diffusion}. If $j$ and $\nabla u$ are parallel, the diffusion\index{diffusion} is isotropic\index{diffusion!isotropic}\index{isotropic diffusion|see {diffusion}}. In the general anisotropic\index{diffusion!anisotropic}\index{anisotropic diffusion|see {diffusion}} case, $j$ and $\nabla u$ are not parallel.

Since we know that a diffusion\index{diffusion} process preserves the total mass, we can express this fact in the \emph{continuity equation}\index{continuity equation}
\begin{equation}
	\partial_t u = - \text{div}\; j,
\end{equation}
where $t$ denotes time.

Plugging Fick's law\index{Fick's law} into the continuity equation\index{continuity equation} we get
\begin{equation}\label{diffusion_equation}
	\partial_t u = \text{div}(D \cdot \nabla u).
\end{equation}
Equation \eqref{diffusion_equation} is called \emph{diffusion equation}\index{diffusion!equation} and in the case of heat transfer, it is called \emph{heat equation}\index{heat equation}.

A detailed investigation of diffusion processes can be found in \cite{weickert98}.

\clearpage
\section{The Structure Tensor\index{structure tensor}}
The \emph{structure tensor}\index{structure tensor} was introduced by F{\"o}rstner and G{\"u}lch \cite{foerstner87} and is a well-known tool for analyzing local orientation. Part of this analysis can be features like corners or edges.

Let us consider a matrix $J$ which results from the tensor product
\begin{equation}
	J(\nabla u_\sigma) := \nabla u_\sigma \nabla u_\sigma^\top
\end{equation}
where $\nabla u_\sigma$ is the gradient of $u$ presmoothed by a Gaussian kernel of size $\sigma$. $\nabla u_\sigma$ is a vector-valued structure descriptor. The matrix $J(\nabla u_\sigma)$ has an orthonormal basis of eigenvectors $v_1$ and $v_2$ with $v_1 \Vert \nabla u_\sigma$ and $v_2 \bot \nabla u_\sigma$. For a grayscale image, $J(\nabla u_\sigma)$ is a matrix of rank 1 which means that the eigenvalue in the direction of the eigenvector $v_2$ is 0.

Convolving $J(\nabla u_\sigma)$ component-wise with a Gaussian kernel $K_\rho$, we obtain the structure tensor\index{structure tensor}
\begin{equation}
	J_\rho(\nabla u_\sigma) := K_\rho * (\nabla u_\sigma \nabla u_\sigma^\top) = \left(\begin{array}{cc}K_\rho * (u_x^2) & K_\rho * (u_x u_y) \\ K_\rho * (u_x u_y) & K_\rho * (u_y^2)\end{array}\right) \quad \text{with } \rho \geq 0.
\end{equation}
The structure tensor\index{structure tensor} $J_\rho$ is symmetric and positive definite. Its orthonormal eigenvectors $v_1$ and $v_2$ specify the preferred structure directions within some \emph{integration scale} $\rho$. The corresponding eigenvalues $\mu_1$ and $\mu_2$ describe the average contrast along the eigendirections and allow for a useful analysis of the local image structure. Thus, the integration scale $\rho$ should reflect the characteristic window size over which the orientation is to be analyzed. Presmoothing in order to obtain $\nabla u_\sigma$ makes the structure tensor insensitive to noise. The parameter $\sigma$ is called \emph{noise scale}.

Let us assume that $\mu_1 \geq \mu_2 \geq 0$. In this case we observe that $v_1$ is the orientation with the highest gray level fluctuations and $v_2$ points us to the preferred local orientation. Furthermore, the eigenvalues tell us useful information for the analysis of the local structure: constant areas are characterized by $\mu_1 = \mu_2 = 0$, straight edges by $\mu_1 \gg \mu_2 = 0$, corners by $\mu_1 \geq \mu_2 \gg 0$ and the expression
\begin{equation*}
	(\mu_1 - \mu_2)^2 = (j_{11} - j_{22})^2 + 4 j_{12}^2
\end{equation*}
becomes large for anisotropic structures, i.e.\ $(\mu_1 - \mu_2)^2$ is a measure of anisotropy.

\clearpage
\section{Edge-Enhancing Diffusion\index{diffusion!edge-enhancing}\index{EED|see{diffusion, edge-enhancing}}}\label{EED}
In the setting of \emph{edge-enhancing diffusion} (EED)\index{diffusion!edge-enhancing} \cite{weickert98, eed-weickert96} we want to prevent the diffusion process from diffusing across edges. Let us assume $\mu_1$ and $\mu_2$ with $\mu_1 \geq \mu_2$ be the eigenvalues of the structure tensor $J_\rho$ and $v_1$ and $v_2$ are the corresponding orthonormal eigenvectors.

Since the diffusion tensor\index{diffusion!tensor} should reflect the local image structure, it should be chosen such that it has the same set of eigenvectors $v_1$ and $v_2$ as the structure tensor\index{structure tensor} $J_\rho$.

To achieve that we do not diffuse across edges, we reduce the diffusivity $\lambda_1$ perpendicular more and more as $\mu_1$ grows. This can be realized by choosing the eigenvalues of the diffusion tensor such that
\begin{eqnarray}\label{eed_eigenvalues}
	\lambda_1(\mu_1) & := & g(\mu_1), \\
	\lambda_2 & := & 1
\end{eqnarray}
with
\begin{equation}\label{eed_g-function}
	g(s) := \left\{\begin{array}{ll}1 & \text{if } s \leq 0 \\ \frac{1}{1 + s / \lambda^2} & \text{if } s > 0\end{array}\right. ,
\end{equation}
where $\mu_1 = |\nabla u_\sigma|^2$ and $\lambda > 0$.

The application of edge-enhancing diffusion\index{diffuion!edge-enhancing} on a fingerprint image is depicted in Figure \ref{diffusion.eed_fingerprint}.

\begin{figure}[htb]
	\centering
	\includegraphics[scale=0.42]{figures/finger.eps}
	\includegraphics[scale=0.42]{figures/finger_eed.eps}
	\caption[Image processed with EED.]{Example for an image processed with edge-enhancing diffusion\index{diffusion!edge-enahncing}. \textbf{Left:} Noisy fingerprint image. \textbf{Right:} Same image processed using edge-enhancing diffusion\index{diffusion!edge-enhancing}: $\lambda = 1$, $\sigma = 0.5$, $\tau = 0.25$, $t = 20$.}
	\label{diffusion.eed_fingerprint}
\end{figure}

\clearpage
\section{Coherence-Enhancing Diffusion\index{diffusion!coherence-enhancing}\index{CED|see{diffusion, coherence-enhancing}}}\label{CED}
In the context of \emph{coherence-enhancing diffusion} (CED)\index{diffusion!coherence-enhancing} \cite{weickert98, ced-weickert99} we want to use the information provided by the structure tensor\index{structure tensor} to design an anisotropic filter that enhances \emph{flow-like structures}. Again, the diffusion tensor should reflect the local image structure.% Thus it has the same eigenvectors $v_1$ and $v_2$ as $J_\rho$.

Let again $\mu_1$ and $\mu_2$ with $\mu_1 \geq \mu_2$ be the eigenvalues of the structure tensor $J_\rho$ and $v_1$ and $v_2$ be the corresponding eigenvectors. If we want to enhance coherent structures, we have to smooth preferably along $v_2$ with a diffusivity of $\lambda_2$. This can be achieved by choosing the eigenvalues of the diffusion tensor\index{diffusion!tensor} in the following way:
\begin{eqnarray}
\label{ced_eigenvalues1}	\lambda_1 & := & \alpha, \\
\label{ced_eigenvalues2}	\lambda_2 & := & \left\{\begin{array}{ll}\alpha & \text{if } \mu_1 = \mu_2 \\ \alpha + (1 - \alpha) \exp\left(\frac{-C}{(\mu_1 - \mu_2)^2}\right) & \text{else}\end{array}\right. ,
\end{eqnarray}
where $C > 0$ acts as a threshold parameter and $\alpha \in (0,1)$ is a small positive parameter. The exponential function and the parameter $\alpha$ are introduced for theoretical reasons: The exponential function ensures the smoothness of the diffusion tensor\index{diffusion!tensor} and the parameter $\alpha$ keeps the diffusion tensor\index{diffusion!tensor} uniformly positive definite. In other words, $\alpha$ acts as a regularization parameter.

The enhancement of a fingerprint image using coherence-enhancing diffusion\index{diffusion!coherence-enhacing} is shown in Figure \ref{diffusion.ced_fingerprint}.

\begin{figure}[htb]
	\centering
	\includegraphics[scale=0.42]{figures/finger.eps}
	\includegraphics[scale=0.42]{figures/finger_ced.eps}
	\caption[Image processed with CED.]{Example for an image processed with coherence-enhancing diffusion\index{diffusion!coherence-enahncing}. \textbf{Left:} Noisy fingerprint image. \textbf{Right:} Same image processed using coherence-enhancing diffusion\index{diffusion!coherence-enhancing}: $C = 1$, $\sigma = 0.5$, $\rho = 4$, $\alpha = 0.1$, $\tau = 0.25$, $t = 20$.}
	\label{diffusion.ced_fingerprint}
\end{figure}

\clearpage
\section{Numerical Method}
We briefly mention how to realize a diffusion process algorithmically. The first step is to rewrite the continuity equation \eqref{diffusion_equation}, resulting in
\begin{equation}
	\partial_t u = \partial_x \left(D \partial_x u\right) + \partial_y \left(D \partial_y u\right).
\end{equation}
We use a \emph{forward difference} for the time derivative $\partial_t u$
\begin{equation}
	\partial_t u \approx \frac{u_{i,j}^{k+1} - u_{i,j}^k}{\tau}.
\end{equation}
%and \emph{central differences} for the spatial derivatives $\partial_x u$ and $\partial_y u$
To compute the structure tensor in each pixel $(i,j)$, we discretize the spatial derivatives $\partial_x u$ and $\partial_y u$ by \emph{central differences}
\begin{eqnarray}
	\label{central_difference1}(\partial_x u)_{i,j}^k & \approx & \frac{u_{i+1,j}^k - u_{i-1,j}^k}{2 h_x}, \\
	\label{central_difference2}(\partial_y u)_{i,j}^k & \approx & \frac{u_{i,j+1}^k - u_{i,j-1}^k}{2 h_y},
\end{eqnarray}
where $h_x$ denotes the spatial grid size in $x$-direction and $h_y$ the spatial grid size in $y$-direction (both usually set to 1). The diffusion tensor
\begin{equation}
	D_{i,j}^k = \left(\begin{array}{cc}a_{i,j}^k & c_{i,j}^k \\[0.3cm] c_{i,j}^k & b_{i,j}^k\end{array}\right)
\end{equation}
is computed using the structure tensor, compare Sections \ref{EED} and \ref{CED}.

Furthermore, we use central differences twice with half the step size for $\partial_x (a u_x)$ and $\partial_y (b u_y)$.
\begin{eqnarray}
	\Big(\partial_x (a u_x)\Big)_{i,j}^k & \approx & \frac{(a u_x)_{i+\frac12,j}^k-(a u_x)_{i-\frac12,j}^k}{h_x} \\
\nonumber	& \approx & \frac{1}{h_x}\left(\frac{a_{i+1,j}^k + a_{i,j}^k}{2} \cdot \frac{u_{i+1,j}^k- u_{i,j}^k}{h_x} \right. \\
\nonumber	&&	\quad \quad \left. - \frac{a_{i,j}^k + a_{i-1,j}^k}{2} \cdot \frac{u_{i,j}^k - u_{i-1,j}^k}{h_x} \right),
\end{eqnarray}
where $(\partial_y (b u_y))_{i,j}^k$ can be computed analogously. Furthermore, we compute the mixed derivatives by
\begin{equation}
	u_{xy} = (u_{\xi\xi} - u_{\eta\eta}) \cdot \frac{h_x^2 + h_y^2}{4h_xh_y}
\end{equation}
using the unit vectors
\begin{eqnarray}
	\xi & = & \frac{1}{\sqrt{h_x^2 + h_y^2}} \left(\begin{array}{c}h_x \\ h_y\end{array}\right), \\
	\eta & = & \frac{1}{\sqrt{h_x^2 + h_y^2}} \left(\begin{array}{r}h_x \\ -h_y\end{array}\right).
\end{eqnarray}
The discretization of $u_{\xi\xi}$ and $u_{\eta\eta}$ works analogously to the discretization above where the step size is $\sqrt{h_x^2 + h_y^2}$ along the diagonal.

The resulting \emph{explicit scheme} reads
\begin{eqnarray}
	u_{i,j}^{k+1} & = & u_{i,j}^k \\
\nonumber		& + &	\tau \left( \frac{(a_{i+1,j}^k+a_{i,j}^k)(u_{i+1,j}^k-u_{i,j}^k)}{2h_x^2} -\frac{(a_{i,j}^k+a_{i-1,j}^k)(u_{i,j}^k-u_{i-1,j}^k)}{2h_x^2} \right.\\
\nonumber			&&	\quad +\frac{(b_{i,j+1}^k+b_{i,j}^k)(u_{i,j+1}^k-u_{i,j}^k)}{2h_y^2} -\frac{(b_{i,j}^k+b_{i,j-1}^k)(u_{i,j}^k-u_{i,j-1}^k)}{2h_y^2} \\
\nonumber			&&	\quad +\frac{(c_{i+1,j+1}^k+c_{i,j}^k)(u_{i+1,j+1}^k-u_{i,j}^k)}{4h_xh_y} -\frac{(c_{i,j}^k+c_{i-1,j-1}^k)(u_{i,j}^k-u_{i-1,j-1}^k)}{4h_xh_y} \\
\nonumber			&&	\quad -\frac{(c_{i+1,j-1}^k+c_{i,j}^k)(u_{i+1,j-1}^k-u_{i,j}^k)}{4h_xh_y} \left. +\frac{(c_{i,j}^k+c_{i-1,j+1}^k)(u_{i,j}^k-u_{i-1,j+1}^k)}{4h_xh_y} \right).
\end{eqnarray}
This scheme is called explicit, since $u_{i,j}^{k+1}$ can be computed explicitly from the known values at the time step $k$.
