From: Bjørn Rustad Date: Thu, 29 Jan 2015 10:42:47 +0000 (+0100) Subject: Refactor files and stuff X-Git-Url: http://git.rustad.me/?a=commitdiff_plain;h=8bf3aa5e9cd67619df0bf229ce3c5cbdf9baa9d6;p=master Refactor files and stuff --- diff --git a/continuous.tex b/continuous.tex new file mode 100644 index 0000000..39fb95d --- /dev/null +++ b/continuous.tex @@ -0,0 +1,1268 @@ +\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. + diff --git a/discrete.tex b/discrete.tex new file mode 100644 index 0000000..ee35cb3 --- /dev/null +++ b/discrete.tex @@ -0,0 +1,722 @@ +\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 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. + diff --git a/results.tex b/results.tex index 10e5e93..c15b23e 100644 --- a/results.tex +++ b/results.tex @@ -513,10 +513,10 @@ approximated is that of an actual continuous circle $u : \mathbb{R}^2 \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} @@ -678,6 +678,6 @@ variation method it is based on, or other methods. In different 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. diff --git a/theory.tex b/theory.tex index 75b8ad1..e69de29 100644 --- a/theory.tex +++ b/theory.tex @@ -1,3033 +0,0 @@ -\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 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. -