--- /dev/null
+\chapter{Continuous formulation}
+
+From the introduction we see that there are many different approaches to
+the image restoration problem, all with their own strengths and
+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 respective weaknesses.
+
+This chapter will be devoted to the continuous formulation of the
+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 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}
+
+\fixme{Rating: 8/10}
+
+The method considered will build on the total variation regularization
+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 dependent. We introduce
+the anisotropic total variation
+\begin{equation}
+ \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)$. 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 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$ 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 = \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 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^*
+ \leq 1 \right\}.
+ \label{eq:extended_tv}
+ \end{equation}
+ \label{def:extended_tv}
+\end{definition}
+\nomenclature{$\TVA(u)$}{Anisotropic total variation of image $u$, with
+anisotropy tensor $A(u)$.}%
+
+With this extended definition, we have arrived at a minimization problem
+where we seek to find a minimizer of the functional
+\begin{equation}
+ F(u) = \int_\Omega (u - f)^2 \, dx + \beta \, \TVA(u).
+ \label{eq:first_anisotropic_functional}
+\end{equation}
+Similar functionals have been considered in
+\cite{grasmair2010anisotropic} and \cite{olsson2009extending}. The
+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}
+
+\fixme{Rating: 6/10}
+
+There are many possible choices for the anisotropy tensor $A(x)$. Our
+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
+none, so we somehow want to be sure about the edges we find.
+
+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$.
+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 \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)(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{}
+ \includegraphics[width=0.4\textwidth]{fig/edges.png}
+ \caption{%
+ The largest eigenvalue of the structure tensor is $\abs{\nabla
+ f_\sigma(x)}^2$ and can as we see function as an edge detector.
+ }
+ \label{fig:edges}
+\end{figure}
+
+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} -
+ s_{22})^2 + 4 s_{12}^2} \right).
+ \label{eq:s_eigenvalues}
+\end{equation}
+
+The vector
+$v_1$ will then indicate the direction of most variation in the
+neighborhood. From this we see that an edge would give $\lambda_1 \gg
+\lambda_2 \approx 0$, while smooth areas would give $\lambda_1 \approx
+\lambda_2 \approx 0$. In corners we have variation in the direction of
+$v_1$ but also perpendicular to $v_1$, so we will have $\lambda_1
+\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.
+
+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
+\begin{equation}
+ \Lambda(x) = \begin{pmatrix}
+ \lambda_1 & 0 \\
+ 0 & \lambda_2
+ \end{pmatrix}
+\end{equation}
+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}.
+ \label{eq:sigma_construction}
+\end{equation}
+and for $\sigma_1$ and $\sigma_2$ we choose
+\begin{equation}
+ \begin{aligned}
+ \sigma_1 &= \left(1 +
+ \frac{(\lambda_1 - \lambda_2)^2}{\omega^2}\right)^{-1}, \\
+ \sigma_2 &= 1.
+ \end{aligned}
+ \label{eq:sigma_construction}
+\end{equation}
+Thus the eigenvectors of $A(x)$ and $S_\rho(x)$ are equal, while the
+eigenvalues are different. A visualization of the two tensors can be
+seen in Figure~\ref{fig:tensor_viz} where the two tensors are shown at
+an edge in the image.
+
+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.
+
+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.
+
+The parameter $\omega$ controls the amount of anisotropy in the method,
+such that if it is very large we are left with the identity matrix. Note
+that changing the parameter $\omega$ implicitly affects how much
+restoration is done. For an image $u$, decreasing $\omega$ will, all
+else being equal, decrease the lowest eigenvalue of $A(x)$ and in turn
+decrease the anisotropic total variation $\TVA(u)$.
+
+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)$.
+
+See \cite{weickert1999coherence} for a different tensor construction,
+made to enhance flow structures in the image, relevant in for example
+fingerprint analysis.
+
+\fixme{argue that $A$ is continuous?}
+
+\begin{figure}
+ \input{fig/tensor_viz}
+\end{figure}
+
+\section{Well-posedness}
+
+\fixme{Rating: 7/10}
+
+The theory of existence and uniqueness for these kinds of variational
+methods is a minefield of more or less subtle problems. Even if we
+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 others will be skipped
+with a reference to further theory.
+
+The basic things we ask of our functional
+\begin{equation}
+ F(u) = \int_\Omega (u - f)^2 + \beta \, \TVA(u)
+\end{equation}
+to have a \emph{well-posed} problem are lower semicontinuity
+and coercivity for existence, convexity for uniqueness and stability. We restrict
+ourself to $L^2(\Omega)$ which makes sense with our fidelity term,
+assuming that $f \in L^2(\Omega)$ initially.
+
+We consider the weak topology, as it will allow us to arrive at an
+existence result relatively easily.
+We say that a sequence $f_n$ in $L^2(\Omega)$ converges weakly to $f$ if
+\begin{equation}
+ \lim_{n \to \infty} \int_\Omega f_n \, \xi \, dx = \int_\Omega f \,
+ \xi \, dx
+\end{equation}
+for all $\xi \in L^2(\Omega)$ and we write $f_n \rightharpoonup f$. A
+weakly convergent sequence is a sequence that converges in the weak topology.
+
+\subsection{Convexity}
+
+We start with convexity as it is the easiest to show. Being quadratic,
+the fidelity term of our functional
+\begin{equation}
+ \int_\Omega (u - f)^2 \, dx
+\end{equation}
+is obviously strictly convex. This can be shown by expanding and
+rearranging the strict convexity condition
+\begin{equation}
+ \int_\Omega (\lambda u_1 + (1-\lambda)u_2 - f)^2 \, dx < \lambda
+ \int_\Omega (u_1 - f)^2 \, dx + (1 - \lambda) \int_\Omega (u_2 - f)^2 \,
+ dx
+\end{equation}
+to obtain that it is equivalent to
+\begin{equation}
+ - \lambda(1 - \lambda) \int_\Omega (u_1 - u_2)^2 \, dx < 0
+\end{equation}
+which is true for $0 < \lambda < 1$ and $u_1 \neq u_2$.
+
+The anisotropic total variation
+\begin{equation}
+ \TVA(u) = \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega u \diver \xi \,
+ dx
+\end{equation}
+can be thought of as---and has the properties of---a semi-norm, and is
+therefore convex. The sum of the two is thus strictly convex, which,
+given the existence of a minimizer, implies uniqueness.
+
+\subsection{Coercivity}
+
+Coercivity relates to how the functional behaves when the norm of the
+image $u$ tends to infinity. What we need in order to conclude with existence
+is sequential coercivity. Thus we need all level sets $F^\alpha = \{ u \in
+L^2(\Omega) : F(u) \leq \alpha \}$ to be \emph{sequentially
+pre-compact}, meaning that all sequences in the set contain a
+subsequence converging to an element of the closure of the set.
+
+It is obvious from the fidelity term that for some fixed $f \in
+L^2(\Omega)$, if $\lnorm{u}_{L^2} \to \infty$ then $F(u) \to \infty$.
+This implies that all the level sets $F^\alpha$ are bounded.
+Since $L^2(\Omega)$ is a Hilbert space, all bounded sequences contain a
+weakly convergent subsequence. Thus all the level sets $F^\alpha$ are
+weakly sequentially pre-compact, since all sequences in these sets have
+a subsequence weakly converging to a point in the closure of the set.
+
+\begin{figure}
+ \input{fig/lower_semicont}
+\end{figure}
+
+\subsection{Lower semi-continuity}
+
+The lower semicontinuity is the most tricky part, and this is where we
+will take some shortcuts. Lower semicontinuity for a functional $F$ at a
+point $u$ means that at points $u_\epsilon$ close to $u$, the functional
+takes values either close to or above $F(u)$. More specifically, for
+every sequence $u_k$ converging to $u$, we have $F(u) \leq \liminf_k
+F(u_k)$. For a function $f : \mathbb{R} \to \mathbb{R}$ this can be
+visualized as in Figure~\ref{fig:lower_semicont}.
+
+Since our space $L^2(\Omega)$ is of infinite dimensions things become a
+little problematic here. The problem lies in the fact that a
+functional which is continuous with respect to sequences is not
+necessarily continuous with respect to the underlying topology. In
+other words, in these spaces, there can be a difference between sequential
+continuity and topological continuity. Topological continuity implies
+sequential continuity, but the converse does not hold. One way to get
+around this would be to consider topological \emph{nets}, an extension
+of sequences, but for simplicity, and because it might not add much to
+the understanding of the restoration method, we will stick to proving
+sequential lower semi-continuity and referring to further theory. For
+further reading on the theory of sequential versus topological
+continuity see for example Megginson's book on Banach space theory
+\cite{megginson}.
+
+The mapping $u
+\mapsto \int_\Omega u \, \xi \, dx$ is weakly continuous for all $\xi
+\in L^2(\Omega)$. Note that when we write weakly continuous it is not a
+weaker version of continuity, but rather continuity in the weak
+topology, and the same goes for weak lower semi-continuity.
+
+Before arguing that our own functional is sequentially weakly lower
+semi-continuous, we present a needed result.
+\begin{lemma}
+ Assume that the functional $F : L^2(\Omega) \to \mathbb{R}$ is
+ defined by
+ \begin{equation}
+ F = \sup_i F_i
+ \end{equation}
+ where all the $F_i$ are sequentially weakly lower semi-continuous, then $F$
+ is sequentially weakly lower semi-continuous, meaning that for any sequence
+ $u_k \rightharpoonup u$ we have $F(u) \leq \liminf_k F(u_k)$.
+ \label{lem:sup_semi_cont}
+\end{lemma}
+\begin{proof}
+ For any sequence $u_k \rightharpoonup u$ in $L^2(\Omega)$ we have
+ \begin{equation}
+ F(u) = \sup_i F_i(u) \leq \sup_i \liminf_{k \to \infty} F_i(u_k)
+ \end{equation}
+ from the sequential weak lower semi-continuity of $F_i$. Using that
+ $\liminf_{k \to \infty} u_k = \sup_k \inf_{l \geq k} u_l$, we obtain
+ \begin{equation}
+ \begin{aligned}
+ F(u) &\leq \sup_i \sup_k \inf_{l \geq k} F_i(u_l) \\
+ &= \sup_k \sup_i \inf_{l \geq k} F_i(u_l) \\
+ &\leq \sup_k \inf_{l \geq k} \sup_i F_i(u_l) \\
+ &= \liminf_{k \to \infty} F(u_k)
+ \end{aligned}
+ \end{equation}
+ which proves that $F$ is sequentially weakly lower semi-continuous.
+\end{proof}
+
+In our functional in \eqref{eq:first_anisotropic_functional}, we first
+consider the fidelity term, and rewrite it as a supremum
+\begin{equation}
+ \int_\Omega (u - f)^2 \, dx
+ %= \sup_{\substack{\xi \in L^2(\Omega) \\ \norm{\xi}_{L^2} \leq
+ %\norm{u-v}_{L^2}}} \int_\Omega (u - v) \, \xi \, dx.
+ = \sup \left\{\int_\Omega (u - f) \, \xi \, dx : \xi \in
+ L^2(\Omega), \abs{\xi(x)} \leq \abs{u(x)-f(x)} \right\}
+\end{equation}
+As the map $u \mapsto \int_\Omega (u - v) \xi\, dx$ is continuous in the
+weak topology, the fidelity term is thus a supremum of weakly continuous
+functionals, and is thus by Lemma~\ref{lem:sup_semi_cont} sequentially
+lower semi-continuous.
+
+For the regularization term the approach is similar. With our extended
+definition from \eqref{eq:extended_tv}, we have
+\begin{equation}
+ %\TVA(u) = \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega u \, \diver \xi \,
+ %dx
+ \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\}
+\end{equation}
+This is again a
+supremum of weakly continuous functionals. Thus the regularization term
+is by Lemma~\ref{lem:sup_semi_cont} also sequentially weakly lower
+semi-continuous.
+
+The sum of the two terms is trivially sequentially weakly lower
+semi-continuous functional since
+\begin{equation}
+ \begin{aligned}
+ F_1(u) + F_2(u) &\leq \liminf_{k \to \infty} F_1(u_k) + \liminf_{k
+ \to \infty} F_2(u_k) \\
+ &= \lim_{k \to \infty} \left( \inf_{l \geq k}
+ F_1(u_l) + \inf_{l \geq k} F_2(u_l) \right) \\
+ &\leq \liminf_{k \to
+ \infty} \left( F_1(u_k) + F_2(u_k) \right),
+ \end{aligned}
+\end{equation}
+and thus our functional is sequentially weakly lower semi-continuous.
+
+The usual ways of going from coercivity and lower semicontinuity to
+existence do not work in infinite dimensions. But with sequential
+coercivity and sequential lower semi-continuity in the weak topology we
+can conclude that we have existence from \cite[Theorem
+5.1]{scherzer2008variational}.
+
+\subsection{Stability}
+
+The last property normally required for well-posedness is stability.
+Thus the question is how the solutions behave when we perturb the
+problem. We stick to $L^2(\Omega)$ and adopt a proof given by Acar and
+Vogel in \cite{acar1994analysis}.
+
+\begin{proposition}
+ Consider a sequence of perturbed problems
+ \begin{equation}
+ \min_{u \in L^2(\Omega)} F_n(u).
+ \end{equation}
+ Assume that each of the functionals $F_n$ is weakly sequentially
+ lower-semicontinuous, has a unique minimizer $u_n$, and in addition that they are
+ uniformly coercive in the sense that for any sequence $v_n \in
+ L^2(\Omega)$
+ \begin{equation}
+ \lim_{n \to \infty} F_n(v_n) = \infty \quad \text{whenever } \, \lim_{n \to
+ \infty} \norm{v_n}_{L^2} = \infty.
+ \label{eq:fn_coercive}
+ \end{equation}
+ Assume also that the problems are consistent in the sense that $F_n \to
+ F$ uniformly, i.e.\ given an $R > 0$ and $\epsilon > 0$ there exists an
+ $N$ such that
+ \begin{equation}
+ \abs{F_n(u) - F(u)} < \epsilon \quad \text{whenever } \, n \geq N \text{ and
+ } \norm{u}_{L^2} \leq R.
+ \label{eq:fn_concistent}
+ \end{equation}
+ The problem is then stable with respect to the perturbations $F_n$ in
+ the sense that if $u^*$ is a minimizer of $F$ we have
+ \begin{equation}
+ u_n \rightharpoonup u^*.
+ \label{eq:un_u_weak}
+ \end{equation}
+\end{proposition}
+\begin{proof}
+ As $u_n$ is a minimizer of $F_n$ we have
+ \begin{equation}
+ F_n(u_n) \leq F_n(u^*)
+ \end{equation}
+ and using \eqref{eq:fn_concistent} we get
+ \begin{equation}
+ %\lim F_n(u_n) \leq \limsup F_n(u_n) \leq \limsup F_n(u^*)
+ %= F(u^*) < \infty.
+ \lim_{n \to \infty} F_n(u_n) \leq \lim_{n \to \infty} F_n(u^*) =
+ F(u^*) < \infty.
+ \end{equation}
+ From the coercivity in \eqref{eq:fn_coercive} we obtain that
+ the sequence $u_n$ is bounded in $L^2(\Omega)$ and thus contains a
+ weakly convergent subsequence.
+
+ Assume that the convergence in \eqref{eq:un_u_weak} does not hold,
+ and denote the weakly convergent subsequence by $u_{n_j}
+ \rightharpoonup \bar{u} \neq u^*$. By the weak sequential lower
+ semicontinuity we have
+ \begin{equation}
+ \begin{aligned}
+ F(\bar{u}) &\leq \liminf_{j \to \infty} F\left(u_{n_j}\right) \\
+ &= \liminf_{j \to \infty} \left(F\left(u_{n_j}\right) -
+ F_{n_j}\left(u_{n_j}\right)\right) + \liminf_{j \to \infty}
+ F_{n_j}\left(u_{n_j}\right) \\
+ &= \liminf_{j \to \infty} F_{n_j}\left(u_{n_j}\right) \\
+ &\leq F(u^*)
+ \end{aligned}
+ \end{equation}
+ which contradicts the uniqueness of $u^*$ and thus we can conclude
+ that the solutions of the perturbed problems converge $u_n
+ \rightharpoonup \bar{u} = u^*$.
+\end{proof}
+
+Note that these are properties of the continuous problem, and that the
+discretizations and numerical methods used for approximating a solution
+can have its own issues with consistency, convergence and stability.
+
+\section{Anisotropic coarea formula}
+
+\fixme{Rating: 6.5/10, especially the part after the proof needs some
+work.}
+
+\fixme{More flow between the sections. Somehow.}
+
+The anisotropic coarea formula we will present here allows us to write
+the anisotropic total variation as an integral over the levels of the
+image. For a similar presentation of the regular coarea formula for all
+$f \in \BV(\Omega)$ see \cite{evans1991measure}.
+
+First we define the thresholded image at level $s$.
+\begin{definition}[Thresholded image]
+ The thresholded image at level $s$ is the function
+ \begin{equation}
+ u^s(x) = \begin{cases}
+ 1 & \text{if } u(x) > s, \\
+ 0 & \text{otherwise.}
+ \end{cases}
+ \end{equation}
+ \label{def:thresholded_image}
+\end{definition}
+This will be used throughout the rest of the thesis. Note that given the
+thresholded image for every level, we are able to reconstruct the image
+as
+\begin{equation}
+ u(x) = \sup \left\{ s : u^s(x) = 1 \right\}.
+\end{equation}
+The thresholded image definition also allows us to write a non-negative
+image $u \geq 0$ as an integral over all the layers
+\begin{equation}
+ u(x) = \int_0^\infty u^s(x) \, ds.
+ \label{eq:positive_int}
+\end{equation}
+Note that \eqref{eq:positive_int} only holds for non-negative images,
+which complicates the following proof a little.
+
+\begin{figure}
+ \input{fig/eta_r}
+\end{figure}
+
+\begin{theorem}[Anisotropic coarea formula]
+ Given an image $u \in \BV(\Omega)$, the anisotropic total variation
+ can be written as an integral over all the levels
+ \begin{equation}
+ \TVA(u) = \int_{-\infty}^\infty \TVA(u^s) \, ds.
+ \label{eq:anisotropic_coarea}
+ \end{equation}
+ \label{thm:anisotropic_coarea}
+\end{theorem}
+For the proof we will avoid measure theory and follow a proof given
+in \cite{olsson2009extending}, but first we will present a necessary
+result from measure theory.
+\begin{theorem}[Lebesgue's Dominated Convergence theorem]
+ Let $\{ f_n \}$ be a sequence of real-valued measurable functions on
+ a space $S$ with measure $d\mu$ which converges almost
+ everywhere to a real-valued measurable function $f$. If there exists
+ an integrable function $g$ such that $\abs{f_n} \leq g$ for all $n$,
+ then $f$ is integrable and
+ \begin{equation}
+ \lim_{n \to \infty} \int_S f_n \, d\mu = \int_S f \, d\mu.
+ \end{equation}
+\end{theorem}
+For a proof and further background on measure theory and Lebesgue
+integration theory see for example \cite{bartle1995elements}.
+\begin{proof}[Proof of the anisotropic coarea formula.]
+ Assume that $u \in C^1(\Omega) \cap \BV(\Omega)$. The extension to
+ all functions $u \in \BV(\Omega)$ will not be considered here, but
+ for the case of regular total variation see \cite[Theorem
+ 5.3.3]{ziemer1989}.
+
+ %\paragraph{First we prove that $\TVA(u) \leq \int_{-\infty}^\infty
+ %\TVA(u^s) \, ds$.}
+ \paragraph{Proof of upper bound.}
+ Assume that $u \geq 0$ such that the integral representation in
+ \eqref{eq:positive_int} holds, then inserting
+ \eqref{eq:positive_int} into the extended
+ total variation definition in \eqref{eq:extended_tv} gives
+ \begin{equation}
+ \begin{aligned}
+ \TVA(u) &= \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega
+ \left( \int_0^\infty u^s ds \right) \diver \xi \, dx
+ = \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega \int_0^\infty u^s
+ \diver \xi \, ds \, dx \\
+ &\leq \int_0^\infty
+ \left( \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega u^s \diver
+ \xi \, dx \right) \, ds
+ = \int_0^\infty \TVA(u^s) \, ds.
+ \end{aligned}
+ \label{eq:tv_u_pos}
+ \end{equation}
+ For $u \leq 0$ we use that $\TVA(-v) = \TVA(v)$ and that $\TVA(c + v) =
+ \TVA(v)$ for any constant $c$. Note that $-u \geq 0$ and that its
+ thresholded image $(-u)^s$ will be exactly the opposite of $u^{-s}$,
+ that is $(-u)^s = 1 - u^{-s}$. This allows us to show that
+ \begin{equation}
+ \begin{aligned}
+ \TVA(u) &= \TVA(-u) \leq \int_0^\infty \TVA \big( (-u)^r
+ \big) \, dr
+ = \int_0^\infty \TVA(1 - u^{-r}) \, dr \\ &= \int_0^\infty
+ \TVA(u^{-r}) \, dr = \int_{-\infty}^0 \TVA(u^s) \, ds.
+ \end{aligned}
+ \label{eq:tv_u_neg}
+ \end{equation}
+ Following from the supremum definition of the anisotropic total
+ variation in \eqref{eq:extended_tv}, we obtain the inequality
+ \begin{equation}
+ \begin{aligned}
+ \TVA(u_1 + u_2) &= \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega (u_1 +
+ u_2) \diver \xi \, dx \\
+ &\leq \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega u_1
+ \diver \xi \, dx + \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega u_2
+ \diver \xi \, dx\\
+ &= \TVA(u_1) + \TVA(u_2).
+ \end{aligned}
+ \label{eq:tv_sum}
+ \end{equation}
+ Next, we write a general $u$ as a difference of two positive
+ functions $u = u_+ - u_-$ where $u_+ = \max\{u,0\}$ and $u_- =
+ -\min\{u,0\}$.
+ Inserting \eqref{eq:tv_u_pos} and
+ \eqref{eq:tv_u_neg} into \eqref{eq:tv_sum} we obtain
+ \begin{equation}
+ \begin{aligned}
+ \TVA(u) &\leq \TVA(u_-) + \TVA(u_+) = \TVA(-u_-) + \TVA(u_+) \\
+ &\leq \int_{-\infty}^0 \TVA\big((-u_-)^s\big) \, ds +
+ \int_0^\infty \TVA(u_+^s) \, ds \\
+ &= \int_{-\infty}^0 \TVA(u^s) \, ds + \int_0^\infty \TVA(u^s) \,
+ ds = \int_{-\infty}^\infty \TVA(u^s) \, ds.
+ \end{aligned}
+ \end{equation}
+ Note that $u_+$ and $u_-$ will not be differentiable everywhere, but
+ we did not use the differentiability of $u$ in this part of the
+ proof.
+
+ %\paragraph{Next we prove that $\TVA(u) \geq \int_{-\infty}^\infty
+ %\TVA(u^s) \, ds$.}
+ \paragraph{Proof of lower bound.}
+ Define the function
+ \begin{equation}
+ m(t) = \int_{\{ x \in \Omega : u(x) \leq t\}} \norm{\nabla u}_A
+ \, dx,
+ \label{eq:mdef}
+ \end{equation}
+ and note that $m(\infty) = \TVA(u)$ and $m(-\infty) = 0$. Since
+ $m(t)$ is non-decreasing with $t$, we can apply the existence
+ theorems of Lebesgue \cite[Thm.\ 17.12, 18.14]{hewstrom} to conclude
+ that $m\prime(t)$ exists almost everywhere and that the following
+ inequality holds:
+ \begin{equation}
+ \int_{-\infty}^\infty m\prime(t)\, dt \leq m(\infty) - m(-\infty) =
+ \TVA(u).
+ \label{eq:tva_geq_mder}
+ \end{equation}
+ Next, fix an $s \in \mathbb{R}$ and define the cut-off function
+ \begin{equation}
+ \begin{aligned}
+ \eta_r(t) = \begin{cases}
+ 0 & \text{if } t < s, \\
+ (t - s)/r & \text{if } s \leq t < s + r, \\
+ 1 & \text{if } t \geq s + r,
+ \end{cases}
+ & \quad
+ \eta_r\prime(t) = \begin{cases}
+ 0 & \text{if } t < s, \\
+ 1 & \text{if } s < t < s + r, \\
+ 0 & \text{if } t > s + r,
+ \end{cases}
+ \end{aligned}
+ \end{equation}
+ visualized in Figure~\ref{fig:eta_r} and \ref{fig:eta_r_diff}. By
+ composing the function $\eta_r$ with our image $u$ and using Green's
+ formula, for example from \cite[Corollary
+ 9.32]{grasmair2010anisotropic} we obtain
+ \begin{equation}
+ \int_\Omega - \eta_r(u) \diver \xi \, dx
+ = \int_\Omega \eta_r\prime(u) \nabla u\cdot \xi \, dx
+ = \frac{1}{r} \int_{\{ s < u < s + r \}} \nabla u\cdot \xi
+ \, dx,
+ \label{eq:eta_greens}
+ \end{equation}
+ for all vector fields $\xi \in C_c^\infty(\Omega, \mathbb{R}^2)$.
+ The measure of $\{ x : u(x) = \lambda \text{ and } \nabla u(x) \neq
+ 0\}$ is zero for all $\lambda$ following from \cite[Corollary \rom{1},
+ Section 3.1.2]{evans1991measure}. Thus we can ignore the sets $\{ u
+ = s\}$ and $\{ u = s + r \}$.
+ Assuming that $\norm{\xi}_A^* \leq 1$ we obtain from \eqref{eq:mdef}
+ and \eqref{eq:eta_greens} that
+ \begin{equation}
+ \begin{aligned}
+ \frac{m(s+r) - m(s)}{r}
+ &= \frac{1}{r} \int_{\{ s < u \leq s+r \}} \norm{\nabla u}_A
+ \, dx \\
+ &\geq \frac{1}{r} \int_{\{ s < u \leq s + r\}} \nabla u \cdot
+ \xi \, dx
+ = \int_\Omega -\eta_r(u) \diver \xi \, dx.
+ \end{aligned}
+ \label{eq:m_ineq_sr}
+ \end{equation}
+ As the limit when $r \to 0$ of the left-hand side exists almost
+ everywhere, suppose it exists at $s \in \mathbb{R}$. The integrand
+ on the right-hand side $-\eta_r(u) \diver \xi$ approaches $-u^s
+ \diver \xi$ pointwise almost everywhere. We apply
+ Lebesgue's dominated convergence theorem
+ using that $\abs{\eta_r(u) \diver \xi} \leq
+ \abs{u^s \diver \xi}$ and that $\xi \in C^\infty_c(\Omega,
+ \mathbb{R}^2)$ is bounded by the extreme value theorem. Thus
+ $\abs{u^s \diver \xi}$ is integrable and
+ \begin{equation}
+ \int_\Omega -\eta_r(u)
+ \diver \xi \, dx \to \int_\Omega -u^s \diver \xi \, dx
+ \end{equation}
+ From \eqref{eq:m_ineq_sr} we then obtain
+ \begin{equation}
+ m\prime(s) \geq - \int_\Omega u^s \diver \xi \, dx.
+ \end{equation}
+ As this holds for any
+ $\norm{\xi}_A^* \leq 1$, we get from the extended total variation
+ definition in \eqref{eq:extended_tv} that $m'(s) \geq \TVA(u_s)$
+ almost everywhere and conclude using \eqref{eq:tva_geq_mder} that
+ \begin{equation}
+ \TVA(u) \geq \int_{-\infty}^\infty m'(t) \, dt \geq
+ \int_{-\infty}^\infty \TVA(u^s) \, ds.
+ \end{equation}
+ Combining the upper and lower bounds just proved, we have equality.
+\end{proof}
+This coarea formula is our first step in transforming the anisotropic
+total variation into an easily discretizable expression.
+
+The anisotropic total variation of the thresholded images occurring in
+the anisotropic coarea formula is very much related to the size of the
+boundary of the level set, as the only variation in a characteristic
+function occurs at the boundary of the set. This is why we introduce
+the following definition of the anisotropic set perimeter.
+\begin{definition}[The anisotropic set perimeter]
+ Given an anisotropy tensor $A$ the anisotropic perimeter of a set
+ $U$ in $\Omega$ is defined as
+ \begin{equation}
+ \PerA(U;\Omega) = \TVA(\idfun_U).
+ \end{equation}
+\end{definition}
+\nomenclature{$\PerA(u;\Omega)$}{Anisotropic perimeter of set $U$ using
+anisotropy tensor $A$.}%
+The anisotropic set perimeter is not like the regular set perimeter and
+does not measure the length of the boundary of the set, but it can for
+sufficiently nice level sets be calculated in the following way
+\begin{equation}
+ \begin{aligned}
+ \PerA(\{ u > s \}; \Omega) &= \TVA(u^s) \\ &= \sup_{\norm{\xi}_A^* \leq 1}
+ \int_\Omega u^s \diver \xi \, dx \\
+ &= \sup_{\norm{\xi}_A^* \leq 1} \int_{\{ u > s \}} \diver \xi \,
+ dx \\
+ &= \sup_{\norm{\xi}_A^* \leq 1} \int_{\partial \{ u > s\} } \nu_s
+ \cdot \xi \, dt \\
+ &= \sup_{\norm{\eta} \leq 1} \int_{\partial \{ u > s\} } \nu_s
+ \cdot \Ahalf \eta \, dt \\
+ &= \int_{\partial \{ u > s\} } \sqrt{\nu_s A \nu_s} \, dt.
+ \end{aligned}
+ \label{eq:perimeter_calc}
+\end{equation}
+Here, $\nu_s$ is the unit exterior normal of the level set $\{ u > s \}$.
+Note that because of the compact support of $\xi$ in
+Definition~\ref{def:extended_tv}, the parts of the boundary of $U$ that
+overlap with the boundary of $\Omega$ will not be included in the
+perimeter.
+
+Exterior normals and perimeters of
+level sets of any function $u \in \BV(\Omega)$ will not be considered
+here, but can for the isotropic case be found in for example
+\cite[Section 5.4 and 5.5]{ziemer1989}.
+
+The anisotropic coarea formula allows us to transform the anisotropic
+total variation such that we are left with minimizing the following
+functional
+\begin{equation}
+ F(u) = \int_\Omega (u - f)^2 \, dx + \beta \, \int_{-\infty}^\infty
+ \PerA(\{u > \lambda\}; \Omega) \, d\lambda.
+ \label{eq:functional_per}
+\end{equation}
+The transformation is motivated by our upcoming Cauchy--Crofton
+integration formula, and discretization, where an approximation of the
+perimeter will be computed using a graph cut machinery.
+
+\section{Cauchy--Crofton formulas}
+
+\fixme{Rating: 7/10, the stuff after the proof could use some
+work.}
+
+\begin{figure}
+ \input{fig/line_param}
+\end{figure}
+
+In the fields of integral geometry and geometric measure theory there are
+a number of interesting integral formulas. Several of them fall in a
+category often referred to as \emph{Cauchy--Crofton style formulas}, and
+give ways to measure geometric objects using the set of all lines in the
+plane. The formulas presented here will give a way to measure the length
+of a curve by counting the times it intersects line in the set of all
+lines. The first formula will be for the isotropic case, and we will
+use this to prove the anisotropic formula following it.
+
+We write $\mathcal{L}$ for the set of all lines in the plane,
+and parametrize them as shown in Figure~\ref{fig:line_param}. A
+line is parametrized by the angle $\phi \in [0, 2\pi)$ of the normal going to the
+origin, and the distance $\rho \in [0, \infty)$ from origin to the line. Sometimes it is
+more convenient to consider a unit vector $\nu$ giving the direction of
+the line instead of the angle parameter $\phi$. We denote a line by
+$\ell_{\phi, \rho} = \ell_{\nu, \rho}$ where $\nu$ is a unit vector
+along the line, i.e.\ $\nu = (-\sin \phi, \cos \phi)^T$. By defining the
+measure on this set $d\mathcal{L} = \dpdr$ we are ready to introduce the
+Cauchy--Crofton formula. Note that the measure $d\mathcal{L}$ is
+invariant under rotations.
+\nomenclature{$\mathcal{L}$}{The set of all straight lines in the
+plane.}%
+\nomenclature{$\ell_{\phi, \rho}$}{A line given by the angle of the
+normal $\phi$ and the distance to origin $\rho$.}%
+\nomenclature{$\ell_{\nu, \rho}$}{A line given by a tangent vector $\nu$
+and the distance to origin $\rho$.}%
+\begin{theorem}[The Euclidean Cauchy--Crofton formula]
+ Given a differentiable curve $C$ in $\mathbb{R}^2$, the length of
+ this curve $\abs{C}$ is related to the set of lines $\mathcal{L}$ as
+ follows
+ \begin{equation}
+ \int_\mathcal{L} \#(\ell_{\phi, \rho} \cap C) \,
+ d\mathcal{L}(\ell_{\phi, \rho}) = 2 \abs{C},
+ \end{equation}
+ where $\#(\ell_{\phi, \rho} \cap C)$ is the number of times the line
+ $\ell_{\phi, \rho}$ intersects the curve $C$.
+ \label{thm:euclidean_cauchy_crofton}
+\end{theorem}
+\begin{proof}
+ See \cite[Theorem 3, Section 1-7.]{do1976differential}.
+\end{proof}
+If our space is equipped with a metric tensor $M(x)$ such that the inner
+product of two vectors $a$ and $b$ in a point $x$ is calculated as
+$\langle a, b\rangle_M = \langle a, M(x) b \rangle$, then the length of
+a curve $\gamma$ parametrized by some parameter $t$ becomes
+\nomenclature{$M(x)$}{A metric tensor.}%
+\nomenclature{$\abs{C}$}{The length of the curve $C$.}%
+\nomenclature{$\abs{C}_M$}{The length of the curve $C$ calculated using
+the metric tensor $M$.}%
+\begin{equation}
+ \abs{\gamma}_M = \int_\gamma \sqrt{\langle \dot{\gamma},
+ M\big(\gamma(t)\big) \, \dot{\gamma} \rangle} \, dt.
+ \label{eq:riemannian_length}
+\end{equation}
+We will now present and prove a Cauchy--Crofton formula in this case
+where our domain is equipped with a metric tensor in each point.
+This elegant formula is very useful when we later will discretize our
+perimeter calculation. The set of lines $\mathcal{L}$ is then discretized
+in a reasonable way, and the length of the curve $C$ can be approximated
+by a sum over all these lines.
+\begin{theorem}[The Riemannian Cauchy--Crofton formula]
+ Assume that our space $\Omega$ is equipped with a continuous
+ positive definite metric
+ tensor $M(x)$, whose eigenvalues are bounded by $0 < k \leq
+ \lambda_2 \leq \lambda_1 \leq K < \infty$ for all $x \in \Omega$.
+ The Cauchy--Crofton formula for a differentiable curve $C$ of finite
+ length then becomes
+ \begin{equation}
+ \abs{C}_M = \int_\mathcal{L} \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}}}
+ \, d\mathcal{L}(\ell_{\nu, \rho}).
+ \label{eq:riemannian_cauchy_crofton}
+ \end{equation}
+ \label{thm:riemannian_cauchy_crofton}
+\end{theorem}
+\begin{proof}[Proof of the Riemannian Cauchy--Crofton formula]
+ Assume first that our space is equipped with a constant metric
+ tensor $M$. The length of a curve in this space can be
+ calculated by transforming the curve and applying the Euclidean
+ Cauchy--Crofton formula
+ \begin{align}
+ \abs{C}_M &= \int_C \sqrt{\langle \dot{C}, M \, \dot{C} \rangle} \, dt
+ = \int_C \sqrt{\langle \Mhalf \dot{C}, \Mhalf
+ \dot{C} \rangle} = \abs{M^{\sfrac{1}{2}}C} \\
+ &= \int_\mathcal{L} \#(\ell_{\phi, \rho} \cap M^{\sfrac{1}{2}}C) \,
+ d\mathcal{L}(\ell_{\phi, \rho}) \\
+ &= \int_\mathcal{L} \#(M^{-\sfrac{1}{2}}\ell_{\phi, \rho} \cap C) \,
+ d\mathcal{L}(\ell_{\phi, \rho}) \\
+ &= \int_\mathcal{L} \#(m_{\phi, \rho} \cap C) \abs{J_M(\ell_{\phi, \rho})}\,
+ d\mathcal{L}(m_{\phi, \rho}).
+ \end{align}
+ Here $J_M(\ell_{\phi, \rho})$ is the Jacobian of the coordinate
+ transformation $F : \mathcal{L} \to \mathcal{L}$, which maps
+ $\ell_{\phi, \rho} \mapsto \Mhalf \ell_{\phi, \rho}$.
+ \nomenclature{$J_M(\ell_{\phi, \rho})$}{Jacobian of the coordinate
+ transformation induced by the metric tensor $M$.}
+
+ We will now compute the Jacobian $J_M(\ell_{\phi, \rho})$. As $M\in
+ \mathbb{R}^{2\times2}$ is symmetric, so is $M^{\sfrac{1}{2}}$, and
+ it admits a decomposition $\Mhalf = U\Sigma U^T$ where the
+ components correspond to the following coordinate transformations
+ \begin{align}
+ %U(\phi,\rho) &= (\phi + \xi, \rho), \\
+ %U^T(\phi,\rho) &= (\phi - \xi, \rho), \\
+ U(\ell_{\nu,\rho}) &= \ell_{\phi + \xi, \rho} = \ell_{U\nu, \rho} \\
+ U^T(\ell_{\nu,\rho}) &= \ell_{\phi - \xi, \rho} = \ell_{U^T\nu, \rho} \\
+ \Sigma &= \begin{pmatrix}
+ \sigma_1 & 0 \\
+ 0 & \sigma_2
+ \end{pmatrix}
+ = \begin{pmatrix}
+ \sqrt{\lambda_1} & 0 \\
+ 0 & \sqrt{\lambda_2}
+ \end{pmatrix}
+ \end{align}
+ As $U$ and $U^T$ correspond to rotations and our measure
+ $\mathcal{L}$ is invariant under rigid motions, $U$ and $U^T$ do not
+ have direct contributions to the Jacobian. They do however affect the
+ input angle of the operator $\Lambda$ such that $J_M(\ell_{\phi,
+ \rho}) = J_{\Sigma^2}(U^T \ell_{\phi, \rho})$. Thus we will now
+ compute $J_{\Sigma^2}(\ell_{\phi, \rho})$. Given a line
+ \begin{equation}
+ \ell_{\phi, \rho} = \begin{pmatrix}
+ \rho\cdot \cos \phi \\
+ \rho\cdot \sin \phi
+ \end{pmatrix}
+ + \mathbb{R} \begin{pmatrix}
+ -\sin \phi \\
+ \cos \phi
+ \end{pmatrix},
+ \end{equation}
+ the operator $\Sigma$ transforms it into
+ \begin{equation}
+ \Sigma \ell_{\phi, \rho} = \begin{pmatrix}
+ \sigma_1 \rho\cdot \cos \phi \\
+ \sigma_2 \rho\cdot \sin \phi
+ \end{pmatrix}
+ + \mathbb{R} \begin{pmatrix}
+ -\sigma_1 \sin \phi \\
+ \sigma_2 \cos \phi
+ \end{pmatrix},
+ \end{equation}
+ which equals the line $\ell_{\theta, \eta}$ with
+ \begin{align}
+ \theta &= \arctan \left( \frac{\sigma_1}{\sigma_2} \tan \phi
+ \right) \label{eq:theta}\\
+ \eta &= \left\langle \begin{pmatrix}
+ \sigma_1 \rho\cdot \cos \phi \\
+ \sigma_2 \rho\cdot \sin \phi
+ \end{pmatrix}, \begin{pmatrix}
+ \cos \theta \\
+ \sin \theta
+ \end{pmatrix}
+ \right\rangle = \sigma_1 \rho \cdot \cos \phi \cdot \cos \theta
+ + \sigma_2 \rho \cdot \sin \phi \cdot \sin \theta.
+ \end{align}
+ As $\partial_\rho \theta = 0$, the Jacobian becomes
+ $\abs{J_{\Sigma^2}(\ell_{\phi, \rho})} = \partial_\phi\theta \cdot
+ \partial_\rho \eta$. Differentiation yields
+ \begin{align}
+ \partial_\phi \theta &= \frac{\frac{\sigma_1}{\sigma_2}
+ \sec^2 \phi}{1 + \frac{\sigma_1^2}{\sigma_2^2}
+ \tan^2 \phi}
+ = \frac{\sigma_1 \sigma_2}{\sigma_1^2 \sin^2 \phi + \sigma_2^2
+ \cos^2 \phi}, \\
+ \partial_\rho \eta &= \sigma_1 \cos \phi \cdot \cos \theta +
+ \sigma_2 \sin \phi \cdot \sin \theta.
+ \end{align}
+ In the expression for $\partial_\rho \eta$ we insert $\theta$ from
+ \eqref{eq:theta} and use that $\sin
+ (\arctan( x )) = x / \sqrt{1 + x^2}$ and that $\cos ( \arctan ( x ))
+ = 1 / \sqrt{1 + x^2}$ to obtain
+ \begin{equation}
+ \partial_\rho \eta = \frac{\sigma_1 \cos \phi + \sigma_2 \sin
+ \phi \, \frac{\sigma_1}{\sigma_2} \tan \phi}
+ {\sqrt{1 + \frac{\sigma_1^2}{\sigma_2^2}
+ \tan^2 \phi}}
+ = \frac{\sigma_1 \sigma_2}{\sqrt{\sigma_1^2 \sin^2 \phi +
+ \sigma_2^2 \cos^2 \phi}}.
+ \end{equation}
+ If $\nu = (\nu_x, \nu_y)^T$ is a unit vector along the line
+ $\ell_{\phi, \rho} = \ell_{\nu, \rho}$ then
+ \begin{equation}
+ \abs{J_{\Sigma^2}(\ell_{\nu, \rho})}
+ = \frac{\sigma_1^2 \sigma_2^2}
+ {\left(\sigma_1^2 \sin^2 \phi + \sigma_2^2 \cos^2 \phi
+ \right)^{\sfrac{3}{2}}}
+ = \frac{\sigma_1^2 \sigma_2^2
+ }{\left(\sigma_1^2 \nu_x^2 + \sigma_2^2
+ \nu_y^2\right)^{\sfrac{3}{2}}}
+ = \frac{\det \Sigma^2}{(\nu^T \cdot \Sigma^2 \cdot
+ \nu)^{\sfrac{3}{2}}}.
+ \end{equation}
+ We are interested in the Jacobian of the whole transformation
+ $J_{\Sigma^2}(U^T \ell_{\theta, \rho})$, so
+ all that is left to do is insert $U^T \ell_{\nu, \rho}$ to obtain
+ \begin{equation}
+ \abs{J_M(\ell_{\nu, \rho})} = \abs{J_{\Sigma^2}(U^T \ell_{\nu, \rho})}
+ = \frac{\det M}{\left(\nu^T U \cdot \Sigma^2 \cdot U^T
+ \nu\right)^{\sfrac{3}{2}}}
+ = \frac{\det M}{\left(\nu^T \cdot M \cdot \nu\right)^{\sfrac{3}{2}}}
+ \end{equation}
+ We have now proved that for a constant metric tensor $M$, the length
+ of the differentiable curve $C$ with regards to this tensor can be
+ calculated as
+ \begin{equation}
+ \abs{C}_M = \int_C \sqrt{\langle \dot{C}, M \dot{C} \rangle} \,
+ dt
+ = \int_\mathcal{L} \# (\ell_{\nu, \rho} \cap C)
+ \frac{\det M}{\left(\nu^T \cdot M \cdot \nu
+ \right)^{\sfrac{3}{2}}} \, d\mathcal{L} (\ell_{\nu, \rho}).
+ \label{eq:riemannian_const_m}
+ \end{equation}
+ We now argue that the similar formula in
+ \eqref{eq:riemannian_cauchy_crofton} holds for a non-constant but
+ continuous metric tensor $M(x)$. By partitioning the domain into
+ disjoint sets $U_i$ such that $\Omega = \cup_i U_i$, we make a
+ piecewise constant approximation $M_\pi(x)$ such that if $x \in U_i$
+ then $M_\pi(x) = M(x_i)$ for some fixed $x_i \in U_i$. We then
+ approximate \eqref{eq:riemannian_const_m} by
+ \begin{equation}
+ \abs{C}_{M_\pi} = \sum_i \int_\mathcal{L} \#(\ell_{\nu, \rho} \cap
+ C \cap U_i) \, w_i(\nu) \, d\mathcal{L}(\ell_{\nu, \rho})
+ \label{eq:mpi_approx}
+ \end{equation}
+ where $w_i$ is the weight-function used in the set $U_i$, that is,
+ \begin{equation}
+ w_i(\nu) = \frac{\det M(x_i)}{\left(\nu^T \cdot M(x_i) \cdot \nu
+ \right)^{\sfrac{3}{2}}}.
+ \end{equation}
+ We further simplify the approximation by introducing the global
+ weight-function $w_\pi(\nu, x)$ which is equal to
+ $w_i(\nu)$ when $x \in U_i$. It can be written as
+ \begin{equation}
+ w_\pi(\nu, x) = \frac{\det M_\pi(x)}{\left(\nu^T \cdot M_\pi(x)
+ \cdot \nu \right)^{\sfrac{3}{2}}}.
+ \label{eq:wpi_def}
+ \end{equation}
+ Using this weight in \eqref{eq:mpi_approx} we can get rid of the sum
+ over the partition $i$ and form a sum of all intersection point of
+ $C$ and the line $\ell_{\nu, \rho}$ currently being integrated over.
+ The approximation becomes
+ \begin{equation}
+ \begin{aligned}
+ \abs{C}_{M_\pi} &=
+ \sum_i \int_\mathcal{L} \sum_{x \in \ell_{\nu, \rho} \cap C \cap
+ U_i} w_\pi(\nu, x) \, d\mathcal{L}(\ell_{\nu, \rho}) \\
+ &= \int_\mathcal{L} \sum_{x \in \ell_{\nu, \rho} \cap C}
+ w_\pi(\nu, x) \, d\mathcal{L}(\ell_{\nu, \rho}).
+ \end{aligned}
+ \label{eq:riemannian_segments}
+ \end{equation}
+ Now it only remains
+ to show that the left- and right-hand side of
+ \eqref{eq:riemannian_segments} converges to the left- and right-hand
+ side of \eqref{eq:riemannian_cauchy_crofton}.
+
+ As our partition $\pi$ is refined, the weight $w_\pi(x)$ converges
+ pointwise to the continuously varying weight
+ \begin{equation}
+ w(\nu, x) = \frac{\det M(x)}{\left( \nu^T \cdot M(x) \cdot \nu
+ \right)^{\sfrac{3}{2}}}
+ \end{equation}
+ found in \eqref{eq:riemannian_cauchy_crofton}.
+
+ Recall from \eqref{eq:riemannian_length} that the left-hand side is
+ calculated as
+ \begin{equation}
+ \abs{C}_{M_\pi}
+ = \int_C \abs{\dot{C}(t)}_{M_\pi} \, dt
+ = \int_C \sqrt{\dot{C}(t)^T M_\pi\big( C(t) \big) \dot{C}(t)}
+ \, dt.
+ \label{eq:piecewise_curve_length}
+ \end{equation}
+ We know that $M_\pi(x)$ converges pointwise to $M(x)$, and thus
+ $\labs{\dot{C}(t)}_{M_\pi}$ converges pointwise to
+ $\labs{\dot{C}(t)}_M$. We have assumed bounds on the eigenvalues
+ of $M(x)$ such that, according to the Rayleigh principle
+ \begin{equation}
+ K \geq \lambda_1 = \max_\xi \frac{\xi^T M_\pi(x) \xi}{\xi^T \xi
+ }
+ \end{equation}
+ and therefore we have the bound
+ \begin{equation}
+ \xi^T M_\pi(x) \xi \leq K \, \norm{\xi}^2, \quad \forall \xi.
+ \end{equation}
+ Thus the integrand of \eqref{eq:piecewise_curve_length} is bounded
+ by $g(t) = \sqrt{K \cdot \dot{C}(t)^T \dot{C}(t)}$. We know that
+ $g(t)$ is integrable as its integral is exactly $\sqrt{K} \abs{C}$
+ and we have assumed that the curve is of finite length. This means
+ we can apply Lebesgue's dominated convergence theorem to see that
+ $\labs{C}_{M_\pi} \to \labs{C}_M$.
+
+ We apply the same theorem to show that the right-hand side of
+ \eqref{eq:riemannian_segments} converges. Recall the definition of
+ $w_\pi$ in \eqref{eq:wpi_def}. The numerator is equal to $\sigma_1^2
+ \sigma_2^2 = \lambda_1 \lambda_2$ and is by assumption bounded from
+ above by $K^2$.
+
+ Next we need to bound $\nu^T M_\pi(x) \nu$ away from zero. According
+ to the Rayleigh principle
+ \begin{equation}
+ \lambda_2 = \min_{\norm{\xi} = 1} \xi^T M_\pi(x) \xi
+ \end{equation}
+ and thus $\nu^T M_\pi(x) \nu \geq \lambda_2 \geq k$. The weight
+ function $w_\pi$ is then bounded such that
+ \begin{equation}
+ \sum_{x \in \ell_{\nu, \rho} \cap C} w_\pi(\nu, x)
+ \leq
+ \sum_{x \in \ell_{\nu, \rho} \cap C}
+ \frac{K^2}{k^{\sfrac{3}{2}}}
+ =
+ \frac{K^2}{k^{\sfrac{3}{2}}} \cdot \# (\ell_{\nu, \rho} \cap C)
+ =: g(\ell_{\nu, \rho}).
+ \end{equation}
+ This is integrable following from the Euclidean Cauchy--Crofton
+ formula in Theorem~\ref{thm:euclidean_cauchy_crofton} and the fact
+ that we assumed $C$ to be of finite length:
+ \begin{equation}
+ \int_\mathcal{L} g(\ell_{\nu, \rho}) \, d\mathcal{L}(\ell_{\nu,
+ \rho}) = \frac{K^2}{k^{\sfrac{3}{2}}} \abs{C} < \infty.
+ \end{equation}
+ Thus we can apply the dominated convergence theorem again and
+ conclude that
+ \begin{equation}
+ \int_\mathcal{L} \sum_{x \in \ell_{\nu, \rho} \cap C}
+ w_\pi(\nu, x) \, d\mathcal{L}(\ell_{\nu, \rho})
+ \to
+ \int_\mathcal{L} \sum_{x \in \ell_{\nu, \rho} \cap C}
+ w(\nu, x) \, d\mathcal{L}(\ell_{\nu, \rho})
+ \end{equation}
+ which---as both sides of the equality in
+ \eqref{eq:riemannian_segments} have been shown to converge---leaves
+ us with what we wanted to prove
+ \begin{equation}
+ \abs{C}_M =
+ \int_\mathcal{L} \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}}} \,
+ d\mathcal{L}(\ell_{\nu, \rho}).
+ \end{equation}
+\end{proof}
+
+With the anisotropic coarea formula in
+Theorem~\ref{thm:anisotropic_coarea} we have a way to calculate the
+anisotropic total variation by integrating the anisotropic perimeter of
+each level set of the image as introduced in \eqref{eq:functional_per}.
+In the Euclidean case, which here would amount to setting the anisotropy
+tensor $A$ equal to the identity matrix $I$, the perimeter coincides
+nicely with the length of the boundary curve, assuming some regularity
+for the boundary.
+
+In the general case we need
+to be more careful. As can be seen in \eqref{eq:perimeter_calc}, the
+anisotropic perimeter is calculated by integrating the norm of the
+\emph{normal vector} around the boundary, while the anisotropic curve
+length in \eqref{eq:riemannian_length} is the integral of the norm of
+the \emph{tangent vector} of the curve. Thus a 90\textdegree{} rotation
+is needed.
+
+If $P$ is a 90\textdegree{} rotation matrix we have
+\begin{equation}
+ \begin{aligned}
+ \PerA(U; \Omega) &= \int_{\partial U} \sqrt{\langle
+ \nu_{\partial U}, A(x) \nu_{\partial U} \rangle} \, dt \\
+ &= \int_{\partial U} \sqrt{\langle P \nu_{\partial U}, P A(x)
+ P^T \, P \nu_{\partial U} \rangle } \, dt.
+ %&\stackrel{\mathclap{\gamma \hskip 0.07em = \hskip 0.07em
+ %\partial U}}{=}
+ %\int_{\gamma} \sqrt{ \langle \dot{\gamma}, M(x) \,
+ %\dot{\gamma} \rangle } \, dt.
+ \end{aligned}
+ \label{eq:per_to_length1}
+\end{equation}
+We simplify the equation by defining the metric tensor $M(x) = P A(x)
+P^T$ and letting $\gamma = \partial U \cap \Omega$ be an arclength parametrization of
+the boundary of $U$ that does not overlap with the boundary of $\Omega$
+\begin{equation}
+ \PerA(U; \Omega) = \int_\gamma \sqrt{\langle \dot{\gamma}, M(x) \,
+ \dot{\gamma} \rangle } \, dt.
+ \label{eq:per_to_length2}
+\end{equation}
+
+Now we make sure that all the assumptions of the Riemannian
+Cauchy--Crofton formula in Theorem~\ref{thm:riemannian_cauchy_crofton}
+are fulfilled so that it can be applied to the curve length integral we
+have constructed in \eqref{eq:per_to_length2}.
+
+The structure tensor is constructed as described in
+Section~\ref{sec:anisotropy_tensor}
+\begin{equation}
+ S_\rho(x) = \left(K_{\rho} * \left( \nabla f_{\sigma}
+ \otimes \nabla f_{\sigma} \right)\right)(x),
+\end{equation}
+Because of the convolutions with the Gaussian function, this is a smooth
+map from $\bar{\Omega}$ to $\mathbb{R}^{2 \times 2}$. As we
+can see in \eqref{eq:s_eigenvalues}, the eigenvalues
+depend continuously on the coefficients of the
+elements in
+the structure tensor $S_\rho(x)$. The extreme value theorem states that
+a continuous real-valued function on a nonempty compact
+space is bounded above. Thus the eigenvalues $s_1 \geq s_2$ of $S_\rho(x)$
+are bounded above and by the construction in
+\eqref{eq:sigma_construction}, there exists uniform bound $k$ such that
+the smallest eigenvalue of the anisotropy tensor $A(x)$ is bounded away
+from zero, as
+\begin{equation}
+ \lambda_2 = \left(1
+ + \frac{(s_1 - s_2)^2}{\omega^2}\right)^{-1} \geq \left(1 +
+ \frac{s_1^2}{\omega^2}\right)^{-1} \geq k > 0.
+\end{equation}
+Hence our metric tensor $M(x) = P A(x) P^T$ is continuous and positive
+definite and thus the curve length calculation in
+\eqref{eq:per_to_length2} fulfills all the assumptions of the
+anisotropic Cauchy--Crofton formula in
+Theorem~\ref{thm:riemannian_cauchy_crofton}.
+Thus we can apply the formula to calculate the
+perimeter in \eqref{eq:per_to_length2} as
+\begin{equation}
+ \PerA(U; \Omega) = \int_\mathcal{L} \sum_{x \in
+ \ell_{\nu, \rho} \cap \gamma } \frac{\det M(x)}
+ {2 \left( \nu^T \cdot M(x) \cdot \nu \right)^{\sfrac{3}{2}}}
+ \, d\mathcal{L}(\ell_{\nu, \rho}) \, ds,
+\end{equation}
+where $\gamma = \partial U \cap \Omega$.
+Note that $P$ does not affect the determinant, i.e.\
+$\det A = \det PAP^T = \det M$, and from our decomposition in
+\eqref{eq:sigma_construction} we see that the transformation $PAP^T \to
+M$ actually amounts to switching the two eigenvalues $\sigma_1$ and
+$\sigma_2$ in $\Sigma$.
+
+This concludes the treatment of the continuous problem. We have seen how
+the anisotropic coarea formula in Theorem~\ref{thm:anisotropic_coarea}
+allows us to calculate the anisotropic total variation as an integral of
+the perimeter of all the level sets. Through the anisotropic
+Cauchy--Crofton formula in Theorem~\ref{thm:riemannian_cauchy_crofton}
+these perimeters are calculated by an integral over the set of all
+lines. We are then left with the functional
+\begin{equation}
+ F(u) = \int_\Omega (u - f)^2 + \beta \, \TVA(u),
+ \label{eq:final_cont_rep}
+\end{equation}
+where
+\begin{equation}
+ \TVA(u) = \int_{-\infty}^\infty \int_{\mathcal{L}} \sum_{x \in
+ \ell_{\nu, \rho} \cap \gamma_s} \frac{\det M(x)}
+ {2 \left( \nu^T \cdot M(x) \cdot \nu \right)^{\sfrac{3}{2}}}
+ \, d\mathcal{L}(\ell_{\nu, \rho}) \, ds,
+\end{equation}
+and $\gamma_s = \partial \{u > s\} \cap \Omega$.
+Within the restrictions that these theorems put on the tensor $M(x)$, we
+have chosen a construction where one eigenvalue is always 1, while the
+other varies from 1 in smooth areas towards 0 around edges, with the
+corresponding eigenvector perpendicular to the edge. \fixme{meh}
+
+As mentioned previously, these transformations are motivated by the
+discretization that follows in the next chapter.
+
--- /dev/null
+\chapter{Discrete formulation}
+\label{chap:discrete}
+
+The whole transformation from the initial functional in
+\eqref{eq:first_anisotropic_functional} through the anisotropic coarea
+formula and the Cauchy--Crofton formula was motivated by the
+discrete formulation which will be described here. After discretizing
+the functional, we will se how a graph cut approach can be used to find
+a global minimizer in polynomial time.
+
+\section{Discretization}
+
+\fixme{Rating: 4/10}
+
+Assume that our discrete images are given on a uniform grid $\mathcal{G}$,
+where each grid point is called a pixel. The image is a function giving
+each pixel a value in the set of levels $\mathcal{P} = \{0, \hdots,
+L-1\}$. This is a reasonable assumption for digital grayscale
+images.
+
+Thus, when discretizing the functional in \eqref{eq:final_cont_rep}, we
+have to consider that our images now have both discrete domain and
+co-domain.
+
+The integrals in \eqref{eq:final_cont_rep} will be approximated by
+discrete sums. First the fidelity term is discretized without too much
+trouble, while with the regularization term, there is more choice as to
+how to discretize the set of lines $\mathcal{L}$. In the end we will
+verify that our discretization is consistent with the continuous
+functional.
+
+\subsection{Fidelity term}
+
+\fixme{Rating: 7/10}
+
+Since it is not affected by our introduction of the anisotropy tensor,
+the fidelity term can be discretized as in my project work
+\cite{project}. For some pixel position $x \in \mathcal{G}$ and some
+level value $k \in
+\mathcal{P}$, we define the following function
+\begin{equation}
+ N_x(k) = \abs{k - f_x}^2
+\end{equation}
+which is the value of the fidelity term if we were to give $u_x$ a value
+of $k$. This allows us write
+\begin{equation}
+ \int_\Omega \abs{u - f}^2\, dx \approx \sum_{x \in \mathcal{G}}
+ \abs{u_x - f_x}^2 \Delta x = \sum_{x \in \mathcal{G}} N_x(u_x)
+ \Delta x.
+ \label{eq:fidelity_approx_1}
+\end{equation}
+The reason we introduce the function $N_x(k)$ is that we want to apply
+the following decomposition formula, which holds for any function $F(k)$
+taking values $k \in \mathcal{P}$:
+\begin{equation}
+ \begin{aligned}
+ F(k) &= \sum_{\lambda = 0}^{k-1} \big( F(\lambda + 1) -
+ F(\lambda) \big) + F(0) \\
+ &= \sum_{\lambda = 0}^{L-2} \big( F(\lambda + 1) - F(\lambda) \big)
+ I( \lambda < k) + F(0),
+ \end{aligned}
+\end{equation}
+where $I(x)$ is the indicator function that takes the value 1 if $x$ is
+true, and 0 if $x$ is false. Since $I(\lambda < u_x) = u_x^\lambda$ we
+rewrite \eqref{eq:fidelity_approx_1} and obtain
+\begin{equation}
+ \sum_{x \in \mathcal{G}} \abs{u_x - f_x}^2 =
+ \sum_{x \in \mathcal{G}} N_x(u_x) = \sum_{\lambda = 0}^{L-2} \sum_{x
+ \in \mathcal{G}}
+ \big( N_x(\lambda + 1) - N_x(\lambda) \big) \, u_x^\lambda + N_x(0).
+\end{equation}
+As our domain is discretized uniformly, we drop the constant
+$\Delta x$, and absorb it into our parameter $\beta$ of
+\eqref{eq:final_cont_rep}. Note that since
+our image takes values in $\mathcal{P} = \{0, \hdots, L-1\}$, the
+thresholded image $u^{L-1}$ is equal to zero everywhere.
+
+\subsection{Regularization term}
+\label{sec:disc_regularization}
+
+\fixme{Rating: 6/10}
+
+Discretizing the regularization term is more challenging. We introduce
+the discrete levels to get
+\begin{equation}
+ \int_{-\infty}^\infty \PerA( \{ u > \lambda \}; \Omega) \, d\lambda
+ \approx \sum_{\lambda = 0}^{L-2} \PerA( \{ u > \lambda \}; \Omega)
+ \, \Delta \lambda.
+ \label{eq:per_approx1}
+\end{equation}
+As with the $\Delta x$ difference, we can absorb the $\Delta \lambda$
+difference into the $\beta$ parameter of
+\eqref{eq:final_cont_rep}. 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
+has zero perimeter.
+
+\subsubsection{Discrete Riemannian Cauchy--Crofton formula}
+
+\begin{figure}
+ \input{fig/line_disc}
+\end{figure}
+By approximating the integral
+Theorem~\ref{thm:riemannian_cauchy_crofton} by a discrete sum we obtain
+the approximation
+\begin{equation}
+ \begin{aligned}
+ \abs{C}_M &= \int_\mathcal{L} \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}}} \,
+ d\mathcal{L}(\ell_{\nu, \rho}) \\
+ &\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\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 \,
+ \Delta\nu.
+ \label{eq:cauchy_crofton_approx1}
+ \end{aligned}
+\end{equation}
+The set of lines $\mathcal{L}$ has been discretized into the set
+$\mathcal{L}_D$. Note that we are approximating the length of the
+\emph{differentiable} curve $C$ in $\Omega$. Being a difference in the
+$\rho$ parameter of our line discretization in
+Figure~\ref{fig:line_param}, the difference $\Delta \rho$ represents the
+distance from one line to the next in a line family as shown in
+Figure~\ref{fig:line_family}, and thus depends on the angle $\phi$
+considered. The difference $\Delta \phi$ is taken to
+be the average of the distance to the two neighboring line families as
+shown in Figure~\ref{fig:line_neigh}, and thus also depends on $\phi$.
+
+In the discrete setting our domain $\Omega$ is discretized as a uniform
+grid $\mathcal{G}$. Our image is then reduced to a function $u :
+\mathcal{G} \to \mathcal{P}$.
+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}.
+
+The choice of our discrete set of lines $\mathcal{L}_D$ is important, as
+it will decide the accuracy of our approximation.
+We need some sensible restrictions on
+the set $\mathcal{L}_D$ to simplify the further discussion. All lines
+intersect at least two grid points, and from the periodicity of our grid
+they thus intersect an infinite number of grid points. This puts some
+restrictions on the angles we can choose. For each angle, we include all
+possible lines of that family, meaning there are no grid points without
+a line of that family intersecting it.
+
+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~\ref{fig:line_family} shows all lines of a given family, i.e.\
+lines having the same angle parameter $\phi$.
+
+Thus not only have we discretized the set of lines, but each line is
+made up of edges going from one grid point to the next. We will
+denote such an edge by $e$ or $e_{ab}$ when its endpoints are $a, b \in
+\mathcal{G}$. Thus we rewrite the discretization of
+\eqref{eq:cauchy_crofton_approx1}, and sum over all the edges in the
+discretization $\mathcal{L}_D$ to obtain
+\begin{equation}
+ \abs{C}_M \approx \sum_{e} \sum_{x \in e \cap C} \, \frac{\det M(x)
+ \norm{e}^3}{2 \left(e^T \cdot M(x) \cdot e\right)^{\sfrac{3}{2}}} \,
+ \Delta\phi \, \Delta\rho.
+ \label{eq:cauchy_crofton_approx2}
+\end{equation}
+This is beginning to look like something we can calculate. One
+difficulty is finding the intersections $e \cap C$. The exact
+calculations of these points will not fit into our graph cut framework
+later, and thus for an edge $e$ we will only consider the question of
+``did $e$ cross $C$ or not?'' This amounts to checking whether the
+terminals of $e$ lie on different sides 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}, be wrong when we have more.
+
+\begin{figure}
+ \input{fig/curve_edge}
+\end{figure}
+
+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 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}
+ M(x) \approx M(e_{ab}) = \frac{M(a) + M(b)}{2},
+ \label{eq:tensor_approx}
+\end{equation}
+the component-wise average of the tensors in the two end points of the
+edge. Recall that we have already done some spatial smoothing of the
+structure tensor in \eqref{eq:s_def} corresponding to the
+\emph{integration scale} $\rho$, and thus we expect the tensors $M(a)$
+and $M(b)$ to be similar for edges $e$ of reasonably short length.
+
+We also remark that using the Rayleigh principle, it is easy to conclude
+that the eigenvalues of the tensor approximation $M(e_{ab})$ are bounded
+below and above by the smallest and largest eigenvalues of $M(a)$ and
+$M(b)$.
+
+\begin{figure}
+ \input{fig/area_proof}
+\end{figure}
+
+We have now almost arrived at our final curve length approximation, but
+we need a way to calculate the inter-line distance $\Delta \rho$ which
+will be provided by the following lemma.
+\begin{lemma}
+ For each family of lines given by an angle parameter $\phi$ in the
+ uniform 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 a line $\ell$ intersecting the point $(p, q) \in
+ \mathbb{Z}^2$ in the grid.
+ The distance $\Delta \rho$ from this line $\ell$ to the neighboring
+ lines can then be calculated as a minimum over the distance to all
+ other grid points.
+
+ The lines are split into edges $e = (\delta s, \delta t)^T$ where
+ $s, t \in \mathbb{Z}$ are coprime such that $e$ does not intersect
+ any other grid points than its two endpoints.
+
+ We then calculate the minimal distance to a grid point not on the
+ line $\ell$ as
+ \begin{equation}
+ \begin{aligned}
+ \Delta \rho &= \min_{(p\prime, q\prime) \in \mathcal{G}
+ \setminus \ell} \left\{
+ \left\langle \delta [p - p\prime, q - q\prime],
+ \frac{e^\perp}{\norm{e^\perp}} \right\rangle \right\} \\
+ &= \min_{(p\prime, q\prime) \in \mathcal{G} \setminus \ell}
+ \left\{\delta^2 \cdot \frac{t(p-p\prime) - s(q -
+ q\prime)}{\norm{e}} \right\}.
+ \end{aligned}
+ \end{equation}
+ Since $s$ and $t$ are coprime, there exists $a, b \in \mathbb{Z}$
+ such that $at - bs = 1$, and since the $\Delta \rho$ cannot be zero,
+ we obtain
+ \begin{equation}
+ \Delta \rho = \frac{\delta^2}{\norm{e}}.
+ \end{equation}
+\end{proof}
+Inserting $\Delta \rho = \delta^2 / \norm{e}$ and the tensor
+approximation of \eqref{eq:tensor_approx} into the curve length
+approximation of \eqref{eq:cauchy_crofton_approx2} 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}}},
+ \label{eq:cauchy_crofton_approx3}
+\end{equation}
+where the sum is over all edges crossing the curve.
+
+The curve length we initially wanted to calculate was the perimeter
+$\PerA(\{u > \lambda\}; \Omega)$ in \eqref{eq:per_approx1}. To find
+edges that crosses the curve, we identify the edges that have one
+terminal inside the level set, and the other outside.
+Thus we rewrite the sum over $e \cap C$ such that
+\begin{equation}
+ \PerA(\{u > \lambda\}; \Omega) \approx \sum_{e_{ab}}
+ \abs{u^\lambda_a - u^\lambda_b} \frac{\det M(e_{ab})
+ \norm{e_{ab}}^2 \, \delta^2 \, \Delta\phi}{2 \left(e_{ab}^T
+ \cdot M(e_{ab}) \cdot e_{ab}\right)^{\sfrac{3}{2}}}.
+ \label{eq:per_approx2}
+\end{equation}
+The absolute value $\abs{u^\lambda_a - u^\lambda_b}$ is 1 if one of
+$a$ and $b$ lie inside the level set and the other lies outside, and
+0 otherwise.
+
+In other words the absolute value is one if $e_{ab}$ crosses the
+perimeter of $\{ u > \lambda\}$ an odd number of times, and zero
+otherwise.
+
+Thus we have arrived at our final discretization, which takes the form
+\begin{gather}
+ F(u) = \sum_\lambda \sum_x F_x^\lambda(u_x^\lambda) + \beta \sum_\lambda
+ \sum_{(x, y)} F_{x,y}^\lambda(u_x^\lambda, u_y^\lambda) =:
+ F^\lambda(u^\lambda), \\
+ \begin{aligned}
+ F_x^\lambda(u_x^\lambda) &= \big(N_x(\lambda + 1) -
+ N_x(\lambda)\big) \cdot u_x^\lambda, \\
+ F_{x,y}^\lambda(u_x^\lambda, u_y^\lambda)
+ &= \abs{u_x^\lambda - u_y^\lambda} \frac{\det M(e_{xy})
+ \norm{e_{xy}}^2 \delta^2 \Delta \phi}{2 \left( e_{xy}^T \cdot
+ M(e_{xy}) \cdot e_{xy} \right)^{\sfrac{3}{2}}}.
+ \end{aligned}
+ \label{eq:final_discretization}
+\end{gather}
+Recall that $N_x(\lambda) = \abs{ \lambda - f_x }^2$.
+
+If we minimize $F_\lambda$ to obtain $u^\lambda$ for each level
+separately, it is obvious that we will also minimize the sum over all
+$F_\lambda$. However, it is not guaranteed that the obtained thresholded
+images $u^\lambda$ can be combined to make an output image $u$. They
+were defined as $u^\lambda = \idfun_{u > \lambda}$, so we need them to
+be monotonically decreasing \fixme{(?)} in increasing level values, i.e.\
+\begin{equation}
+ u_x^\lambda \geq u_x^\mu, \quad \forall \lambda \leq \mu, \quad
+ \forall x \in \mathcal{G}.
+ \label{eq:decreasing_sets}
+\end{equation}
+Later we will present a graph cut algorithm that find the thresholded
+images minimizing each level, \emph{while guaranteeing that they meet
+this requirement.}
+
+\fixme{maybe with a $w_{xy}$ definition}
+
+\subsubsection{Consistency}
+
+\fixme{Rating: 6/10}
+
+Consistency relates to whether a solution to the continuous problem
+fits in the discretized equation, in other words, whether 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. The sum is a midpoint rule
+approximation of the integral. As the grid is refined and $\delta \to 0$
+the sum will converge to 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 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}.
+
+If we consider the edges $e$ of each family separately, the curve length
+approximation in \eqref{eq:cauchy_crofton_approx3} can be written
+\begin{equation}
+ \begin{aligned}
+ \abs{C}_M &=
+ \int_\nu \int_\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}}} \, d\rho \, d\nu \\
+ &\approx \sum_\nu \sum_\rho \sum_{e_{\nu, \rho} \cap C}
+ \, \frac{\det M(e_{\nu, \rho}) \norm{e_{\nu, \rho}}^3}{2\left(
+ e_{\nu, \rho}^T \cdot M(e_{\nu, \rho}) \cdot e_{\nu, \rho}
+ \right)^{\sfrac{3}{2}}} \, \Delta \rho \, \Delta \nu.
+ \end{aligned}
+\end{equation}
+
+As described in the construction of this formula, there are four main
+approximations used. Firstly there is the fact that we do not consider
+the actual intersection points, but only whether an edge crosses the
+curve or not. Secondly we have the tensor which is averaged as in
+\eqref{eq:tensor_approx}. And then we have the discretizations of our
+two line parameters $\nu$ and $\rho$.
+
+It is intuitive that if $\sup \norm{e} \to 0$, the number of times the
+differentiable curve $C$ can cross a given edge decreases. We will not
+prove convergence, but rather assume that the special cases where it
+might not work, are negligible.
+
+Further, if $\sup \norm{e} \to 0$ it is obvious that the tensor average
+in \eqref{eq:tensor_approx} converges to the tensor in the intersection
+point.
+
+For each $\phi$ parameter, our discretization in the $\rho$ dimension
+can be regarded as a midpoint rule as shown in
+Figure~\ref{fig:line_midpoint}. Thus if $\sup \Delta \rho \to 0$, this
+part of the discretization is fine.
+
+\begin{figure}
+ \input{fig/line_midpoint}
+\end{figure}
+
+The discretization in the $\phi$ dimension can also be regarded as a
+version of the \emph{rectangle method}, although not the midpoint rule.
+As shown in Figure~\ref{fig:circ_rule}, the circle is split into
+intervals
+\begin{equation}
+ \left[
+ \frac{\phi_{k-1} + \phi_{k}}{2}, \frac{\phi_{k} + \phi_{k+1}}{2}
+ \right]
+\end{equation}
+of length $\Delta \phi_k = (\phi_{k+1} + \phi_{k-1}) / 2$. The summand
+is evaluated at $\phi_k$, somewhere inside the interval.
+Thus if
+$\sup \Delta \phi_k \to 0$, this discretization is also consistent.
+
+\begin{figure}
+ \input{fig/circ_rule}
+\end{figure}
+
+To show that all these properties can be fulfilled, we look at a
+particular neighborhood stencil construction.
+Consider a square centered around a grid point with side lengths
+$\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 $n^2 = \lfloor 1 / \sqrt{\delta} \rfloor^2$ grid
+points. This means that the number of grid points along the outer edge
+of this square $n$ goes to infinity.
+
+\begin{figure}
+ \input{fig/square_cons}
+\end{figure}
+
+We include all grid points inside the square in our neighborhood, except
+for multiple points that lie on the same line from the origin. Thus if
+two or more grid points lie on the same line, we include only the one
+closest to the origin. This implies that 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 construction can be seen in Figure~\ref{fig:square_cons} for $n =
+5$.
+The maximal angle between two lines $\phi_k - \phi_{k-1}$ will always be
+when one of $\phi_k$ or $\phi_{k-1}$ is horizontal or vertical, shown in
+Figure~\ref{fig:square_cons} as the angle $a$.
+Thus the
+largest $\Delta \phi_k$ will then be around the horizontal or
+vertical edge.
+These angles can be calculated to be
+\begin{equation}
+ \sup \Delta \phi_k = \frac{2}{2} \arctan \frac{1/n}{n/2} = \arctan
+ \frac{2}{n^2} \to 0.
+\end{equation}
+
+Further we see that the edge length will be bounded by half of the
+diagonal of the square such that
+\begin{equation}
+ \norm{e} \leq \sqrt{\delta / 2} \to 0.
+\end{equation}
+And finally we know from Lemma~\ref{lem:delta_rho} that for each line family
+$\delta^2 = \Delta \rho \norm{e}$ and the fact that $\norm{e} \geq
+\delta$. Thus for the inter-line distance $\Delta \rho$ we have
+\begin{equation}
+ \sup \Delta \rho = \sup \frac{\delta^2}{\norm{e}} \leq
+ \frac{\delta^2}{\delta} = \delta \to 0.
+\end{equation}
+
+And thus the approximation has been shown to be equivalent to
+well-known, and consistent integral approximations, where the
+summand converges to the integrand, and the differences $\Delta \phi$
+and $\Delta \rho$ go to zero. Thus the perimeter approximation in
+\eqref{eq:per_approx1} is consistent with the continuous formulation in
+Theorem~\ref{thm:riemannian_cauchy_crofton}.
+
+Note that as we will work with digital images with fixed resolutions, we
+do not really have the chance to refine our discretization. We do
+however have to take these things into account when creating our
+neighborhood stencil, to make sure that we get a reasonable
+approximation of the perimeter lengths.
+
+\section{Graph cut formulation}
+
+The discretization we arrived at in \eqref{eq:final_discretization} can
+be minimized using graph cuts. For each level $\lambda$, a minimum graph
+cut is found to produce the corresponding level set $\{ u > \lambda \}$. These
+are then combined to form the final restored image $u$.
+
+In this section we will look at how these graphs are constructed such
+that their minimum cuts correspond to the minimizers of the functional
+$F^\lambda$. The description is taken from my project work
+\cite{project}, and is included here for completeness.
+
+\subsection{Graphs}
+
+\fixme{Maybe just call them graphs this time? from $E$ to $F$!}
+
+Using the notation of \cite{cormen2009introduction} we will denote a
+directed graph as $G = (V, E)$ where $V$ is a finite set of vertices, and
+$E$ is a binary relation on $V$. If $(u, v) \in E$ we say that there is
+an edge from $u$ to $v$ in the graph $G$.
+
+We introduce the non-negative capacity function $c : V \times V \to
+\left[0, \infty \right)$. Only edges $(u, v) \in E$ can have a positive
+capacity $c(u, v) = q > 0$ and it means that it is possible to send a
+\emph{flow} of maximum $q$ units from $u$ to $v$. For convenience we
+will let $c(u, v) = 0$ for any pair $(u, v) \not\in E$, and we do not
+allow self-loops in our graph. When a directed graph $G$ is equipped
+with capacity function $c$, one might call it a capacitated graph or a
+graph, but as all our graphs will be capacitated from this point, we
+will just call them graphs and we write $G = (V, E, c)$.
+
+There are two special vertices in the graph, the source $s$ and the
+sink $t$. Contrary to other vertices, which can neither produce nor
+receive excess flow, the source can produce and the sink can receive an
+unlimited amount of flow. The most basic problem in graph flow theory
+is the question of how much flow it is possible to send through the
+graph from the source to the sink.
+
+What we seek in our final graph is a minimum $s$-$t$-cut, a ``minimal''
+line through the graph that cuts a set of edges and
+divides the vertex set in two, separating the source from the sink.
+\begin{definition}[$s$-$t$-cut]
+ Given a graph $G = (V, E, c)$, an $s$-$t$-cut $(S, T)$ of $G$ is a
+ partition of $V$ into $S$ and $T = V - S$ such that $s \in S$ and $t
+ \in T$. The capacity of the cut is
+ \begin{equation}
+ c(S, T) = \sum_{u \in S} \sum_{v \in T} c(u, v),
+ \end{equation}
+ and a minimum $s$-$t$-cut is a cut whose capacity is minimum over
+ all $s$-$t$-cuts.
+ \label{def:s_t_cut}
+\end{definition}
+Note that there might exist several minimum $s$-$t$-cuts in a graph,
+resulting in different partitions of $V$. This is why we need to verify
+later that the cuts we obtain result in stackable thresholded images.
+
+\subsection{Graph representable functionals}
+\label{sec:graph_representable}
+
+The next step is to find a way to construct a graph such that we can
+minimize the functional in \eqref{eq:final_discretization} by finding a minimum
+$s$-$t$-cut. We will do this by creating small and simple graphs
+representing the separate summands of the functional. For these small
+graphs it will be easy to verify that the minimal cut also minimizes
+the corresponding part of the functional, and they can then be
+merged giving a graph representing the complete functional in
+\eqref{eq:final_discretization}.
+
+First we need to establish the definition of a graph representable
+function, presented by Kolmogorov and Zabih in
+\cite{kolmogorov2002energy}.
+\begin{definition}[Graph representable functions]
+ A function $\mathcal{E}(x_1, \hdots, x_n)$ of $n$ binary
+ variables is graph-representable if
+ there exists a graph $G = (V, E, c)$ with terminals $s$ and $t$,
+ and a subset of vertices $V_0 = \{v_1, \hdots, v_n\} \subseteq V -
+ \{s, t\}$ such that, for any configuration $(x_1, \hdots, x_n) \in
+ \{0, 1\}^n$, the value of the energy $\mathcal{E}(x_1, \hdots, x_n)$
+ is equal to a constant plus the cost of the minimum $s$-$t$-cut
+ among all cuts $C = (S, T)$ where $x_i = 0 \Leftrightarrow v_i \in
+ S$ and $x_i = 1 \Leftrightarrow v_i \in T$, $\forall\, 1 \leq i \leq
+ n$.
+ \label{def:graph_representable}
+\end{definition}
+From this definition we see that if we have a graph-representable
+function $\mathcal{E}$ it is possible to find an exact global minimum
+of $\mathcal{E}$ by finding a minimal $s$-$t$-cut in a graph
+representing $\mathcal{E}$.
+
+Furthermore Kolmogorov and Zabih present an important result concerning
+what kinds of functions are graph-representable.
+\begin{theorem}[Identification of graph representable functions]
+ Given an energy function $\mathcal{E}$ of $n$ binary variables of
+ the form
+ \begin{equation}
+ \mathcal{E}(x_1, \hdots, x_n) = \sum_i \mathcal{E}^i(x_i)
+ + \sum_{i<j} \mathcal{E}^{i,j}(x_i, x_j),
+ \end{equation}
+ then $\mathcal{E}$ is graph representable if and only if each term
+ $\mathcal{E}^{i,j}$ satisfies the inequality
+ \begin{equation}
+ \mathcal{E}^{i,j}(0,0) + \mathcal{E}^{i,j}(1,1) \leq
+ \mathcal{E}^{i,j}(0,1) + \mathcal{E}^{i,j}(1,0).
+ \end{equation}
+ \label{thm:graph_rep_id}
+\end{theorem}
+This theorem will allow us to verify that our functional actually
+is graph representable.
+
+Finally, the following theorem, proved by Kolmogorov and Zabih in
+\cite{kolmogorov2002energy}, will be crucial in our graph
+construction.
+\begin{theorem}[Additivity]
+ The sum of a finite number of graph-representable functions
+ \begin{equation}
+ \mathcal{E}(x_1, \hdots, x_n) = \sum_k \mathcal{E}^k(x_1,
+ \hdots, x_n),
+ \end{equation}
+ each represented by a graph $G^k = (V^k, E^k, c^k)$, is
+ graph-representable by $G = (V, E, c)$ where $V = \cup_k V^k$,
+ $E = \cup_k E^k$ and $c(u,v) = \sum_k c^k(u,v)$.
+\end{theorem}
+It allows us to construct small graphs representing the different
+summands of our functional \eqref{eq:final_discretization}, before adding
+them together to create a final graph representing complete functional.
+
+Note that when we apply this theorem later, we will assume that all the
+summands of \eqref{eq:final_discretization} have the whole picture as their
+domain. It is unproblematic to extend $F_\lambda^x(u^\lambda_x)$ and
+$F^{x,y}(u^\lambda_x, u^\lambda_y)$ such that they take all the pixels
+their argument and then ignore all pixels except the ones they actually
+depend on.
+
+\subsection{Graph construction}
+
+\begin{figure}[t]
+ \input{fig/norm_subgraph}
+\end{figure}
+
+We will construct a graph in such a way that if a vertex
+$u^\lambda_x$ ends up on the source side of the cut we set $u^\lambda_x
+= 0$, and if it ends up on the sink side we set $u^\lambda_x = 1$, as in
+Definition~\ref{def:graph_representable}. This is an arbitrary choice,
+but still something we have to keep in mind through the rest of the
+section.
+
+We will now consider the two kinds of summands in the discrete
+functional
+\eqref{eq:final_discretization}. The fidelity term coming from our aim to keep
+the output image close to the original image, and the regularization
+term coming from our aim to minimize the total variation.
+
+\begin{table}[b]
+ \centering
+ \caption{Each row represents one of the two possible values of
+ $u^\lambda_x \in \{0,1\}$. The functional $F_\lambda^x(u^\lambda_x)$
+ and minimum cut obtaining this configuration is shown. The last
+ two columns show the capacities of the cut for each of the two
+ graph constructions in Figure~\ref{fig:norm_subgraph}. We
+ verify that for each of the two graph constructions, the cut
+ capacities are equal to the functional value, plus a constant.
+ }
+ \begin{tabular}{ccccc}
+ \hline
+ $u^\lambda_x$ & $F_\lambda^x(u^\lambda_x)$ & Min.\ cut $(S, T)$
+ & Graph \subref{fig:norm_subgraph_pos} cut cap. & Graph
+ \subref{fig:norm_subgraph_neg} cut cap. \\
+ \hline
+ $0$ & $0$ & $(\{s, u^\lambda_x\}, \{t\})$ & $0$ & $-F_\lambda^x(1)$ \\
+ $1$ & $F_\lambda^x(1)$ & $(\{s\}, \{u^\lambda_x, t\})$ & $F_\lambda^x(1)$ & $0$
+ \end{tabular}
+ \label{tab:fid_energy}
+\end{table}
+
+\subsubsection{Fidelity term}
+
+\begin{figure}[t]
+ \input{fig/neigh_subgraph}
+\end{figure}
+
+The fidelity term of our discrete functional
+\eqref{eq:final_discretization}
+simplifies to
+\begin{align}
+ F_\lambda^x(0) &= 0
+ \label{eq:fidelity_energy0} \\
+ F_\lambda^x(1) &=
+ N_x(\lambda + 1) -
+ N_x(\lambda)
+ \label{eq:fidelity_energy1}
+\end{align}
+where $F_\lambda^x(1)$ might be positive or negative depending on
+$\lambda$ and the pixel value $v_x$.
+
+Figure~\ref{fig:norm_subgraph} shows how graphs can be constructed to
+represent this part of functional. The construction differs
+depending on whether $F_\lambda^x(1)$ is positive or negative.
+Table~\ref{tab:fid_energy} shows how the cuts correspond to the values
+of $u^\lambda_x$ and we can easily verify that the constructed graph
+actually represents the fidelity term.
+
+\subsubsection{Regularization term}
+
+For our regularization term in \eqref{eq:final_discretization} on the form
+\begin{equation}
+ F^{x,y}(u^\lambda_x, u^\lambda_y) =
+ w_{xy} \abs{u^\lambda_x - u^\lambda_y}
+\end{equation}
+we have
+\begin{equation}
+\begin{aligned}
+ F^{x,y}(0, 0) &= 0, \\
+ F^{x,y}(0, 1) &= w_{xy}, \\
+ F^{x,y}(1, 0) &= w_{xy}, \\
+ F^{x,y}(1, 1) &= 0,
+\end{aligned}
+\label{eq:neigh_energies}
+\end{equation}
+and by Theorem~\ref{thm:graph_rep_id} our functional is graph
+representable. In \cite{kolmogorov2002energy}, Kolmogorov and Zabih
+presents a way to construct a graph for any graph representable function
+on the form shown in Theorem~\ref{thm:graph_rep_id}. Since the
+energies in \eqref{eq:neigh_energies} are especially simple, the
+construction and presentation is simplified.
+
+Figure~\ref{fig:neigh_subgraph} shows two different ways of how a graph
+can be constructed to represent the regularization term.
+See Table~\ref{tab:neigh_energy} for an overview of how the two values
+of $u^\lambda_x$ corresponds to cuts in the graph.
+
+\begin{table}[t]
+ \centering
+ \caption{An overview of the possible configurations of the variables
+ in the term $F^{x,y}(u^\lambda_x, u^\lambda_y)$. For each
+ configuration the corresponding functional value and the cut yielding this
+ configuration is shown. The last two columns show the capacities
+ of the cut in the two alternative graph constructions shown in
+ Figure~\ref{fig:neigh_subgraph}. We verify that for each of the
+ two graph constructions, the cut capacities are equal to the
+ functional value, plus a constant.
+ }
+ \begin{tabular}{ccccc}
+ \hline
+ $(u^\lambda_x, u^\lambda_y)$
+ & $F^{x,y}(u^\lambda_x, u^\lambda_y)$
+ & Min.\ cut $(S, T)$
+ & Alt.\ \subref{fig:neigh_subgraph_alt1} cut cap. & Alt.\
+ \subref{fig:neigh_subgraph_alt2} cut cap. \\
+ \hline
+ $(0, 0)$ & $0$ & $(\{s, u^\lambda_x, u^\lambda_y\}, \{t\})$ & $w_{xy}$ & $0$ \\
+ $(0, 1)$ & $w_{xy}$ & $(\{s, u^\lambda_x\}, \{u^\lambda_y, t\})$ & $2w_{xy}$ & $w_{xy}$ \\
+ $(1, 0)$ & $w_{xy}$ & $(\{s, u^\lambda_y\}, \{u^\lambda_x, t\})$ & $2w_{xy}$ & $w_{xy}$ \\
+ $(1, 1)$ & $0$ & $(\{s\}, \{u^\lambda_x, u^\lambda_y, t\})$ & $w_{xy}$ & $0$
+ \end{tabular}
+ \label{tab:neigh_energy}
+\end{table}
+
+Figure~\ref{fig:big_graph} shows a visualization of how the final
+graph might look with all its edges. The source will have a lot of
+outgoing edges, one for each pixel, while the sink has one incoming edge
+from each pixel. The vertices corresponding to the pixels are only
+connected to the source, the sink, and their neighboring pixels, so
+their edge degree is much smaller than for $s$ and $t$.
+
+\begin{figure}[bh]
+ \input{fig/big_graph}
+\end{figure}
+
\path[edge] (s) -- node[weight noslope] {$\max\{F_{L-1}^x(1), 0\}$} (u);
\path[edge] (u) -- node[weight noslope] {$\max\{F_{L-1}^x(1), 0\} - F_\lambda^x(1)$} (t);
\end{tikzpicture}
-\caption{%
- \fixme{boop}
+\caption[Modified fidelity subgraph for the Boykov--Kolmogorov
+algorithm]{%
+ A modified version of the subgraph corresponding to the fidelity
+ term $F_\lambda^x(u^\lambda_x)$ that has positive and non-decreasing
+ edges with decreasing $\lambda$.
}
\label{fig:bk_norm_subgraph}
% Inkluder alle de andre kildefilene
\input{introduction}
-\input{theory}
+\input{continuous}
+\input{discrete}
+\input{maxflow}
\input{results}
\cleardoublepage
--- /dev/null
+\chapter{Maximum flow approach}
+
+In the previous section we have seen how finding the minimum cut of
+carefully constructed graph can give us the thresholded image
+minimizing the functional for one level value $\lambda$. We will in
+this and the next section see how such a minimum cut can be found by
+sending flow through the graph and trying to identify the
+``bottleneck''. This chapter, except for the description of the
+Boykov--Kolmogorov algorithm is taken from my project work
+\cite{project} and is included here for completeness.
+
+\section{Flow graphs}
+
+We have already introduced capacities, and briefly mentioned the notion
+of flow as something limited by the capacity. In other words, flow is
+something we can send through our graph, but the capacity limits how
+much we can send along each edge. It is useful to imagine a water supply
+graph with pipes of different sizes and a water source and sink.
+
+Formally, we introduce the \emph{flow} as the function $f : V \times V
+\to [0, \infty)$. This function keeps count of how much flow we are
+sending through each edge of our graph and must satisfy the following
+two constraints
+\begin{description}
+ \item[Capacity constraint:] For all $u, v \in V, 0 \leq f(u, v) \leq c(u,
+ v)$, i.e., for every pair of vertices, the flow is less than or
+ equal to the capacity.
+ \item[Flow conservation:] For all $u \in V - \{s, t\}$
+ \begin{equation}
+ \sum_{v \in V} f(v, u) = \sum_{v \in V} f(u, v),
+ \end{equation}
+ i.e., for any vertex except the source and the sink, the
+ flow into the vertex must be equal to the flow out of the
+ vertex.
+\end{description}
+Note that we have defined $f$ with all pairs of vertices as its domain,
+even though it is only non-zero on edges $(u, v) \in E$. This makes it
+easier to write sums as in the flow conservation constraint. We say that
+an edge $(u, v)$ is \emph{saturated} if $f(u, v) = c(u, v)$.
+
+We define $\abs{f}$ as the net amount of flow from the source to the
+sink in the graph. Because of the flow conservation constraint, this
+can be calculated as the net amount of flow going out of the source
+\begin{equation}
+ \abs{f} = \sum_{v \in V} f(s, v) - \sum_{u \in V} f(u, s).
+ \label{eq:total_flow}
+\end{equation}
+Furthermore we denote the \emph{net flow} across an $s$-$t$-cut $C = (S,
+T)$ as
+\begin{equation}
+ f(S, T) = \sum_{u \in S} \sum_{v \in T} f(u, v) - \sum_{v \in T}
+ \sum_{u \in S} f(v, u).
+ \label{eq:net_flow_cut}
+\end{equation}
+Note how this definition differs from the capacity of a cut $c(S,T)$ in
+Definition~\ref{def:s_t_cut}. While the capacity of a cut represents how
+much flow it is maximally possible to send from $S$ to $T$, the net flow
+across a cut represents the net amount of flow going across the cut,
+counting negatively the flow that goes back from $T$ to $S$.
+
+For any $s$-$t$-cut we have that
+\begin{equation}
+ \abs{f} = f(S, T).
+ \label{eq:total_flow_eq_cut}
+\end{equation}
+This is quite intuitive given the flow conservation constraint, and a
+full chain of arguments can be found in \cite{cormen2009introduction}.
+
+Following from \eqref{eq:net_flow_cut} and \eqref{eq:total_flow_eq_cut},
+we find that
+\begin{equation}
+ \begin{aligned}
+ \abs{f} &= \sum_{u \in S} \sum_{v \in T} f(u, v)
+ - \sum_{v \in T} \sum_{u \in S} f(v, u) \\
+ &\leq \sum_{u \in S} \sum_{v \in T} f(u, v) \\
+ &\leq \sum_{u \in S} \sum_{v \in T} c(u, v) \\
+ &\leq c(S, T)
+ \label{eq:flow_leq_cut}
+ \end{aligned}
+\end{equation}
+for any $s$-$t$-cut $C = (S, T)$. A very central result in graph flow
+theory called the max-flow min-cut theorem will be presented later. It
+states that the inequality of \eqref{eq:flow_leq_cut} becomes an equality
+when $f$ is a maximum flow for some cut $C = (S, T)$, and that all such
+cuts are minimum cuts.
+
+%Recall that we want to find the minimum $s$-$t$ cut in our graph. When
+%finding this minimum cut, we make use of an important duality theorem in
+%graph flow theory, stating that the capacity of a minimum $s$-$t$-cut
+%in a graph, is equal to the maximum flow from the source to the sink.
+
+But how does this help us? We know that if we know the maximum flow
+value, and we have an $s$-$t$-cut with capacity equal to the maximum
+flow, we actually have a minimum cut. The question is then, how do we
+find a maximum flow, and how do we find a minimum cut?
+
+\section{Augmenting path algorithms}
+The family of augmenting flow algorithms represent a popular approach to
+the maximum flow problem. The idea is simply to look for paths from the
+source to the sink that can carry additional flow, so-called augmenting
+paths, and then send the maximum possible amount of flow along such a
+path. When no such path exists anymore, no more flow can be sent from
+the source to the sink, and a maximum flow has been reached.
+
+\subsection{Residual graph}
+When further discussing approaches to solving the maximum flow problem
+we will need the notion of a residual graph $G_f = (V_f, E_f, c_f)$,
+which is derived from the original graph $G$ and contains the edges
+along which it is possible to send additional flow. This means that
+$E_f$ contains the edges $(u, v)$ from $E$ where $f(u,v) < c(u,v)$. But
+that is not all; an important realisation is that it is also possible to
+push flow \emph{back} along an edge where the flow is already positive.
+In other words, sending flow from $v$ to $u$ by cancelling some or all
+of the flow that is already going from $u$ to $v$.
+
+Thus the capacity function $c_f$ of our residual graph becomes
+\begin{equation}
+ c_f(u, v) =
+ \begin{cases}
+ c(u, v) - f(u, v) & \text{if } (u, v) \in E, \\
+ f(v, u) & \text{if } (v, u) \in E, \\
+ 0 & \text{otherwise.}
+ \end{cases}
+\end{equation}
+The vertices $V_f$ of $G_f$ are the same as the original graph $G$,
+while the edges $E_f$ are taken to be all pairs of vertices $(u, v)$
+with $c_f(u, v) > 0$. Since we can in $E_f$ at most have all the
+original edges, and their reversals, we have $\abs{E_f} \leq 2 \abs{E}$.
+
+Note that there is ambiguity in the definition of the residual graph
+in the case where the original graph contains anti-parallel edges. One
+could avoid this by defining $c_f(u, v) = f(v, u) + c(u, v) - f(u, v)$
+instead, or as they do in \cite{cormen2009introduction}, disallow
+anti-parallel edges altogether. In any case it is not something we will
+have to think about in the implementation, since we will not actually
+construct the residual graph.
+
+With the residual graph defined, we are ready to formally present the
+max-flow min-cut theorem.
+\begin{theorem}[Max-flow min-cut theorem]
+ If $f$ is a flow in a graph $G = (V, E, c)$ with source $s$ and
+ sink $t$, then the following statements are equivalent:
+ \begin{enumerate}
+ \item $f$ is a maximum flow in $G$.
+ \item The residual graph $G_f$ contains no augmenting paths.
+ \item $\abs{f} = c(S, T)$ for some cut $(S, T)$ of $G$.
+ \label{it:max_flow_min_cut}
+ \end{enumerate}
+\end{theorem}
+See \cite{cormen2009introduction} for a proof, and remark that because
+of the inequality in \eqref{eq:flow_leq_cut}, the cut in
+Statement~\ref{it:max_flow_min_cut} is a minimum cut. The theorem does
+not tell us how to find such a cut, and there are multiple ways. One
+possibility is to take $S$ to be all vertices reachable from the source
+in the residual graph and $T = V - S$.
+
+Figure~\ref{fig:aug_flow} shows a simple graph which already has five
+units flowing from $s$ to $t$. The marked path is a possible augmenting
+path, and note that it follows an edge in $E$ in the reverse direction,
+made possible by the construction of the residual graph just
+described.
+
+\begin{figure}
+ \input{fig/aug_flow}
+\end{figure}
+
+\subsection{Ford-Fulkerson}
+The Ford-Fulkerson algorithm is the most basic augmenting path
+algorithm, which can be extended to more advanced algorithms. It is
+stated in pseudocode in Algorithm~\ref{alg:ford_fulkerson}, and the idea
+is to augment the flow along paths from $s$ to $t$ until it is no longer
+possible.
+\begin{algorithm}
+\begin{algorithmic}
+ \Function{Ford-Fulkerson}{$G$, $s$, $t$}
+ \While{there exists a path $p$ from $s$ to $t$ in the residual
+ graph $G_f$}
+ \State $\alpha \gets \min\{c_f(u, v) : (u, v) \in p \}$
+ \ForAll{$(u, v) \in p$}
+ \If{$(u, v) \in E$}
+ \State $f(u, v) \mathrel{+}= \Delta f$
+ \Else
+ \State $f(v, u) \mathrel{-}= \Delta f$
+ \Comment{Push flow back}
+ \EndIf
+ \EndFor
+ \EndWhile
+ \EndFunction
+\end{algorithmic}
+\caption{\sf The Ford-Fulkerson max-flow algorithm}
+\label{alg:ford_fulkerson}
+\end{algorithm}
+
+There are different ways to find augmenting paths, and a common choice
+is to do a breadth-first search from the source until the sink is found,
+as this will yield the shortest possible augmenting path. This version
+of the algorithm is called Edmonds-Karp and has a running time of
+$O(\abs{V}\abs{E}^2)$. See for example \cite{cormen2009introduction}
+for a description of the breadth-first search, and a formal proof of the
+running time of the Edmonds-Karp algorithm.
+
+\section{Other algorithms}
+
+There are many different maximum flow algorithms that fall into the
+augmenting path category, although we will see a different approach in
+the next section.
+
+The algorithm of Dinitz, originally published in 1970, later improved
+on, and described by the original author in \cite{dinitz2006dinitz}, is
+a variant of the augmenting path algorithm. It maintains a distance
+labeling $d(u)$ of the vertices $u \in V$ in the graph, where $d(u)$
+is the shortest path from the $s$ to $u$ in the residual graph. This
+can be computed with a simple breadth-first search. The next step is to
+construct a \emph{blocking flow} $f'$, using only edges in $E_f' = \{
+(u,v) \in E_f : d(u) + 1 = d(v) \}$. The blocking flow is such that if
+we augment the flow $f$ by $f'$, there is no longer any paths from $s$
+to $t$ following edges in $E_f'$. After the blocking flow $f'$ has been
+found and added to $f$, the distance labels are recalculated, and the
+label of the sink will be increased by at least one.
+
+Boykov and Kolmogorov present a variant of the augmenting path algorithm
+in \cite{boykov2004experimental}, specialized for the kinds of graphs
+occuring in graphical applications. It will be presented in
+Section~\ref{sec:bk}
+
+\section{Push--relabel algorithm}
+
+The push-relabel algorithm is a different approach to the maximum flow
+problem, presented by Goldberg and Tarjan in \cite{goldberg1988new}.
+Unlike the augmenting flow algorithms, it does not maintain a valid flow
+$f$ in the graph at all times, but still obtains a valid maximum flow
+when the algorithm terminates.
+
+\subsection{Preflow}
+Instead of maintaining a valid flow, we introduce the concept of a
+\emph{preflow} by relaxing the flow conservation constraint from
+earlier. We allow positive excess in the vertices and the flow
+conservation constraint from before then becomes
+\begin{description}
+ \item[Preflow conservation:] For all $u \in V - \{s, t\}$
+ \begin{equation}
+ \sum_{v \in V} f(v, u) \geq \sum_{v \in V} f(u, v),
+ \end{equation}
+ i.e., for any vertex except the source and the sink, the
+ flow into the vertex must greater or equal to the flow out of
+ the vertex.
+\end{description}
+
+As in most of the cited push-relabel literature, we define $N =
+\abs{V}$, and for all vertices $u \in V$ we define the excess
+\begin{equation}
+ e(u) = \sum_{v \in V} f(v, u) - \sum_{v \in V} f(u, v),
+\end{equation}
+which represents the amount of flow which \emph{disappears} in vertex
+$u$. Equivalent to the preflow conservation constraint is stating that
+$e(u) \geq 0$ for all vertices $u \in V - \{s,t\}$.
+
+%The idea of the algorithm is to maintain a height map of the
+%vertices in the graph where connected vertices can not have a large
+%height difference. Then we ``lift'' the source vertex to let as
+%many units as possible flow through the edges of the graph towards the
+%sink. When a maximum preflow is reached, there will normally be excess
+%flow in some of the vertices, which has to be pushed back towards the
+%source in order to obtain a valid flow.
+
+In addition to the flow, we maintain a height map $d : V \to \mathbb{N}$
+that satisfies $d(t) = 0$, and for every edge $(u, v)$ in the residual
+graph, i.e.\ every edge with $c_f(u, v) > 0$, we require that $d(u) \leq
+d(v) + 1$. For all vertices $u$, the label $d(u)$ will be a lower bound
+on the length from $u$ to $t$ in $G_f$ which is why it is also often
+called a distance labeling.
+
+A vertex $u$ is \emph{active} if $u \in V - \{s,t\}$, it has positive
+excess ($e(u) > 0$) and $d(u) < N$. These are the vertices we want to
+operate on to increase the preflow.
+
+\subsection{Basic operations}
+
+The algorithm performs two basic operations, the \emph{push} and
+\emph{relabel} operations, while always maintaining a valid preflow $f$
+and a valid distance labeling $d$.
+
+\subsubsection{The push procedure}
+
+The push procedure moves excess flow from an active vertex along an edge
+$(u, v) \in E_f$ for which $d(u) = d(v) + 1$, i.e.\ to a vertex with a
+smaller distance label. We call such edges \emph{admissible}. See
+Algorithm~\ref{alg:push} for a pseudocode implementation of the push
+operation.
+
+Assuming that $f$ is a valid preflow, it is easy to verify that the
+preflow $f$ and labeling $d$ remain valid after running the push
+procedure on some admissible edge $(u, v)$.
+
+%The capacity constraint is fulfilled since we
+%at most increase the flow along $(u, v)$ with the residual capacity
+%$c_f(u,v)$. The preflow constraint is fulfilled since the excess $e$
+%increases for $v$, remains non-negative for $u$ and remains the same for
+%all other vertices.
+
+%The residual graph might change during the push procedure. The edge
+%$(v, u)$ will appear if it does not already exist. Initially $d(u) =
+%d(v) + 1$, so for the new edge $(v, u)$ we have
+%\begin{equation}
+% d(v) = d(u) - 1 \leq d(u) + 1
+%\end{equation}
+%and $d$ remains a valid labeling after the push procedure is finished.
+
+\begin{algorithm}
+\begin{algorithmic}
+ \Function{Push}{$u$, $v$}
+ \State $\Delta f \gets \min(c_f(u, v), e(u))$
+ \If{$(u, v) \in E$}
+ \State $f(u, v) \mathrel{+}= \Delta f$
+ \Else
+ \State $f(v, u) \mathrel{-}= \Delta f$
+ \Comment{Push flow back}
+ \EndIf
+ \Comment{Excess $e(u)$ and $e(v)$ will also change}
+ %\State $f(v, u) \mathrel{-}= f_\text{aug}$
+ %\State $e(u) \mathrel{-}= f_\text{aug}$
+ %\State $e(v) \mathrel{+}= f_\text{aug}$
+ \EndFunction
+\end{algorithmic}
+\caption{%
+ \sf The push procedure of the Push-Relabel algorithm.
+}
+\label{alg:push}
+\end{algorithm}
+
+\subsubsection{The relabel procedure}
+
+The relabel procedure is our tool for changing the distance labels of
+the vertices. It changes the label of a vertex to the greatest possible
+value, which is one more than the lowest label among its neighbors in
+the residual graph. See Algorithm~\ref{alg:relabel} for a pseudocode
+implementation.
+\begin{algorithm}
+\begin{algorithmic}
+ \Function{Relabel}{$u$}
+% \If{$u$ is only vertex at its height}
+% \Call{Gap}{$u$}
+% \Else
+ \If{there is a $v \in V$ such that $(u, v) \in E_f$}
+ \State $d(u) \gets \min\{d(v), \; \forall v \in V : (u,v) \in E_f\} + 1$
+ \Else
+ \State $d(u) \gets N$
+ \Comment{$u$ becomes inactive}
+ \EndIf
+% \EndIf
+ \EndFunction
+\end{algorithmic}
+\caption{\sf The relabel procedure of the Push-Relabel algorithm}
+\label{alg:relabel}
+\end{algorithm}
+
+If $d$ was a valid labeling before running the relabel procedure, then
+we still have $d(u) \leq d(v) + 1$ for all neighbors $v$ of $u$ in the
+residual graph, and $d$ remains a valid labeling. The capacity
+constraint and preflow constraint remain satisfied assuming they were
+satisfied before the procedure was started.
+
+\subsection{Putting it all together}
+These basic procedures are then applied to active vertices and admissible
+edges until we obtain our minimum cut. We will see later that when there
+are no more active vertices, we can extract the minimum cut from the
+graph.
+
+In the first phase of the algorithm we initialize a valid preflow and
+distance labeling by saturating all edges out of the source $s$, and
+then setting its distance label $d(s) = N$. We then apply the push and
+relabel procedures where applicable until there are no more active
+vertices and a maximum preflow is obtained.
+
+A vertex $u$ can only be successfully relabeled to obtain a new label if
+the outgoing edges of $u$ in the residual graph have changed since the
+previous relabeling. This is why the push and relabel procedures often
+are combined into a \emph{discharge} procedure as shown in
+Algorithm~\ref{alg:discharge}. When it is run on an active vertex $u$,
+we push as much as possible of the excess flow to other vertices before
+the vertex is relabeled.
+\begin{algorithm}
+\begin{algorithmic}
+ \Function{Discharge}{$u$}
+ \ForAll{$v \in V$ such that $(u, v) \in E_f$}
+ \If{$c_f(u, v) > 0$ and $d(u) = d(v) + 1$}
+ \Call{Push}{$u$, $v$}
+ \EndIf
+ \EndFor
+
+ \If{$e(u) > 0$}
+ \Call{Relabel}{$u$}
+ \EndIf
+ \EndFunction
+\end{algorithmic}
+\caption{\sf The discharge procedure of the Push-Relabel algorithm}
+\label{alg:discharge}
+\end{algorithm}
+
+In the second phase of the algorithm this maximum preflow is turned into a
+maximum flow by returning excess flow which did not reach the sink from
+inside the graph back to the source. We can skip this part of the
+algorithm, as it is possible to identify a minimum cut as soon as the
+first phase is finished, and the following theorem allows us to do that.
+\begin{theorem}[Cut identification]
+ Given a graph $G = (V, E, c)$, assume that the first phase of
+ the push-relabel algorithm has terminated so that no more active
+ vertices remain. Then there exists a $k \in \{1, \ldots, N-1\}$
+ such that there is no vertex with label $k$. For every
+ such $k$ the vertex sets $S = \{ u \in V : d(u) > k\}$ and $T = \{ u
+ \in V : d(u) < k \}$ define a minimum cut $C = (S, T)$ in $G$.
+ \label{thm:cut_identification}
+\end{theorem}
+\begin{proof}
+ There are $N$ vertices, the source has label $N$ and the sink has
+ label $0$, and the $N - 2$ remaining vertices can not occupy all the
+ $N-1$ labels in $\{1, \ldots, N-1\}$, so there must exist an $k$ as
+ described.
+
+ There can be no edge $(u, v) \in E_f$ such that $u \in S$ and $v \in
+ T$, as this would imply $k \leq d(u) - 1 \leq d(v) < k$. From the
+ construction of $E_f$ we now know that all edges in $E$ from $S$ to
+ $T$ are saturated, and all edges from $T$ to $S$ carry no flow. This
+ implies the capacity of the cut is equal to the flow through the
+ cut, i.e.\ $c(S, T) = f(S, T)$.
+
+ Since the first phase of the algorithm has terminated, there can be
+ no active vertices, and therefore no excess in $T$, except for the
+ sink. If all flow excess in vertices in $S$ is returned to the
+ source, we can apply the max-flow min-cut theorem to conclude that
+ $C = (S, T)$ is a minimum $s$-$t$-cut, since the cut capacity is
+ equal to the flow.
+\end{proof}
+
+We will see later that with the gap relabeling heuristic, there will
+always be a gap at label $k = N - 1$ such that we can construct our cut
+by taking $S = \{ u \in V : d(u) \geq N \}$.
+
+Note that the vertices in $S$ are vertices earlier described as being on
+the source side of the cut, as no additional flow can go from these
+vertices to the sink.
+
+\subsection{Complexity}
+
+In their original article \cite{goldberg1988new}, Goldberg and Tarjan
+analyze the complexity of the push-relabel algoritm by considering the
+maximum number of basic operations we can possibly do before the
+algorithm terminates.
+
+The number of relabelings is in $O(\abs{V}^2)$ since every time the
+procedure is applicable to a vertex $u$, the label $d(u)$ increases by
+at least one.
+
+The number of saturating pushes is in $O(\abs{V}\abs{E})$. When a push
+along $(u,v)$ is saturating, the label of $v$ has to increase with at
+least 2 before a push can saturate the same edge (in the opposite
+direction). Since the number of relabelings of a single vertex is
+bounded by $\abs{V}$, and we have $\abs{E}$ edges, this gives the stated
+number of saturating pushes.
+
+The number of non-saturating pushes is the most complicated to bound,
+and will also make up the asymptotic running time of the algorithm. The
+idea is to define
+\begin{equation}
+ \phi = \sum_{\mathclap{u \text{ active}}} d(u),
+\end{equation}
+and look at how much this number changes throughout the algorithm. It
+starts at zero and ends at zero. Every non-saturating push from $u$ to
+$v$ makes $\phi$ decrease by at least one since it makes $u$ inactive
+(but might activate $v$). The total increase in $\phi$ due to
+relabelings is less than $\abs{V}^2$. A saturating push from $u$ to $v$
+increases $\phi$ by at most $\abs{V}$, since $v$ might become active.
+
+Even if $\phi$ is always increased by relabelings and saturating pushes,
+we can bound the number of non-saturating pushes by
+\begin{equation}
+ \abs{V}^2 + \abs{V} \underbrace{c \abs{V}
+ \abs{E}}_{\mathclap{\#(\text{saturating pushes})}}
+\end{equation}
+which means that in the general case, the algorithm has a complexity of
+$O(\abs{E}\abs{V}^2)$.
+
+\subsection{Vertex selection rules}
+Until now we have stated that the discharge procedure is run on active
+vertices until there are no more active vertices left. The choice of the
+order in which to discharge these active vertices remain, and multiple
+possibilities exist.
+
+The First In, First Out (FIFO) approach is to always maintain a queue of
+active vertices. When the vertex from the beginning of the queue is
+discharged, other vertices might become active, and these are added at
+the end of the queue. The original article of Goldberg and Tarjan
+\cite{goldberg1988new} contains a proof that this selection rule gives a
+complexity of $O(\abs{V}^3)$, which can be very good if you have a dense
+graph.
+
+Using highest level selection rule one always discharges the active
+vertex with the largest distance label. Goldberg and Tarjan state that
+this rule also gives a complexity of $O(\abs{V}^3)$ while this bound is
+improved to $O(\abs{V}^2 \sqrt{\abs{E}})$ in an article by Cheriyan and
+Maheshwari \cite{cheriyan1989analysis}.
+
+These are complexity bounds, and the actual running time of the
+algorithm, which can only be determined by implementing it and running
+it, varies a lot with the structure of the input graph.
+
+Cherkassky and Goldberg describe the algorithm along with different
+selection rules, heuristics and their implementation in
+\cite{cherkassky1997implementing}.
+
+\subsection{Heuristics}
+Different heuristics exist that can speed up the algorithm
+considerably. Being heuristics, they are not guaranteed to work, and
+might perform differently on different kinds of graphs. The most used
+heuristics are the gap and global relabeling heuristics, both aiming to
+reduce the total number of relabeling steps.
+
+The gap relabeling heuristic aims to find a label $k$ as in
+Theorem~\ref{thm:cut_identification} such that no vertex has that label.
+From vertices $u$ with $d(u) > k$, there are no unsaturated edges going
+to vertices with smaller distance labels, so no more flow can ever find
+its way from these vertices to the sink. These vertices are therefore
+given the label $N$ and never considered again as they will never become
+active. Algorithm~\ref{alg:gap} shows a pseudocode representation of
+what is done once a gap $k$ is found.
+
+\begin{algorithm}
+\begin{algorithmic}
+ \Function{Gap}{$k$}
+ \ForAll{$u \in V$ such that $d(u) > k$}
+ \State $d(u) \gets N$
+ \EndFor
+ \EndFunction
+\end{algorithmic}
+\caption{\sf The gap procedure of the Push-Relabel algorithm}
+\label{alg:gap}
+\end{algorithm}
+
+The preflow and capacity constraints are still valid after the gap
+relabeling procedure, as only the labels $d$ are changed. For the
+labeling $d$ one has to verify that $d(u) \leq d(v) + 1$ for every edge
+$(u,v) \in E_f$ in the residual network. If none, both, or only $v$ is
+relabeled, this is trivial. It is not possible that only $v$ would be
+relabeled, as this would imply that $d(u) \geq d(v) + 2$ which is not a
+valid initial labeling.
+
+%Before integrating the gap relabeling procedure into our algorithm we
+%need to verify that it does not invalidate our preflow $f$ or distance
+%labeling $d$.
+
+%But why does this work? The only thing we need to verify is that given a
+%graph with a valid preflow and a valid labeling, the gap relabeling
+%procedure will not change the validity of these two things.
+
+%\begin{lemma}[Gap relabeling validity]
+% Given a graph $G = (V, E, c)$, a distance labeling $d$ and a
+% preflow $f$, assume there exists a gap $k$ such that no vertex has
+% label $k$. Running the gap relabeling procedure on label $k$ will
+% yield a valid distance labeling and an unchanged and valid preflow
+% $f$.
+%\end{lemma}
+%\begin{proof}
+% No new edges are created, no edges disappear, the preflow is
+% unchanged, so the preflow and capacity constraint remain satisfied
+% after the gap relabeling.
+%
+% Define $S = \{ u \in V : d(u) > k\}$ and $T = V - S$. Initially
+% $d(u) \leq d(v) + 1$ for every edge $(u,v) \in E_f$. These
+% inequalities have to hold after the gap procedure is finished, when
+% $d(u) = N$ for all $u \in S$.
+%
+% For $(u, v) \in E_f$ we have four possibilities
+% \begin{description}
+% \item[$u, v \in T$\textmd{:}]
+% The labels $d(u)$ and $d(v)$ remain unchanged and the
+% inequality still holds.
+% \item[$u, v \in S$\textmd{:}]
+% After the gap procedure we have $d(u) = d(v)$ so the
+% inequality still holds.
+% \item[$u \in S, v \in T$\textmd{:}]
+% This is not possible as it would imply $d(u) \geq d(v) + 2$
+% and we have assumed an initial valid labeling.
+% \item[$u \in T, v \in S$\textmd{:}]
+% After relabeling we have $d(u) < k < N < d(v) + 1$.
+% \end{description}
+%
+% Hence, both the preflow $f$ and distance labeling $d$ are valid.
+%\end{proof}
+
+When running the push-relabel algorithm with the gap heuristic, we can
+be sure that there will never be a vertex $u$ with label $d(u) = N-1$ at
+the end of the algorithm, i.e.\ we know that there will always be a gap
+at label $N-1$. This can be seen using the same reasoning as in
+Theorem~\ref{thm:cut_identification}, because if there was a vertex with label
+$N-1$, there would only be $N-3$ vertices possibly having labels in
+$\{1, \ldots, N-2\}$, so a gap must exist somewhere in that interval.
+When using the gap relabeling heuristic, such a gap can not exist, so we
+can conclude that there is no vertex with label $N-1$.
+
+Using Theorem~\ref{thm:cut_identification} we can then conclude that the
+sets $S = \{ u \in V : d(u) \geq N\}$ and $T = V - S$ form a minimum cut
+of the graph.
+
+\subsection{Parametric push-relabel algorithm}
+
+Now we have an algorithm for finding a minimum $s$-$t$-cut in a graph,
+so let's return to the graph constructed in
+Section~\ref{sec:graph_construction}. For every level $\lambda \in \{0,
+\ldots, L\}$ we want to find a minimum $s$-$t$-cut which gives us the
+thresholded image $u^\lambda$. These can then hopefully be stacked
+together to form the final image $u$.
+
+\subsubsection{Graph reuse}
+
+Solving $L$ separate minimum cut problems seems like a lot of work, but
+when using the push-relabel algorithm we will, if we do things in the
+right order, be able to reuse the graph when going from one label to
+the next.
+
+\begin{figure}
+ \input{fig/norm_evolution}
+\end{figure}
+
+Going back to the graph representations in
+Figure~\ref{fig:norm_subgraph} and Figure~\ref{fig:neigh_subgraph} we
+know that only the capacity of edges from sub-graphs representing the
+fidelity term depend on our level parameter $\lambda$. From the
+expression in \eqref{eq:fidelity_energy0}, visualized in
+Figure~\ref{fig:norm_evolution}, we see that the energy term
+$F_\lambda^x(1)$ increases monotonically with increasing $\lambda$
+parameter. Let $u, v \in V - \{s, t\}$. Since the edges in
+Figure~\ref{fig:norm_subgraph} now are the only ones depending on
+$\lambda$, the following is true for \emph{decreasing} values of
+$\lambda$
+\begin{description}
+ \item[Edges from $s$ to $u$\textmd{:}]
+ As seen in Figure~\ref{fig:norm_subgraph_neg} the capacity of
+ these edges will increase monotonically with decreasing
+ $\lambda$ parameter.
+ \item[Edges from $u$ to $v$\textmd{:}]
+ These edges have no $\lambda$-dependence and will remain
+ unchanged.
+ \item[Edges from $v$ to $t$\textmd{:}]
+ As seen in Figure~\ref{fig:norm_subgraph_pos} the capacity of
+ these edges will decrease monotonically with decreasing
+ $\lambda$ parameter.
+\end{description}
+
+After running the push-relabel algorithm for $\lambda = k$, we are left
+with a graph $G = (V, E, c)$, a preflow $f$ and a labeling $d$. To
+obtain the graph for $\lambda = k-1$ we have to change the capacity of
+two different kinds of edges, and this is done in the following way to
+keep the capacity, preflow and labeling constraints satisfied.
+\begin{description}
+ \item[Edges from $s$ to $u$\textmd{:}]
+ The capacity $c(s, u)$ is increased, and the flow is set to be
+ equal to the capacity $f(s, u) = c(s, u)$. The vertex $u$ might
+ have an increased excess $e(u)$, which might in turn make it
+ active.
+ \item[Edges from $v$ to $t$\textmd{:}]
+ The capacity $c(v, t)$ is decreased. If it is decreased to a
+ value below the current flow value, we set $f(v, t) = c(v, t)$
+ which will decrease the excess of the sink $t$, and increase the
+ excess of $v$.
+\end{description}
+
+None of these actions will create new edges in the residual graph, and
+we do not change the labeling $d$, so the labeling constraints are also
+satisfied in the new graph.
+
+Through this procedure we have easily created the graph for $\lambda =
+k-1$, and the distance labels remain the same. As these labels always
+increase monotonically, we have a head start compared to the case where
+we reset the flow and labels.
+
+\subsubsection{Output image construction}
+
+We mentioned already in Chapter~\ref{chap:discrete} that in order to
+be able to construct our output image $u$, the thresholded images
+$u^\lambda$ would have to stack one on top of the other.
+Because of the reuse of the distance labels
+between the iterations of the push-relabel algorithm, we can guarantee
+that this is possible.
+
+Consider two subsequent runs of the push-relabel algorithm, for labels
+$\lambda$ and $\lambda - 1$ ending with distance labels $d^\lambda$ and
+$d^{\lambda-1}$ respectively. We already know that the distance labels
+$d$ are monotonically increasing. This means that the set $S = \{ u \in
+V : d(u) \geq N \}$ is increasing in size, more precisely, we have the
+inclusion
+\begin{equation}
+ \{ u \in V : d^\lambda(u) \geq N \} \subseteq \{ u \in V :
+ d^{\lambda-1}(u) \geq N \}.
+\end{equation}
+For a pixel $x \in S$ we will set $u^\lambda_x = 0$, which together with
+the previous inclusion property implies
+\begin{equation}
+ u^\lambda_x \geq u^{\lambda-1}_x
+\end{equation}
+for all $x \in \mathcal{G}$. Being equivalent with the inequality in
+\eqref{eq:decreasing_sets}, this means our algorithm produces stackable
+thresholded images $u^\lambda$.
+
+We then construct our output image $u$ by giving each pixel the value
+\begin{equation}
+ u_x = \min \{ \lambda \in \{0, \ldots, L-1\} : u^\lambda_x = 1 \}.
+\end{equation}
+This marks the end of the description of the implemented algorithm, but
+we will further discuss some possible improvements, and also look at the
+results when using the method on different kinds of noisy images.
+
+\subsection{Divide and conquer}
+The possibility of re-using the graph between separate level
+is a very nice property of the push-relabel algorithm, but there are
+further room for improvements. Consider one pixel $x$ with value $u_x$,
+and imagine we only wanted to find the value of this pixel. One could go
+through all pixel values $\lambda \in (L-1, \ldots, 0)$, and see when
+$u^\lambda_x$ changes from $1$ to $0$, just as we do for all the pixels
+in the algorithm above. Ignoring graph re-use this would have us solve
+$O(L)$ maximum flow problems.
+
+Improving on this we could employ the idea of binary search to find the
+value of $u_x$ in only $O(\log_2 L)$ time. After finding one cut, we
+know whether $u_x$ is above or below the current $\lambda$ value, and by
+choosing $\lambda$ as the midpoint of the current possible range of
+$u_x$, we can cut the search space in half for each iteration of the
+algorithm.
+
+We can extend this idea to the problem of finding all pixel values.
+Instead of running the algorithm for successively decreasing values of
+$\lambda$, we choose some $\lambda$ in the middle of the range $\{0,
+\ldots, L-1\}$. The cut we obtain consists of two sets $S = \{ u \in V :
+d(u) \geq N\}$ and $T = V - S$. We know that no more flow can be sent
+from $S$ to $T$, even if we decrease the value of $\lambda$ and adjust
+the capacities accordingly.
+
+The idea is now that we have halved the possible $\lambda$ interval for
+\emph{all} pixels. We continue by considering the two sets $S$ and $T$
+separately, and applying the algorithm recursively, at each time halving
+the $\lambda$ interval until we have the value of every pixel.
+
+Combining the divide and conquer approach with the parametric
+push-relabel algorithm is a bit problematic. For all pixels $x \in T$ we
+know that $u_x \leq \lambda$, and we can reuse the graph when
+decreasing $\lambda$. However for pixels $x \in S$, we seek to find $u_x
+> \lambda$, meaning we have to increase $\lambda$ which does not allow
+graph re-use.
+
+Goldfarb and Yin \cite{goldfarb2009parametric} have found that the
+divide and conquer approach only yields improved performance when using
+the $L^2$ norm in the fidelity term. This has to do with the fact that
+the fidelity term for the $L^1$ norm only changes once for each pixel,
+as we can see in Figure~\ref{fig:norm_evolution}, reducing the amount of
+work that has to be done in each iteration of the regular parametric
+push-relabel algorithm.
+
+See \cite{gallo1989fast}, \cite{hochbaum2001efficient} and
+\cite{goldfarb2009parametric} for more information.
+
+\section{Boykov--Kolmogorov algorithm}
+\label{sec:bk}
+
+\begin{algorithm}
+ \begin{algorithmic}
+ \Function{BoykovKolmogorov}{$G$, $s$, $t$}
+ \State $A := {s, t}$, $O := \emptyset$, $S = {s}$, $T = {t}$
+ \State $s\text{.color} = S$, $t\text{.color} = T$
+ \While{\texttt{True}}
+ \State $e \gets$ \Call{Grow}{$G$, $A$}
+ \If{not $e$}
+ \State break
+ \EndIf
+ \State \Call{Augment}{$G$, $s$, $t$, $e$, $O$}
+ \State \Call{Adopt}{$G$, $s$, $t$, $O$}
+ \EndWhile
+ \EndFunction
+ \end{algorithmic}
+ \caption{\sf The Boykov--Kolmogorov maximum flow algorithm}
+ \label{alg:boykov_kolmogorov}
+\end{algorithm}
+
+A maximum flow algorithm specialized for the type of graphs found in
+imaging applications is described by Boykov and Kolmogorov in
+\cite{boykov2004experimental}. Their algorithm is an augmenting path
+algorithm where the paths are found using trees that keep track of
+possible paths.
+
+There is one tree denoted $S$ with the source $s$ as its root, and
+another tree denoted $T$ with the sink $t$ as its root. These trees are
+non-overlapping, and all edges in the tree $S$ can carry flow towards
+the leafs of the tree, while all edges in $T$ can carry flow towards the
+sink $t$.
+
+A vertex can either be part of these trees, or be a free vertex, and we
+write
+\begin{equation}
+ u\text{.color} = \begin{cases}
+ S & \text{ if } u \in S \\
+ T & \text{ if } u \in T \\
+ \texttt{Free} & \text{ otherwise.}
+ \end{cases}
+\end{equation}
+The idea of the algorithm is then to grow these trees towards eachother,
+and when they touch, we have found an augmenting path from $s$ to $t$.
+As much flow as possible is then pushed along this path, and some edges
+will become saturated. If an edge in the tree $S$ or $T$ becomes
+saturated, the tree is split up into multiple trees, also called a
+forest. The parts no longer connected to the root are called orphans,
+and we write $u \in O$ for an orphan vertex $u$. These vertices are---if
+possible---joined back into the tree, or become free nodes.
+
+The vertices in the trees $S$ and $T$ are either active or passive, and
+we write $u \in A$ if $u$ is active. The active vertices are those at
+the boundary of the tree, which can possibly connect to other vertices
+to grow the tree. The passive vertices are internal in the tree, and
+edges to their neighbors are either a part of the tree, or completely
+saturated.
+
+The main loop of the algorithm is then as shown in
+Algorithm~\ref{alg:boykov_kolmogorov}. In the grow procedure, the trees are grown
+from their active vertices until a new path is found. The flow
+is then augmented along this path in the augment procedure, and orphan
+vertices might be created. These orphan vertices are then either joined
+back into their respective trees, or become free, in the adopt
+procedure.
+
+\begin{algorithm}
+\begin{algorithmic}
+ \Function{Grow}{$G$, $A$}
+ \While{$\abs{A} \neq 0$}
+ \State $u \gets $ one node from $A$
+ \ForAll{$v$ such that \Call{treeCap}{$u$, $v$} > 0}
+ \If{$v\text{.color} = \texttt{Free}$}
+ \State $v\text{.color} \gets u\text{.color}$
+ \State $v\text{.parent} \gets u$
+ \State $A := A \cup \{v\}$
+ \ElsIf{$v\text{.color} \neq u\text{.color}$}
+ \State \Return $e := (u, v)$
+ \EndIf
+ \EndFor
+ \State Remove $u$ from $A$
+ \EndWhile
+ \State \Return \texttt{False}
+ \EndFunction
+\end{algorithmic}
+\caption{\sf The grow procedure of the Boykov--Kolmogorov maximum flow
+algorithm}
+\label{alg:grow}
+\end{algorithm}
+
+The grow procedure is shown in Algorithm~\ref{alg:grow}. It goes through
+the set of active vertices $A$ to try to expand the trees $S$ and $T$.
+When considering an active vertex $u \in S$, we want to grow by finding
+vertices $v$ such that $c_f(u, v) > 0$, while when considering an active
+vertex $u \in T$, we want to grow by finding vertices $v$ such that
+$c_f(v, u) > 0$. This is why the \textsc{treeCap} is introduced which is
+defined as
+\begin{equation}
+ \textsc{treeCap}(u, v) = \begin{cases}
+ c_f(u, v) & \text{ if } u \in S, \\
+ c_f(v, u) & \text{ if } u \in T.
+ \end{cases}
+\end{equation}
+Thus if a non-saturated edge is found from an active vertex $u$ to a free
+vertex $v$, then $v$ is added to the tree of $u$. The tree structure is
+stored by keeping a parent attribute in each non-free node, such that
+for example $v\text{.parent} = u$.
+
+If a non-saturated edge is found from an active vertex $u$ to a vertex
+$v$ in the other tree, the two trees connect, and we can return the edge
+$e$ since we have an augmenting path. Note that if this happens, the
+vertex $u$ is still active. It only becomes passive when all neighbors
+are considered without yielding an augmenting path.
+
+\begin{algorithm}
+ \begin{algorithmic}
+ \Function{Augment}{$G$, $s$, $t$, $e$, $O$}
+ \State $p \gets $ path from $s$ to $t$ through $e$
+ \Comment through the trees $S$ and $T$
+ \State $\Delta f \gets \min \{c_f(u, v) : (u,v) \in p\}$
+ \ForAll{$(u, v) \in p$}
+ \State $f(u, v) \mathrel{+}= \Delta f$
+ \State $f(v, u) \mathrel{-}= \Delta f$
+ \If{$c_f(u, v) = 0$}
+ \If{$u\text{.color} = v\text{.color} = S$}
+ \State $v\text{.parent} = \texttt{Null}$
+ \State $O := O \cup \{v\}$
+ \ElsIf{$u\text{.color} = v\text{.color} = T$}
+ \State $u\text{.parent} = \texttt{Null}$
+ \State $O := O \cup \{u\}$
+ \EndIf
+ \EndIf
+ \EndFor
+ \EndFunction
+ \end{algorithmic}
+ \caption{\sf The augment procedure of the Boykov--Kolmogorov maximum
+ flow algorithm}
+ \label{alg:augment}
+\end{algorithm}
+
+The augment procedure finds the augmenting path going through the tree
+$S$, the edge $e$ and the tree $T$. The maximal possible flow is then
+sent along this path. At least one edge will then become saturated. If a
+saturated edge occurs in the trees $S$ or $T$, the edge terminal
+farthest from the tree root is marked as an orphan. Note that even if a
+whole subtree is disconnected from the main tree, only the root of this
+subtree is marked as an orphan. The adopt procedure will take care of
+either reconnecting, or freeing all the vertices in the disconnected
+subtree.
+
+\begin{algorithm}
+\begin{algorithmic}
+ \Function{Adopt}{$G$, $s$, $t$, $O$}
+ \While{$\abs{O} \neq 0$}
+ \State $u \gets $ one node from $O$
+ \State Remove $u$ from $O$
+ \State found $\gets \texttt{False}$
+ \ForAll{$v$ such that \Call{treeCap}{$v$, $u$} > 0}
+ \If{$u\text{.color} \neq v\text{.color}$}
+ \State continue
+ \EndIf
+ \State orig $\gets$ \Call{treeOrigin}{$v$}
+ \If{orig $\neq s$ and orig $\neq t$}
+ \State continue
+ \EndIf
+ \State found $\gets \texttt{True}$
+ \State $u\text{.parent} \gets v$
+ \State break
+ \EndFor
+ \If{found $\neq \texttt{True}$}
+ \ForAll{$v$ such that $u\text{.color} = v\text{.color}$}
+ \If{\Call{treeCap}{$v$, $u$} > 0}
+ \State $A := A \cup \{ v \}$
+ \EndIf
+ \If{$v\text{.parent} = u$}
+ \State $O := O \cup \{ v \}$
+ \State $v\text{.parent} \gets \texttt{Null}$
+ \EndIf
+ \State $u\text{.color} \gets \texttt{Free}$
+ \State $A := A - \{ u \}$
+ \EndFor
+ \EndIf
+ \EndWhile
+ \EndFunction
+\end{algorithmic}
+\caption{\sf The adopt procedure of the Boykov--Kolmogorov maximum flow
+algorithm}
+\label{alg:adopt}
+\end{algorithm}
+
+The adopt procedure processes all the vertices in the set of orphans
+$O$. These vertices are either single vertices, or roots of disconnected
+subtrees. For an orphaned vertex $u$ we look through its neighbors to
+find a possible parent vertex $v$. It should belong to the same tree as
+$u$, and also satisfy $\textsc{treeCap}(v, u) > 0$. In addition, $v$
+should be connected to one of the tree roots $s$ or $t$. This keeps us
+from connecting to other orphaned vertices or subtrees, and is checked
+by the procedure $\textsc{treeOrigin}(v)$, which follows the parent
+information until reaching either $s$, $t$ or an orphaned vertex.
+
+If a potential parent $v$ of $u$ is \emph{not} found, $u$ becomes a free
+vertex. All vertices that had $u$ as its parent are orphaned, and are thus
+treeted by the adopt procedure later. Vertices that are in one of the
+trees $S$ or $T$, and have a non-saturated edge to this newly freed
+vertex $u$, become active.
+
+When no more active vertices remain, there are no more possible
+augmenting paths, and the algorithms terminates with a maximum flow.
+A proof of correctness can be found in \fixme{kolmogorov-phd}.
+
+\subsection{Graph reuse}
+
+As in the push-relabel algorithm described earlier, it is not necessary
+to completely restart the Boykov--Kolmogorov algorithm for every level
+of the image, but we can reuse the graph and the trees $S$ and $T$ in
+successive runs.
+
+\begin{figure}
+ \input{fig/bk_norm_subgraph}
+\end{figure}
+
+Since this is algorithm does not deal with the relaxed preflow concept,
+the updates from one level to the next has to be done in a different
+way. We can no longer decrease the capacity of edges, as this could
+break the flow conservation constraint. However, the graph construction
+in Section~\ref{sec:graph_representable} did allow for the addition of an arbitrary constant to
+all the edges of the sub-graph. Thus the construction in
+Figure~\ref{fig:bk_norm_subgraph} is also valid. We have added the
+non-negative constant $\max\{F_{L-1}^x(1), 0\}$ to all the edges of the
+graph in Figure~\ref{fig:norm_subgraph_neg}. As $\lambda$ goes from $L-1$ to $0$ both edges
+stay non-negative. The edge from $u_\lambda^x$ to $t$ is non-decreasing
+with decreasing $\lambda$ parameter.
+
+This construction allows us to update the capacities of the edges, while
+retaining the flow and the trees $S$ and $T$. A nice property that
+follows is that the partition $(V - T, T)$ is after each run a minimum
+cut since there are is no more residual capacity between the trees $S$
+and $T$. Further, the set of vertices $T$ is non-decreasing, or more
+importantly, no vertex already
+in the tree $T$ will leave $T$ for succeeding lower $\lambda$ values.
+This can be seen from the fact that if $u$ is in $T$ and the algorithm
+has terminated, all paths from $s$ to $t$ going
+through $u$ has a saturated edge somewhere before $u$. For all these
+paths, the capacity will only change for edges after $u$, as we only
+change the capacity of edges $(v, t)$. Thus these saturated edges will
+stay saturated, and $u$ will stay in $T$. \fixme{shaky}
+
+\subsection{Performance improvements}
+
+There are several open choices in the implementation of the algorithm,
+for example the order in which active vertices are processed. As
+recommended in \cite{boykov2004experimental}, we have implemented a
+«First-In-First-Out» queue. This ensures that at least the first
+augmenting path found is a shortest path, although later the distance
+information is lost in the adoption stage.
+
+In the adoption stage it is possible, and perhaps preferrable to seek
+possible parent that is closest to the root of the tree, and adopt that
+vertex as parent, instead of the first possible.
+
+Because of the particular graph construction, all vertices except $s$
+and $t$ are connected to $s$ and $t$. Thus there are as many two-edged
+paths from $s$ to $t$ as there are pixels in the image. When increasing
+the capacity of edges $(v, t)$, a quick sweep over these two-edged paths
+to send any possible flow may speed up the algorithm.
+
+\section{Implementation}
+
+A \cpp{} implementation can be found in
+Appendix~\ref{app:c++implementation}. It uses the open computer vision
+library OpenCV \cite{opencv_library} to load and save image files.
+
+Note that when implementing maximum flow algorithms it is not a good
+idea, memory- and performance-wise, to actually construct the residual
+graph $G_f$. Instead, every time we update the flow $f(u,v)$ we set
+the flow in the opposite direction to its negative value $f(v,u) =
+-f(u,v)$. Then we can at any time, consider the value $c(u,v) - f(u,v)$
+in the place of the residual capacity $c_f(u,v)$.
+
+For the gap relabeling heuristic, we need to have a easy way of finding
+when a gap occurs. This is done by keeping track of how many vertices
+exist with each label.
+
+\fixme{something on bk?}
+\fixme{move this to the appendix possibly}
+
+\section{Performance comparison}
+
+Maybe.
+
\label{fig:finger_contrast_atv2}
\end{subfigure}
\caption[Comparison contrast loss in regular an anisotropic restoration]{%
- \fixme{nope}
A noisy fingerprint restored using both isotropic and
- anisotropic total variation. All parameters are kept constant
- except for $\omega$, the anisotropy parameter.
+ anisotropic total variation to see how the contrast loss
+ compares. Parameters: $\abs{\mathcal{N}} = 32$, $\sigma = 3$,
+ $\rho = 10$.
}
\label{fig:finger_contrast}
\end{figure}
applications the input images have different properties, and one might
also have different hopes for, and restrictions on the output image.
Thus I hope a thorough description in addition to an inspection of
-restoration results can give insight into the good and not-so-good sides
+restoration results has given insight into the good and not-so-good sides
of this particular method.
-\chapter{Continuous formulation}
-
-From the introduction we see that there are many different approaches to
-the image restoration problem, all with their own strengths and
-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 respective weaknesses.
-
-This chapter will be devoted to the continuous formulation of the
-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 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}
-
-\fixme{Rating: 8/10}
-
-The method considered will build on the total variation regularization
-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 dependent. We introduce
-the anisotropic total variation
-\begin{equation}
- \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)$. 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 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$ 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 = \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 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^*
- \leq 1 \right\}.
- \label{eq:extended_tv}
- \end{equation}
- \label{def:extended_tv}
-\end{definition}
-\nomenclature{$\TVA(u)$}{Anisotropic total variation of image $u$, with
-anisotropy tensor $A(u)$.}%
-
-With this extended definition, we have arrived at a minimization problem
-where we seek to find a minimizer of the functional
-\begin{equation}
- F(u) = \int_\Omega (u - f)^2 \, dx + \beta \, \TVA(u).
- \label{eq:first_anisotropic_functional}
-\end{equation}
-Similar functionals have been considered in
-\cite{grasmair2010anisotropic} and \cite{olsson2009extending}. The
-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}
-
-\fixme{Rating: 6/10}
-
-There are many possible choices for the anisotropy tensor $A(x)$. Our
-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
-none, so we somehow want to be sure about the edges we find.
-
-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$.
-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 \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)(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{}
- \includegraphics[width=0.4\textwidth]{fig/edges.png}
- \caption{%
- The largest eigenvalue of the structure tensor is $\abs{\nabla
- f_\sigma(x)}^2$ and can as we see function as an edge detector.
- }
- \label{fig:edges}
-\end{figure}
-
-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} -
- s_{22})^2 + 4 s_{12}^2} \right).
- \label{eq:s_eigenvalues}
-\end{equation}
-
-The vector
-$v_1$ will then indicate the direction of most variation in the
-neighborhood. From this we see that an edge would give $\lambda_1 \gg
-\lambda_2 \approx 0$, while smooth areas would give $\lambda_1 \approx
-\lambda_2 \approx 0$. In corners we have variation in the direction of
-$v_1$ but also perpendicular to $v_1$, so we will have $\lambda_1
-\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.
-
-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
-\begin{equation}
- \Lambda(x) = \begin{pmatrix}
- \lambda_1 & 0 \\
- 0 & \lambda_2
- \end{pmatrix}
-\end{equation}
-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}.
- \label{eq:sigma_construction}
-\end{equation}
-and for $\sigma_1$ and $\sigma_2$ we choose
-\begin{equation}
- \begin{aligned}
- \sigma_1 &= \left(1 +
- \frac{(\lambda_1 - \lambda_2)^2}{\omega^2}\right)^{-1}, \\
- \sigma_2 &= 1.
- \end{aligned}
- \label{eq:sigma_construction}
-\end{equation}
-Thus the eigenvectors of $A(x)$ and $S_\rho(x)$ are equal, while the
-eigenvalues are different. A visualization of the two tensors can be
-seen in Figure~\ref{fig:tensor_viz} where the two tensors are shown at
-an edge in the image.
-
-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.
-
-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.
-
-The parameter $\omega$ controls the amount of anisotropy in the method,
-such that if it is very large we are left with the identity matrix. Note
-that changing the parameter $\omega$ implicitly affects how much
-restoration is done. For an image $u$, decreasing $\omega$ will, all
-else being equal, decrease the lowest eigenvalue of $A(x)$ and in turn
-decrease the anisotropic total variation $\TVA(u)$.
-
-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)$.
-
-See \cite{weickert1999coherence} for a different tensor construction,
-made to enhance flow structures in the image, relevant in for example
-fingerprint analysis.
-
-\fixme{argue that $A$ is continuous?}
-
-\begin{figure}
- \input{fig/tensor_viz}
-\end{figure}
-
-\section{Well-posedness}
-
-\fixme{Rating: 7/10}
-
-The theory of existence and uniqueness for these kinds of variational
-methods is a minefield of more or less subtle problems. Even if we
-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 others will be skipped
-with a reference to further theory.
-
-The basic things we ask of our functional
-\begin{equation}
- F(u) = \int_\Omega (u - f)^2 + \beta \, \TVA(u)
-\end{equation}
-to have a \emph{well-posed} problem are lower semicontinuity
-and coercivity for existence, convexity for uniqueness and stability. We restrict
-ourself to $L^2(\Omega)$ which makes sense with our fidelity term,
-assuming that $f \in L^2(\Omega)$ initially.
-
-We consider the weak topology, as it will allow us to arrive at an
-existence result relatively easily.
-We say that a sequence $f_n$ in $L^2(\Omega)$ converges weakly to $f$ if
-\begin{equation}
- \lim_{n \to \infty} \int_\Omega f_n \, \xi \, dx = \int_\Omega f \,
- \xi \, dx
-\end{equation}
-for all $\xi \in L^2(\Omega)$ and we write $f_n \rightharpoonup f$. A
-weakly convergent sequence is a sequence that converges in the weak topology.
-
-\subsection{Convexity}
-
-We start with convexity as it is the easiest to show. Being quadratic,
-the fidelity term of our functional
-\begin{equation}
- \int_\Omega (u - f)^2 \, dx
-\end{equation}
-is obviously strictly convex. This can be shown by expanding and
-rearranging the strict convexity condition
-\begin{equation}
- \int_\Omega (\lambda u_1 + (1-\lambda)u_2 - f)^2 \, dx < \lambda
- \int_\Omega (u_1 - f)^2 \, dx + (1 - \lambda) \int_\Omega (u_2 - f)^2 \,
- dx
-\end{equation}
-to obtain that it is equivalent to
-\begin{equation}
- - \lambda(1 - \lambda) \int_\Omega (u_1 - u_2)^2 \, dx < 0
-\end{equation}
-which is true for $0 < \lambda < 1$ and $u_1 \neq u_2$.
-
-The anisotropic total variation
-\begin{equation}
- \TVA(u) = \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega u \diver \xi \,
- dx
-\end{equation}
-can be thought of as---and has the properties of---a semi-norm, and is
-therefore convex. The sum of the two is thus strictly convex, which,
-given the existence of a minimizer, implies uniqueness.
-
-\subsection{Coercivity}
-
-Coercivity relates to how the functional behaves when the norm of the
-image $u$ tends to infinity. What we need in order to conclude with existence
-is sequential coercivity. Thus we need all level sets $F^\alpha = \{ u \in
-L^2(\Omega) : F(u) \leq \alpha \}$ to be \emph{sequentially
-pre-compact}, meaning that all sequences in the set contain a
-subsequence converging to an element of the closure of the set.
-
-It is obvious from the fidelity term that for some fixed $f \in
-L^2(\Omega)$, if $\lnorm{u}_{L^2} \to \infty$ then $F(u) \to \infty$.
-This implies that all the level sets $F^\alpha$ are bounded.
-Since $L^2(\Omega)$ is a Hilbert space, all bounded sequences contain a
-weakly convergent subsequence. Thus all the level sets $F^\alpha$ are
-weakly sequentially pre-compact, since all sequences in these sets have
-a subsequence weakly converging to a point in the closure of the set.
-
-\begin{figure}
- \input{fig/lower_semicont}
-\end{figure}
-
-\subsection{Lower semi-continuity}
-
-The lower semicontinuity is the most tricky part, and this is where we
-will take some shortcuts. Lower semicontinuity for a functional $F$ at a
-point $u$ means that at points $u_\epsilon$ close to $u$, the functional
-takes values either close to or above $F(u)$. More specifically, for
-every sequence $u_k$ converging to $u$, we have $F(u) \leq \liminf_k
-F(u_k)$. For a function $f : \mathbb{R} \to \mathbb{R}$ this can be
-visualized as in Figure~\ref{fig:lower_semicont}.
-
-Since our space $L^2(\Omega)$ is of infinite dimensions things become a
-little problematic here. The problem lies in the fact that a
-functional which is continuous with respect to sequences is not
-necessarily continuous with respect to the underlying topology. In
-other words, in these spaces, there can be a difference between sequential
-continuity and topological continuity. Topological continuity implies
-sequential continuity, but the converse does not hold. One way to get
-around this would be to consider topological \emph{nets}, an extension
-of sequences, but for simplicity, and because it might not add much to
-the understanding of the restoration method, we will stick to proving
-sequential lower semi-continuity and referring to further theory. For
-further reading on the theory of sequential versus topological
-continuity see for example Megginson's book on Banach space theory
-\cite{megginson}.
-
-The mapping $u
-\mapsto \int_\Omega u \, \xi \, dx$ is weakly continuous for all $\xi
-\in L^2(\Omega)$. Note that when we write weakly continuous it is not a
-weaker version of continuity, but rather continuity in the weak
-topology, and the same goes for weak lower semi-continuity.
-
-Before arguing that our own functional is sequentially weakly lower
-semi-continuous, we present a needed result.
-\begin{lemma}
- Assume that the functional $F : L^2(\Omega) \to \mathbb{R}$ is
- defined by
- \begin{equation}
- F = \sup_i F_i
- \end{equation}
- where all the $F_i$ are sequentially weakly lower semi-continuous, then $F$
- is sequentially weakly lower semi-continuous, meaning that for any sequence
- $u_k \rightharpoonup u$ we have $F(u) \leq \liminf_k F(u_k)$.
- \label{lem:sup_semi_cont}
-\end{lemma}
-\begin{proof}
- For any sequence $u_k \rightharpoonup u$ in $L^2(\Omega)$ we have
- \begin{equation}
- F(u) = \sup_i F_i(u) \leq \sup_i \liminf_{k \to \infty} F_i(u_k)
- \end{equation}
- from the sequential weak lower semi-continuity of $F_i$. Using that
- $\liminf_{k \to \infty} u_k = \sup_k \inf_{l \geq k} u_l$, we obtain
- \begin{equation}
- \begin{aligned}
- F(u) &\leq \sup_i \sup_k \inf_{l \geq k} F_i(u_l) \\
- &= \sup_k \sup_i \inf_{l \geq k} F_i(u_l) \\
- &\leq \sup_k \inf_{l \geq k} \sup_i F_i(u_l) \\
- &= \liminf_{k \to \infty} F(u_k)
- \end{aligned}
- \end{equation}
- which proves that $F$ is sequentially weakly lower semi-continuous.
-\end{proof}
-
-In our functional in \eqref{eq:first_anisotropic_functional}, we first
-consider the fidelity term, and rewrite it as a supremum
-\begin{equation}
- \int_\Omega (u - f)^2 \, dx
- %= \sup_{\substack{\xi \in L^2(\Omega) \\ \norm{\xi}_{L^2} \leq
- %\norm{u-v}_{L^2}}} \int_\Omega (u - v) \, \xi \, dx.
- = \sup \left\{\int_\Omega (u - f) \, \xi \, dx : \xi \in
- L^2(\Omega), \abs{\xi(x)} \leq \abs{u(x)-f(x)} \right\}
-\end{equation}
-As the map $u \mapsto \int_\Omega (u - v) \xi\, dx$ is continuous in the
-weak topology, the fidelity term is thus a supremum of weakly continuous
-functionals, and is thus by Lemma~\ref{lem:sup_semi_cont} sequentially
-lower semi-continuous.
-
-For the regularization term the approach is similar. With our extended
-definition from \eqref{eq:extended_tv}, we have
-\begin{equation}
- %\TVA(u) = \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega u \, \diver \xi \,
- %dx
- \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\}
-\end{equation}
-This is again a
-supremum of weakly continuous functionals. Thus the regularization term
-is by Lemma~\ref{lem:sup_semi_cont} also sequentially weakly lower
-semi-continuous.
-
-The sum of the two terms is trivially sequentially weakly lower
-semi-continuous functional since
-\begin{equation}
- \begin{aligned}
- F_1(u) + F_2(u) &\leq \liminf_{k \to \infty} F_1(u_k) + \liminf_{k
- \to \infty} F_2(u_k) \\
- &= \lim_{k \to \infty} \left( \inf_{l \geq k}
- F_1(u_l) + \inf_{l \geq k} F_2(u_l) \right) \\
- &\leq \liminf_{k \to
- \infty} \left( F_1(u_k) + F_2(u_k) \right),
- \end{aligned}
-\end{equation}
-and thus our functional is sequentially weakly lower semi-continuous.
-
-The usual ways of going from coercivity and lower semicontinuity to
-existence do not work in infinite dimensions. But with sequential
-coercivity and sequential lower semi-continuity in the weak topology we
-can conclude that we have existence from \cite[Theorem
-5.1]{scherzer2008variational}.
-
-\subsection{Stability}
-
-The last property normally required for well-posedness is stability.
-Thus the question is how the solutions behave when we perturb the
-problem. We stick to $L^2(\Omega)$ and adopt a proof given by Acar and
-Vogel in \cite{acar1994analysis}.
-
-\begin{proposition}
- Consider a sequence of perturbed problems
- \begin{equation}
- \min_{u \in L^2(\Omega)} F_n(u).
- \end{equation}
- Assume that each of the functionals $F_n$ is weakly sequentially
- lower-semicontinuous, has a unique minimizer $u_n$, and in addition that they are
- uniformly coercive in the sense that for any sequence $v_n \in
- L^2(\Omega)$
- \begin{equation}
- \lim_{n \to \infty} F_n(v_n) = \infty \quad \text{whenever } \, \lim_{n \to
- \infty} \norm{v_n}_{L^2} = \infty.
- \label{eq:fn_coercive}
- \end{equation}
- Assume also that the problems are consistent in the sense that $F_n \to
- F$ uniformly, i.e.\ given an $R > 0$ and $\epsilon > 0$ there exists an
- $N$ such that
- \begin{equation}
- \abs{F_n(u) - F(u)} < \epsilon \quad \text{whenever } \, n \geq N \text{ and
- } \norm{u}_{L^2} \leq R.
- \label{eq:fn_concistent}
- \end{equation}
- The problem is then stable with respect to the perturbations $F_n$ in
- the sense that if $u^*$ is a minimizer of $F$ we have
- \begin{equation}
- u_n \rightharpoonup u^*.
- \label{eq:un_u_weak}
- \end{equation}
-\end{proposition}
-\begin{proof}
- As $u_n$ is a minimizer of $F_n$ we have
- \begin{equation}
- F_n(u_n) \leq F_n(u^*)
- \end{equation}
- and using \eqref{eq:fn_concistent} we get
- \begin{equation}
- %\lim F_n(u_n) \leq \limsup F_n(u_n) \leq \limsup F_n(u^*)
- %= F(u^*) < \infty.
- \lim_{n \to \infty} F_n(u_n) \leq \lim_{n \to \infty} F_n(u^*) =
- F(u^*) < \infty.
- \end{equation}
- From the coercivity in \eqref{eq:fn_coercive} we obtain that
- the sequence $u_n$ is bounded in $L^2(\Omega)$ and thus contains a
- weakly convergent subsequence.
-
- Assume that the convergence in \eqref{eq:un_u_weak} does not hold,
- and denote the weakly convergent subsequence by $u_{n_j}
- \rightharpoonup \bar{u} \neq u^*$. By the weak sequential lower
- semicontinuity we have
- \begin{equation}
- \begin{aligned}
- F(\bar{u}) &\leq \liminf_{j \to \infty} F\left(u_{n_j}\right) \\
- &= \liminf_{j \to \infty} \left(F\left(u_{n_j}\right) -
- F_{n_j}\left(u_{n_j}\right)\right) + \liminf_{j \to \infty}
- F_{n_j}\left(u_{n_j}\right) \\
- &= \liminf_{j \to \infty} F_{n_j}\left(u_{n_j}\right) \\
- &\leq F(u^*)
- \end{aligned}
- \end{equation}
- which contradicts the uniqueness of $u^*$ and thus we can conclude
- that the solutions of the perturbed problems converge $u_n
- \rightharpoonup \bar{u} = u^*$.
-\end{proof}
-
-Note that these are properties of the continuous problem, and that the
-discretizations and numerical methods used for approximating a solution
-can have its own issues with consistency, convergence and stability.
-
-\section{Anisotropic coarea formula}
-
-\fixme{Rating: 6.5/10, especially the part after the proof needs some
-work.}
-
-\fixme{More flow between the sections. Somehow.}
-
-The anisotropic coarea formula we will present here allows us to write
-the anisotropic total variation as an integral over the levels of the
-image. For a similar presentation of the regular coarea formula for all
-$f \in \BV(\Omega)$ see \cite{evans1991measure}.
-
-First we define the thresholded image at level $s$.
-\begin{definition}[Thresholded image]
- The thresholded image at level $s$ is the function
- \begin{equation}
- u^s(x) = \begin{cases}
- 1 & \text{if } u(x) > s, \\
- 0 & \text{otherwise.}
- \end{cases}
- \end{equation}
- \label{def:thresholded_image}
-\end{definition}
-This will be used throughout the rest of the thesis. Note that given the
-thresholded image for every level, we are able to reconstruct the image
-as
-\begin{equation}
- u(x) = \sup \left\{ s : u^s(x) = 1 \right\}.
-\end{equation}
-The thresholded image definition also allows us to write a non-negative
-image $u \geq 0$ as an integral over all the layers
-\begin{equation}
- u(x) = \int_0^\infty u^s(x) \, ds.
- \label{eq:positive_int}
-\end{equation}
-Note that \eqref{eq:positive_int} only holds for non-negative images,
-which complicates the following proof a little.
-
-\begin{figure}
- \input{fig/eta_r}
-\end{figure}
-
-\begin{theorem}[Anisotropic coarea formula]
- Given an image $u \in \BV(\Omega)$, the anisotropic total variation
- can be written as an integral over all the levels
- \begin{equation}
- \TVA(u) = \int_{-\infty}^\infty \TVA(u^s) \, ds.
- \label{eq:anisotropic_coarea}
- \end{equation}
- \label{thm:anisotropic_coarea}
-\end{theorem}
-For the proof we will avoid measure theory and follow a proof given
-in \cite{olsson2009extending}, but first we will present a necessary
-result from measure theory.
-\begin{theorem}[Lebesgue's Dominated Convergence theorem]
- Let $\{ f_n \}$ be a sequence of real-valued measurable functions on
- a space $S$ with measure $d\mu$ which converges almost
- everywhere to a real-valued measurable function $f$. If there exists
- an integrable function $g$ such that $\abs{f_n} \leq g$ for all $n$,
- then $f$ is integrable and
- \begin{equation}
- \lim_{n \to \infty} \int_S f_n \, d\mu = \int_S f \, d\mu.
- \end{equation}
-\end{theorem}
-For a proof and further background on measure theory and Lebesgue
-integration theory see for example \cite{bartle1995elements}.
-\begin{proof}[Proof of the anisotropic coarea formula.]
- Assume that $u \in C^1(\Omega) \cap \BV(\Omega)$. The extension to
- all functions $u \in \BV(\Omega)$ will not be considered here, but
- for the case of regular total variation see \cite[Theorem
- 5.3.3]{ziemer1989}.
-
- %\paragraph{First we prove that $\TVA(u) \leq \int_{-\infty}^\infty
- %\TVA(u^s) \, ds$.}
- \paragraph{Proof of upper bound.}
- Assume that $u \geq 0$ such that the integral representation in
- \eqref{eq:positive_int} holds, then inserting
- \eqref{eq:positive_int} into the extended
- total variation definition in \eqref{eq:extended_tv} gives
- \begin{equation}
- \begin{aligned}
- \TVA(u) &= \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega
- \left( \int_0^\infty u^s ds \right) \diver \xi \, dx
- = \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega \int_0^\infty u^s
- \diver \xi \, ds \, dx \\
- &\leq \int_0^\infty
- \left( \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega u^s \diver
- \xi \, dx \right) \, ds
- = \int_0^\infty \TVA(u^s) \, ds.
- \end{aligned}
- \label{eq:tv_u_pos}
- \end{equation}
- For $u \leq 0$ we use that $\TVA(-v) = \TVA(v)$ and that $\TVA(c + v) =
- \TVA(v)$ for any constant $c$. Note that $-u \geq 0$ and that its
- thresholded image $(-u)^s$ will be exactly the opposite of $u^{-s}$,
- that is $(-u)^s = 1 - u^{-s}$. This allows us to show that
- \begin{equation}
- \begin{aligned}
- \TVA(u) &= \TVA(-u) \leq \int_0^\infty \TVA \big( (-u)^r
- \big) \, dr
- = \int_0^\infty \TVA(1 - u^{-r}) \, dr \\ &= \int_0^\infty
- \TVA(u^{-r}) \, dr = \int_{-\infty}^0 \TVA(u^s) \, ds.
- \end{aligned}
- \label{eq:tv_u_neg}
- \end{equation}
- Following from the supremum definition of the anisotropic total
- variation in \eqref{eq:extended_tv}, we obtain the inequality
- \begin{equation}
- \begin{aligned}
- \TVA(u_1 + u_2) &= \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega (u_1 +
- u_2) \diver \xi \, dx \\
- &\leq \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega u_1
- \diver \xi \, dx + \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega u_2
- \diver \xi \, dx\\
- &= \TVA(u_1) + \TVA(u_2).
- \end{aligned}
- \label{eq:tv_sum}
- \end{equation}
- Next, we write a general $u$ as a difference of two positive
- functions $u = u_+ - u_-$ where $u_+ = \max\{u,0\}$ and $u_- =
- -\min\{u,0\}$.
- Inserting \eqref{eq:tv_u_pos} and
- \eqref{eq:tv_u_neg} into \eqref{eq:tv_sum} we obtain
- \begin{equation}
- \begin{aligned}
- \TVA(u) &\leq \TVA(u_-) + \TVA(u_+) = \TVA(-u_-) + \TVA(u_+) \\
- &\leq \int_{-\infty}^0 \TVA\big((-u_-)^s\big) \, ds +
- \int_0^\infty \TVA(u_+^s) \, ds \\
- &= \int_{-\infty}^0 \TVA(u^s) \, ds + \int_0^\infty \TVA(u^s) \,
- ds = \int_{-\infty}^\infty \TVA(u^s) \, ds.
- \end{aligned}
- \end{equation}
- Note that $u_+$ and $u_-$ will not be differentiable everywhere, but
- we did not use the differentiability of $u$ in this part of the
- proof.
-
- %\paragraph{Next we prove that $\TVA(u) \geq \int_{-\infty}^\infty
- %\TVA(u^s) \, ds$.}
- \paragraph{Proof of lower bound.}
- Define the function
- \begin{equation}
- m(t) = \int_{\{ x \in \Omega : u(x) \leq t\}} \norm{\nabla u}_A
- \, dx,
- \label{eq:mdef}
- \end{equation}
- and note that $m(\infty) = \TVA(u)$ and $m(-\infty) = 0$. Since
- $m(t)$ is non-decreasing with $t$, we can apply the existence
- theorems of Lebesgue \cite[Thm.\ 17.12, 18.14]{hewstrom} to conclude
- that $m\prime(t)$ exists almost everywhere and that the following
- inequality holds:
- \begin{equation}
- \int_{-\infty}^\infty m\prime(t)\, dt \leq m(\infty) - m(-\infty) =
- \TVA(u).
- \label{eq:tva_geq_mder}
- \end{equation}
- Next, fix an $s \in \mathbb{R}$ and define the cut-off function
- \begin{equation}
- \begin{aligned}
- \eta_r(t) = \begin{cases}
- 0 & \text{if } t < s, \\
- (t - s)/r & \text{if } s \leq t < s + r, \\
- 1 & \text{if } t \geq s + r,
- \end{cases}
- & \quad
- \eta_r\prime(t) = \begin{cases}
- 0 & \text{if } t < s, \\
- 1 & \text{if } s < t < s + r, \\
- 0 & \text{if } t > s + r,
- \end{cases}
- \end{aligned}
- \end{equation}
- visualized in Figure~\ref{fig:eta_r} and \ref{fig:eta_r_diff}. By
- composing the function $\eta_r$ with our image $u$ and using Green's
- formula, for example from \cite[Corollary
- 9.32]{grasmair2010anisotropic} we obtain
- \begin{equation}
- \int_\Omega - \eta_r(u) \diver \xi \, dx
- = \int_\Omega \eta_r\prime(u) \nabla u\cdot \xi \, dx
- = \frac{1}{r} \int_{\{ s < u < s + r \}} \nabla u\cdot \xi
- \, dx,
- \label{eq:eta_greens}
- \end{equation}
- for all vector fields $\xi \in C_c^\infty(\Omega, \mathbb{R}^2)$.
- The measure of $\{ x : u(x) = \lambda \text{ and } \nabla u(x) \neq
- 0\}$ is zero for all $\lambda$ following from \cite[Corollary \rom{1},
- Section 3.1.2]{evans1991measure}. Thus we can ignore the sets $\{ u
- = s\}$ and $\{ u = s + r \}$.
- Assuming that $\norm{\xi}_A^* \leq 1$ we obtain from \eqref{eq:mdef}
- and \eqref{eq:eta_greens} that
- \begin{equation}
- \begin{aligned}
- \frac{m(s+r) - m(s)}{r}
- &= \frac{1}{r} \int_{\{ s < u \leq s+r \}} \norm{\nabla u}_A
- \, dx \\
- &\geq \frac{1}{r} \int_{\{ s < u \leq s + r\}} \nabla u \cdot
- \xi \, dx
- = \int_\Omega -\eta_r(u) \diver \xi \, dx.
- \end{aligned}
- \label{eq:m_ineq_sr}
- \end{equation}
- As the limit when $r \to 0$ of the left-hand side exists almost
- everywhere, suppose it exists at $s \in \mathbb{R}$. The integrand
- on the right-hand side $-\eta_r(u) \diver \xi$ approaches $-u^s
- \diver \xi$ pointwise almost everywhere. We apply
- Lebesgue's dominated convergence theorem
- using that $\abs{\eta_r(u) \diver \xi} \leq
- \abs{u^s \diver \xi}$ and that $\xi \in C^\infty_c(\Omega,
- \mathbb{R}^2)$ is bounded by the extreme value theorem. Thus
- $\abs{u^s \diver \xi}$ is integrable and
- \begin{equation}
- \int_\Omega -\eta_r(u)
- \diver \xi \, dx \to \int_\Omega -u^s \diver \xi \, dx
- \end{equation}
- From \eqref{eq:m_ineq_sr} we then obtain
- \begin{equation}
- m\prime(s) \geq - \int_\Omega u^s \diver \xi \, dx.
- \end{equation}
- As this holds for any
- $\norm{\xi}_A^* \leq 1$, we get from the extended total variation
- definition in \eqref{eq:extended_tv} that $m'(s) \geq \TVA(u_s)$
- almost everywhere and conclude using \eqref{eq:tva_geq_mder} that
- \begin{equation}
- \TVA(u) \geq \int_{-\infty}^\infty m'(t) \, dt \geq
- \int_{-\infty}^\infty \TVA(u^s) \, ds.
- \end{equation}
- Combining the upper and lower bounds just proved, we have equality.
-\end{proof}
-This coarea formula is our first step in transforming the anisotropic
-total variation into an easily discretizable expression.
-
-The anisotropic total variation of the thresholded images occurring in
-the anisotropic coarea formula is very much related to the size of the
-boundary of the level set, as the only variation in a characteristic
-function occurs at the boundary of the set. This is why we introduce
-the following definition of the anisotropic set perimeter.
-\begin{definition}[The anisotropic set perimeter]
- Given an anisotropy tensor $A$ the anisotropic perimeter of a set
- $U$ in $\Omega$ is defined as
- \begin{equation}
- \PerA(U;\Omega) = \TVA(\idfun_U).
- \end{equation}
-\end{definition}
-\nomenclature{$\PerA(u;\Omega)$}{Anisotropic perimeter of set $U$ using
-anisotropy tensor $A$.}%
-The anisotropic set perimeter is not like the regular set perimeter and
-does not measure the length of the boundary of the set, but it can for
-sufficiently nice level sets be calculated in the following way
-\begin{equation}
- \begin{aligned}
- \PerA(\{ u > s \}; \Omega) &= \TVA(u^s) \\ &= \sup_{\norm{\xi}_A^* \leq 1}
- \int_\Omega u^s \diver \xi \, dx \\
- &= \sup_{\norm{\xi}_A^* \leq 1} \int_{\{ u > s \}} \diver \xi \,
- dx \\
- &= \sup_{\norm{\xi}_A^* \leq 1} \int_{\partial \{ u > s\} } \nu_s
- \cdot \xi \, dt \\
- &= \sup_{\norm{\eta} \leq 1} \int_{\partial \{ u > s\} } \nu_s
- \cdot \Ahalf \eta \, dt \\
- &= \int_{\partial \{ u > s\} } \sqrt{\nu_s A \nu_s} \, dt.
- \end{aligned}
- \label{eq:perimeter_calc}
-\end{equation}
-Here, $\nu_s$ is the unit exterior normal of the level set $\{ u > s \}$.
-Note that because of the compact support of $\xi$ in
-Definition~\ref{def:extended_tv}, the parts of the boundary of $U$ that
-overlap with the boundary of $\Omega$ will not be included in the
-perimeter.
-
-Exterior normals and perimeters of
-level sets of any function $u \in \BV(\Omega)$ will not be considered
-here, but can for the isotropic case be found in for example
-\cite[Section 5.4 and 5.5]{ziemer1989}.
-
-The anisotropic coarea formula allows us to transform the anisotropic
-total variation such that we are left with minimizing the following
-functional
-\begin{equation}
- F(u) = \int_\Omega (u - f)^2 \, dx + \beta \, \int_{-\infty}^\infty
- \PerA(\{u > \lambda\}; \Omega) \, d\lambda.
- \label{eq:functional_per}
-\end{equation}
-The transformation is motivated by our upcoming Cauchy--Crofton
-integration formula, and discretization, where an approximation of the
-perimeter will be computed using a graph cut machinery.
-
-\section{Cauchy--Crofton formulas}
-
-\fixme{Rating: 7/10, the stuff after the proof could use some
-work.}
-
-\begin{figure}
- \input{fig/line_param}
-\end{figure}
-
-In the fields of integral geometry and geometric measure theory there are
-a number of interesting integral formulas. Several of them fall in a
-category often referred to as \emph{Cauchy--Crofton style formulas}, and
-give ways to measure geometric objects using the set of all lines in the
-plane. The formulas presented here will give a way to measure the length
-of a curve by counting the times it intersects line in the set of all
-lines. The first formula will be for the isotropic case, and we will
-use this to prove the anisotropic formula following it.
-
-We write $\mathcal{L}$ for the set of all lines in the plane,
-and parametrize them as shown in Figure~\ref{fig:line_param}. A
-line is parametrized by the angle $\phi \in [0, 2\pi)$ of the normal going to the
-origin, and the distance $\rho \in [0, \infty)$ from origin to the line. Sometimes it is
-more convenient to consider a unit vector $\nu$ giving the direction of
-the line instead of the angle parameter $\phi$. We denote a line by
-$\ell_{\phi, \rho} = \ell_{\nu, \rho}$ where $\nu$ is a unit vector
-along the line, i.e.\ $\nu = (-\sin \phi, \cos \phi)^T$. By defining the
-measure on this set $d\mathcal{L} = \dpdr$ we are ready to introduce the
-Cauchy--Crofton formula. Note that the measure $d\mathcal{L}$ is
-invariant under rotations.
-\nomenclature{$\mathcal{L}$}{The set of all straight lines in the
-plane.}%
-\nomenclature{$\ell_{\phi, \rho}$}{A line given by the angle of the
-normal $\phi$ and the distance to origin $\rho$.}%
-\nomenclature{$\ell_{\nu, \rho}$}{A line given by a tangent vector $\nu$
-and the distance to origin $\rho$.}%
-\begin{theorem}[The Euclidean Cauchy--Crofton formula]
- Given a differentiable curve $C$ in $\mathbb{R}^2$, the length of
- this curve $\abs{C}$ is related to the set of lines $\mathcal{L}$ as
- follows
- \begin{equation}
- \int_\mathcal{L} \#(\ell_{\phi, \rho} \cap C) \,
- d\mathcal{L}(\ell_{\phi, \rho}) = 2 \abs{C},
- \end{equation}
- where $\#(\ell_{\phi, \rho} \cap C)$ is the number of times the line
- $\ell_{\phi, \rho}$ intersects the curve $C$.
- \label{thm:euclidean_cauchy_crofton}
-\end{theorem}
-\begin{proof}
- See \cite[Theorem 3, Section 1-7.]{do1976differential}.
-\end{proof}
-If our space is equipped with a metric tensor $M(x)$ such that the inner
-product of two vectors $a$ and $b$ in a point $x$ is calculated as
-$\langle a, b\rangle_M = \langle a, M(x) b \rangle$, then the length of
-a curve $\gamma$ parametrized by some parameter $t$ becomes
-\nomenclature{$M(x)$}{A metric tensor.}%
-\nomenclature{$\abs{C}$}{The length of the curve $C$.}%
-\nomenclature{$\abs{C}_M$}{The length of the curve $C$ calculated using
-the metric tensor $M$.}%
-\begin{equation}
- \abs{\gamma}_M = \int_\gamma \sqrt{\langle \dot{\gamma},
- M\big(\gamma(t)\big) \, \dot{\gamma} \rangle} \, dt.
- \label{eq:riemannian_length}
-\end{equation}
-We will now present and prove a Cauchy--Crofton formula in this case
-where our domain is equipped with a metric tensor in each point.
-This elegant formula is very useful when we later will discretize our
-perimeter calculation. The set of lines $\mathcal{L}$ is then discretized
-in a reasonable way, and the length of the curve $C$ can be approximated
-by a sum over all these lines.
-\begin{theorem}[The Riemannian Cauchy--Crofton formula]
- Assume that our space $\Omega$ is equipped with a continuous
- positive definite metric
- tensor $M(x)$, whose eigenvalues are bounded by $0 < k \leq
- \lambda_2 \leq \lambda_1 \leq K < \infty$ for all $x \in \Omega$.
- The Cauchy--Crofton formula for a differentiable curve $C$ of finite
- length then becomes
- \begin{equation}
- \abs{C}_M = \int_\mathcal{L} \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}}}
- \, d\mathcal{L}(\ell_{\nu, \rho}).
- \label{eq:riemannian_cauchy_crofton}
- \end{equation}
- \label{thm:riemannian_cauchy_crofton}
-\end{theorem}
-\begin{proof}[Proof of the Riemannian Cauchy--Crofton formula]
- Assume first that our space is equipped with a constant metric
- tensor $M$. The length of a curve in this space can be
- calculated by transforming the curve and applying the Euclidean
- Cauchy--Crofton formula
- \begin{align}
- \abs{C}_M &= \int_C \sqrt{\langle \dot{C}, M \, \dot{C} \rangle} \, dt
- = \int_C \sqrt{\langle \Mhalf \dot{C}, \Mhalf
- \dot{C} \rangle} = \abs{M^{\sfrac{1}{2}}C} \\
- &= \int_\mathcal{L} \#(\ell_{\phi, \rho} \cap M^{\sfrac{1}{2}}C) \,
- d\mathcal{L}(\ell_{\phi, \rho}) \\
- &= \int_\mathcal{L} \#(M^{-\sfrac{1}{2}}\ell_{\phi, \rho} \cap C) \,
- d\mathcal{L}(\ell_{\phi, \rho}) \\
- &= \int_\mathcal{L} \#(m_{\phi, \rho} \cap C) \abs{J_M(\ell_{\phi, \rho})}\,
- d\mathcal{L}(m_{\phi, \rho}).
- \end{align}
- Here $J_M(\ell_{\phi, \rho})$ is the Jacobian of the coordinate
- transformation $F : \mathcal{L} \to \mathcal{L}$, which maps
- $\ell_{\phi, \rho} \mapsto \Mhalf \ell_{\phi, \rho}$.
- \nomenclature{$J_M(\ell_{\phi, \rho})$}{Jacobian of the coordinate
- transformation induced by the metric tensor $M$.}
-
- We will now compute the Jacobian $J_M(\ell_{\phi, \rho})$. As $M\in
- \mathbb{R}^{2\times2}$ is symmetric, so is $M^{\sfrac{1}{2}}$, and
- it admits a decomposition $\Mhalf = U\Sigma U^T$ where the
- components correspond to the following coordinate transformations
- \begin{align}
- %U(\phi,\rho) &= (\phi + \xi, \rho), \\
- %U^T(\phi,\rho) &= (\phi - \xi, \rho), \\
- U(\ell_{\nu,\rho}) &= \ell_{\phi + \xi, \rho} = \ell_{U\nu, \rho} \\
- U^T(\ell_{\nu,\rho}) &= \ell_{\phi - \xi, \rho} = \ell_{U^T\nu, \rho} \\
- \Sigma &= \begin{pmatrix}
- \sigma_1 & 0 \\
- 0 & \sigma_2
- \end{pmatrix}
- = \begin{pmatrix}
- \sqrt{\lambda_1} & 0 \\
- 0 & \sqrt{\lambda_2}
- \end{pmatrix}
- \end{align}
- As $U$ and $U^T$ correspond to rotations and our measure
- $\mathcal{L}$ is invariant under rigid motions, $U$ and $U^T$ do not
- have direct contributions to the Jacobian. They do however affect the
- input angle of the operator $\Lambda$ such that $J_M(\ell_{\phi,
- \rho}) = J_{\Sigma^2}(U^T \ell_{\phi, \rho})$. Thus we will now
- compute $J_{\Sigma^2}(\ell_{\phi, \rho})$. Given a line
- \begin{equation}
- \ell_{\phi, \rho} = \begin{pmatrix}
- \rho\cdot \cos \phi \\
- \rho\cdot \sin \phi
- \end{pmatrix}
- + \mathbb{R} \begin{pmatrix}
- -\sin \phi \\
- \cos \phi
- \end{pmatrix},
- \end{equation}
- the operator $\Sigma$ transforms it into
- \begin{equation}
- \Sigma \ell_{\phi, \rho} = \begin{pmatrix}
- \sigma_1 \rho\cdot \cos \phi \\
- \sigma_2 \rho\cdot \sin \phi
- \end{pmatrix}
- + \mathbb{R} \begin{pmatrix}
- -\sigma_1 \sin \phi \\
- \sigma_2 \cos \phi
- \end{pmatrix},
- \end{equation}
- which equals the line $\ell_{\theta, \eta}$ with
- \begin{align}
- \theta &= \arctan \left( \frac{\sigma_1}{\sigma_2} \tan \phi
- \right) \label{eq:theta}\\
- \eta &= \left\langle \begin{pmatrix}
- \sigma_1 \rho\cdot \cos \phi \\
- \sigma_2 \rho\cdot \sin \phi
- \end{pmatrix}, \begin{pmatrix}
- \cos \theta \\
- \sin \theta
- \end{pmatrix}
- \right\rangle = \sigma_1 \rho \cdot \cos \phi \cdot \cos \theta
- + \sigma_2 \rho \cdot \sin \phi \cdot \sin \theta.
- \end{align}
- As $\partial_\rho \theta = 0$, the Jacobian becomes
- $\abs{J_{\Sigma^2}(\ell_{\phi, \rho})} = \partial_\phi\theta \cdot
- \partial_\rho \eta$. Differentiation yields
- \begin{align}
- \partial_\phi \theta &= \frac{\frac{\sigma_1}{\sigma_2}
- \sec^2 \phi}{1 + \frac{\sigma_1^2}{\sigma_2^2}
- \tan^2 \phi}
- = \frac{\sigma_1 \sigma_2}{\sigma_1^2 \sin^2 \phi + \sigma_2^2
- \cos^2 \phi}, \\
- \partial_\rho \eta &= \sigma_1 \cos \phi \cdot \cos \theta +
- \sigma_2 \sin \phi \cdot \sin \theta.
- \end{align}
- In the expression for $\partial_\rho \eta$ we insert $\theta$ from
- \eqref{eq:theta} and use that $\sin
- (\arctan( x )) = x / \sqrt{1 + x^2}$ and that $\cos ( \arctan ( x ))
- = 1 / \sqrt{1 + x^2}$ to obtain
- \begin{equation}
- \partial_\rho \eta = \frac{\sigma_1 \cos \phi + \sigma_2 \sin
- \phi \, \frac{\sigma_1}{\sigma_2} \tan \phi}
- {\sqrt{1 + \frac{\sigma_1^2}{\sigma_2^2}
- \tan^2 \phi}}
- = \frac{\sigma_1 \sigma_2}{\sqrt{\sigma_1^2 \sin^2 \phi +
- \sigma_2^2 \cos^2 \phi}}.
- \end{equation}
- If $\nu = (\nu_x, \nu_y)^T$ is a unit vector along the line
- $\ell_{\phi, \rho} = \ell_{\nu, \rho}$ then
- \begin{equation}
- \abs{J_{\Sigma^2}(\ell_{\nu, \rho})}
- = \frac{\sigma_1^2 \sigma_2^2}
- {\left(\sigma_1^2 \sin^2 \phi + \sigma_2^2 \cos^2 \phi
- \right)^{\sfrac{3}{2}}}
- = \frac{\sigma_1^2 \sigma_2^2
- }{\left(\sigma_1^2 \nu_x^2 + \sigma_2^2
- \nu_y^2\right)^{\sfrac{3}{2}}}
- = \frac{\det \Sigma^2}{(\nu^T \cdot \Sigma^2 \cdot
- \nu)^{\sfrac{3}{2}}}.
- \end{equation}
- We are interested in the Jacobian of the whole transformation
- $J_{\Sigma^2}(U^T \ell_{\theta, \rho})$, so
- all that is left to do is insert $U^T \ell_{\nu, \rho}$ to obtain
- \begin{equation}
- \abs{J_M(\ell_{\nu, \rho})} = \abs{J_{\Sigma^2}(U^T \ell_{\nu, \rho})}
- = \frac{\det M}{\left(\nu^T U \cdot \Sigma^2 \cdot U^T
- \nu\right)^{\sfrac{3}{2}}}
- = \frac{\det M}{\left(\nu^T \cdot M \cdot \nu\right)^{\sfrac{3}{2}}}
- \end{equation}
- We have now proved that for a constant metric tensor $M$, the length
- of the differentiable curve $C$ with regards to this tensor can be
- calculated as
- \begin{equation}
- \abs{C}_M = \int_C \sqrt{\langle \dot{C}, M \dot{C} \rangle} \,
- dt
- = \int_\mathcal{L} \# (\ell_{\nu, \rho} \cap C)
- \frac{\det M}{\left(\nu^T \cdot M \cdot \nu
- \right)^{\sfrac{3}{2}}} \, d\mathcal{L} (\ell_{\nu, \rho}).
- \label{eq:riemannian_const_m}
- \end{equation}
- We now argue that the similar formula in
- \eqref{eq:riemannian_cauchy_crofton} holds for a non-constant but
- continuous metric tensor $M(x)$. By partitioning the domain into
- disjoint sets $U_i$ such that $\Omega = \cup_i U_i$, we make a
- piecewise constant approximation $M_\pi(x)$ such that if $x \in U_i$
- then $M_\pi(x) = M(x_i)$ for some fixed $x_i \in U_i$. We then
- approximate \eqref{eq:riemannian_const_m} by
- \begin{equation}
- \abs{C}_{M_\pi} = \sum_i \int_\mathcal{L} \#(\ell_{\nu, \rho} \cap
- C \cap U_i) \, w_i(\nu) \, d\mathcal{L}(\ell_{\nu, \rho})
- \label{eq:mpi_approx}
- \end{equation}
- where $w_i$ is the weight-function used in the set $U_i$, that is,
- \begin{equation}
- w_i(\nu) = \frac{\det M(x_i)}{\left(\nu^T \cdot M(x_i) \cdot \nu
- \right)^{\sfrac{3}{2}}}.
- \end{equation}
- We further simplify the approximation by introducing the global
- weight-function $w_\pi(\nu, x)$ which is equal to
- $w_i(\nu)$ when $x \in U_i$. It can be written as
- \begin{equation}
- w_\pi(\nu, x) = \frac{\det M_\pi(x)}{\left(\nu^T \cdot M_\pi(x)
- \cdot \nu \right)^{\sfrac{3}{2}}}.
- \label{eq:wpi_def}
- \end{equation}
- Using this weight in \eqref{eq:mpi_approx} we can get rid of the sum
- over the partition $i$ and form a sum of all intersection point of
- $C$ and the line $\ell_{\nu, \rho}$ currently being integrated over.
- The approximation becomes
- \begin{equation}
- \begin{aligned}
- \abs{C}_{M_\pi} &=
- \sum_i \int_\mathcal{L} \sum_{x \in \ell_{\nu, \rho} \cap C \cap
- U_i} w_\pi(\nu, x) \, d\mathcal{L}(\ell_{\nu, \rho}) \\
- &= \int_\mathcal{L} \sum_{x \in \ell_{\nu, \rho} \cap C}
- w_\pi(\nu, x) \, d\mathcal{L}(\ell_{\nu, \rho}).
- \end{aligned}
- \label{eq:riemannian_segments}
- \end{equation}
- Now it only remains
- to show that the left- and right-hand side of
- \eqref{eq:riemannian_segments} converges to the left- and right-hand
- side of \eqref{eq:riemannian_cauchy_crofton}.
-
- As our partition $\pi$ is refined, the weight $w_\pi(x)$ converges
- pointwise to the continuously varying weight
- \begin{equation}
- w(\nu, x) = \frac{\det M(x)}{\left( \nu^T \cdot M(x) \cdot \nu
- \right)^{\sfrac{3}{2}}}
- \end{equation}
- found in \eqref{eq:riemannian_cauchy_crofton}.
-
- Recall from \eqref{eq:riemannian_length} that the left-hand side is
- calculated as
- \begin{equation}
- \abs{C}_{M_\pi}
- = \int_C \abs{\dot{C}(t)}_{M_\pi} \, dt
- = \int_C \sqrt{\dot{C}(t)^T M_\pi\big( C(t) \big) \dot{C}(t)}
- \, dt.
- \label{eq:piecewise_curve_length}
- \end{equation}
- We know that $M_\pi(x)$ converges pointwise to $M(x)$, and thus
- $\labs{\dot{C}(t)}_{M_\pi}$ converges pointwise to
- $\labs{\dot{C}(t)}_M$. We have assumed bounds on the eigenvalues
- of $M(x)$ such that, according to the Rayleigh principle
- \begin{equation}
- K \geq \lambda_1 = \max_\xi \frac{\xi^T M_\pi(x) \xi}{\xi^T \xi
- }
- \end{equation}
- and therefore we have the bound
- \begin{equation}
- \xi^T M_\pi(x) \xi \leq K \, \norm{\xi}^2, \quad \forall \xi.
- \end{equation}
- Thus the integrand of \eqref{eq:piecewise_curve_length} is bounded
- by $g(t) = \sqrt{K \cdot \dot{C}(t)^T \dot{C}(t)}$. We know that
- $g(t)$ is integrable as its integral is exactly $\sqrt{K} \abs{C}$
- and we have assumed that the curve is of finite length. This means
- we can apply Lebesgue's dominated convergence theorem to see that
- $\labs{C}_{M_\pi} \to \labs{C}_M$.
-
- We apply the same theorem to show that the right-hand side of
- \eqref{eq:riemannian_segments} converges. Recall the definition of
- $w_\pi$ in \eqref{eq:wpi_def}. The numerator is equal to $\sigma_1^2
- \sigma_2^2 = \lambda_1 \lambda_2$ and is by assumption bounded from
- above by $K^2$.
-
- Next we need to bound $\nu^T M_\pi(x) \nu$ away from zero. According
- to the Rayleigh principle
- \begin{equation}
- \lambda_2 = \min_{\norm{\xi} = 1} \xi^T M_\pi(x) \xi
- \end{equation}
- and thus $\nu^T M_\pi(x) \nu \geq \lambda_2 \geq k$. The weight
- function $w_\pi$ is then bounded such that
- \begin{equation}
- \sum_{x \in \ell_{\nu, \rho} \cap C} w_\pi(\nu, x)
- \leq
- \sum_{x \in \ell_{\nu, \rho} \cap C}
- \frac{K^2}{k^{\sfrac{3}{2}}}
- =
- \frac{K^2}{k^{\sfrac{3}{2}}} \cdot \# (\ell_{\nu, \rho} \cap C)
- =: g(\ell_{\nu, \rho}).
- \end{equation}
- This is integrable following from the Euclidean Cauchy--Crofton
- formula in Theorem~\ref{thm:euclidean_cauchy_crofton} and the fact
- that we assumed $C$ to be of finite length:
- \begin{equation}
- \int_\mathcal{L} g(\ell_{\nu, \rho}) \, d\mathcal{L}(\ell_{\nu,
- \rho}) = \frac{K^2}{k^{\sfrac{3}{2}}} \abs{C} < \infty.
- \end{equation}
- Thus we can apply the dominated convergence theorem again and
- conclude that
- \begin{equation}
- \int_\mathcal{L} \sum_{x \in \ell_{\nu, \rho} \cap C}
- w_\pi(\nu, x) \, d\mathcal{L}(\ell_{\nu, \rho})
- \to
- \int_\mathcal{L} \sum_{x \in \ell_{\nu, \rho} \cap C}
- w(\nu, x) \, d\mathcal{L}(\ell_{\nu, \rho})
- \end{equation}
- which---as both sides of the equality in
- \eqref{eq:riemannian_segments} have been shown to converge---leaves
- us with what we wanted to prove
- \begin{equation}
- \abs{C}_M =
- \int_\mathcal{L} \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}}} \,
- d\mathcal{L}(\ell_{\nu, \rho}).
- \end{equation}
-\end{proof}
-
-With the anisotropic coarea formula in
-Theorem~\ref{thm:anisotropic_coarea} we have a way to calculate the
-anisotropic total variation by integrating the anisotropic perimeter of
-each level set of the image as introduced in \eqref{eq:functional_per}.
-In the Euclidean case, which here would amount to setting the anisotropy
-tensor $A$ equal to the identity matrix $I$, the perimeter coincides
-nicely with the length of the boundary curve, assuming some regularity
-for the boundary.
-
-In the general case we need
-to be more careful. As can be seen in \eqref{eq:perimeter_calc}, the
-anisotropic perimeter is calculated by integrating the norm of the
-\emph{normal vector} around the boundary, while the anisotropic curve
-length in \eqref{eq:riemannian_length} is the integral of the norm of
-the \emph{tangent vector} of the curve. Thus a 90\textdegree{} rotation
-is needed.
-
-If $P$ is a 90\textdegree{} rotation matrix we have
-\begin{equation}
- \begin{aligned}
- \PerA(U; \Omega) &= \int_{\partial U} \sqrt{\langle
- \nu_{\partial U}, A(x) \nu_{\partial U} \rangle} \, dt \\
- &= \int_{\partial U} \sqrt{\langle P \nu_{\partial U}, P A(x)
- P^T \, P \nu_{\partial U} \rangle } \, dt.
- %&\stackrel{\mathclap{\gamma \hskip 0.07em = \hskip 0.07em
- %\partial U}}{=}
- %\int_{\gamma} \sqrt{ \langle \dot{\gamma}, M(x) \,
- %\dot{\gamma} \rangle } \, dt.
- \end{aligned}
- \label{eq:per_to_length1}
-\end{equation}
-We simplify the equation by defining the metric tensor $M(x) = P A(x)
-P^T$ and letting $\gamma = \partial U \cap \Omega$ be an arclength parametrization of
-the boundary of $U$ that does not overlap with the boundary of $\Omega$
-\begin{equation}
- \PerA(U; \Omega) = \int_\gamma \sqrt{\langle \dot{\gamma}, M(x) \,
- \dot{\gamma} \rangle } \, dt.
- \label{eq:per_to_length2}
-\end{equation}
-
-Now we make sure that all the assumptions of the Riemannian
-Cauchy--Crofton formula in Theorem~\ref{thm:riemannian_cauchy_crofton}
-are fulfilled so that it can be applied to the curve length integral we
-have constructed in \eqref{eq:per_to_length2}.
-
-The structure tensor is constructed as described in
-Section~\ref{sec:anisotropy_tensor}
-\begin{equation}
- S_\rho(x) = \left(K_{\rho} * \left( \nabla f_{\sigma}
- \otimes \nabla f_{\sigma} \right)\right)(x),
-\end{equation}
-Because of the convolutions with the Gaussian function, this is a smooth
-map from $\bar{\Omega}$ to $\mathbb{R}^{2 \times 2}$. As we
-can see in \eqref{eq:s_eigenvalues}, the eigenvalues
-depend continuously on the coefficients of the
-elements in
-the structure tensor $S_\rho(x)$. The extreme value theorem states that
-a continuous real-valued function on a nonempty compact
-space is bounded above. Thus the eigenvalues $s_1 \geq s_2$ of $S_\rho(x)$
-are bounded above and by the construction in
-\eqref{eq:sigma_construction}, there exists uniform bound $k$ such that
-the smallest eigenvalue of the anisotropy tensor $A(x)$ is bounded away
-from zero, as
-\begin{equation}
- \lambda_2 = \left(1
- + \frac{(s_1 - s_2)^2}{\omega^2}\right)^{-1} \geq \left(1 +
- \frac{s_1^2}{\omega^2}\right)^{-1} \geq k > 0.
-\end{equation}
-Hence our metric tensor $M(x) = P A(x) P^T$ is continuous and positive
-definite and thus the curve length calculation in
-\eqref{eq:per_to_length2} fulfills all the assumptions of the
-anisotropic Cauchy--Crofton formula in
-Theorem~\ref{thm:riemannian_cauchy_crofton}.
-Thus we can apply the formula to calculate the
-perimeter in \eqref{eq:per_to_length2} as
-\begin{equation}
- \PerA(U; \Omega) = \int_\mathcal{L} \sum_{x \in
- \ell_{\nu, \rho} \cap \gamma } \frac{\det M(x)}
- {2 \left( \nu^T \cdot M(x) \cdot \nu \right)^{\sfrac{3}{2}}}
- \, d\mathcal{L}(\ell_{\nu, \rho}) \, ds,
-\end{equation}
-where $\gamma = \partial U \cap \Omega$.
-Note that $P$ does not affect the determinant, i.e.\
-$\det A = \det PAP^T = \det M$, and from our decomposition in
-\eqref{eq:sigma_construction} we see that the transformation $PAP^T \to
-M$ actually amounts to switching the two eigenvalues $\sigma_1$ and
-$\sigma_2$ in $\Sigma$.
-
-This concludes the treatment of the continuous problem. We have seen how
-the anisotropic coarea formula in Theorem~\ref{thm:anisotropic_coarea}
-allows us to calculate the anisotropic total variation as an integral of
-the perimeter of all the level sets. Through the anisotropic
-Cauchy--Crofton formula in Theorem~\ref{thm:riemannian_cauchy_crofton}
-these perimeters are calculated by an integral over the set of all
-lines. We are then left with the functional
-\begin{equation}
- F(u) = \int_\Omega (u - f)^2 + \beta \, \TVA(u),
- \label{eq:final_cont_rep}
-\end{equation}
-where
-\begin{equation}
- \TVA(u) = \int_{-\infty}^\infty \int_{\mathcal{L}} \sum_{x \in
- \ell_{\nu, \rho} \cap \gamma_s} \frac{\det M(x)}
- {2 \left( \nu^T \cdot M(x) \cdot \nu \right)^{\sfrac{3}{2}}}
- \, d\mathcal{L}(\ell_{\nu, \rho}) \, ds,
-\end{equation}
-and $\gamma_s = \partial \{u > s\} \cap \Omega$.
-Within the restrictions that these theorems put on the tensor $M(x)$, we
-have chosen a construction where one eigenvalue is always 1, while the
-other varies from 1 in smooth areas towards 0 around edges, with the
-corresponding eigenvector perpendicular to the edge. \fixme{meh}
-
-As mentioned previously, these transformations are motivated by the
-discretization that follows in the next chapter.
-
-\chapter{Discrete formulation}
-
-The whole transformation from the initial functional in
-\eqref{eq:first_anisotropic_functional} through the anisotropic coarea
-formula and the Cauchy--Crofton formula was motivated by the
-discrete formulation which will be described here. After discretizing
-the functional, we will se how a graph cut approach can be used to find
-a global minimizer in polynomial time.
-
-\section{Discretization}
-
-\fixme{Rating: 4/10}
-
-Assume that our discrete images are given on a uniform grid $\mathcal{G}$,
-where each grid point is called a pixel. The image is a function giving
-each pixel a value in the set of levels $\mathcal{P} = \{0, \hdots,
-L-1\}$. This is a reasonable assumption for digital grayscale
-images.
-
-Thus, when discretizing the functional in \eqref{eq:final_cont_rep}, we
-have to consider that our images now have both discrete domain and
-co-domain.
-
-The integrals in \eqref{eq:final_cont_rep} will be approximated by
-discrete sums. First the fidelity term is discretized without too much
-trouble, while with the regularization term, there is more choice as to
-how to discretize the set of lines $\mathcal{L}$. In the end we will
-verify that our discretization is consistent with the continuous
-functional.
-
-\subsection{Fidelity term}
-
-\fixme{Rating: 7/10}
-
-Since it is not affected by our introduction of the anisotropy tensor,
-the fidelity term can be discretized as in my project work
-\cite{project}. For some pixel position $x \in \mathcal{G}$ and some
-level value $k \in
-\mathcal{P}$, we define the following function
-\begin{equation}
- N_x(k) = \abs{k - f_x}^2
-\end{equation}
-which is the value of the fidelity term if we were to give $u_x$ a value
-of $k$. This allows us write
-\begin{equation}
- \int_\Omega \abs{u - f}^2\, dx \approx \sum_{x \in \mathcal{G}}
- \abs{u_x - f_x}^2 \Delta x = \sum_{x \in \mathcal{G}} N_x(u_x)
- \Delta x.
- \label{eq:fidelity_approx_1}
-\end{equation}
-The reason we introduce the function $N_x(k)$ is that we want to apply
-the following decomposition formula, which holds for any function $F(k)$
-taking values $k \in \mathcal{P}$:
-\begin{equation}
- \begin{aligned}
- F(k) &= \sum_{\lambda = 0}^{k-1} \big( F(\lambda + 1) -
- F(\lambda) \big) + F(0) \\
- &= \sum_{\lambda = 0}^{L-2} \big( F(\lambda + 1) - F(\lambda) \big)
- I( \lambda < k) + F(0),
- \end{aligned}
-\end{equation}
-where $I(x)$ is the indicator function that takes the value 1 if $x$ is
-true, and 0 if $x$ is false. Since $I(\lambda < u_x) = u_x^\lambda$ we
-rewrite \eqref{eq:fidelity_approx_1} and obtain
-\begin{equation}
- \sum_{x \in \mathcal{G}} \abs{u_x - f_x}^2 =
- \sum_{x \in \mathcal{G}} N_x(u_x) = \sum_{\lambda = 0}^{L-2} \sum_{x
- \in \mathcal{G}}
- \big( N_x(\lambda + 1) - N_x(\lambda) \big) \, u_x^\lambda + N_x(0).
-\end{equation}
-As our domain is discretized uniformly, we drop the constant
-$\Delta x$, and absorb it into our parameter $\beta$ of
-\eqref{eq:final_cont_rep}. Note that since
-our image takes values in $\mathcal{P} = \{0, \hdots, L-1\}$, the
-thresholded image $u^{L-1}$ is equal to zero everywhere.
-
-\subsection{Regularization term}
-\label{sec:disc_regularization}
-
-\fixme{Rating: 6/10}
-
-Discretizing the regularization term is more challenging. We introduce
-the discrete levels to get
-\begin{equation}
- \int_{-\infty}^\infty \PerA( \{ u > \lambda \}; \Omega) \, d\lambda
- \approx \sum_{\lambda = 0}^{L-2} \PerA( \{ u > \lambda \}; \Omega)
- \, \Delta \lambda.
- \label{eq:per_approx1}
-\end{equation}
-As with the $\Delta x$ difference, we can absorb the $\Delta \lambda$
-difference into the $\beta$ parameter of
-\eqref{eq:final_cont_rep}. 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
-has zero perimeter.
-
-\subsubsection{Discrete Riemannian Cauchy--Crofton formula}
-
-\begin{figure}
- \input{fig/line_disc}
-\end{figure}
-By approximating the integral
-Theorem~\ref{thm:riemannian_cauchy_crofton} by a discrete sum we obtain
-the approximation
-\begin{equation}
- \begin{aligned}
- \abs{C}_M &= \int_\mathcal{L} \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}}} \,
- d\mathcal{L}(\ell_{\nu, \rho}) \\
- &\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\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 \,
- \Delta\nu.
- \label{eq:cauchy_crofton_approx1}
- \end{aligned}
-\end{equation}
-The set of lines $\mathcal{L}$ has been discretized into the set
-$\mathcal{L}_D$. Note that we are approximating the length of the
-\emph{differentiable} curve $C$ in $\Omega$. Being a difference in the
-$\rho$ parameter of our line discretization in
-Figure~\ref{fig:line_param}, the difference $\Delta \rho$ represents the
-distance from one line to the next in a line family as shown in
-Figure~\ref{fig:line_family}, and thus depends on the angle $\phi$
-considered. The difference $\Delta \phi$ is taken to
-be the average of the distance to the two neighboring line families as
-shown in Figure~\ref{fig:line_neigh}, and thus also depends on $\phi$.
-
-In the discrete setting our domain $\Omega$ is discretized as a uniform
-grid $\mathcal{G}$. Our image is then reduced to a function $u :
-\mathcal{G} \to \mathcal{P}$.
-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}.
-
-The choice of our discrete set of lines $\mathcal{L}_D$ is important, as
-it will decide the accuracy of our approximation.
-We need some sensible restrictions on
-the set $\mathcal{L}_D$ to simplify the further discussion. All lines
-intersect at least two grid points, and from the periodicity of our grid
-they thus intersect an infinite number of grid points. This puts some
-restrictions on the angles we can choose. For each angle, we include all
-possible lines of that family, meaning there are no grid points without
-a line of that family intersecting it.
-
-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~\ref{fig:line_family} shows all lines of a given family, i.e.\
-lines having the same angle parameter $\phi$.
-
-Thus not only have we discretized the set of lines, but each line is
-made up of edges going from one grid point to the next. We will
-denote such an edge by $e$ or $e_{ab}$ when its endpoints are $a, b \in
-\mathcal{G}$. Thus we rewrite the discretization of
-\eqref{eq:cauchy_crofton_approx1}, and sum over all the edges in the
-discretization $\mathcal{L}_D$ to obtain
-\begin{equation}
- \abs{C}_M \approx \sum_{e} \sum_{x \in e \cap C} \, \frac{\det M(x)
- \norm{e}^3}{2 \left(e^T \cdot M(x) \cdot e\right)^{\sfrac{3}{2}}} \,
- \Delta\phi \, \Delta\rho.
- \label{eq:cauchy_crofton_approx2}
-\end{equation}
-This is beginning to look like something we can calculate. One
-difficulty is finding the intersections $e \cap C$. The exact
-calculations of these points will not fit into our graph cut framework
-later, and thus for an edge $e$ we will only consider the question of
-``did $e$ cross $C$ or not?'' This amounts to checking whether the
-terminals of $e$ lie on different sides 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}, be wrong when we have more.
-
-\begin{figure}
- \input{fig/curve_edge}
-\end{figure}
-
-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 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}
- M(x) \approx M(e_{ab}) = \frac{M(a) + M(b)}{2},
- \label{eq:tensor_approx}
-\end{equation}
-the component-wise average of the tensors in the two end points of the
-edge. Recall that we have already done some spatial smoothing of the
-structure tensor in \eqref{eq:s_def} corresponding to the
-\emph{integration scale} $\rho$, and thus we expect the tensors $M(a)$
-and $M(b)$ to be similar for edges $e$ of reasonably short length.
-
-We also remark that using the Rayleigh principle, it is easy to conclude
-that the eigenvalues of the tensor approximation $M(e_{ab})$ are bounded
-below and above by the smallest and largest eigenvalues of $M(a)$ and
-$M(b)$.
-
-\begin{figure}
- \input{fig/area_proof}
-\end{figure}
-
-We have now almost arrived at our final curve length approximation, but
-we need a way to calculate the inter-line distance $\Delta \rho$ which
-will be provided by the following lemma.
-\begin{lemma}
- For each family of lines given by an angle parameter $\phi$ in the
- uniform 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 a line $\ell$ intersecting the point $(p, q) \in
- \mathbb{Z}^2$ in the grid.
- The distance $\Delta \rho$ from this line $\ell$ to the neighboring
- lines can then be calculated as a minimum over the distance to all
- other grid points.
-
- The lines are split into edges $e = (\delta s, \delta t)^T$ where
- $s, t \in \mathbb{Z}$ are coprime such that $e$ does not intersect
- any other grid points than its two endpoints.
-
- We then calculate the minimal distance to a grid point not on the
- line $\ell$ as
- \begin{equation}
- \begin{aligned}
- \Delta \rho &= \min_{(p\prime, q\prime) \in \mathcal{G}
- \setminus \ell} \left\{
- \left\langle \delta [p - p\prime, q - q\prime],
- \frac{e^\perp}{\norm{e^\perp}} \right\rangle \right\} \\
- &= \min_{(p\prime, q\prime) \in \mathcal{G} \setminus \ell}
- \left\{\delta^2 \cdot \frac{t(p-p\prime) - s(q -
- q\prime)}{\norm{e}} \right\}.
- \end{aligned}
- \end{equation}
- Since $s$ and $t$ are coprime, there exists $a, b \in \mathbb{Z}$
- such that $at - bs = 1$, and since the $\Delta \rho$ cannot be zero,
- we obtain
- \begin{equation}
- \Delta \rho = \frac{\delta^2}{\norm{e}}.
- \end{equation}
-\end{proof}
-Inserting $\Delta \rho = \delta^2 / \norm{e}$ and the tensor
-approximation of \eqref{eq:tensor_approx} into the curve length
-approximation of \eqref{eq:cauchy_crofton_approx2} 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}}},
- \label{eq:cauchy_crofton_approx3}
-\end{equation}
-where the sum is over all edges crossing the curve.
-
-The curve length we initially wanted to calculate was the perimeter
-$\PerA(\{u > \lambda\}; \Omega)$ in \eqref{eq:per_approx1}. To find
-edges that crosses the curve, we identify the edges that have one
-terminal inside the level set, and the other outside.
-Thus we rewrite the sum over $e \cap C$ such that
-\begin{equation}
- \PerA(\{u > \lambda\}; \Omega) \approx \sum_{e_{ab}}
- \abs{u^\lambda_a - u^\lambda_b} \frac{\det M(e_{ab})
- \norm{e_{ab}}^2 \, \delta^2 \, \Delta\phi}{2 \left(e_{ab}^T
- \cdot M(e_{ab}) \cdot e_{ab}\right)^{\sfrac{3}{2}}}.
- \label{eq:per_approx2}
-\end{equation}
-The absolute value $\abs{u^\lambda_a - u^\lambda_b}$ is 1 if one of
-$a$ and $b$ lie inside the level set and the other lies outside, and
-0 otherwise.
-
-In other words the absolute value is one if $e_{ab}$ crosses the
-perimeter of $\{ u > \lambda\}$ an odd number of times, and zero
-otherwise.
-
-Thus we have arrived at our final discretization, which takes the form
-\begin{gather}
- F(u) = \sum_\lambda \sum_x F_x^\lambda(u_x^\lambda) + \beta \sum_\lambda
- \sum_{(x, y)} F_{x,y}^\lambda(u_x^\lambda, u_y^\lambda) =:
- F^\lambda(u^\lambda), \\
- \begin{aligned}
- F_x^\lambda(u_x^\lambda) &= \big(N_x(\lambda + 1) -
- N_x(\lambda)\big) \cdot u_x^\lambda, \\
- F_{x,y}^\lambda(u_x^\lambda, u_y^\lambda)
- &= \abs{u_x^\lambda - u_y^\lambda} \frac{\det M(e_{xy})
- \norm{e_{xy}}^2 \delta^2 \Delta \phi}{2 \left( e_{xy}^T \cdot
- M(e_{xy}) \cdot e_{xy} \right)^{\sfrac{3}{2}}}.
- \end{aligned}
- \label{eq:final_discretization}
-\end{gather}
-Recall that $N_x(\lambda) = \abs{ \lambda - f_x }^2$.
-
-If we minimize $F_\lambda$ to obtain $u^\lambda$ for each level
-separately, it is obvious that we will also minimize the sum over all
-$F_\lambda$. However, it is not guaranteed that the obtained thresholded
-images $u^\lambda$ can be combined to make an output image $u$. They
-were defined as $u^\lambda = \idfun_{u > \lambda}$, so we need them to
-be monotonically decreasing \fixme{(?)} in increasing level values, i.e.\
-\begin{equation}
- u_x^\lambda \geq u_x^\mu, \quad \forall \lambda \leq \mu, \quad
- \forall x \in \mathcal{G}.
-\end{equation}
-Later we will present a graph cut algorithm that find the thresholded
-images minimizing each level, \emph{while guaranteeing that they meet
-this requirement.}
-
-\fixme{maybe with a $w_{xy}$ definition}
-
-\subsubsection{Consistency}
-
-\fixme{Rating: 6/10}
-
-Consistency relates to whether a solution to the continuous problem
-fits in the discretized equation, in other words, whether 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. The sum is a midpoint rule
-approximation of the integral. As the grid is refined and $\delta \to 0$
-the sum will converge to 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 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}.
-
-If we consider the edges $e$ of each family separately, the curve length
-approximation in \eqref{eq:cauchy_crofton_approx3} can be written
-\begin{equation}
- \begin{aligned}
- \abs{C}_M &=
- \int_\nu \int_\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}}} \, d\rho \, d\nu \\
- &\approx \sum_\nu \sum_\rho \sum_{e_{\nu, \rho} \cap C}
- \, \frac{\det M(e_{\nu, \rho}) \norm{e_{\nu, \rho}}^3}{2\left(
- e_{\nu, \rho}^T \cdot M(e_{\nu, \rho}) \cdot e_{\nu, \rho}
- \right)^{\sfrac{3}{2}}} \, \Delta \rho \, \Delta \nu.
- \end{aligned}
-\end{equation}
-
-As described in the construction of this formula, there are four main
-approximations used. Firstly there is the fact that we do not consider
-the actual intersection points, but only whether an edge crosses the
-curve or not. Secondly we have the tensor which is averaged as in
-\eqref{eq:tensor_approx}. And then we have the discretizations of our
-two line parameters $\nu$ and $\rho$.
-
-It is intuitive that if $\sup \norm{e} \to 0$, the number of times the
-differentiable curve $C$ can cross a given edge decreases. We will not
-prove convergence, but rather assume that the special cases where it
-might not work, are negligible.
-
-Further, if $\sup \norm{e} \to 0$ it is obvious that the tensor average
-in \eqref{eq:tensor_approx} converges to the tensor in the intersection
-point.
-
-For each $\phi$ parameter, our discretization in the $\rho$ dimension
-can be regarded as a midpoint rule as shown in
-Figure~\ref{fig:line_midpoint}. Thus if $\sup \Delta \rho \to 0$, this
-part of the discretization is fine.
-
-\begin{figure}
- \input{fig/line_midpoint}
-\end{figure}
-
-The discretization in the $\phi$ dimension can also be regarded as a
-version of the \emph{rectangle method}, although not the midpoint rule.
-As shown in Figure~\ref{fig:circ_rule}, the circle is split into
-intervals
-\begin{equation}
- \left[
- \frac{\phi_{k-1} + \phi_{k}}{2}, \frac{\phi_{k} + \phi_{k+1}}{2}
- \right]
-\end{equation}
-of length $\Delta \phi_k = (\phi_{k+1} + \phi_{k-1}) / 2$. The summand
-is evaluated at $\phi_k$, somewhere inside the interval.
-Thus if
-$\sup \Delta \phi_k \to 0$, this discretization is also consistent.
-
-\begin{figure}
- \input{fig/circ_rule}
-\end{figure}
-
-To show that all these properties can be fulfilled, we look at a
-particular neighborhood stencil construction.
-Consider a square centered around a grid point with side lengths
-$\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 $n^2 = \lfloor 1 / \sqrt{\delta} \rfloor^2$ grid
-points. This means that the number of grid points along the outer edge
-of this square $n$ goes to infinity.
-
-\begin{figure}
- \input{fig/square_cons}
-\end{figure}
-
-We include all grid points inside the square in our neighborhood, except
-for multiple points that lie on the same line from the origin. Thus if
-two or more grid points lie on the same line, we include only the one
-closest to the origin. This implies that 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 construction can be seen in Figure~\ref{fig:square_cons} for $n =
-5$.
-The maximal angle between two lines $\phi_k - \phi_{k-1}$ will always be
-when one of $\phi_k$ or $\phi_{k-1}$ is horizontal or vertical, shown in
-Figure~\ref{fig:square_cons} as the angle $a$.
-Thus the
-largest $\Delta \phi_k$ will then be around the horizontal or
-vertical edge.
-These angles can be calculated to be
-\begin{equation}
- \sup \Delta \phi_k = \frac{2}{2} \arctan \frac{1/n}{n/2} = \arctan
- \frac{2}{n^2} \to 0.
-\end{equation}
-
-Further we see that the edge length will be bounded by half of the
-diagonal of the square such that
-\begin{equation}
- \norm{e} \leq \sqrt{\delta / 2} \to 0.
-\end{equation}
-And finally we know from Lemma~\ref{lem:delta_rho} that for each line family
-$\delta^2 = \Delta \rho \norm{e}$ and the fact that $\norm{e} \geq
-\delta$. Thus for the inter-line distance $\Delta \rho$ we have
-\begin{equation}
- \sup \Delta \rho = \sup \frac{\delta^2}{\norm{e}} \leq
- \frac{\delta^2}{\delta} = \delta \to 0.
-\end{equation}
-
-And thus the approximation has been shown to be equivalent to
-well-known, and consistent integral approximations, where the
-summand converges to the integrand, and the differences $\Delta \phi$
-and $\Delta \rho$ go to zero. Thus the perimeter approximation in
-\eqref{eq:per_approx1} is consistent with the continuous formulation in
-Theorem~\ref{thm:riemannian_cauchy_crofton}.
-
-Note that as we will work with digital images with fixed resolutions, we
-do not really have the chance to refine our discretization. We do
-however have to take these things into account when creating our
-neighborhood stencil, to make sure that we get a reasonable
-approximation of the perimeter lengths.
-
-\section{Graph cut formulation}
-
-The discretization we arrived at in \eqref{eq:final_discretization} can
-be minimized using graph cuts. For each level $\lambda$, a minimum graph
-cut is found to produce the corresponding level set $\{ u > \lambda \}$. These
-are then combined to form the final restored image $u$.
-
-In this section we will look at how these graphs are constructed such
-that their minimum cuts correspond to the minimizers of the functional
-$F^\lambda$. The description is taken from my project work
-\cite{project}, and is included here for completeness.
-
-\subsection{Graphs}
-
-\fixme{Maybe just call them graphs this time? from $E$ to $F$!}
-
-Using the notation of \cite{cormen2009introduction} we will denote a
-directed graph as $G = (V, E)$ where $V$ is a finite set of vertices, and
-$E$ is a binary relation on $V$. If $(u, v) \in E$ we say that there is
-an edge from $u$ to $v$ in the graph $G$.
-
-We introduce the non-negative capacity function $c : V \times V \to
-\left[0, \infty \right)$. Only edges $(u, v) \in E$ can have a positive
-capacity $c(u, v) = q > 0$ and it means that it is possible to send a
-\emph{flow} of maximum $q$ units from $u$ to $v$. For convenience we
-will let $c(u, v) = 0$ for any pair $(u, v) \not\in E$, and we do not
-allow self-loops in our graph. When a directed graph $G$ is equipped
-with capacity function $c$, one might call it a capacitated graph or a
-graph, but as all our graphs will be capacitated from this point, we
-will just call them graphs and we write $G = (V, E, c)$.
-
-There are two special vertices in the graph, the source $s$ and the
-sink $t$. Contrary to other vertices, which can neither produce nor
-receive excess flow, the source can produce and the sink can receive an
-unlimited amount of flow. The most basic problem in graph flow theory
-is the question of how much flow it is possible to send through the
-graph from the source to the sink.
-
-What we seek in our final graph is a minimum $s$-$t$-cut, a ``minimal''
-line through the graph that cuts a set of edges and
-divides the vertex set in two, separating the source from the sink.
-\begin{definition}[$s$-$t$-cut]
- Given a graph $G = (V, E, c)$, an $s$-$t$-cut $(S, T)$ of $G$ is a
- partition of $V$ into $S$ and $T = V - S$ such that $s \in S$ and $t
- \in T$. The capacity of the cut is
- \begin{equation}
- c(S, T) = \sum_{u \in S} \sum_{v \in T} c(u, v),
- \end{equation}
- and a minimum $s$-$t$-cut is a cut whose capacity is minimum over
- all $s$-$t$-cuts.
- \label{def:s_t_cut}
-\end{definition}
-Note that there might exist several minimum $s$-$t$-cuts in a graph,
-resulting in different partitions of $V$. This is why we need to verify
-later that the cuts we obtain result in stackable thresholded images.
-
-\subsection{Graph representable energy functions}
-
-The next step is to find a way to construct a graph such that we can
-minimize the energy in \eqref{eq:total_energy} by finding a minimum
-$s$-$t$-cut. We will do this by creating small and simple graphs
-representing the separate summands of the energy. For these small
-graphs it will be easy to verify that the minimal cut also minimizes
-the corresponding part of the energy function, and they can then be
-merged giving a graph representing the complete energy function in
-\eqref{eq:total_energy}.
-
-First we need to establish the definition of a graph representable
-function, presented by Kolmogorov and Zabih in
-\cite{kolmogorov2002energy}.
-\begin{definition}[Graph representable functions]
- A function $\mathcal{E}(x_1, \hdots, x_n)$ of $n$ binary
- variables is graph-representable if
- there exists a graph $G = (V, E, c)$ with terminals $s$ and $t$,
- and a subset of vertices $V_0 = \{v_1, \hdots, v_n\} \subseteq V -
- \{s, t\}$ such that, for any configuration $(x_1, \hdots, x_n) \in
- \{0, 1\}^n$, the value of the energy $\mathcal{E}(x_1, \hdots, x_n)$
- is equal to a constant plus the cost of the minimum $s$-$t$-cut
- among all cuts $C = (S, T)$ where $x_i = 0 \Leftrightarrow v_i \in
- S$ and $x_i = 1 \Leftrightarrow v_i \in T$, $\forall\, 1 \leq i \leq
- n$.
- \label{def:graph_representable}
-\end{definition}
-From this definition we see that if we have a graph-representable
-function $\mathcal{E}$ it is possible to find an exact global minimum
-of $\mathcal{E}$ by finding a minimal $s$-$t$-cut in a graph
-representing $\mathcal{E}$.
-
-Furthermore Kolmogorov and Zabih present an important result concerning
-what kinds of functions are graph-representable.
-\begin{theorem}[Identification of graph representable functions]
- Given an energy function $\mathcal{E}$ of $n$ binary variables of
- the form
- \begin{equation}
- \mathcal{E}(x_1, \hdots, x_n) = \sum_i \mathcal{E}^i(x_i)
- + \sum_{i<j} \mathcal{E}^{i,j}(x_i, x_j),
- \end{equation}
- then $\mathcal{E}$ is graph representable if and only if each term
- $\mathcal{E}^{i,j}$ satisfies the inequality
- \begin{equation}
- \mathcal{E}^{i,j}(0,0) + \mathcal{E}^{i,j}(1,1) \leq
- \mathcal{E}^{i,j}(0,1) + \mathcal{E}^{i,j}(1,0).
- \end{equation}
- \label{thm:graph_rep_id}
-\end{theorem}
-This theorem will allow us to verify that our energy function actually
-is graph representable.
-
-Finally, the following theorem, proved by Kolmogorov and Zabih in
-\cite{kolmogorov2002energy}, will be crucial in our graph
-construction.
-\begin{theorem}[Additivity]
- The sum of a finite number of graph-representable functions
- \begin{equation}
- \mathcal{E}(x_1, \hdots, x_n) = \sum_k \mathcal{E}^k(x_1,
- \hdots, x_n),
- \end{equation}
- each represented by a graph $G^k = (V^k, E^k, c^k)$, is
- graph-representable by $G = (V, E, c)$ where $V = \cup_k V^k$,
- $E = \cup_k E^k$ and $c(u,v) = \sum_k c^k(u,v)$.
-\end{theorem}
-It allows us to construct small graphs representing the different
-summands of our energy function \eqref{eq:total_energy}, before adding
-them together to create a final graph representing the total energy.
-
-Note that when we apply this theorem later, we will assume that all the
-summands of \eqref{eq:total_energy} have the whole picture as their
-domain. It is unproblematic to extend $E_\lambda^x(u^\lambda_x)$ and
-$E^{x,y}(u^\lambda_x, u^\lambda_y)$ such that they take all the pixels
-their argument and then ignore all pixels except the ones they actually
-depend on.
-
-\subsection{Graph construction}
-
-\begin{figure}[t]
- \input{fig/norm_subgraph}
-\end{figure}
-
-We will construct a graph in such a way that if a vertex
-$u^\lambda_x$ ends up on the source side of the cut we set $u^\lambda_x
-= 0$, and if it ends up on the sink side we set $u^\lambda_x = 1$, as in
-Definition~\ref{def:graph_representable}. This is an arbitrary choice,
-but still something we have to keep in mind through the rest of the
-section.
-
-We will now consider the two kinds of summands in the energy function in
-\eqref{eq:total_energy}. The fidelity term coming from our aim to keep
-the output image close to the original image, and the regularization
-term coming from our aim to minimize the total variation.
-
-\begin{table}[b]
- \centering
- \caption{Each row represents one of the two possible values of
- $u^\lambda_x \in \{0,1\}$. The energy $E_\lambda^x(u^\lambda_x)$
- and minimum cut obtaining this configuration is shown. The last
- two columns show the capacities of the cut for each of the two
- graph constructions in Figure~\ref{fig:norm_subgraph}. We
- verify that for each of the two graph constructions, the cut
- capacities are equal to the energies, plus a constant.
- }
- \begin{tabular}{ccccc}
- \hline
- $u^\lambda_x$ & $E_\lambda^x(u^\lambda_x)$ & Min.\ cut $(S, T)$
- & Graph \subref{fig:norm_subgraph_pos} cut cap. & Graph
- \subref{fig:norm_subgraph_neg} cut cap. \\
- \hline
- $0$ & $0$ & $(\{s, u^\lambda_x\}, \{t\})$ & $0$ & $-E_\lambda^x(1)$ \\
- $1$ & $E_\lambda^x(1)$ & $(\{s\}, \{u^\lambda_x, t\})$ & $E_\lambda^x(1)$ & $0$
- \end{tabular}
- \label{tab:fid_energy}
-\end{table}
-
-\subsubsection{Fidelity term}
-
-\begin{figure}[t]
- \input{fig/neigh_subgraph}
-\end{figure}
-
-The fidelity term of our energy function in \eqref{eq:total_energy}
-simplifies to
-\begin{align}
- E_\lambda^x(0) &= 0
- \label{eq:fidelity_energy0} \\
- E_\lambda^x(1) &=
- N_x(\lambda + 1) -
- N_x(\lambda)
- \label{eq:fidelity_energy1}
-\end{align}
-where $E_\lambda^x(1)$ might be positive or negative depending on
-$\lambda$ and the pixel value $v_x$.
-
-Figure~\ref{fig:norm_subgraph} shows how graphs can be constructed to
-represent this part of the total energy. The construction differs
-depending on whether $E_\lambda^x(1)$ is positive or negative.
-Table~\ref{tab:fid_energy} shows how the cuts correspond to the values
-of $u^\lambda_x$ and we can easily verify that the constructed graph
-actually represents the fidelity term in the energy function.
-
-\subsubsection{Regularization term}
-
-For our neighboring relation in \eqref{eq:total_energy} of the form
-\begin{equation}
- E^{x,y}(u^\lambda_x, u^\lambda_y) =
- w_{xy} \abs{u^\lambda_x - u^\lambda_y}
-\end{equation}
-we have
-\begin{equation}
-\begin{aligned}
- E^{x,y}(0, 0) &= 0, \\
- E^{x,y}(0, 1) &= w_{xy}, \\
- E^{x,y}(1, 0) &= w_{xy}, \\
- E^{x,y}(1, 1) &= 0,
-\end{aligned}
-\label{eq:neigh_energies}
-\end{equation}
-and by Theorem~\ref{thm:graph_rep_id} our energy function is graph
-representable. In \cite{kolmogorov2002energy}, Kolmogorov and Zabih
-presents a way to construct a graph for any graph representable function
-on the form shown in Theorem~\ref{thm:graph_rep_id}. Since the
-energies in \eqref{eq:neigh_energies} are especially simple, the
-construction and presentation is simplified.
-
-
-Figure~\ref{fig:neigh_subgraph} shows two different ways of how a graph
-can be constructed to represent the regularization term.
-See Table~\ref{tab:neigh_energy} for an overview of how the two values
-of $u^\lambda_x$ corresponds to cuts in the graph.
-
-\begin{table}[t]
- \centering
- \caption{An overview of the possible configurations of the variables
- in the term $E^{x,y}(u^\lambda_x, u^\lambda_y)$. For each
- configuration the corresponding energy and the cut yielding this
- configuration is shown. The last two columns show the capacities
- of the cut in the two alternative graph constructions shown in
- Figure~\ref{fig:neigh_subgraph}. We verify that for each of the
- two graph constructions, the cut capacities are equal to the
- energies, plus a constant.
- }
- \begin{tabular}{ccccc}
- \hline
- $(u^\lambda_x, u^\lambda_y)$
- & $E^{x,y}(u^\lambda_x, u^\lambda_y)$
- & Min.\ cut $(S, T)$
- & Alt.\ \subref{fig:neigh_subgraph_alt1} cut cap. & Alt.\
- \subref{fig:neigh_subgraph_alt2} cut cap. \\
- \hline
- $(0, 0)$ & $0$ & $(\{s, u^\lambda_x, u^\lambda_y\}, \{t\})$ & $w_{xy}$ & $0$ \\
- $(0, 1)$ & $w_{xy}$ & $(\{s, u^\lambda_x\}, \{u^\lambda_y, t\})$ & $2w_{xy}$ & $w_{xy}$ \\
- $(1, 0)$ & $w_{xy}$ & $(\{s, u^\lambda_y\}, \{u^\lambda_x, t\})$ & $2w_{xy}$ & $w_{xy}$ \\
- $(1, 1)$ & $0$ & $(\{s\}, \{u^\lambda_x, u^\lambda_y, t\})$ & $w_{xy}$ & $0$
- \end{tabular}
- \label{tab:neigh_energy}
-\end{table}
-
-Figure~\ref{fig:big_graph} shows a visualization of how the final
-graph might look with all its edges. The source will have a lot of
-outgoing edges, one for each pixel, while the sink has one incoming edge
-for each pixel. The vertices corresponding to the pixels are only
-connected to the source, the sink, and their neighboring pixels, so
-their edge degree is much smaller than for $s$ and $t$.
-
-\begin{figure}[bh]
- \input{fig/big_graph}
-\end{figure}
-
-\chapter{Maximum flow approach}
-
-In the previous section we have seen how finding the minimum cut of
-carefully constructed graph can give us the thresholded image
-minimizing the energy function for one level value $\lambda$. We will in
-this and the next section see how such a minimum cut can be found by
-sending flow through the graph and trying to identify the
-``bottleneck''. This chapter, except for the description of the
-Boykov--Kolmogorov algorithm is taken from my project work
-\cite{project} and is included here for completeness.
-
-\section{Flow graphs}
-
-We have already introduced capacities, and briefly mentioned the notion
-of flow as something limited by the capacity. In other words, flow is
-something we can send through our graph, but the capacity limits how
-much we can send along each edge. It is useful to imagine a water supply
-graph with pipes of different sizes and a water source and sink.
-
-Formally, we introduce the \emph{flow} as the function $f : V \times V
-\to [0, \infty)$. This function keeps count of how much flow we are
-sending through each edge of our graph and must satisfy the following
-two constraints
-\begin{description}
- \item[Capacity constraint:] For all $u, v \in V, 0 \leq f(u, v) \leq c(u,
- v)$, i.e., for every pair of vertices, the flow is less than or
- equal to the capacity.
- \item[Flow conservation:] For all $u \in V - \{s, t\}$
- \begin{equation}
- \sum_{v \in V} f(v, u) = \sum_{v \in V} f(u, v),
- \end{equation}
- i.e., for any vertex except the source and the sink, the
- flow into the vertex must be equal to the flow out of the
- vertex.
-\end{description}
-Note that we have defined $f$ with all pairs of vertices as its domain,
-even though it is only non-zero on edges $(u, v) \in E$. This makes it
-easier to write sums as in the flow conservation constraint. We say that
-an edge $(u, v)$ is \emph{saturated} if $f(u, v) = c(u, v)$.
-
-We define $\abs{f}$ as the net amount of flow from the source to the
-sink in the graph. Because of the flow conservation constraint, this
-can be calculated as the net amount of flow going out of the source
-\begin{equation}
- \abs{f} = \sum_{v \in V} f(s, v) - \sum_{u \in V} f(u, s).
- \label{eq:total_flow}
-\end{equation}
-Furthermore we denote the \emph{net flow} across an $s$-$t$-cut $C = (S,
-T)$ as
-\begin{equation}
- f(S, T) = \sum_{u \in S} \sum_{v \in T} f(u, v) - \sum_{v \in T}
- \sum_{u \in S} f(v, u).
- \label{eq:net_flow_cut}
-\end{equation}
-Note how this definition differs from the capacity of a cut $c(S,T)$ in
-Definition~\ref{def:s_t_cut}. While the capacity of a cut represents how
-much flow it is maximally possible to send from $S$ to $T$, the net flow
-across a cut represents the net amount of flow going across the cut,
-counting negatively the flow that goes back from $T$ to $S$.
-
-For any $s$-$t$-cut we have that
-\begin{equation}
- \abs{f} = f(S, T).
- \label{eq:total_flow_eq_cut}
-\end{equation}
-This is quite intuitive given the flow conservation constraint, and a
-full chain of arguments can be found in \cite{cormen2009introduction}.
-
-Following from \eqref{eq:net_flow_cut} and \eqref{eq:total_flow_eq_cut},
-we find that
-\begin{equation}
- \begin{aligned}
- \abs{f} &= \sum_{u \in S} \sum_{v \in T} f(u, v)
- - \sum_{v \in T} \sum_{u \in S} f(v, u) \\
- &\leq \sum_{u \in S} \sum_{v \in T} f(u, v) \\
- &\leq \sum_{u \in S} \sum_{v \in T} c(u, v) \\
- &\leq c(S, T)
- \label{eq:flow_leq_cut}
- \end{aligned}
-\end{equation}
-for any $s$-$t$-cut $C = (S, T)$. A very central result in graph flow
-theory called the max-flow min-cut theorem will be presented later. It
-states that the inequality of \eqref{eq:flow_leq_cut} becomes an equality
-when $f$ is a maximum flow for some cut $C = (S, T)$, and that all such
-cuts are minimum cuts.
-
-%Recall that we want to find the minimum $s$-$t$ cut in our graph. When
-%finding this minimum cut, we make use of an important duality theorem in
-%graph flow theory, stating that the capacity of a minimum $s$-$t$-cut
-%in a graph, is equal to the maximum flow from the source to the sink.
-
-But how does this help us? We know that if we know the maximum flow
-value, and we have an $s$-$t$-cut with capacity equal to the maximum
-flow, we actually have a minimum cut. The question is then, how do we
-find a maximum flow, and how do we find a minimum cut?
-
-\section{Augmenting path algorithms}
-The family of augmenting flow algorithms represent a popular approach to
-the maximum flow problem. The idea is simply to look for paths from the
-source to the sink that can carry additional flow, so-called augmenting
-paths, and then send the maximum possible amount of flow along such a
-path. When no such path exists anymore, no more flow can be sent from
-the source to the sink, and a maximum flow has been reached.
-
-\subsection{Residual graph}
-When further discussing approaches to solving the maximum flow problem
-we will need the notion of a residual graph $G_f = (V_f, E_f, c_f)$,
-which is derived from the original graph $G$ and contains the edges
-along which it is possible to send additional flow. This means that
-$E_f$ contains the edges $(u, v)$ from $E$ where $f(u,v) < c(u,v)$. But
-that is not all; an important realisation is that it is also possible to
-push flow \emph{back} along an edge where the flow is already positive.
-In other words, sending flow from $v$ to $u$ by cancelling some or all
-of the flow that is already going from $u$ to $v$.
-
-Thus the capacity function $c_f$ of our residual graph becomes
-\begin{equation}
- c_f(u, v) =
- \begin{cases}
- c(u, v) - f(u, v) & \text{if } (u, v) \in E, \\
- f(v, u) & \text{if } (v, u) \in E, \\
- 0 & \text{otherwise.}
- \end{cases}
-\end{equation}
-The vertices $V_f$ of $G_f$ are the same as the original graph $G$,
-while the edges $E_f$ are taken to be all pairs of vertices $(u, v)$
-with $c_f(u, v) > 0$. Since we can in $E_f$ at most have all the
-original edges, and their reversals, we have $\abs{E_f} \leq 2 \abs{E}$.
-
-Note that there is ambiguity in the definition of the residual graph
-in the case where the original graph contains anti-parallel edges. One
-could avoid this by defining $c_f(u, v) = f(v, u) + c(u, v) - f(u, v)$
-instead, or as they do in \cite{cormen2009introduction}, disallow
-anti-parallel edges altogether. In any case it is not something we will
-have to think about in the implementation, since we will not actually
-construct the residual graph.
-
-With the residual graph defined, we are ready to formally present the
-max-flow min-cut theorem.
-\begin{theorem}[Max-flow min-cut theorem]
- If $f$ is a flow in a graph $G = (V, E, c)$ with source $s$ and
- sink $t$, then the following statements are equivalent:
- \begin{enumerate}
- \item $f$ is a maximum flow in $G$.
- \item The residual graph $G_f$ contains no augmenting paths.
- \item $\abs{f} = c(S, T)$ for some cut $(S, T)$ of $G$.
- \label{it:max_flow_min_cut}
- \end{enumerate}
-\end{theorem}
-See \cite{cormen2009introduction} for a proof, and remark that because
-of the inequality in \eqref{eq:flow_leq_cut}, the cut in
-Statement~\ref{it:max_flow_min_cut} is a minimum cut. The theorem does
-not tell us how to find such a cut, and there are multiple ways. One
-possibility is to take $S$ to be all vertices reachable from the source
-in the residual graph and $T = V - S$.
-
-Figure~\ref{fig:aug_flow} shows a simple graph which already has five
-units flowing from $s$ to $t$. The marked path is a possible augmenting
-path, and note that it follows an edge in $E$ in the reverse direction,
-made possible by the construction of the residual graph just
-described.
-
-\begin{figure}
- \input{fig/aug_flow}
-\end{figure}
-
-\subsection{Ford-Fulkerson}
-The Ford-Fulkerson algorithm is the most basic augmenting path
-algorithm, which can be extended to more advanced algorithms. It is
-stated in pseudocode in Algorithm~\ref{alg:ford_fulkerson}, and the idea
-is to augment the flow along paths from $s$ to $t$ until it is no longer
-possible.
-\begin{algorithm}
-\begin{algorithmic}
- \Function{Ford-Fulkerson}{$G$, $s$, $t$}
- \While{there exists a path $p$ from $s$ to $t$ in the residual
- graph $G_f$}
- \State $\alpha \gets \min\{c_f(u, v) : (u, v) \in p \}$
- \ForAll{$(u, v) \in p$}
- \If{$(u, v) \in E$}
- \State $f(u, v) \mathrel{+}= \Delta f$
- \Else
- \State $f(v, u) \mathrel{-}= \Delta f$
- \Comment{Push flow back}
- \EndIf
- \EndFor
- \EndWhile
- \EndFunction
-\end{algorithmic}
-\caption{\sf The Ford-Fulkerson max-flow algorithm}
-\label{alg:ford_fulkerson}
-\end{algorithm}
-
-There are different ways to find augmenting paths, and a common choice
-is to do a breadth-first search from the source until the sink is found,
-as this will yield the shortest possible augmenting path. This version
-of the algorithm is called Edmonds-Karp and has a running time of
-$O(\abs{V}\abs{E}^2)$. See for example \cite{cormen2009introduction}
-for a description of the breadth-first search, and a formal proof of the
-running time of the Edmonds-Karp algorithm.
-
-\section{Other algorithms}
-
-There are many different maximum flow algorithms that fall into the
-augmenting path category, although we will see a different approach in
-the next section.
-
-The algorithm of Dinitz, originally published in 1970, later improved
-on, and described by the original author in \cite{dinitz2006dinitz}, is
-a variant of the augmenting path algorithm. It maintains a distance
-labeling $d(u)$ of the vertices $u \in V$ in the graph, where $d(u)$
-is the shortest path from the $s$ to $u$ in the residual graph. This
-can be computed with a simple breadth-first search. The next step is to
-construct a \emph{blocking flow} $f'$, using only edges in $E_f' = \{
-(u,v) \in E_f : d(u) + 1 = d(v) \}$. The blocking flow is such that if
-we augment the flow $f$ by $f'$, there is no longer any paths from $s$
-to $t$ following edges in $E_f'$. After the blocking flow $f'$ has been
-found and added to $f$, the distance labels are recalculated, and the
-label of the sink will be increased by at least one.
-
-Boykov and Kolmogorov present a variant of the augmenting path algorithm
-in \cite{boykov2004experimental}, specialized for the kinds of graphs
-occuring in graphical applications. It will be presented in
-Section~\ref{sec:bk}
-
-\section{Push--relabel algorithm}
-
-The push-relabel algorithm is a different approach to the maximum flow
-problem, presented by Goldberg and Tarjan in \cite{goldberg1988new}.
-Unlike the augmenting flow algorithms, it does not maintain a valid flow
-$f$ in the graph at all times, but still obtains a valid maximum flow
-when the algorithm terminates.
-
-\subsection{Preflow}
-Instead of maintaining a valid flow, we introduce the concept of a
-\emph{preflow} by relaxing the flow conservation constraint from
-earlier. We allow positive excess in the vertices and the flow
-conservation constraint from before then becomes
-\begin{description}
- \item[Preflow conservation:] For all $u \in V - \{s, t\}$
- \begin{equation}
- \sum_{v \in V} f(v, u) \geq \sum_{v \in V} f(u, v),
- \end{equation}
- i.e., for any vertex except the source and the sink, the
- flow into the vertex must greater or equal to the flow out of
- the vertex.
-\end{description}
-
-As in most of the cited push-relabel literature, we define $N =
-\abs{V}$, and for all vertices $u \in V$ we define the excess
-\begin{equation}
- e(u) = \sum_{v \in V} f(v, u) - \sum_{v \in V} f(u, v),
-\end{equation}
-which represents the amount of flow which \emph{disappears} in vertex
-$u$. Equivalent to the preflow conservation constraint is stating that
-$e(u) \geq 0$ for all vertices $u \in V - \{s,t\}$.
-
-%The idea of the algorithm is to maintain a height map of the
-%vertices in the graph where connected vertices can not have a large
-%height difference. Then we ``lift'' the source vertex to let as
-%many units as possible flow through the edges of the graph towards the
-%sink. When a maximum preflow is reached, there will normally be excess
-%flow in some of the vertices, which has to be pushed back towards the
-%source in order to obtain a valid flow.
-
-In addition to the flow, we maintain a height map $d : V \to \mathbb{N}$
-that satisfies $d(t) = 0$, and for every edge $(u, v)$ in the residual
-graph, i.e.\ every edge with $c_f(u, v) > 0$, we require that $d(u) \leq
-d(v) + 1$. For all vertices $u$, the label $d(u)$ will be a lower bound
-on the length from $u$ to $t$ in $G_f$ which is why it is also often
-called a distance labeling.
-
-A vertex $u$ is \emph{active} if $u \in V - \{s,t\}$, it has positive
-excess ($e(u) > 0$) and $d(u) < N$. These are the vertices we want to
-operate on to increase the preflow.
-
-\subsection{Basic operations}
-
-The algorithm performs two basic operations, the \emph{push} and
-\emph{relabel} operations, while always maintaining a valid preflow $f$
-and a valid distance labeling $d$.
-
-\subsubsection{The push procedure}
-
-The push procedure moves excess flow from an active vertex along an edge
-$(u, v) \in E_f$ for which $d(u) = d(v) + 1$, i.e.\ to a vertex with a
-smaller distance label. We call such edges \emph{admissible}. See
-Algorithm~\ref{alg:push} for a pseudocode implementation of the push
-operation.
-
-Assuming that $f$ is a valid preflow, it is easy to verify that the
-preflow $f$ and labeling $d$ remain valid after running the push
-procedure on some admissible edge $(u, v)$.
-
-%The capacity constraint is fulfilled since we
-%at most increase the flow along $(u, v)$ with the residual capacity
-%$c_f(u,v)$. The preflow constraint is fulfilled since the excess $e$
-%increases for $v$, remains non-negative for $u$ and remains the same for
-%all other vertices.
-
-%The residual graph might change during the push procedure. The edge
-%$(v, u)$ will appear if it does not already exist. Initially $d(u) =
-%d(v) + 1$, so for the new edge $(v, u)$ we have
-%\begin{equation}
-% d(v) = d(u) - 1 \leq d(u) + 1
-%\end{equation}
-%and $d$ remains a valid labeling after the push procedure is finished.
-
-\begin{algorithm}
-\begin{algorithmic}
- \Function{Push}{$u$, $v$}
- \State $\Delta f \gets \min(c_f(u, v), e(u))$
- \If{$(u, v) \in E$}
- \State $f(u, v) \mathrel{+}= \Delta f$
- \Else
- \State $f(v, u) \mathrel{-}= \Delta f$
- \Comment{Push flow back}
- \EndIf
- \Comment{Excess $e(u)$ and $e(v)$ will also change}
- %\State $f(v, u) \mathrel{-}= f_\text{aug}$
- %\State $e(u) \mathrel{-}= f_\text{aug}$
- %\State $e(v) \mathrel{+}= f_\text{aug}$
- \EndFunction
-\end{algorithmic}
-\caption{%
- \sf The push procedure of the Push-Relabel algorithm.
-}
-\label{alg:push}
-\end{algorithm}
-
-\subsubsection{The relabel procedure}
-
-The relabel procedure is our tool for changing the distance labels of
-the vertices. It changes the label of a vertex to the greatest possible
-value, which is one more than the lowest label among its neighbors in
-the residual graph. See Algorithm~\ref{alg:relabel} for a pseudocode
-implementation.
-\begin{algorithm}
-\begin{algorithmic}
- \Function{Relabel}{$u$}
-% \If{$u$ is only vertex at its height}
-% \Call{Gap}{$u$}
-% \Else
- \If{there is a $v \in V$ such that $(u, v) \in E_f$}
- \State $d(u) \gets \min\{d(v), \; \forall v \in V : (u,v) \in E_f\} + 1$
- \Else
- \State $d(u) \gets N$
- \Comment{$u$ becomes inactive}
- \EndIf
-% \EndIf
- \EndFunction
-\end{algorithmic}
-\caption{\sf The relabel procedure of the Push-Relabel algorithm}
-\label{alg:relabel}
-\end{algorithm}
-
-If $d$ was a valid labeling before running the relabel procedure, then
-we still have $d(u) \leq d(v) + 1$ for all neighbors $v$ of $u$ in the
-residual graph, and $d$ remains a valid labeling. The capacity
-constraint and preflow constraint remain satisfied assuming they were
-satisfied before the procedure was started.
-
-\subsection{Putting it all together}
-These basic procedures are then applied to active vertices and admissible
-edges until we obtain our minimum cut. We will see later that when there
-are no more active vertices, we can extract the minimum cut from the
-graph.
-
-In the first phase of the algorithm we initialize a valid preflow and
-distance labeling by saturating all edges out of the source $s$, and
-then setting its distance label $d(s) = N$. We then apply the push and
-relabel procedures where applicable until there are no more active
-vertices and a maximum preflow is obtained.
-
-A vertex $u$ can only be successfully relabeled to obtain a new label if
-the outgoing edges of $u$ in the residual graph have changed since the
-previous relabeling. This is why the push and relabel procedures often
-are combined into a \emph{discharge} procedure as shown in
-Algorithm~\ref{alg:discharge}. When it is run on an active vertex $u$,
-we push as much as possible of the excess flow to other vertices before
-the vertex is relabeled.
-\begin{algorithm}
-\begin{algorithmic}
- \Function{Discharge}{$u$}
- \ForAll{$v \in V$ such that $(u, v) \in E_f$}
- \If{$c_f(u, v) > 0$ and $d(u) = d(v) + 1$}
- \Call{Push}{$u$, $v$}
- \EndIf
- \EndFor
-
- \If{$e(u) > 0$}
- \Call{Relabel}{$u$}
- \EndIf
- \EndFunction
-\end{algorithmic}
-\caption{\sf The discharge procedure of the Push-Relabel algorithm}
-\label{alg:discharge}
-\end{algorithm}
-
-In the second phase of the algorithm this maximum preflow is turned into a
-maximum flow by returning excess flow which did not reach the sink from
-inside the graph back to the source. We can skip this part of the
-algorithm, as it is possible to identify a minimum cut as soon as the
-first phase is finished, and the following theorem allows us to do that.
-\begin{theorem}[Cut identification]
- Given a graph $G = (V, E, c)$, assume that the first phase of
- the push-relabel algorithm has terminated so that no more active
- vertices remain. Then there exists a $k \in \{1, \ldots, N-1\}$
- such that there is no vertex with label $k$. For every
- such $k$ the vertex sets $S = \{ u \in V : d(u) > k\}$ and $T = \{ u
- \in V : d(u) < k \}$ define a minimum cut $C = (S, T)$ in $G$.
- \label{thm:cut_identification}
-\end{theorem}
-\begin{proof}
- There are $N$ vertices, the source has label $N$ and the sink has
- label $0$, and the $N - 2$ remaining vertices can not occupy all the
- $N-1$ labels in $\{1, \ldots, N-1\}$, so there must exist an $k$ as
- described.
-
- There can be no edge $(u, v) \in E_f$ such that $u \in S$ and $v \in
- T$, as this would imply $k \leq d(u) - 1 \leq d(v) < k$. From the
- construction of $E_f$ we now know that all edges in $E$ from $S$ to
- $T$ are saturated, and all edges from $T$ to $S$ carry no flow. This
- implies the capacity of the cut is equal to the flow through the
- cut, i.e.\ $c(S, T) = f(S, T)$.
-
- Since the first phase of the algorithm has terminated, there can be
- no active vertices, and therefore no excess in $T$, except for the
- sink. If all flow excess in vertices in $S$ is returned to the
- source, we can apply the max-flow min-cut theorem to conclude that
- $C = (S, T)$ is a minimum $s$-$t$-cut, since the cut capacity is
- equal to the flow.
-\end{proof}
-
-We will see later that with the gap relabeling heuristic, there will
-always be a gap at label $k = N - 1$ such that we can construct our cut
-by taking $S = \{ u \in V : d(u) \geq N \}$.
-
-Note that the vertices in $S$ are vertices earlier described as being on
-the source side of the cut, as no additional flow can go from these
-vertices to the sink.
-
-\subsection{Complexity}
-
-In their original article \cite{goldberg1988new}, Goldberg and Tarjan
-analyze the complexity of the push-relabel algoritm by considering the
-maximum number of basic operations we can possibly do before the
-algorithm terminates.
-
-The number of relabelings is in $O(\abs{V}^2)$ since every time the
-procedure is applicable to a vertex $u$, the label $d(u)$ increases by
-at least one.
-
-The number of saturating pushes is in $O(\abs{V}\abs{E})$. When a push
-along $(u,v)$ is saturating, the label of $v$ has to increase with at
-least 2 before a push can saturate the same edge (in the opposite
-direction). Since the number of relabelings of a single vertex is
-bounded by $\abs{V}$, and we have $\abs{E}$ edges, this gives the stated
-number of saturating pushes.
-
-The number of non-saturating pushes is the most complicated to bound,
-and will also make up the asymptotic running time of the algorithm. The
-idea is to define
-\begin{equation}
- \phi = \sum_{\mathclap{u \text{ active}}} d(u),
-\end{equation}
-and look at how much this number changes throughout the algorithm. It
-starts at zero and ends at zero. Every non-saturating push from $u$ to
-$v$ makes $\phi$ decrease by at least one since it makes $u$ inactive
-(but might activate $v$). The total increase in $\phi$ due to
-relabelings is less than $\abs{V}^2$. A saturating push from $u$ to $v$
-increases $\phi$ by at most $\abs{V}$, since $v$ might become active.
-
-Even if $\phi$ is always increased by relabelings and saturating pushes,
-we can bound the number of non-saturating pushes by
-\begin{equation}
- \abs{V}^2 + \abs{V} \underbrace{c \abs{V}
- \abs{E}}_{\mathclap{\#(\text{saturating pushes})}}
-\end{equation}
-which means that in the general case, the algorithm has a complexity of
-$O(\abs{E}\abs{V}^2)$.
-
-\subsection{Vertex selection rules}
-Until now we have stated that the discharge procedure is run on active
-vertices until there are no more active vertices left. The choice of the
-order in which to discharge these active vertices remain, and multiple
-possibilities exist.
-
-The First In, First Out (FIFO) approach is to always maintain a queue of
-active vertices. When the vertex from the beginning of the queue is
-discharged, other vertices might become active, and these are added at
-the end of the queue. The original article of Goldberg and Tarjan
-\cite{goldberg1988new} contains a proof that this selection rule gives a
-complexity of $O(\abs{V}^3)$, which can be very good if you have a dense
-graph.
-
-Using highest level selection rule one always discharges the active
-vertex with the largest distance label. Goldberg and Tarjan state that
-this rule also gives a complexity of $O(\abs{V}^3)$ while this bound is
-improved to $O(\abs{V}^2 \sqrt{\abs{E}})$ in an article by Cheriyan and
-Maheshwari \cite{cheriyan1989analysis}.
-
-These are complexity bounds, and the actual running time of the
-algorithm, which can only be determined by implementing it and running
-it, varies a lot with the structure of the input graph.
-
-Cherkassky and Goldberg describe the algorithm along with different
-selection rules, heuristics and their implementation in
-\cite{cherkassky1997implementing}.
-
-\subsection{Heuristics}
-Different heuristics exist that can speed up the algorithm
-considerably. Being heuristics, they are not guaranteed to work, and
-might perform differently on different kinds of graphs. The most used
-heuristics are the gap and global relabeling heuristics, both aiming to
-reduce the total number of relabeling steps.
-
-The gap relabeling heuristic aims to find a label $k$ as in
-Theorem~\ref{thm:cut_identification} such that no vertex has that label.
-From vertices $u$ with $d(u) > k$, there are no unsaturated edges going
-to vertices with smaller distance labels, so no more flow can ever find
-its way from these vertices to the sink. These vertices are therefore
-given the label $N$ and never considered again as they will never become
-active. Algorithm~\ref{alg:gap} shows a pseudocode representation of
-what is done once a gap $k$ is found.
-
-\begin{algorithm}
-\begin{algorithmic}
- \Function{Gap}{$k$}
- \ForAll{$u \in V$ such that $d(u) > k$}
- \State $d(u) \gets N$
- \EndFor
- \EndFunction
-\end{algorithmic}
-\caption{\sf The gap procedure of the Push-Relabel algorithm}
-\label{alg:gap}
-\end{algorithm}
-
-The preflow and capacity constraints are still valid after the gap
-relabeling procedure, as only the labels $d$ are changed. For the
-labeling $d$ one has to verify that $d(u) \leq d(v) + 1$ for every edge
-$(u,v) \in E_f$ in the residual network. If none, both, or only $v$ is
-relabeled, this is trivial. It is not possible that only $v$ would be
-relabeled, as this would imply that $d(u) \geq d(v) + 2$ which is not a
-valid initial labeling.
-
-%Before integrating the gap relabeling procedure into our algorithm we
-%need to verify that it does not invalidate our preflow $f$ or distance
-%labeling $d$.
-
-%But why does this work? The only thing we need to verify is that given a
-%graph with a valid preflow and a valid labeling, the gap relabeling
-%procedure will not change the validity of these two things.
-
-%\begin{lemma}[Gap relabeling validity]
-% Given a graph $G = (V, E, c)$, a distance labeling $d$ and a
-% preflow $f$, assume there exists a gap $k$ such that no vertex has
-% label $k$. Running the gap relabeling procedure on label $k$ will
-% yield a valid distance labeling and an unchanged and valid preflow
-% $f$.
-%\end{lemma}
-%\begin{proof}
-% No new edges are created, no edges disappear, the preflow is
-% unchanged, so the preflow and capacity constraint remain satisfied
-% after the gap relabeling.
-%
-% Define $S = \{ u \in V : d(u) > k\}$ and $T = V - S$. Initially
-% $d(u) \leq d(v) + 1$ for every edge $(u,v) \in E_f$. These
-% inequalities have to hold after the gap procedure is finished, when
-% $d(u) = N$ for all $u \in S$.
-%
-% For $(u, v) \in E_f$ we have four possibilities
-% \begin{description}
-% \item[$u, v \in T$\textmd{:}]
-% The labels $d(u)$ and $d(v)$ remain unchanged and the
-% inequality still holds.
-% \item[$u, v \in S$\textmd{:}]
-% After the gap procedure we have $d(u) = d(v)$ so the
-% inequality still holds.
-% \item[$u \in S, v \in T$\textmd{:}]
-% This is not possible as it would imply $d(u) \geq d(v) + 2$
-% and we have assumed an initial valid labeling.
-% \item[$u \in T, v \in S$\textmd{:}]
-% After relabeling we have $d(u) < k < N < d(v) + 1$.
-% \end{description}
-%
-% Hence, both the preflow $f$ and distance labeling $d$ are valid.
-%\end{proof}
-
-When running the push-relabel algorithm with the gap heuristic, we can
-be sure that there will never be a vertex $u$ with label $d(u) = N-1$ at
-the end of the algorithm, i.e.\ we know that there will always be a gap
-at label $N-1$. This can be seen using the same reasoning as in
-Theorem~\ref{thm:cut_identification}, because if there was a vertex with label
-$N-1$, there would only be $N-3$ vertices possibly having labels in
-$\{1, \ldots, N-2\}$, so a gap must exist somewhere in that interval.
-When using the gap relabeling heuristic, such a gap can not exist, so we
-can conclude that there is no vertex with label $N-1$.
-
-Using Theorem~\ref{thm:cut_identification} we can then conclude that the
-sets $S = \{ u \in V : d(u) \geq N\}$ and $T = V - S$ form a minimum cut
-of the graph.
-
-\subsection{Parametric push-relabel algorithm}
-
-Now we have an algorithm for finding a minimum $s$-$t$-cut in a graph,
-so let's return to the graph constructed in
-Section~\ref{sec:graph_construction}. For every level $\lambda \in \{0,
-\ldots, L\}$ we want to find a minimum $s$-$t$-cut which gives us the
-thresholded image $u^\lambda$. These can then hopefully be stacked
-together to form the final image $u$.
-
-\subsubsection{Graph reuse}
-
-Solving $L$ separate minimum cut problems seems like a lot of work, but
-when using the push-relabel algorithm we will, if we do things in the
-right order, be able to reuse the graph when going from one label to
-the next.
-
-\begin{figure}
- \input{fig/norm_evolution}
-\end{figure}
-
-Going back to the graph representations in
-Figure~\ref{fig:norm_subgraph} and Figure~\ref{fig:neigh_subgraph} we
-know that only the capacity of edges from sub-graphs representing the
-fidelity term depend on our level parameter $\lambda$. From the
-expression in \eqref{eq:fidelity_energy0}, visualized in
-Figure~\ref{fig:norm_evolution}, we see that the energy term
-$E_\lambda^x(1)$ increases monotonically with increasing $\lambda$
-parameter. Let $u, v \in V - \{s, t\}$. Since the edges in
-Figure~\ref{fig:norm_subgraph} now are the only ones depending on
-$\lambda$, the following is true for \emph{decreasing} values of
-$\lambda$
-\begin{description}
- \item[Edges from $s$ to $u$\textmd{:}]
- As seen in Figure~\ref{fig:norm_subgraph_neg} the capacity of
- these edges will increase monotonically with decreasing
- $\lambda$ parameter.
- \item[Edges from $u$ to $v$\textmd{:}]
- These edges have no $\lambda$-dependence and will remain
- unchanged.
- \item[Edges from $v$ to $t$\textmd{:}]
- As seen in Figure~\ref{fig:norm_subgraph_pos} the capacity of
- these edges will decrease monotonically with decreasing
- $\lambda$ parameter.
-\end{description}
-
-After running the push-relabel algorithm for $\lambda = k$, we are left
-with a graph $G = (V, E, c)$, a preflow $f$ and a labeling $d$. To
-obtain the graph for $\lambda = k-1$ we have to change the capacity of
-two different kinds of edges, and this is done in the following way to
-keep the capacity, preflow and labeling constraints satisfied.
-\begin{description}
- \item[Edges from $s$ to $u$\textmd{:}]
- The capacity $c(s, u)$ is increased, and the flow is set to be
- equal to the capacity $f(s, u) = c(s, u)$. The vertex $u$ might
- have an increased excess $e(u)$, which might in turn make it
- active.
- \item[Edges from $v$ to $t$\textmd{:}]
- The capacity $c(v, t)$ is decreased. If it is decreased to a
- value below the current flow value, we set $f(v, t) = c(v, t)$
- which will decrease the excess of the sink $t$, and increase the
- excess of $v$.
-\end{description}
-
-None of these actions will create new edges in the residual graph, and
-we do not change the labeling $d$, so the labeling constraints are also
-satisfied in the new graph.
-
-Through this procedure we have easily created the graph for $\lambda =
-k-1$, and the distance labels remain the same. As these labels always
-increase monotonically, we have a head start compared to the case where
-we reset the flow and labels.
-
-\subsubsection{Output image construction}
-
-We mentioned already in Section~\ref{sec:total_energy} that in order to
-be able to construct our output image $u$, the thresholded images
-$u^\lambda$ would have to stack one on top of the other as shown in
-Figure~\ref{fig:img_decomp}. Because of the reuse of the distance labels
-between the iterations of the push-relabel algorithm, we can guarantee
-that this is possible.
-
-Consider two subsequent runs of the push-relabel algorithm, for labels
-$\lambda$ and $\lambda - 1$ ending with distance labels $d^\lambda$ and
-$d^{\lambda-1}$ respectively. We already know that the distance labels
-$d$ are monotonically increasing. This means that the set $S = \{ u \in
-V : d(u) \geq N \}$ is increasing in size, more precisely, we have the
-inclusion
-\begin{equation}
- \{ u \in V : d^\lambda(u) \geq N \} \subseteq \{ u \in V :
- d^{\lambda-1}(u) \geq N \}.
-\end{equation}
-For a pixel $x \in S$ we will set $u^\lambda_x = 0$, which together with
-the previous inclusion property implies
-\begin{equation}
- u^\lambda_x \geq u^{\lambda-1}_x
-\end{equation}
-for all $x \in \mathcal{G}$. Being equivalent with the inequality in
-\eqref{eq:stackable}, this means our algorithm produces stackable
-thresholded images $u^\lambda$.
-
-We then construct our output image $u$ by giving each pixel the value
-\begin{equation}
- u_x = \min \{ \lambda \in \{0, \ldots, L-1\} : u^\lambda_x = 1 \}.
-\end{equation}
-This marks the end of the description of the implemented algorithm, but
-we will further discuss some possible improvements, and also look at the
-results when using the method on different kinds of noisy images.
-
-\subsection{Divide and conquer}
-The possibility of re-using the graph between separate level
-is a very nice property of the push-relabel algorithm, but there are
-further room for improvements. Consider one pixel $x$ with value $u_x$,
-and imagine we only wanted to find the value of this pixel. One could go
-through all pixel values $\lambda \in (L-1, \ldots, 0)$, and see when
-$u^\lambda_x$ changes from $1$ to $0$, just as we do for all the pixels
-in the algorithm above. Ignoring graph re-use this would have us solve
-$O(L)$ maximum flow problems.
-
-Improving on this we could employ the idea of binary search to find the
-value of $u_x$ in only $O(\log_2 L)$ time. After finding one cut, we
-know whether $u_x$ is above or below the current $\lambda$ value, and by
-choosing $\lambda$ as the midpoint of the current possible range of
-$u_x$, we can cut the search space in half for each iteration of the
-algorithm.
-
-We can extend this idea to the problem of finding all pixel values.
-Instead of running the algorithm for successively decreasing values of
-$\lambda$, we choose some $\lambda$ in the middle of the range $\{0,
-\ldots, L-1\}$. The cut we obtain consists of two sets $S = \{ u \in V :
-d(u) \geq N\}$ and $T = V - S$. We know that no more flow can be sent
-from $S$ to $T$, even if we decrease the value of $\lambda$ and adjust
-the capacities accordingly.
-
-The idea is now that we have halved the possible $\lambda$ interval for
-\emph{all} pixels. We continue by considering the two sets $S$ and $T$
-separately, and applying the algorithm recursively, at each time halving
-the $\lambda$ interval until we have the value of every pixel.
-
-Combining the divide and conquer approach with the parametric
-push-relabel algorithm is a bit problematic. For all pixels $x \in T$ we
-know that $u_x \leq \lambda$, and we can reuse the graph when
-decreasing $\lambda$. However for pixels $x \in S$, we seek to find $u_x
-> \lambda$, meaning we have to increase $\lambda$ which does not allow
-graph re-use.
-
-Goldfarb and Yin \cite{goldfarb2009parametric} have found that the
-divide and conquer approach only yields improved performance when using
-the $L^2$ norm in the fidelity term. This has to do with the fact that
-the fidelity term for the $L^1$ norm only changes once for each pixel,
-as we can see in Figure~\ref{fig:norm_evolution}, reducing the amount of
-work that has to be done in each iteration of the regular parametric
-push-relabel algorithm.
-
-See \cite{gallo1989fast}, \cite{hochbaum2001efficient} and
-\cite{goldfarb2009parametric} for more information.
-
-\section{Boykov--Kolmogorov algorithm}
-\label{sec:bk}
-
-\begin{algorithm}
- \begin{algorithmic}
- \Function{BoykovKolmogorov}{$G$, $s$, $t$}
- \State $A := {s, t}$, $O := \emptyset$, $S = {s}$, $T = {t}$
- \State $s\text{.color} = S$, $t\text{.color} = T$
- \While{\texttt{True}}
- \State $e \gets$ \Call{Grow}{$G$, $A$}
- \If{not $e$}
- \State break
- \EndIf
- \State \Call{Augment}{$G$, $s$, $t$, $e$, $O$}
- \State \Call{Adopt}{$G$, $s$, $t$, $O$}
- \EndWhile
- \EndFunction
- \end{algorithmic}
- \caption{\sf The Boykov--Kolmogorov maximum flow algorithm}
- \label{alg:boykov_kolmogorov}
-\end{algorithm}
-
-A maximum flow algorithm specialized for the type of graphs found in
-imaging applications is described by Boykov and Kolmogorov in
-\cite{boykov2004experimental}. Their algorithm is an augmenting path
-algorithm where the paths are found using trees that keep track of
-possible paths.
-
-There is one tree denoted $S$ with the source $s$ as its root, and
-another tree denoted $T$ with the sink $t$ as its root. These trees are
-non-overlapping, and all edges in the tree $S$ can carry flow towards
-the leafs of the tree, while all edges in $T$ can carry flow towards the
-sink $t$.
-
-A vertex can either be part of these trees, or be a free vertex, and we
-write
-\begin{equation}
- u\text{.color} = \begin{cases}
- S & \text{ if } u \in S \\
- T & \text{ if } u \in T \\
- \texttt{Free} & \text{ otherwise.}
- \end{cases}
-\end{equation}
-The idea of the algorithm is then to grow these trees towards eachother,
-and when they touch, we have found an augmenting path from $s$ to $t$.
-As much flow as possible is then pushed along this path, and some edges
-will become saturated. If an edge in the tree $S$ or $T$ becomes
-saturated, the tree is split up into multiple trees, also called a
-forest. The parts no longer connected to the root are called orphans,
-and we write $u \in O$ for an orphan vertex $u$. These vertices are---if
-possible---joined back into the tree, or become free nodes.
-
-The vertices in the trees $S$ and $T$ are either active or passive, and
-we write $u \in A$ if $u$ is active. The active vertices are those at
-the boundary of the tree, which can possibly connect to other vertices
-to grow the tree. The passive vertices are internal in the tree, and
-edges to their neighbors are either a part of the tree, or completely
-saturated.
-
-The main loop of the algorithm is then as shown in
-Algorithm~\ref{alg:boykov_kolmogorov}. In the grow procedure, the trees are grown
-from their active vertices until a new path is found. The flow
-is then augmented along this path in the augment procedure, and orphan
-vertices might be created. These orphan vertices are then either joined
-back into their respective trees, or become free, in the adopt
-procedure.
-
-\begin{algorithm}
-\begin{algorithmic}
- \Function{Grow}{$G$, $A$}
- \While{$\abs{A} \neq 0$}
- \State $u \gets $ one node from $A$
- \ForAll{$v$ such that \Call{treeCap}{$u$, $v$} > 0}
- \If{$v\text{.color} = \texttt{Free}$}
- \State $v\text{.color} \gets u\text{.color}$
- \State $v\text{.parent} \gets u$
- \State $A := A \cup \{v\}$
- \ElsIf{$v\text{.color} \neq u\text{.color}$}
- \State \Return $e := (u, v)$
- \EndIf
- \EndFor
- \State Remove $u$ from $A$
- \EndWhile
- \State \Return \texttt{False}
- \EndFunction
-\end{algorithmic}
-\caption{\sf The grow procedure of the Boykov--Kolmogorov maximum flow
-algorithm}
-\label{alg:grow}
-\end{algorithm}
-
-The grow procedure is shown in Algorithm~\ref{alg:grow}. It goes through
-the set of active vertices $A$ to try to expand the trees $S$ and $T$.
-When considering an active vertex $u \in S$, we want to grow by finding
-vertices $v$ such that $c_f(u, v) > 0$, while when considering an active
-vertex $u \in T$, we want to grow by finding vertices $v$ such that
-$c_f(v, u) > 0$. This is why the \textsc{treeCap} is introduced which is
-defined as
-\begin{equation}
- \textsc{treeCap}(u, v) = \begin{cases}
- c_f(u, v) & \text{ if } u \in S, \\
- c_f(v, u) & \text{ if } u \in T.
- \end{cases}
-\end{equation}
-Thus if a non-saturated edge is found from an active vertex $u$ to a free
-vertex $v$, then $v$ is added to the tree of $u$. The tree structure is
-stored by keeping a parent attribute in each non-free node, such that
-for example $v\text{.parent} = u$.
-
-If a non-saturated edge is found from an active vertex $u$ to a vertex
-$v$ in the other tree, the two trees connect, and we can return the edge
-$e$ since we have an augmenting path. Note that if this happens, the
-vertex $u$ is still active. It only becomes passive when all neighbors
-are considered without yielding an augmenting path.
-
-\begin{algorithm}
- \begin{algorithmic}
- \Function{Augment}{$G$, $s$, $t$, $e$, $O$}
- \State $p \gets $ path from $s$ to $t$ through $e$
- \Comment through the trees $S$ and $T$
- \State $\Delta f \gets \min \{c_f(u, v) : (u,v) \in p\}$
- \ForAll{$(u, v) \in p$}
- \State $f(u, v) \mathrel{+}= \Delta f$
- \State $f(v, u) \mathrel{-}= \Delta f$
- \If{$c_f(u, v) = 0$}
- \If{$u\text{.color} = v\text{.color} = S$}
- \State $v\text{.parent} = \texttt{Null}$
- \State $O := O \cup \{v\}$
- \ElsIf{$u\text{.color} = v\text{.color} = T$}
- \State $u\text{.parent} = \texttt{Null}$
- \State $O := O \cup \{u\}$
- \EndIf
- \EndIf
- \EndFor
- \EndFunction
- \end{algorithmic}
- \caption{\sf The augment procedure of the Boykov--Kolmogorov maximum
- flow algorithm}
- \label{alg:augment}
-\end{algorithm}
-
-The augment procedure finds the augmenting path going through the tree
-$S$, the edge $e$ and the tree $T$. The maximal possible flow is then
-sent along this path. At least one edge will then become saturated. If a
-saturated edge occurs in the trees $S$ or $T$, the edge terminal
-farthest from the tree root is marked as an orphan. Note that even if a
-whole subtree is disconnected from the main tree, only the root of this
-subtree is marked as an orphan. The adopt procedure will take care of
-either reconnecting, or freeing all the vertices in the disconnected
-subtree.
-
-\begin{algorithm}
-\begin{algorithmic}
- \Function{Adopt}{$G$, $s$, $t$, $O$}
- \While{$\abs{O} \neq 0$}
- \State $u \gets $ one node from $O$
- \State Remove $u$ from $O$
- \State found $\gets \texttt{False}$
- \ForAll{$v$ such that \Call{treeCap}{$v$, $u$} > 0}
- \If{$u\text{.color} \neq v\text{.color}$}
- \State continue
- \EndIf
- \State orig $\gets$ \Call{treeOrigin}{$v$}
- \If{orig $\neq s$ and orig $\neq t$}
- \State continue
- \EndIf
- \State found $\gets \texttt{True}$
- \State $u\text{.parent} \gets v$
- \State break
- \EndFor
- \If{found $\neq \texttt{True}$}
- \ForAll{$v$ such that $u\text{.color} = v\text{.color}$}
- \If{\Call{treeCap}{$v$, $u$} > 0}
- \State $A := A \cup \{ v \}$
- \EndIf
- \If{$v\text{.parent} = u$}
- \State $O := O \cup \{ v \}$
- \State $v\text{.parent} \gets \texttt{Null}$
- \EndIf
- \State $u\text{.color} \gets \texttt{Free}$
- \State $A := A - \{ u \}$
- \EndFor
- \EndIf
- \EndWhile
- \EndFunction
-\end{algorithmic}
-\caption{\sf The adopt procedure of the Boykov--Kolmogorov maximum flow
-algorithm}
-\label{alg:adopt}
-\end{algorithm}
-
-The adopt procedure processes all the vertices in the set of orphans
-$O$. These vertices are either single vertices, or roots of disconnected
-subtrees. For an orphaned vertex $u$ we look through its neighbors to
-find a possible parent vertex $v$. It should belong to the same tree as
-$u$, and also satisfy $\textsc{treeCap}(v, u) > 0$. In addition, $v$
-should be connected to one of the tree roots $s$ or $t$. This keeps us
-from connecting to other orphaned vertices or subtrees, and is checked
-by the procedure $\textsc{treeOrigin}(v)$, which follows the parent
-information until reaching either $s$, $t$ or an orphaned vertex.
-
-If a potential parent $v$ of $u$ is \emph{not} found, $u$ becomes a free
-vertex. All vertices that had $u$ as its parent are orphaned, and are thus
-treeted by the adopt procedure later. Vertices that are in one of the
-trees $S$ or $T$, and have a non-saturated edge to this newly freed
-vertex $u$, become active.
-
-When no more active vertices remain, there are no more possible
-augmenting paths, and the algorithms terminates with a maximum flow.
-A proof of correctness can be found in \fixme{kolmogorov-phd}.
-
-\subsection{Graph reuse}
-
-As in the push-relabel algorithm described earlier, it is not necessary
-to completely restart the Boykov--Kolmogorov algorithm for every level
-of the image, but we can reuse the graph and the trees $S$ and $T$ in
-successive runs.
-
-\begin{figure}
- \input{fig/bk_norm_subgraph}
-\end{figure}
-
-Since this is algorithm does not deal with the relaxed preflow concept,
-the updates from one level to the next has to be done in a different
-way. We can no longer decrease the capacity of edges, as this could
-break the flow conservation constraint. However, the graph construction
-in \fixme{ref} did allow for the addition of an arbitrary constant to
-all the edges of the sub-graph. Thus the construction in
-Figure~\ref{bk_norm_subgraph} is also valid. We have added the
-non-negative constant $\max\{E_{L-1}^x(1), 0\}$ to all the edges of the
-graph in \fixme{ref}. As $\lambda$ goes from $L-1$ to $0$ both edges
-stay non-negative. The edge from $u_\lambda^x$ to $t$ is non-decreasing
-with decreasing $\lambda$ parameter.
-
-This construction allows us to update the capacities of the edges, while
-retaining the flow and the trees $S$ and $T$. A nice property that
-follows is that the partition $(V - T, T)$ is a minimum cut, and the set
-of vertices $T$ is non-decreasing. More importantly, no vertex already
-in the tree $T$ will leave $T$ for succeeding lower $\lambda$ values.
-This can be seen from the fact that if $u$ is in $T$ and the algorithm
-has terminated, all paths from $s$ to $t$ going
-through $u$ has a saturated edge somewhere before $u$. For all these
-paths, the capacity will only change for edges after $u$, as we only
-change the capacity of edges $(v, t)$. Thus these saturated edges will
-stay saturated, and $u$ will stay in $T$.
-
-\subsection{Performance improvements}
-
-There are several open choices in the implementation of the algorithm,
-for example the order in which active vertices are processed. As
-recommended in \cite{boykov2004experimental}, we have implemented a
-«First-In-First-Out» queue. This ensures that at least the first
-augmenting path found is a shortest path, although later the distance
-information is lost in the adoption stage.
-
-In the adoption stage it is possible, and perhaps preferrable to seek
-possible parent that is closest to the root of the tree, and adopt that
-vertex as parent, instead of the first possible.
-
-Because of the particular graph construction, all vertices except $s$
-and $t$ are connected to $s$ and $t$. Thus there are as many two-edged
-paths from $s$ to $t$ as there are pixels in the image. When increasing
-the capacity of edges $(v, t)$, a quick sweep over these two-edged paths
-to send any possible flow may speed up the algorithm.
-
-\section{Implementation}
-
-A \cpp{} implementation can be found in
-Appendix~\ref{app:c++implementation}. It uses the open computer vision
-library OpenCV \cite{opencv_library} to load and save image files.
-
-Note that when implementing maximum flow algorithms it is not a good
-idea, memory- and performance-wise, to actually construct the residual
-graph $G_f$. Instead, every time we update the flow $f(u,v)$ we set
-the flow in the opposite direction to its negative value $f(v,u) =
--f(u,v)$. Then we can at any time, consider the value $c(u,v) - f(u,v)$
-in the place of the residual capacity $c_f(u,v)$.
-
-For the gap relabeling heuristic, we need to have a easy way of finding
-when a gap occurs. This is done by keeping track of how many vertices
-exist with each label.
-
-\fixme{something on bk?}
-
-\section{Performance comparison}
-
-Maybe.
-
projects / master / commitdiff