\section{Diffusion filtering}
-First we say something about normal Gaussian filtering, how it is
-related to solving the heat equation etc. Weickert
-\cite{weickert1998anisotropic} is a good reference here, and in many of
-the other sections. End with the motivation of controlling the diffusion
-coefficient.
+Diffusion filtering is a broad group of filtering or restoration methods
+which take the original image as initial value of some diffusion
+process. The most well-known method being the Gaussian filter or
+Gaussian blur, which convolves the image with the Gaussian function
+\begin{equation}
+ G_\sigma(x,y) = \frac{1}{2\pi\sigma^2} \exp \left( - \frac{x^2 +
+ y^2}{2\sigma^2} \right).
+ \label{eq:gaussian_function}
+\end{equation}
+In the discrete setting where the image consists of a grid of pixels,
+the Gaussian blur amounts to setting each pixel in the output image to a
+weighted average of its neighboring pixels in the original image.
+
+The Gaussian function happens to be the fundamental solution of the heat
+equation $\partial_t u = \Delta u$. Convolving it with the original
+image $v$ is therefore equivalent to solving the heat equation with $v$
+as initial value, until some time $T > 0$ depending on $\sigma$. Care
+must be taken on the boundary, and one possibility, which can be argued
+to be sensible in some imaging contexts, is to symmetrically extend the
+image in all directions.
+
+By basic Fourier analysis it is possible to show that the Gaussian
+filter is a low-pass filter which attenuates high frequencies.
+\fixme{ref weickert?}
+
+This method will, in addition to smoothing out possible noise, remove
+details from the image, which motivates the next set of methods, where
+the amount of diffusion can vary for different parts of the image.
\subsection{Non-linear diffusion}
-This section is for methods that control the diffusion coefficient
-$\alpha$. Here I guess Perona-Malik is the most noteworthy. We write
-something about the motivation, its well-posedness (or lack thereof) and
-other alternatives. But we can control the diffusion even more!
+\fixme{Some note on the different names used in different places. We
+stick to Weickert.}
-Perona--Malik
+In the theory of the heat equation there is a thermal diffusivity
+$\alpha$ that can be introduced such that
\begin{equation}
- \partial_t u = \mathrm{div} \left( \frac{\nabla u}{1 +
+ \begin{cases}
+ \partial_t u &= \diver \big( \alpha(u) \nabla u\big) \\
+ u |_{t=0} &= v
+ \end{cases}
+\end{equation}
+The $\alpha$ is dependant on the material, and can also vary throughout
+the object. We can make use of this in the image restoration context by
+controlling the diffusivity in different parts of the image. Optimally,
+we would like there to be a lot of diffusion in smooth parts of the
+image, and not so much in areas with a lot of details.
+
+One much-studied non-linear diffusion equation is the Perona--Malik
+equation
+\begin{equation}
+ \partial_t u = \diver \left( \frac{\nabla u}{1 +
\frac{\abs{\nabla u}^2}{\lambda^2}} \right).
\label{eq:perona_malik}
\end{equation}
-Total variation flow
+which has a thermal diffusivity that varies from $1$ in smooth areas to
+$\infty$ as the norm of the gradient $\abs{\nabla u}$ grows. The special
+form of the thermal diffusivity has been shown to be related to how
+brightness is percepted by the human visual system \fixme{ref}. The
+model has some theoretical problems related to well-posedness, for more
+information see \cite{weickert1998anisotropic}.
+
+A different kind of non-linear diffusion model is the total variation
+flow which can be formulated as
\begin{equation}
- \partial_t u = \mathrm{div} \frac{\nabla u}{\abs{\nabla u}}.
+ \partial_t u = \diver \frac{\nabla u}{\abs{\nabla u}},
\end{equation}
+where the diffusivity has a similar effect of reducing the diffusion in
+areas of high variation. As the name suggests this model can be related
+to the \fixme{variational} total variation formulation presented later.
+One discrete time-step in the solution of this PDE corresponds to the
+Euler--Lagrange equation of the variational formulation.
+
+Common for all of these are that we control the amount of diffusion done
+in each point of the image, based somehow on the image itself. But we
+can go further by controlling not only the amount, but also the
+direction of the diffusion, such that we might have more diffusion in
+directions with little variation. This is the so-called anisotropic
+diffusion, meaning not uniform in all directions.
\subsection{Anisotropic diffusion}
-Section for methods introducing a diffusion \emph{tensor} $A$ that
+\fixme{Section for methods introducing a diffusion \emph{tensor} $A$ that
varies around the image. This means we can control diffusion strength in
different directions. More specifically we can make the diffusion strong
along edges, and weak across edges. We should give a real theoretical
-background here, so we have something to base the next chapter on.
+background here, so we have something to base the next chapter on.}
-We must also say something about what this method does well, and where
-it fails. Anisotropic filtering might introduce structures that do not
-exist in the original image. Some example images are possible here.
+\fixme{Where are riemannian metrics introduced? Here? Probably later,
+since the tensor does not really relate to a metric here.}
-We consider the initial boundary value problem
+The diffusion is made directionally dependant by introducing a diffusion
+\emph{tensor} $A(u)$ that varies throughout the image such that the
+initial boundary value problem becomes
\begin{align}
\begin{cases}
- \partial_t u &= \mathrm{div} \big(A(u) \nabla u\big)
+ \partial_t u &= \diver \big(A(u) \nabla u\big)
\text{ on } \Omega \times (0, \infty),\\
u(x, 0) &= f(x)
\text{ on } \Omega,\\
A(u) \nabla u \cdot \nu &= 0
- \text{ on } \Gamma \times (0, \infty).
+ \text{ on } \partial \Omega \times (0, \infty).
\end{cases}
\end{align}
-where $\nu$ is the outer normal of $\Omega$ and $A(u)$ is our image
-dependent diffusion tensor. Say something about $\Gamma$, if it hasn't
-been described earlier.
\fixme{Fix the alignment.}
-
-We define the smoothed image
+where $\nu$ is the outer normal of $\Omega$. The tensor $A(u)$ would
+normally diminish the effect of $\nabla u$ across what we believe to be
+edges in the image. Weickert suggests constructing $A(u)$ based on the
+edge estimator $\nabla u_\sigma$ where
\begin{equation}
- u_\sigma(x,t) := \big(K_\sigma * \tilde{u}(\cdot,t)\big)(x)
+ u_\sigma := K_\sigma * \tilde{u}
\end{equation}
-(after we introduce the Gaussian kernel $K_\sigma$ of course), and we
-also have to say something about $\tilde{u}$ being the extension of $u$
-in $\mathbb{R}^2$.
-
-We then convolve again, componentwise, to obtain the structure tensor
+and $\tilde{u}$ is an extension of $u$ from $\Omega$ to $\mathbb{R}^2$.
+This can be done by considering the structure tensor
\begin{equation}
- J_\rho(\nabla u_\sigma) := K_\rho * \left( \nabla u_\sigma \otimes
- \nabla u_\sigma \right).
+ S(x) := K_\rho * (\nabla u_\sigma \otimes \nabla u_\sigma)
\end{equation}
-We then have to say something about what the $\sigma$ and $\rho$
-actually control, and how the eigenvectors of this tensor behaves. Maybe
-say something about existence and uniqueness as well. Weickert does not
-say much about how to proceed from the structure tensor, at least not in
-Chapter 2, but we should definitely say something about that.
+and its eigenvectors and eigenvalues. This will be further discussed
+later. Assuming some smoothness, symmetry and uniform positive
+definiteness on $A(u)$ one can prove well-posedness, regularity and a
+extremum principle as done in \cite{weickert1998anisotropic}.
+
+\fixme{Main source here is Weickert.}
-Main source here is Weickert.
+\fixme{
+ We must also say something about what this method does well, and
+ where it fails. Anisotropic filtering might introduce structures
+ that do not exist in the original image. Some example images are
+ possible here.
+}
\section{Total variation filtering}
-Describe briefly the functional we minimize in this method, and
+\fixme{
+ Describe briefly the functional we minimize in this method, and
intuitively how this leads to a filtered image. But how much do we write
about how we find our minimizer? We have to at least mention that there
are different ways to do it, and that one of them involves finding a
minimum cut for each level.
+}
+\fixme{
We could introduce the normal coarea and perimeter theory here so that
we can refer to it when we introduce the anisotropic versions in the
next chapter?
+}
+\fixme{
We also write something about the limitations of the method. It is good
at preserving sharp edges, but suffers from the so-called stair-casing
effect. Thin objects and corners might gradually lose contrast depending
on the strength of the regularization. Example images here as well.
+}
+\fixme{
Then finish the chapter in some way saying that we will try to combine
these two methods, retaining their strenghts but eliminating their
weaknesses.
+}
-\fixme{This is taken from the project.}
-
-First we introduce this.
+Although related to the total variation flow diffusion equation in
+\fixme{ref}, the total variation image restoration method is usually
+formulated as a minimization problem
\begin{equation}
\min_u \int_\Omega \abs{u - v}^p + \beta \int_\Omega \abs{\nabla u}.
\label{eq:first_min_presentation}
\end{equation}
-Then we say something about difficulties calculating the gradient and
-introduce the total variation.
+The first term penalizes images $u$ that are far from the original
+image, while the second term is the total variation, and minimizing it
+will smooth out and regularize the image. The $\beta$ parameter controls
+the strength of the regularization.
+
+Since we do not only want to consider images $u \in C^1(\Omega)$ for
+which the gradient exists, we introduce the total variation using the
+distributional derivative.
\begin{definition}[Total variation]
Given a function $u \in L^1(\Omega)$, the total variation of $u$,
often written $\int_\Omega \abs{Du}$, where the $D$ is the gradient
Note that since $\Omega$ is open and bounded, the test functions
$\varphi$ vanish on the boundary of $\Omega$ \fixme{which means we
-assume some boundary condition}. Next we introduce the space of
-functions with bounded variation.
+assume some boundary condition}. As this restoration method is the one
+which is extended in this thesis, we will look a little bit more deeply
+into the background and the numerical methods relating to it.
+
+As we are trying to minimize the total variation of the image, it is
+useful to introduce the space of functions of bounded variation.
\begin{definition}[Functions of bounded variation]
The space of functions of bounded variation $\BV(\Omega)$ is the
space of functions $u \in L^1(\Omega)$ for which the total variation
\right\}.
\end{equation}
\end{definition}
+Our optimization problem has thus become
+\begin{equation}
+ \min_{u \in \BV(\Omega)} \int_\Omega \abs{u - v}^p + \beta \, \TV(u)
+\end{equation}
+\fixme{say something about what $p$ can be. is it problematic if $p=2$?}
+
+\subsection{Chambolle's dual approach}
+
+Maybe, maybe not.
+
+\subsection{Graph cut approach}
-Before introducing the coarea formula, we need the formal definition of
-a set perimeter. Here, $\idfun_E$ signifies the characteristic function
-of the set $E$, the function which is equal to one in every point in
-$E$, and zero elsewhere.
+Using graph cuts is also the approach we will be taking later when
+considering the anisotropic total variation minimization, and it is
+therefore valuable to briefly look into how graph cuts are used in the
+case of regular total variation.
+
+The idea of the graph cut approach is to decompose the minimization
+problem into one minimization problem for each \emph{level} of the
+image, before solving them separately and combining the results.
+
+For an image $u$ and a level $\lambda$ we denote the \emph{level set} by
+$\{ u > \lambda\}$ defined as the set $\{ x \in \Omega : u_x > \lambda
+\}$. The thresholded image $u^\lambda$ is then defined as
+\begin{equation}
+ u^\lambda = \idfun_{u > \lambda}.
+\end{equation}
+Here, $\idfun_E$ signifies the characteristic function of the set $E$,
+the function which is equal to one in every point in $E$, and zero
+elsewhere. The graph cut we find later will for each level $\lambda$
+give us the thresholded image $u^\lambda$ which can then be combined to
+form the complete image $u$.
+\fixme{Decide upon which way this inequality should go, as it varies in
+the literature.}
+
+Before introducing the coarea formula, which will allow us to write the
+total variation as an integral over all the level values, we need the
+formal definition of a set perimeter.
\begin{definition}[Set perimeter]
The perimeter of a set $E$ is defined as
\begin{equation}
- P(E,\Omega) = \int_\Omega \abs{D \idfun_E}.
+ \Per(E,\Omega) = \int_\Omega \abs{D \idfun_E}.
\end{equation}
A measurable set $E \subset \Omega$ is of finite perimeter in
$\Omega$ if $\idfun_E \in \BV(\Omega)$.
$\Omega$, then $P(E,\Omega)$ will not include the overlapping parts of
the boundary.
-For an image $u$ and a level $\lambda$ we denote the \emph{level set} by
-$\{u \leq \lambda\}$ defined as the set $\{ x \in \Omega : u_x \leq
-\lambda\}$. The thresholded images are defined as
-\begin{equation}
- u^\lambda = \idfun_{u \leq \lambda}.
- \label{eq:level_set}
-\end{equation}
-\fixme{Decide upon which way this inequality should go, as it varies in
-the literature.}
With all these definitions in place, we are ready to introduce the
coarea formula, which will let us write the total variation as an
-integral over the range of the image levels $\lambda$.
+integral over the range of the image levels $\lambda$. \fixme{repeated}
\begin{theorem}[The coarea formula]
Let $u \in BV(\Omega)$. Then for almost every $\lambda$ the set
- $\{u \leq \lambda\}$ is of finite perimeter, and one has the coarea
+ $\{u > \lambda\}$ is of finite perimeter, and one has the coarea
formula
\begin{equation}
\TV(u)
= \int_\Omega \abs{D u}
- = \int_\mathbb{R} P(\{u \leq \lambda\},\Omega) \, d\lambda
+ = \int_\mathbb{R} \Per(\{u > \lambda\},\Omega) \, d\lambda
\label{eq:coarea_formula}
\end{equation}
\end{theorem}
\begin{equation}
E_v(u) =
\int_\Omega \abs{u - v}^p
- + \beta \int_\mathbb{R} P(\{u \leq \lambda\},\Omega) \, d\lambda.
+ + \beta \int_\mathbb{R} \Per(\{u > \lambda\},\Omega) \, d\lambda.
\label{eq:energy_perimeter}
\end{equation}
\fixme{Should this be a minimization problem or an energy function.}
but we need to fill in some text at least), and add some general
references.}
-Now write a little bit about how one might solve this? Just as a small
-introduction for the things to come.
-
-\subsection{Primal--dual approach}
-
-Maybe, maybe not.
-
-\subsection{Graph cut approach}
+The energy function is now ready to be discretized. The image is now
+defined on a grid $\mathcal{G}$ and each pixel $x \in \mathcal{G}$ can
+take a value in the set of levels $\mathcal{L} = \{0, \hdots, L-1\}$.
+This is a reasonable assumption for digital grayscale images. The energy
+functional in \fixme{ref} is then discretized as
+\begin{equation}
+ E_v(u) =
+ \sum_{\lambda=0}^{L-2} \sum_x E^x_\lambda(u^\lambda_x)
+ + \beta \sum_{\lambda = 0}^{L-2} \sum_{(x, y)} E^{x,y}(u^\lambda_x,
+ u^\lambda_y)
+ =: \sum_{\lambda=0}^{L-2} F_\lambda(u^\lambda)
+ \label{eq:total_energy}
+\end{equation}
+where
+\begin{align}
+ E^x_\lambda(u^\lambda_x) &=
+ \big(
+ N_x(\lambda + 1) -
+ N_x(\lambda)
+ \big)
+ \, u^\lambda_x
+ \label{eq:fidelity_energy} \\
+ E^{x,y}(u^\lambda_x, u^\lambda_y) &=
+ w_{xy} \abs{u_x^\lambda - u_y^\lambda}.
+ \label{eq:neigh_energy}
+\end{align}
+and
+\begin{equation}
+ N_x(k) = \abs{k - v_x}^p.
+\end{equation}
-Just a little bit to show how it works maybe. Describe briefly the
-intuition of it all, how we split the problem into one problem on each
-level, and how these are solved using graph cuts. Describe what the
-graph cuts give us, and how it is combined back into an image.
+Here the sum over $(x, y)$ is over all pixel pairs $x$ and $y$ such
+that $y \in \mathcal{N}(x)$. \fixme{meh}. The intuition is that
+the perimeter of the set $\{ u > \lambda \}$ is related to the number of
+neighborhood relations crossing the perimeter, and if the weights
+$w_{xy}$ are chosen carefully, we even have concistency.
\section{Anisotropic total variation}
-Here we start out with the normal total variation (?), give some
-motivation and introduce the anisotropy. We probably need the definition
-of total variation, and the space of functions of bounded variation. The
-order is not clear to me, yet!
-
+From anisotropic diffusion in \fixme{ref} we borrow the idea of
+making the regularization in each point directionally dependant. We
+introduce the anisotropic total variation
\begin{equation}
- \int_\Omega \left( u - v \right)^2
- + \alpha \int_\Omega \nabla u(x)^T A(u) \nabla u(x)
+ \TVM(u) = \int_\Omega \sqrt{\nabla u(x)^T M(x) \nabla u(x)} \, dx
\end{equation}
\fixme{To differential or not to differential.}
-Describe how this is similar to a Euler time step in a related equation
-\cite{grasmair2010anisotropic}. But how do we calculate the gradients in
-this case? According to Grasmair the $u$ lives in $W^{1,1}(\Omega)$ (at
-least when the square root is involved), and the functional has to be
-extended to $\BV(\Omega)$. Maybe this $dDu$ thing from the other paper
-is important, have to understand what's happening there.
+for all $u \in C^1(\Omega)$. We see that if $M(x)$ is the identity
+matrix we get the normal total variation from \fixme{ref}, but our goal
+is to control $M(x)$ such that $\nabla u(x)$ is weighted down across
+edges, and thus edges will be regularized less.
+\fixme{Energy function here?}
+\fixme{Describe how this is similar to a Euler time step in a related
+equation \cite{grasmair2010anisotropic}. This is only without the square
+root.}
+\fixme{Extend the functional}
\section{Anisotropy tensor}
-We discuss what kind of properties we want from the anisotropy tensor,
-and also introduce the structure tensor (?) described by Weickert. This
-tensor can be visualized, and can be used for edge detection. There is
-probably a bit to be said about the smoothing we do. We also discuss
-whether to use the noisy image (yes, probably), or the smoothed image
-(implicitly, complicated, possible iteratively) in the structure tensor.
-We then describe how the structure tensor is decomposed and transformed
-into our anisotropy tensor. Then there are some problems related to
-stability which we can discuss here, or maybe we should leave it for a
-later, more implementation-focused chapter.
-
-Again, we define the structure tensor
+\fixme{
+ We discuss what kind of properties we want from the anisotropy
+ tensor, and also introduce the structure tensor (?) described by
+ Weickert. This tensor can be visualized, and can be used for edge
+ detection. There is probably a bit to be said about the smoothing we
+ do. We also discuss whether to use the noisy image (yes, probably),
+ or the smoothed image (implicitly, complicated, possible
+ iteratively) in the structure tensor. We then describe how the
+ structure tensor is decomposed and transformed into our anisotropy
+ tensor. Then there are some problems related to stability which we
+ can discuss here, or maybe we should leave it for a later, more
+ implementation-focused chapter.
+}
+
+There are many possible choices for the \fixme{metric} tensor $M(x)$, as
+long as one fulfills the assumptions of \fixme{ref}. Additionally we
+will later need the eigenvalues to be finite and bounded away from zero
+(?).
+
+Our wish for the metric tensor is that it down-weighs the variation
+across edges. We do not want the tensor to be sensitive to noise in the
+image. \fixme{Maybe refer to some figure}
+
+A good starting point is the structure tensor
\begin{equation}
- J_\rho(\nabla u_\sigma) := K_\rho * \left( \nabla u_\sigma \otimes
- \nabla u_\sigma \right).
+ S(x) := K_\rho * \left( \nabla u_\sigma(x) \otimes
+ \nabla u_\sigma(x) \right).
\end{equation}
This is then decomposed
\begin{equation}
\begin{theorem}[Anisotropic coarea formula]
Given an image $u \in \BV(\Omega)$, the anisotropic total variation
- can be calculated written as an integral over all the levels
+ can be written as an integral over all the levels
\begin{equation}
J_A(u) = \int_{-\infty}^\infty J_A(u^s) \, ds.
\end{equation}
to conclude that $m'(t)$ exists a.e.\ (w.r.t.\ Lebesgue measure) and
that Lebesgue's inequality holds:
\begin{equation}
- \int_{-infty}^\infty m'(t)\, dt \leq m(\infty) - m(-\infty) =
+ \int_{-\infty}^\infty m'(t)\, dt \leq m(\infty) - m(-\infty) =
J_A(u).
\end{equation}
Next, fix an $s \in \mathbb{R}$ and define the function
\fixme{something something set of measure zero, lebesgue integral
something something}
\end{proof}
+This is where we magically introduce the perimeter, or more specifically
+\begin{equation}
+ \Per(U;\Omega) = J_A(\idfun_U).
+\end{equation}
+This is just a definition, the magic lies in the fact that we can
+measure this perimeter as usual, just making sure to take the metric
+tensor into account
+\begin{equation}
+ \Per(U;\Omega) = \int_{\partial^*U \cap \Omega} \sqrt{\nu_U(x)^T
+ A(x) \nu_U(x)} \, d \mathcal{H}^{n-1}.
+\end{equation}
+And then we say a little bit about this reduced boundary. It is a bit
+wierd that we here consider the normal vector $\nu_U$ while when
+calculating the length of a curve in a Riemannian metric space we use
+the tangent vector
+\begin{equation}
+ L_a^b(c) = \int_a^b \norm{c'(t)}_A \, dt
+\end{equation}
+where $c: [a,b] \to \Omega$ is a parametrized curve, such that $c'(t)$
+is a tangent vector.
\section{Cauchy--Crofton formulas}
\emph{edges} going between the points in the lattice. We will denote an
edge by $e$. The curve length can then be approximated by
\begin{equation}
- \abs{C}_R \approx \sum_{e} n_C(e) \, \frac{\norm{e}^3 \det M(e)}{2
+ \abs{C}_R \approx \sum_{e} n_C(e) \, \frac{\det M(e) \norm{e}^3}{2
\left(e^T \cdot M(e) \cdot e\right)^{\sfrac{3}{2}}} \, \Delta\phi \,
\Delta\rho.
\end{equation}
+As our metric tensor $M$ is only known in the grid points, we
+approximate the tensor in the intersection point $p$ (between the curve
+$C$ and the edge $e$) by
+\begin{equation}
+ M(p) \approx M(e) = \frac{M\big(p_1(e)\big) + M\big(p_2(e)\big)}{2},
+\end{equation}
+the component-wise average of the tensors in the two end points of the
+edge.
\begin{figure}
\centering
\includegraphics[width=0.7\textwidth]{fig/area_argument.png}
family $k$ is $\Delta \rho_k = \delta^2 / \norm{e_k}$. Thus the curve
length is approximated by
\begin{equation}
- \abs{C}_R \approx \sum_{e} n_C(e) \, \frac{\norm{e}^2 \det M(e)
+ \abs{C}_R \approx \sum_{e\cap C} \frac{\det M(e) \norm{e}^2
\, \delta^2 \, \Delta\phi}{2 \left(e^T \cdot M(e) \cdot
e\right)^{\sfrac{3}{2}}}.
\end{equation}
\fixme{We now have the problem that $C$ is a curve, but later we want it
to be a cut...}
+\subsubsection{Consistency}
+Consider a square centered around grid point $p$ with side lengths
+$\sqrt{\delta}$. The size of this square will go to zero, while the
+number of points along the outer edge $\lfloor 1 / \delta \rfloor$ of
+the square, goes to infinity. Also, the number of points in the square
+(and thus the possible number of edges) goes to infinity so the
+discretization is finer and finer. So, as $\delta \to 0$, the length of
+the edges are bounded by $\sqrt{\delta} \to 0$, and $\Delta \phi \to 0$
+because the number of points around the edge of the square goes to
+infinity! \fixme{waving with hands} Also, since the curve is somewhat
+regular, when $\norm{e}$ goes to zero, the number of times the curve can
+cross $e$ is maximum $1$ so that's great. Also, $M(e) \to M(p)$
+obviously. \fixme{Some figure here aswell. And maybe some consideration
+of $\sup \Delta \phi$ etc.}
+
\section{Graph cut formulation}
Maybe this is more tightly connected with the previous section, but the
projects / master / commitdiff