weaknesses. The method considered in this thesis is an anisotropic total
variation formulation, and the hope is that we keep the strengths of the
anisotropic diffusion and total variation methods, while eliminating
-some of their weaknesses.
+some of their respective weaknesses.
This chapter will be devoted to the continuous formulation of the
-method. We will look at the functional we want to minimize, its
-different forms and briefly discuss its well-posedness. Through the
-coarea formula, the anisotropic total variation is rewritten as an
-integral of the perimeter of all the level sets of the image.
+method. We will look at the functional we want to minimize and its
+different forms, and briefly discuss its well-posedness. Through the
+anisotropic coarea formula, the anisotropic total variation is rewritten
+as an integral of the perimeter of all the level sets of the image.
-Later, the Cauchy--Crofton formula is introduced to make it feasible to
-calculate the perimeter of these level sets. All of this leads up to the
-discretization of our functional in the following chapter.
+Later, the anisotropic Cauchy--Crofton formula is introduced to make it
+feasible to calculate the perimeter of these level sets. All of this
+leads up to the discretization of our functional in the next chapter.
\section{Anisotropic total variation}
\TVA(u) = \int_\Omega \sqrt{\nabla u(x)^T A(x) \nabla u(x)} \, dx
\label{eq:aniso_tv_sqrt}
\end{equation}
-for all $u \in C^1(\Omega)$. If $A(x)$ is the identity matrix we get the
-regular total variation found in \eqref{eq: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.
+for all $u \in C^1(\Omega)$. We assume here that $A(x)$ is continuous
+and positive definite,
+and we will later need the eigenvalues of $A(x)$ to be uniformly bounded
+below and above. If $A(x)$ is the identity matrix we get the
+regular total variation found in \eqref{eq:first_min_presentation}.
+When reducing the regular total variation, we will also try to reduce
+the variation over known edges in the image. This can lead to unwanted
+contrast loss, expecially in fine details.
+By controlling $A(x)$ such that the contribution of $\nabla u(x)$ is
+reduced across known edges, we hope to retain the regularization
+properties of the original method while reducing some of the negative
+effects. If the variation across an edge is ``ignored'' by the
+functional, there is no gain in reducing the height of the edge as
+before.
Note that $u(x)$ and $A(x)$ are always dependent on the position in the
-image $x$, but we will sometimes drop writing the $x$, when no
-misunderstandings are possible.
+image $x$, but we will sometimes drop the $x$, when no misunderstandings
+are possible.
As we will not always be working with differentiable images, we extend
the definition of the total variation functional. 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)}
- \leq 1} \int_\Omega (\Ahalf \nabla u)^T \xi \, dx \\
- &= \sup_{\abs{\xi(x)}\leq 1} \int_\Omega \nabla u \cdot \Ahalf\xi \, dx
- \\
- &= \sup_{\abs{\xi(x)}\leq 1} \int_\Omega u \diver (\Ahalf\xi) \, dx \\
- &= \sup_{\eta^T A^{-1} \eta \leq 1} \int_\Omega u \diver \eta \,
- dx,
-\end{align}
+positive definite, the matrix $A$ can be factored into two symmetric
+matrices as $A = \Ahalf \Ahalf$. We can then write
+\begin{equation}
+ \begin{aligned}
+ \TVA(u) &= \int_\Omega \abs{\Ahalf \nabla u} \, dx \\
+ &= \sup_{\abs{\xi(x)}
+ \leq 1} \int_\Omega (\Ahalf \nabla u)^T \xi \, dx \\
+ &= \sup_{\abs{\xi(x)}\leq 1} \int_\Omega \nabla u \cdot
+ \Ahalf\xi \, dx \\
+ &= \sup_{\abs{\xi(x)}\leq 1} \int_\Omega u \diver (\Ahalf\xi) \,
+ dx \\
+ &= \sup_{\eta^T A^{-1} \eta \leq 1} \int_\Omega u \diver \eta \,
+ dx,
+ \end{aligned}
+\end{equation}
where $\xi$ and $\eta = \Ahalf \xi$ are in $C_c^\infty(\Omega,
\mathbb{R}^2)$, the space of $C^\infty$ vector fields with compact support.
-In the following we let $\norm{\xi}_A = \sqrt{\xi^T A \xi}$ and
-$\norm{\eta}_A^* = \sqrt{\eta^T A^{-1} \eta}$, and with that we present
-the formal definition of the anisotropic total variation.
+In the following we let $\norm{\xi}_A = \sup_x \sqrt{\xi^T A \xi}$ and
+$\norm{\eta}_A^* = \sup_x \sqrt{\eta^T A^{-1} \eta}$, and with that we
+present the formal definition of the anisotropic total variation.
\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
+ For a function $u \in BV(\Omega)$ and a continuous 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^*
\end{equation}
Similar functionals have been considered in
\cite{grasmair2010anisotropic} and \cite{olsson2009extending}. The
-question is now how to construct this anisotropy tensor $A(x)$ to get
-the improvements we hope for, and how this affects our numerical
-solution methods.
+question is now how to construct the anisotropy tensor $A(x)$ to get
+the improvements we hope for, and how the introduction of the tensor
+affects our numerical solution method.
\subsection{Anisotropy tensor}
\label{sec:anisotropy_tensor}
initial total variation $\abs{\nabla u}$.}
There are many possible choices for the anisotropy tensor $A(x)$. Our
-constraints are that we have assumed it to be symmetric positive
-definite, and we have some wishes for its properties. We would first and
-foremost like it to down-weight $\nabla u$ in \eqref{eq:aniso_tv_sqrt}
-across true edges, while maintaining normal regularization properties in
-smooth sections.
+constraints are that we have assumed it to be continuous and symmetric
+positive definite, and we have some wishes for its properties. We would
+first and foremost like it to down-weight $\nabla u$ in
+\eqref{eq:aniso_tv_sqrt} across true edges, while maintaining normal
+regularization properties in smooth sections.
By true edges we mean that that we do not want the tensor to be
sensitive to noise in the image, and thus find edges where there are
Edges can be found in many different ways, but as suggested by Weickert
in his book on Anisotropic Diffusion \cite{weickert1998anisotropic}, and
briefly mentioned in Section~\ref{sec:anisotropic_diffusion}, a good
-starting point is the \emph{edge detector} $\nabla f_\sigma$,
-where $f_\sigma = K_\sigma * \tilde{f}$ and $\tilde{f}$ is the symmetric
-extension of the initial image $f$ in $\mathbb{R}^2$. The smoothing
-parameter $\sigma$ is called the \emph{noise scale}, and with this we
-aim to avoid being too sensitive to noise in $f$.
+starting point is the \emph{edge detector} $\nabla f_\sigma$.
+The image is smoothed by a Gaussian filter as described in
+Section~\ref{sec:diffusion_filtering}: $f_\sigma = K_\sigma * \tilde{f}$,
+where $\tilde{f}$ is the symmetric extension of the initial image $f$ in
+$\mathbb{R}^2$. The smoothing parameter $\sigma$ is called the
+\emph{noise scale}, and with this we aim to avoid being too sensitive to
+noise in $f$.
As seen in Figure~\ref{fig:edges}, the edge detector is fine for
detecting edges, but it can not give us information about larger
structures, like corners and textures, which is why we introduce
-the structure tensor
+the \emph{structure tensor} $S_\rho(x)$.
+First consider the tensor $S_0(x) = \nabla f_\sigma(x) \otimes
+\nabla f_\sigma(x)$. It is symmetric positive semi-definite, and
+obviously contains no more information than the edge detector itself.
+Its eigenvalues are $\lambda_1 = \labs{\nabla f_\sigma(x)}^2$ and
+$\lambda_2 = 0$ with corresponding eigenvectors $v_1$ and $v_2$ parallel
+and perpendicular to $\nabla f_\sigma(x)$ respectively.
+
+We are also interested in identifying and being sensitive to features in
+a neighborhood around the point $x$, such as corners or curved edges and
+\fixme{coherent structures}. This is why we introduce the component-wise
+convolution with $K_\rho$ such that
\begin{equation}
S_\rho(x) := K_\rho * \big( \nabla f_\sigma(x) \otimes
- \nabla f_\sigma(x) \big).
+ \nabla f_\sigma(x) \big)(x).
\label{eq:s_def}
\end{equation}
+The parameter $\rho$, called
+the \emph{integration scale}, controls the size of the neighborhood
+which affects the structure tensor. It controls the size of the
+structures we want our anisotropy tensor to be sensitive to.
\begin{figure}
\centering{}
\label{fig:edges}
\end{figure}
-First consider the tensor $S_0(x) = \nabla f_\sigma(x) \otimes
-\nabla f_\sigma(x)$. It obviously contains no more information
-than the edge detector itself. Its eigenvalues will be 0 and
-$\labs{\nabla f_\sigma(x)}^2$ with corresponding eigenvectors $v_1$ and
-$v_2$ perpendicular and parallel to $\nabla f_\sigma(x)$ respectively.
-
-We are also interested in identifying and being sensitive to features in
-a neighborhood around the point $x$, such as corners or curved edges and
-\fixme{coherent structures}. This is why we introduce the convolution
-with $K_\rho$, which is done component-wise. The parameter $\rho$, called
-the \emph{integration scale} thus controls the size of the neighborhood
-which affects the structure tensor. It controls the size of the
-structures we want our anisotropy tensor to be sensitive to.
-
-Being a smoothed combination of symmetric positive semi-definite
-matrices, $S_\rho(x)$ can easily be verified to be symmetric positive
-semi-definite itself. We order its two real eigenvalues $\lambda_1 \geq
-\lambda_2$ with corresponding eigenvectors $v_1$ and $v_2$. From the
-characteristic polynomial of $S_\rho(x) = (\begin{smallmatrix} s_{11} &
-s_{12} \\ s_{12} & s_{22} \end{smallmatrix})$ we obtain a closed form
+The smoothed tensor $S_\rho(x)$ can easily be verified to be symmetric
+positive semi-definite, just like $S_0(x)$. We order the two real
+eigenvalues $\lambda_1 \geq \lambda_2$ with corresponding eigenvectors
+$v_1$ and $v_2$. From the characteristic polynomial of $S_\rho(x) =
+(\begin{smallmatrix} s_{11} & s_{12} \\ s_{12} & s_{22}
+\end{smallmatrix})$ we obtain a closed form
expression for the eigenvalues
\begin{equation}
\lambda = \frac{1}{2} \left( s_{11} + s_{22} \pm \sqrt{(s_{11} -
\approx \lambda_2 \gg 0$. Thus the quantity $(\lambda_1 - \lambda_2)^2$
will be large around edges and small in smooth or non-coherent areas.
-We begin our anisotropy tensor construction by eigendecomposing the
-structure tensor as
+To extract this information from the structure tensor, we decompose it
+as
\begin{equation}
S_\rho(x) = U(x) \Lambda(x) U(x)^T,
\end{equation}
-where $\Lambda(x)$ is a matrix with the eigenvalues of $S_\rho(x)$ on the
-diagonal, while $U(x)$ is a rotation matrix and has the eigenvectors of
-$S_\rho(x)$ as its columns.
-
-The eigen-decomposition give us
+where
\begin{equation}
\Lambda(x) = \begin{pmatrix}
\lambda_1 & 0 \\
0 & \lambda_2
\end{pmatrix}
\end{equation}
-where $\lambda_1 \geq \lambda_2$, then we construct a new matrix $A(x) =
-U(x) \Sigma(x) U(x)^T$ where
+has the eigenvalues $\lambda_1 \geq \lambda_2$ on its diagonal, while
+$U(x)$ is a rotation matrix and has the eigenvectors of $S_\rho(x)$ as
+its columns. From this we construct a new matrix $A(x) = U(x) \Sigma(x)
+U(x)^T$ where
\begin{equation}
\Sigma(x) = \begin{pmatrix}
\sigma_1 & 0 \\
0 & \sigma_2
\end{pmatrix}.
\end{equation}
-To reverse the effect of $S_\rho(x)$ such that $A(x)$ actually weighs down
-the norm across edges, we need $\sigma_1 \leq \sigma_2$, so we construct
+and for $\sigma_1$ and $\sigma_2$ we choose
\begin{equation}
\begin{aligned}
\sigma_1 &= \left(1 +
\end{aligned}
\label{eq:sigma_construction}
\end{equation}
-This way, $\sigma_1 \in (0, 1]$. The rotation is kept, while the size of
-the eigenvalues are changed.
+Thus the eigenvectors of $A(x)$ and $S_\rho(x)$ are equal, while the
+eigenvalues are different. In smooth areas, $\sigma_1 \approx 1$ and
+$A(x)$ will be close to the identity matrix, while $A(x)$ will reduce
+the effect of $\nabla u$ across edges.
+
+The parameter $\omega$ controls the amount of anisotropy in the method.
+
+Around corners $A(x)$ will be close to the identity matrix, which gives
+regularization similar to smooth areas. This is one possible down-side
+of this tensor choice.
+
+For the case where $\lambda_1 = \lambda_2$, the $U(x)$ in our
+decomposition is not well-defined. This is not a problem though, since
+$\Sigma(x)$ will be the identity matrix, so any orthogonal matrix will
+suffice for $U(x)$.
\fixme{%
Discuss $\omega$. Discuss other tensor choices, the coherency thing
from Weickert. Discuss that this is not optimal in corners.
}
-Note that in smooth parts of the image we have $\sigma_1 \approx
-\sigma_2 = 1$ and the anisotropic total variation is close to the
-regular total variation.
-
\fixme{%
This should be visualized, with a figure showing the length and
direction of the eigenvalues in the area around an edge.
restrict ourself to a nice space such as $L^2(\Omega)$ we will at some
point run into problems. The discussion here is not meant to give the
most rigorous background, but rather an overview of what needs to be
-shown. Some problems will be worked around, while other will be skipped
+shown. Some problems will be worked around, while others will be skipped
with a reference to further theory.
The basic things we ask of our functional
F(u) = \int_\Omega (u - f)^2 + \beta \, \TVA(u)
\end{equation}
are lower semicontinuity
-and coercivity for existence and convexity for uniqueness. We restrict
-ourself to $L^2(\Omega)$, and leave the extension to $\BV(\Omega)$ to
-someone else \fixme{Ref, and is it really an extension? More like a
-constriction?}
+and coercivity for existence, and convexity for uniqueness. We restrict
+ourself to $L^2(\Omega)$ which makes sense with our fidelity term,
+assuming that $f \in L^2(\Omega)$ initially.
\subsection{Convexity}
\begin{equation}
\int_\Omega (u - f)^2 \, dx
\end{equation}
-is obviously strictly convex. It can be shown by expanding and
+is obviously strictly convex. This can be shown by expanding and
rearranging the strict convexity condition
\begin{equation}
\int_\Omega (\lambda x + (1-\lambda)y - f)^2 \, dx < \lambda
projects / master / commitdiff