\section{Diffusion filtering}
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
+based on physical diffusion processes. The basic idea is to take the
+noisy image as initial value of some diffusion process, and then let it
+evolve for some time. The most well-known method is probably the
+Gaussian filter or Gaussian blur, in which one convolves the image with
+the Gaussian function
\begin{equation}
- G_\sigma(x,y) = \frac{1}{2\pi\sigma^2} \exp \left( - \frac{x^2 +
+ K_\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 blur amounts to calculating each pixel in the output image
+as a weighted average of its neighboring pixels in the input 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.
+equation $\partial_t u = \Delta u$. Convolving $K_\sigma(x,y)$ 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$. Boundary conditions have to be specified of course, and one
+common choice is to symmetrically extend the image in $x$ and $y$
+directions, which corresponds to zero flux boundary conditions.
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?}
+filter is a low-pass filter which attenuates high frequencies. This and
+some further theory can be found in Weickert's book on anisotropic
+diffusion \cite{weickert1998anisotropic}.
-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.
+The main concern with the Gaussian filter is that it 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}
-\fixme{Some note on the different names used in different places. We
-stick to Weickert.}
-
-In the theory of the heat equation there is a thermal diffusivity
-$\alpha$ that can be introduced such that
+In the theory of the heat equation one can introduce a \emph{thermal
+diffusivity}
+$\alpha$ such that
\begin{equation}
\begin{cases}
\partial_t u &= \diver \big( \alpha(u) \nabla u\big) \\
- u |_{t=0} &= v
+ u |_{t=0} &= f
\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,
+The thermal diffusivity $\alpha$ is material dependant, and can also vary throughout
+the object. It specifies how well heat travels through the specific
+point in the object. We can make use of this in the image restoration context by
+controlling the diffusivity in different parts of the image, where our
+goal is to reduce noise without loosing image detail. 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.
\frac{\abs{\nabla u}^2}{\lambda^2}} \right).
\label{eq:perona_malik}
\end{equation}
-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}.
+The thermal diffusivity $\alpha(u) = (1 + \sfrac{\abs{\nabla u}^2}{
+\lambda^2})^{-1}$ varies from $1$ in smooth areas to $0$ as the
+norm of the gradient $\abs{\nabla u}$ grows. Although noise will also
+contribute to an increase in the size of $\abs{\nabla u}$, the general
+assumption is that this is a good measure of where the details are in
+the image. This particular form of the thermal diffusivity has been
+shown to be related to how brightness is percepted by the human visual
+system. The model has some theoretical problems related to
+well-posedness, for more information see \cite{weickert1998anisotropic}.
+
+Note that we follow Weickert's terminology when it comes to the
+distinction between non-linear and anisotropic methods. The
+Perona--Malik method, and other methods with varying scalar
+diffusivities, are often called anisotropic, but we will call them
+non-linear, and spare the anisotropic term for methods where the
+diffusivity is a tensor, and thus both location and direction dependant.
A different kind of non-linear diffusion model is the total variation
flow which can be formulated as
\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
+to the 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.
+can go further by not only making the diffusivity location dependant,
+but also direction dependant. This allows us to have have more
+diffusion in directions with little variation. As mentioned, we will
+use the term anisotropic diffusion, meaning ``not uniform in all
+directions.''
\subsection{Anisotropic diffusion}
+\label{sec:anisotropic_diffusion}
-\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.}
-
-\fixme{Where are riemannian metrics introduced? Here? Probably later,
-since the tensor does not really relate to a metric here.}
-
-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}
+The diffusivity is made directionally dependant by introducing a diffusion
+\emph{tensor} $A(u)$ such that the initial boundary value problem
+becomes
+\begin{equation}
\begin{cases}
\partial_t u &= \diver \big(A(u) \nabla u\big)
\text{ on } \Omega \times (0, \infty),\\
- u(x, 0) &= f(x)
+ u|_{t=0} &= f
\text{ on } \Omega,\\
A(u) \nabla u \cdot \nu &= 0
\text{ on } \partial \Omega \times (0, \infty).
\end{cases}
-\end{align}
-\fixme{Fix the alignment.}
+ \label{eq:aniso_diff}
+\end{equation}
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
+normally be constructed such as to diminish the effect of $\nabla
+u$ across what we believe to be edges in the image. This way, there will
+also be less diffusion through these edges. Weickert suggests in
+\cite{weickert1998anisotropic} constructing $A(u)$ based on the
edge estimator $\nabla u_\sigma$ where
\begin{equation}
u_\sigma := K_\sigma * \tilde{u}
\end{equation}
-and $\tilde{u}$ is an extension of $u$ from $\Omega$ to $\mathbb{R}^2$.
-This can be done by considering the structure tensor
+and $\tilde{u}$ is an extension of $u$ from $\Omega$ to $\mathbb{R}^2$
+made by for example symmetrically extending $u$ across the edges.
+Assuming we are at an edge in the image, the
+direction of $\nabla u_\sigma$ should be perpendicular to the edge,
+while its magnitude will provide information on the steepness of the
+edge.
+
+To extract this information, and also to identify features on a larger
+scale, we introduce the \emph{structure tensor}
\begin{equation}
- S(x) := K_\rho * (\nabla u_\sigma \otimes \nabla u_\sigma)
+ S_\rho(x) := K_\rho * (\nabla u_\sigma \otimes \nabla u_\sigma),
\end{equation}
-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.}
-
-\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}
-
-\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.
-}
-
-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
+where the convolution with the Gaussian function $K_\rho$ is done
+componentwise. The anisotropy tensor $A(u)$ can then be constructed
+based on the eigenvectors and eigenvalues of $S_\rho(x)$.
+
+The structure tensor and its properties will be discussed further when
+we introduce our anisotropic total variation functional. 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}.
+
+Even if the diffusivity tensor was introduced to make our method more
+conserving of edges, the solution of \eqref{eq:aniso_diff} will still be
+infinitely differentiable \cite{weickert1998anisotropic}, i.e.\ $u(T) \in
+C^\infty(\Omega)$ for $T > 0$. Thus we have no real discontinuities, and
+no real edges in our solution.
+
+Further, just like its isotropic version, the anisotropic diffusion
+may introduce structure based on noise, when there really was no
+structure to begin with. This is a problem we hope to avoid in our
+anisotropic total variation method.
+
+\begin{figure}
+ \centering
+ \begin{subfigure}[t]{0.4\textwidth}
+ \centering
+ \includegraphics[width=\textwidth]{fig/fnoise.png}
+ \caption{Noisy, $\sigma > 9000$}
+ \end{subfigure}
+ ~
+ \begin{subfigure}[t]{0.4\textwidth}
+ \centering
+ \includegraphics[width=\textwidth]{fig/fanisodiff.png}
+ \caption{Regularized, parameters: bleep bloop}
+ \end{subfigure}
+ \caption{\fixme{UPDATE THIS}}
+ \label{fig:ansi_diff}
+\end{figure}
+
+\section{Total variation}
+\label{sec:total_variation}
+
+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}
+ \begin{gathered}
+ \min_{u \in L^p(\Omega)} F(u) \\
+ F(u) = \underbrace{\int_\Omega \abs{u - f}^p \, dx}_{\text{fidelity
+ term}} + \beta \underbrace{\int_\Omega
+ \abs{\nabla u} \, dx}_{\mathclap{\text{regularization
+ term}}},
+ \label{eq:first_min_presentation}
+ \end{gathered}
\end{equation}
-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.
+where $p$ is normally taken to be 1 or 2.
+The fidelity term penalizes images $u$ that are far from the original
+image $f$, and thus controls the \emph{fidelity} of our solution. The
+regularization term is the total variation, and minimizing it will
+reduce the variation and regularize the image. The $\beta$ parameter
+controls the strength of the regularization.
+
+Since we do not only want to consider differentiable 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
- taken in the distributional sense, is
+ often written $\int_\Omega \abs{Du}\, dx$, where the $D$ is the
+ gradient taken in the distributional sense, is
\begin{equation}
\TV(u)
- = \int_\Omega \abs{Du}
- = \sup \left\{ \int_\Omega u \cdot \mathrm{div} \, \varphi :
+ = \int_\Omega \abs{Du} \, dx
+ = \sup \left\{ \int_\Omega u \cdot \diver \varphi \, dx :
\varphi \in C^\infty_c\left(\Omega, \mathbb{R}^2\right),
\norm{\varphi}_{L^\infty(\Omega)} \leq 1
\right\}.
\end{definition}
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}. 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.
+$\varphi$ vanish on the boundary of $\Omega$. This means we assume
+no-flux conditions on the boundary. As this restoration method is the
+one which will be extended later in this thesis, we will look a little
+bit more deeply into the background and the numerical methods relating
+to it.
+
+As we are searching for an image with low total variation, 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
\end{definition}
Our optimization problem has thus become
\begin{equation}
- \min_{u \in \BV(\Omega)} \int_\Omega \abs{u - v}^p + \beta \, \TV(u)
+ \min_{u \in \BV(\Omega)} \int_\Omega \abs{u - v}^p \, dx + \beta \,
+ \TV(u),
+ \label{eq:second_min_presentation}
\end{equation}
-\fixme{say something about what $p$ can be. is it problematic if $p=2$?}
+and there are many different numerical methods to find or approximate a
+solution.
+
+As with anisotropic diffusion, the total variation method also has its
+strengths and weaknesses. Its main strenght is its ability to recover
+edges in the input image. The total variation of a section only takes
+the absolute change into account, and does not favour gradual changes
+over sharp discontinuities like the diffusion methods. This comes at a
+cost though, and in some cases the method may introduce new edges that
+were not in the original image, since flat sections of zero variation
+are encouraged by the functional, an effect called the
+\emph{stair-casing effect}. Thin objects and corners may also suffer
+from contrast loss since bringing them closer to their surroundings in
+level value reduces the total variation.
+
+\begin{figure}
+ \centering
+ \begin{subfigure}[t]{0.4\textwidth}
+ \centering
+ \includegraphics[width=\textwidth]{fig/finger.png}
+ \caption{FINGER}
+ \end{subfigure}
+ ~
+ \begin{subfigure}[t]{0.4\textwidth}
+ \centering
+ \includegraphics[width=\textwidth]{fig/ftv.png}
+ \caption{TV}
+ \end{subfigure}
+ \caption{
+ Example of how the contrast of thin stuff might be reduced when
+ regularization is too high. \fixme{FIX}
+ }
+ \label{fig:tv_example}
+\end{figure}
\subsection{Chambolle's dual approach}
-Maybe, maybe not.
+A popular approach for solving the minimization problem, is the dual
+algorithm of Chambolle described in \cite{chambolle2004algorithm}. It
+has been improved on by others, and also extended to primal-dual
+methods. An overview of the different methods and further references can
+be found in \cite{caselles2011total}.
\subsection{Graph cut approach}
+In the theory of graphs, a graph cut is a set of edges that when
+removed will separate the graph into two disconnected parts. A minimum
+cut is a set of edges such that the sum of their weight is minimal. It
+has been shown that using special graph constructions, the graph cuts
+can be used to minimize a special type of energy functions.
+
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
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
+\}$. The thresholded image $u^\lambda$, an indicator function, is then
+defined as
\begin{equation}
u^\lambda = \idfun_{u > \lambda}.
\end{equation}
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}
- \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)$.
-\end{definition}
-\fixme{Something about reduced boundary and Hausdorff measures here?}
-The set perimeter behaves mostly as one might expect for ``nice'' sets,
-just note that if the boundary of $E$ overlaps with the boundary of
-$\Omega$, then $P(E,\Omega)$ will not include the overlapping parts of
-the boundary.
-
-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$. \fixme{repeated}
-\begin{theorem}[The coarea formula]
- Let $u \in BV(\Omega)$. Then for almost every $\lambda$ the set
- $\{u > \lambda\}$ is of finite perimeter, and one has the coarea
- formula
- \begin{equation}
- \TV(u)
- = \int_\Omega \abs{D u}
- = \int_\mathbb{R} \Per(\{u > \lambda\},\Omega) \, d\lambda
- \label{eq:coarea_formula}
- \end{equation}
-\end{theorem}
-Note that the equality in \eqref{eq:coarea_formula} holds for all $u \in
-L^1(\Omega)$, not only when the total variation is finite. We can now
-rewrite \eqref{eq:first_min_presentation} as
+
+Using the \emph{coarea formula} which we will spend more space on later,
+the total variation can then be rewritten into an integral over the
+perimeter of the level sets
\begin{equation}
- E_v(u) =
- \int_\Omega \abs{u - v}^p
- + \beta \int_\mathbb{R} \Per(\{u > \lambda\},\Omega) \, d\lambda.
- \label{eq:energy_perimeter}
+ \TV(u)
+ = \int_\Omega \abs{D u} \, dx
+ = \int_\mathbb{R} \int_\Omega \abs{D u^\lambda} \, dx \, d\lambda
+ = \int_\mathbb{R} \Per(\{u > \lambda\},\Omega) \, d\lambda.
+ \label{eq:coarea_formula}
\end{equation}
-\fixme{Should this be a minimization problem or an energy function.}
-
-\fixme{Shorten this section (or maybe it's short enough, only two pages,
-but we need to fill in some text at least), and add some general
-references.}
-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
+Our input images will be digital images given as pixels on a grid, and
+we will represent them as functions $u : \mathcal{G} \to \mathcal{L}$,
+where $\mathcal{G}$ is our regular discrete pixel grid, and
+$\mathcal{L} = \{0, \hdots, L-1\}$ is the set of discrete level values.
+Through careful manipulation of the functional in
+\eqref{eq:second_min_presentation} we obtain a discrete energy
+functional decomposed as a sum over all the level values
\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,
+ F(u) =
+ \sum_{\lambda=0}^{L-2} \sum_x F^x_\lambda(u^\lambda_x)
+ + \beta \sum_{\lambda = 0}^{L-2} \sum_{(x, y)} F^{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}
+where the sum over $(x, y)$ is over all pixel pairs $x$ and $y$ such
+that are in a neighbor relation, i.e.\ are ``close'' to each other.
+
+The actual form of the functional, and the steps to construct it will be
+presented later, but it can be useful to consider the intuition behind
+the construction, and how it relates to a graph cut.
+
+The energy $F^{x,y}$ represents how strongly the pixels at $x$ and $y$
+are connected, and will depend on their spatial relation. Pixels close
+to each other are more strongly related, and thus more likely to have
+similar values. The energy $F^x_\lambda$ represents how strongly we
+believe that the pixel at $x$ is connected to the level $\lambda$, and
+it is related to the level value of the input image $f$ in that pixel.
+
+By finding the minimum cut, we could say that we find the minimum set of
+relations that has to be broken in order to create a thresholded image
+$u^\lambda$.
-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.
+\fixme{wow. such handwavy.}
\section{Anisotropic total variation}
-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
+The method considered will build on the total variation regularisation
+method of Section \ref{sec:total_variation}. From anisotropic diffusion
+in Section \ref{sec:anisotropic_diffusion} we borrow the idea of making
+the regularization in each point directionally dependant. We introduce
+the anisotropic total variation
\begin{equation}
\TVA(u) = \int_\Omega \sqrt{\nabla u(x)^T A(x) \nabla u(x)} \, dx
\end{equation}
-for all $u \in C^1(\Omega)$. We see that if $A(x)$ is the identity
-matrix we get the usual total variation from \fixme{ref}. Our goal
-is to control $A(x)$ such that $\nabla u(x)$ is weighted down across
-edges, and thus the smoothing across the edges will be reduced.
+for all $u \in C^1(\Omega)$. If $A(x)$ is the identity
+matrix we get the regular total variation found in
+\eqref{first_min_presentation}. One problem with the regular total
+variation method is that when reducing the total variation, it will also
+try to reduce the variation over known edges in the image, which can
+lead to contrast loss, especially in thin details. By controlling $A(x)$
+such that $\nabla u(x)$ is weighted down across known edges, we hope to
+retain the regularization properties of the method while reducing some
+of the negative effects. If the variation across an edge is ``ignored''
+by the energy functional, there is no gain in reducing the height of the
+edge as before.
As we will not always be working with differentiable images, we extend
-the definition of the total variation functional to the space
-$\BV(\Omega)$. Being symmetric positive definite, the matrix $A(x)$ can
-be factored into two symmetric matrices as $A(x) = \Ahalf(x) \Ahalf(x)$. We can
-then write
+the definition of the total variation functional to the space of
+functions of bounded variation $\BV(\Omega)$. Being symmetric positive
+definite, the matrix $A(x)$ can be factored into two symmetric matrices
+as $A(x) = \Ahalf(x) \Ahalf(x)$. We can then write
\begin{align}
\TVA(u) &= \int_\Omega \abs{\Ahalf \nabla u} \, dx \\
&= \sup_{\abs{\xi(x)}
where $\xi$ and $\eta = \Ahalf \xi$ are in $C_c^\infty(\Omega,
\mathbb{R}^2)$, the space of $C^1$ vector fields with compact support.
If we let $\norm{\xi}_A = \sqrt{\xi^T A \xi}$ and $\norm{\eta}_A^* =
-\sqrt{\eta^T A^{-1} \eta}$, then our new extended definition reads
-\begin{equation}
- \TVA(u) = \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega u \diver \xi \,
- dx.
- \label{eq:extended_tv}
-\end{equation}
+\sqrt{\eta^T A^{-1} \eta}$.
+\begin{definition}[Anisotropic total variation]
+ For a function $u \in BV(\Omega)$ and a symmetric positive definite
+ tensor $A : \Omega \to \mathbb{R}^{2 \times 2}$ we define the
+ anisotropic total variation
+ \begin{equation}
+ \TVA(u) = \sup \left\{ \int_\Omega u \diver \xi \,
+ dx : \xi \in C_c^\infty(\Omega, \mathbb{R}^2), \norm{\xi}_A^*
+ \leq 1 \right\}.
+ \label{eq:extended_tv}
+ \end{equation}
+\end{definition}
+
+The question is then how to construct this anisotropy tensor $A(x)$ to
+get the improvements we hope for, and how this affects our numerical
+solution methods.
\fixme{Include the whole functional somewhere.}
\fixme{Describe how this is similar to a Euler time step in a related
\, \Delta \lambda.
\end{equation}
Note that we will later ignore the $\Delta \lambda$ difference, as we can
-just absorb this into the $\beta$ parameter. The perimeter is then
+just absorb this into the $\beta$ parameter of \fixme{ref}. The perimeter is then
calculated using a discretized version of the Cauchy--Crofton formula
introduced in Theorem \ref{thm:riemannian_cauchy_crofton}. Again, we
stop the sum at $L-2$ since the level set $\{ u > L - 1\}$ is empty and
&\approx \sum_{\ell_{\nu, \rho} \in \mathcal{L}_D}
\sum_{x \in \ell_{\nu, \rho} \cap C} \,
\frac{\det M(x)}{2\left(\nu^T
- \cdot M(x) \cdot \nu\right)^{\sfrac{3}{2}}} \, \Delta\rho \,
- \Delta\phi \\
+ \cdot M(x) \cdot \nu\right)^{\sfrac{3}{2}}} \, \Delta\ell_{\nu,
+ \rho} \\
&= \sum_\nu \sum_\rho \sum_{x \in \ell_{\nu, \rho} \cap C} \,
\frac{\det M(x)}{2\left(\nu^T
\cdot M(x) \cdot \nu\right)^{\sfrac{3}{2}}} \, \Delta\rho \,
\end{equation}
The set of lines $\mathcal{L}$ has been discretized to the lines
$\mathcal{L}_D$. Note that we are approximating the length of the
-\emph{differentiable} curve $C$ in $\Omega$. Further we need to
-discretize the domain $\Omega$, which will lead to discrete level sets
-$\{ u > \lambda\}$ and ``discrete'' curves.
+\emph{differentiable} curve $C$ in $\Omega$.
Since the final goal is to work with digital images, it makes sense to
discretize our domain $\Omega$ as a regular grid $\mathcal{G}$. Our
image is then reduced to a function $u : \mathcal{G} \to \mathcal{L}$.
\fixme{mathcal L is now two things.} Moreover, the level sets $\{ u >
\lambda\}$ will be functions taking the value of 0 or 1 on this grid, as
-shown in Figure \fixme{ref}. Thus the perimeter of these sets will be
-piecewise linear lines following the boundaries of the pixels in
-$\mathcal{G}$.
+shown in Figure \fixme{ref}.
The choice of our discrete set of lines $\mathcal{L}_D$ is important, as
it will decide the accuracy of our approximation in
\eqref{eq:cauchy_crofton_approx1}. We will only consider lines going
-through points in in our grid $\mathcal{G}$, and for now we will
-consider a discretization which is uniform throughout the domain. The
+through points at least two in in our grid $\mathcal{G}$, and for now we will
+consider a discretization which is uniform throughout the domain,
+meaning that $\Delta \rho$ is constant for each line family, and that in
+each grid point, there is a line from each family. \fixme{moar
+explanation.} The
set of lines can then be represented by the neighborhood of a pixel as
shown in Figure \ref{fig:line_neigh}. Extending the edges shown in the
-Figure gives all lines going through the point considered. Figure
+figure gives all lines going through the point considered. Figure
\ref{fig:line_family} shows all lines of a given family, i.e.\ lines
having the same angle parameter $\phi$.
calculations of these points will not fit into our graph cut framework
later, and thus for an edge $e$ we will consider only the question of
``did $e$ cross $C$ or not?'' This amounts to checking whether the
-terminals of $e$ lie on each side of the perimeter $C$, and the
+terminals of $e$ lie on each side of the curve $C$, and the
approximation is exact for zero or one intersection points, but will, as
we see in Figure \ref{fig:curve_edge}, not be entirely correct when we
have more.
The second difficulty is that in the discrete setting, we will only have
an approximation of the metric tensor $M(x)$ for each point $x \in
\mathcal{G}$, and it is thus not available for arbitrary intersection
-points in $\Omega$. For an edge $e$ we will utilize the average of the
+points in $\Omega$. For an intersection of edge $e$ we will utilize the average of the
tensor in the two endpoints of the edge. Thus for an intersection point
$x$ somewhere on the edge $e_{ab}$, we approximate the metric tensor by
\begin{equation}
\fixme{
we must define what we mean by a reasonable line family. meaning
each line goes through more than one grid point. and there are no
- grid points withoug a line through it
+ grid points without a line through it
}
\begin{figure}
\end{figure}
\begin{lemma}
-For each family of lines given by an angle parameter $\phi$ we have the
-relation
-\begin{equation}
- \delta^2 = \norm{e} \Delta \rho
-\end{equation}
+ For each family of lines given by an angle parameter $\phi$ in the
+ regular grid of size $\delta$ we have the relation
+ \begin{equation}
+ \delta^2 = \norm{e} \Delta \rho
+ \end{equation}
+ \label{lem:delta_rho}
\end{lemma}
\begin{proof}
- Consider two partitions of the plane $\mathbb{R}^2$ as shown in
- Figure \ref{fig:area_proof}. The blue squares area $\delta^2$,
+ This is easy to verify for simple angles $\phi = \{0,
+ \sfrac{\pi}{4}, \sfrac{\pi}{2}, \hdots\}$. For a general argument
+ consider the two partitions of the plane $\mathbb{R}^2$ as shown in
+ Figure \ref{fig:area_proof}. The blue squares have area $\delta^2$,
while the red rectangles have length $\norm{e}$ and width $\Delta
- rho$. The blue grid repeats with a period of $\delta$ in both
- directions. Each grid point is also the terminal of two edges going
+ \rho$. The blue partition repeats with a period of $\delta$ in both
+ directions.
+
+ Each grid point is also the terminal of two edges going
each in the directions $\phi$ and $-\phi$. And also from each grid
point we draw a line of length $\Delta \rho$ to the next line. From
this construction, the red partition is also periodic with period
$\delta$ in both directions.
- Consider a \fixme{donut} tile of size $\delta \cdot \delta$ where
- the left side is connected to the right, and the top to the bottom.
+ Now consider a square of size $n\delta \cdot n\delta$, where we
+ connect the right side to the left side, and the top to the bottom
+ so that it \fixme{topologically} is shaped like a donut. Equipped
+ with a grid structure as before, this surface will contain $n \cdot
+ n$ distinct grid points. Since our blue and red partitions are
+ $\delta$-periodic, we can use them to partition this surface,
+ without trouble.
+
+ We have one blue square for each grid point, and we also have one
+ red rectangle for each grid point, thus we have $n^2$ of each, and
+ their areas must be equal
+ \begin{equation}
+ \delta^2 = \norm{e} \Delta \rho.
+ \end{equation}
- We can use a tile of size $\delta \cdot \delta$ of the original red
- and blue partition to partition this \fixme{donut} tile as well.
-\end{proof}
-\begin{proof}
- This is easy to verify for simple angles $\phi = \{0,
- \sfrac{\pi}{4}, \sfrac{\pi}{2}, \hdots\}$. For a general argument
-consider the Figure \ref{fig:area_proof}, and call the blue squares
-$\delta$-squares, and the red rectangles $\rho$-rectangles. The
-$\delta$-squares have area $\delta^2$, while the $\rho$-rectangles have
-area $\Delta \rho \cdot \abs{e}$. Our goal is to show that these are equal.
-First note that we have a periodicity in both directions with period
-$\delta$, and each $\delta$-square looks exactly the same. Further,
-each $\rho$-rectangle is partitioned by the blue grid into the
-subsets $A_i$ such that
-\begin{equation}
- A = \bigcup_{i=1}^{i \leq m} A_i,
-\end{equation}
-where $A$ is the set of points in a $\rho$-rectangle. All
-$\rho$-rectangles start in a grid point, so the will decompose in the
-same way. Since all $\delta$-squares look the same, we can find all
-$A_i$ in a given $\delta$-square such that
-\begin{equation}
- \delta^2 \geq \abs{\bigcup_{i=1}^{i \leq m} A_i} = \Delta \rho \abs{e}.
-\end{equation}
-Further the $\rho$-rectangles partition the plane and thus each point in
-a $\delta$-square is also in a $\rho$-rectangle, and each
-$\rho$-rectangle consists of the same parts $A_i$, and thus
-\begin{equation}
- \delta^2 \leq \abs{\bigcup_{i=1}^{i \leq m} A_i} = \Delta \rho \abs{e}.
-\end{equation}
-\fixme{we are missing some argument that we don't have more than one
-$A_i$ in the square}
-We can then conclude that $\delta^2 = \Delta \rho \abs{e}$.
+ %Our goal is to show that the areas are equal. Each $\rho$-rectangle
+ %is partitioned by the blue grid into the subsets $A_i$ such that
+ %\begin{equation}
+ % A = \bigcup_{i=1}^{i \leq m} A_i,
+ %\end{equation}
+ %where $A$ is the set of points in a $\rho$-rectangle. All
+ %$\rho$-rectangles start in a grid point, so the will decompose in the
+ %same way. Since all $\delta$-squares look the same, we can find all
+ %$A_i$ in a given $\delta$-square such that
+ %\begin{equation}
+ % \delta^2 \geq \abs{\bigcup_{i=1}^{i \leq m} A_i} = \Delta \rho \abs{e}.
+ %\end{equation}
+ %Further the $\rho$-rectangles partition the plane and thus each point in
+ %a $\delta$-square is also in a $\rho$-rectangle, and each
+ %$\rho$-rectangle consists of the same parts $A_i$, and thus
+ %\begin{equation}
+ % \delta^2 \leq \abs{\bigcup_{i=1}^{i \leq m} A_i} = \Delta \rho \abs{e}.
+ %\end{equation}
+ %\fixme{we are missing some argument that we don't have more than one
+ %$A_i$ in the square}
+ %We can then conclude that $\delta^2 = \Delta \rho \abs{e}$.
\end{proof}
Inserting
-this into the curve length approximation we obtain
+this \fixme{and the tensor approx} into the curve length approximation
+of \fixme{ref} we obtain
\begin{equation}
\abs{C}_M \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}}}.
+ e\right)^{\sfrac{3}{2}}},
\end{equation}
+where the sum is over all edges crossing the curve an odd number of
+times.
Going back to the perimeter of the level set $\{ u > \lambda \}$ we can
rewrite the sum to be a sum over all edges that goes from one side of
the set to the other such that
\subsubsection{Consistency}
-As for any discretization these are important properties to consider.
Consistency relates to how well a solution to the continuous problem
fits in the discretized equation, in other words, how well the
discretized equation approximates the continuous one.
It is obvious that the discretization of the fidelity term in
\eqref{eq:fidelity_approx_1} is consistent. We have left out the pixel
size $\Delta x$ in the sum, and absorbed it into the regularization
-parameter $\beta$, but apart from that, the sum is the midpoint rule
+parameter $\beta$, but apart from that, the sum is a midpoint rule
approximation of the integral.
For the regularization term we will argue that for a differentiable
curve $C$, the discretization of our domain $\Omega$ and the set of
lines $\mathcal{L}$ gives a discrete Cauchy--Crofton formula that is
consistent with the continuous one. We will show that for an
-increasingly refined domain $\mathcal{G}$, there exists a choice for
+increasingly refined discrete domain $\mathcal{G}$, there exists a choice for
$\mathcal{L}_D$ that leads to a consistent Cauchy--Crofton formula. For
convenience we will use a neighborhood representation of $\mathcal{L}_D$
similar to the one in Figure \ref{fig:line_neigh}.
\end{figure}
Consider a square centered around a grid point with side lengths
-$\sqrt{\delta}$. As $\delta$ goes to zero, the size of this square will
-go to zero. Inside this square we can fit a square of $\lfloor 1 /
-\sqrt{\delta} \rfloor^2$ grid points. This means that the number of
-grid points along the outer edge of this square $\lfloor 1 /
-\sqrt{\delta} \rfloor$ goes to infinity. For each grid point along
-the outer edge of this square, we include in our neighborhood a grid point
-having the same angle $\phi$ to the $x$-axis. This means either
-including the actual grid point at the outer edge, or one having the
-same angle, just closer to the center. This construction can be
-seen in Figure \ref{fig:square_cons}.
+$\sqrt{\delta}$ as shown in Figure \ref{fig:square_cons}. As $\delta$
+goes to zero, the size of this square will go to zero. Inside this
+square we can fit a square of $\lfloor 1 / \sqrt{\delta} \rfloor^2$ grid
+points. This means that the number of grid points along the outer edge
+of this square $\lfloor 1 / \sqrt{\delta} \rfloor$ goes to infinity. For
+each grid point along the outer edge of this square, we include in our
+neighborhood a grid point having the same angle $\phi$ to the $x$-axis.
+This means either including the actual grid point at the outer edge, or
+one having the same angle, just closer to the center. This construction
+can be seen in Figure \ref{fig:square_cons}.
The maximal $\Delta \phi$ will then be between the
horizontal or vertical edge and its neighbors, shown in Figure
\ref{fig:square_cons} as angle $a$. These angles can be calculated to be
\rfloor}{2}} \leq \arctan \frac{\delta}{\sfrac{\sqrt{\delta}}{2} -
\delta} = \arctan \frac{1}{\frac{1}{2\sqrt{\delta}} - 1} \to 0.
\end{equation}
+\fixme{the flooring here is not right}
-Further we know that for each line family
+Further we know from Lemma \ref{lem:delta_rho} that for each line family
\begin{equation}
\delta^2 = \Delta \rho \norm{e},
\end{equation}
\sup \Delta \rho \leq \frac{\delta^2}{\delta} = \delta \to 0.
\end{equation}
-Since all edges are inside our square, the edge length is bounded from
-above by $\sqrt{2} \delta$ and will also go to zero. Thus for a given
+Since all edges are inside our square,
+we have $\norm{e} \leq \sqrt{2} \delta \to 0$. Thus for a given
curve $C$, our simplification of only counting an intersection between
$e$ and $C$ when $C$ crosses $e$ an odd number of times becomes
increasingly correct.
projects / master / commitdiff
