\TVM(u) = \int_\Omega \sqrt{\nabla u(x)^T M(x) \nabla u(x)} \, dx
\end{equation}
\fixme{To differential or not to differential.}
-for all $u \in C^1(\Omega)$. We see that if $M(x)$ is the identity
+for all $u \in C_c^1(\Omega)$. We see that if $M(x)$ is the identity
matrix we get the normal total variation from \fixme{ref}, but our goal
is to control $M(x)$ such that $\nabla u(x)$ is weighted down across
edges, and thus edges will be regularized less.
+\fixme{Why compact support?}
+
+Before introducing the anisotropic coarea formula, we extend the
+definition of this functional to the space $\BV(\Omega)$. Being
+symmetric positive definite, the matrix $A(x)$ can be factored into two
+symmetric matrices as $A(x) = Q(x) Q(x)$. We can then write
+\begin{align}
+ \TVM(u) &= \int_\Omega \abs{Q \nabla u} \, dx \\
+ &= \sup_{\abs{\xi(x)}
+ \leq 1} \int_\Omega (Q \nabla u)^T \xi \, dx \\
+ &= \sup_{\abs{\xi(x)}\leq 1} \int_\Omega \nabla u \cdot Q\xi \, dx
+ \\
+ &= \sup_{\abs{\xi(x)}\leq 1} \int_\Omega -u \diver (Q\xi) \, dx \\
+ &= \sup_{\eta^T A^{-1} \eta \leq 1} \int_\Omega -u \diver \eta \,
+ dx.
+\end{align}
+If we let $\norm{\xi}_A = \sqrt{\xi^T A \xi}$ and $\norm{\eta}_A^* =
+\sqrt{\eta^T A^{-1} \eta}$, then our new extended definition reads
+\begin{equation}
+ \TVM(u) = \sup_{\substack{\xi \in C_c^1(\Omega; \mathbb{R}^N) \\
+ \norm{\xi(x)}_A^* \leq 1}} \int_\Omega -u \diver \xi \, dx.
+\end{equation}
+
\fixme{Energy function here?}
\fixme{Describe how this is similar to a Euler time step in a related
equation \cite{grasmair2010anisotropic}. This is only without the square
across edges. We do not want the tensor to be sensitive to noise in the
image. \fixme{Maybe refer to some figure}
+\fixme{A little more motivation.}
A good starting point is the structure tensor
\begin{equation}
S(x) := K_\rho * \left( \nabla u_\sigma(x) \otimes
- \nabla u_\sigma(x) \right).
+ \nabla u_\sigma(x) \right),
+ \label{eq:s_def}
\end{equation}
-This is then decomposed
+where the convolution is done componentwise. As this matrix is real and
+symmetric, it can be decomposed into
\begin{equation}
- J_\rho(x) = U(x)^T \Sigma(x) U(x)
+ S(x) = U(x) \Lambda(x) U(x)^T,
\end{equation}
-and say something about the number of dimensions. The numerical problems
-should be discussed somewhere but maybe not here. The eigenvalues are
-extracted such that
+where $\Lambda(x)$ is a matrix with the eigenvalues of $S(x)$ on the
+diagonal, while $U(x)$ has the eigenvectors of $S(x)$ as its columns.
+
+If we consider the structure tensor of \eqref{eq:s_def} before
+convolving with the Gaussian function, it is clear that there will be
+maximum one non-zero eigenvalue, whose eigenvector is the vector
+$\nabla u_\sigma(x)$ itself. If $\nabla u_\sigma(x) = (u_x, u_y)^T$ the
+eigenvalue equation becomes
\begin{equation}
- \Sigma(x) = \begin{bmatrix}
- \sigma_1(x) & 0 \\
- 0 & \sigma_2(x)
- \end{bmatrix}
+ \begin{pmatrix}
+ u_x^2 & u_x u_y \\
+ u_x u_y & u_y^2
+ \end{pmatrix}
+ \begin{pmatrix}
+ u_x \\
+ u_y
+ \end{pmatrix} =
+ (u_x^2 + u_y^2)
+ \begin{pmatrix}
+ u_x \\
+ u_y
+ \end{pmatrix}
\end{equation}
-where $\sigma_1 \gg \sigma_2$. This is then inserted into a new matrix
-$A = U^T \Lambda U$ where
+and unless the smoothing parameter $\rho$ is very large, the eigenvalues
+of $S(x)$ will always be dominated by the direction of $\nabla
+u_\sigma(x)$. This means that $S(x)$ is not at all suitable for being
+our metric tensor, since it actually weights the direction of $\nabla
+u_\sigma(x)$ more. We only need to change around the eigenvalues a bit.
+If the eigen-decomposition gave us
\begin{equation}
- \Lambda = \begin{bmatrix}
+ \Lambda(x) = \begin{pmatrix}
\lambda_1 & 0 \\
0 & \lambda_2
- \end{bmatrix}
+ \end{pmatrix}
\end{equation}
-and
+where $\lambda_1 \geq \lambda_2$, then we construct a new matrix $M(x) =
+U(x) \Sigma(x) U(x)^T$ where
+\begin{equation}
+ \Sigma(x) = \begin{pmatrix}
+ \sigma_1(x) & 0 \\
+ 0 & \sigma_2(x)
+ \end{pmatrix}.
+\end{equation}
+To reverse the effect of $S(x)$ such that $M(x)$ actually weighs down
+the norm across edges, we need $\sigma_1 \leq \sigma_2$, so we construct
\begin{align}
- \lambda_1 &= \frac{1}{1 + \frac{(\sigma_1 - \sigma_2)^2}{\gamma^2}},
- \\
- \lambda_2 &= 1.
+ \sigma_1 &= \left(1 + \frac{(\lambda_1 -
+ \lambda_2)^2}{\gamma^2}\right)^{-1}, \\
+ \sigma_2 &= 1.
\end{align}
-So the rotation is kept, while the size of the eigenvalues are changed.
+This way, $\sigma_1 \in (0, 1]$.
+The rotation is kept, while the size of the eigenvalues are changed.
This should be visualized, with a figure showing the length and
-direction of the eigenvalues in the area around an edge. But why is one
-of them always 1? Oh, it is not always the shortest, it varies?
+direction of the eigenvalues in the area around an edge.
+
+Note that in smooth parts of the image we have $\sigma_1 \approx
+\sigma_2 = 1$ and the anisotropic total variation is close to the
+regular total variation.
+
+\fixme{
+ The numerical problems should be discussed somewhere but maybe not
+ here. The eigenvalues are extracted such that
+}
\section{Analysis}
\fixme{Section heading!}
-In spirit of my project I have included a section with this title. It
-will be a theory heavy chapter, probably including the anisotropic
-coarea formula and the perimeters and all that. Should probably include
-some proofs and derivations as it is not easy to find this in the
-literature. This is one of the most important parts of the theory
-chapter, I guess.
-
-Existence and uniqueness should probably be mentioned here as well? But
-how much do we say about it? Grasmair is an OK reference, but he has the
-square root in there, what does that change? OK, talked to Markus about
-it, it is very technical and does not give much to do it properly, but
-should be discussed. Did he say that it could be done in some smaller
-space $C^1$ or $C^2$ and then we could leave the extension to someone
-else? For existence we need coercivity (that we can't go infinitely far
-to get a better solution) and weak lower semi-continuity, which means we
-consider the weak topology induced by the weak convergence (?). We do
-some kind of extension with $+\infty$ for functions outside our space,
-but that means we lose coercivity, which we have to fix again.
-Uniqueness is related to convexity. When the anisotropy tensor uses the
-smoothed initial value, this is trivial, but not so much if we use the
-output image in the anisotropy tensor. For us though, this will not be
-the case.
+\fixme{
+ Existence and uniqueness should probably be mentioned here as well?
+ But how much do we say about it? Grasmair is an OK reference. OK,
+ talked to Markus about it, it is very technical and does not give
+ much to do it properly, but should be discussed. Did he say that it
+ could be done in some smaller space $C^1$ or $C^2$ and then we could
+ leave the extension to someone else? For existence we need
+ coercivity (that we can't go infinitely far to get a better
+ solution) and weak lower semi-continuity (the lowest point always
+ exists at discontinuities), which means we consider the weak
+ topology induced by the weak convergence (?). We do some kind of
+ extension with $+\infty$ for functions outside our space, but that
+ means we lose coercivity, which we have to fix again. Uniqueness is
+ related to convexity. When the anisotropy tensor uses the smoothed
+ initial value, this is trivial (?), but not so much if we use the
+ output image in the anisotropy tensor. For us though, this will not
+ be the case.
+}
Earlier we had
\begin{equation}
u = \int_0^\infty u^s \, ds.
\end{equation}
-For any image $u \in C^1_c(\Omega)$ define the anisotropic total
-variation functional
-\begin{equation}
- J_A(u) = \int_\Omega \sqrt{\nabla u(x)^T A(x) \nabla u(x)} \, dx.
-\end{equation}
-\fixme{Why compact support?}
-Before introducing the anisotropic coarea formula, we extend the
-definition of this functional to the space $\BV(\Omega)$. Being
-symmetric positive definite, the matrix $A(x)$ can be factored into two
-symmetric matrices as $A(x) = Q(x) Q(x)$. We can then write
-\begin{align}
- J_A(u) &= \int_\Omega \abs{Q \nabla u} \, dx \\
- &= \sup_{\abs{\xi(x)}
- \leq 1} \int_\Omega (Q \nabla u)^T \xi \, dx \\
- &= \sup_{\abs{\xi(x)}\leq 1} \int_\Omega \nabla u \cdot Q\xi \, dx
- \\
- &= \sup_{\abs{\xi(x)}\leq 1} \int_\Omega -u \diver (Q\xi) \, dx \\
- &= \sup_{\eta^T A^{-1} \eta \leq 1} \int_\Omega -u \diver \eta \,
- dx.
-\end{align}
-If we let $\norm{\xi}_A = \sqrt{\xi^T A \xi}$ and $\norm{\eta}_A^* =
-\sqrt{\eta^T A^{-1} \eta}$, then our new extended definition reads
-\begin{equation}
- J_A(u) = \sup_{\substack{\xi \in C_c^1(\Omega; \mathbb{R}^N) \\
- \norm{\xi(x)}_A^* \leq 1}} \int_\Omega -u \diver \xi \, dx.
-\end{equation}
-
\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}
- J_A(u) = \int_{-\infty}^\infty J_A(u^s) \, ds.
+ \TVM(u) = \int_{-\infty}^\infty \TVM(u^s) \, ds.
\end{equation}
\end{theorem}
\begin{proof}
Assume that $u \in C^1(\Omega) \cap \BV(\Omega)$. \fixme{Something
about the approximation arguments.}
- \paragraph{First we prove that $J_A(u) \leq \int_{-\infty}^\infty
- J_A(u^s) \, ds$.}
+ \paragraph{First we prove that $\TVM(u) \leq \int_{-\infty}^\infty
+ \TVM(u^s) \, ds$.}
Assume that $u \geq 0$ such that the integral in \fixme{ref} holds,
then inserting \fixme{ref} into \fixme{ref} gives
\begin{align}
- J_A(u) &= \sup_{\norm{\xi}_A^* \leq 1} \int_\Omega -\left(
+ \TVM(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^*} \int_\Omega
-u^s \diver \xi \, dx \right) \, ds
- = \int_0^\infty J_A(u^s) \, ds.
+ = \int_0^\infty \TVM(u^s) \, ds.
\end{align}
- For $u \leq 0$ we use that $J_A(-v) = J_A(v)$ and that $J_A(c + v) =
- J_A(v)$ for any constant $c$. This allows us to show that
+ For $u \leq 0$ we use that $\TVM(-v) = \TVM(v)$ and that $\TVM(c + v) =
+ \TVM(v)$ for any constant $c$. This allows us to show that
\fixme{also something about $u^{-r}$}
\begin{align}
- J_A(u) &= J_A(-u) \leq \int_0^\infty J_A \big( (-u)^r \big) \,
+ \TVM(u) &= \TVM(-u) \leq \int_0^\infty \TVM \big( (-u)^r \big) \,
dr \\
- &= \int_0^\infty J_A(1 - u^{-r}) \, dr = \int_0^\infty
- J_A(u^{-r}) \, dr = \int_{-\infty}^0 J_A(u^s) \, ds.
+ &= \int_0^\infty \TVM(1 - u^{-r}) \, dr = \int_0^\infty
+ \TVM(u^{-r}) \, dr = \int_{-\infty}^0 \TVM(u^s) \, ds.
\end{align}
Next, we write a general $u$ as a difference between two positive
function $u = u_+ - u_-$ \fixme{They will not be differentiable
everywhere, but a.e.\ ?} and conclude that
\begin{align}
- J_A(u) &\leq J_A(u_-) + J_A(u_+) \\
- &\leq \int_{-\infty}^0 J_A\big((-u_-)^s\big) \, ds +
- \int_0^\infty J_A(u_+^s) \, ds\\
- &= \int_{-\infty}^0 J_A(u^s) \, ds + \int_0^\infty J_A(u^s) \,
- ds = \int_{-\infty}^\infty J_A(u^s) \, ds
+ \TVM(u) &\leq \TVM(u_-) + \TVM(u_+) \\
+ &\leq \int_{-\infty}^0 \TVM\big((-u_-)^s\big) \, ds +
+ \int_0^\infty \TVM(u_+^s) \, ds\\
+ &= \int_{-\infty}^0 \TVM(u^s) \, ds + \int_0^\infty \TVM(u^s) \,
+ ds = \int_{-\infty}^\infty \TVM(u^s) \, ds
\end{align}
- \paragraph{Then we prove that $J_A(u) \geq \int_{-\infty}^\infty
- J_A(u^s) \, ds$.}
+ \paragraph{Then we prove that $\TVM(u) \geq \int_{-\infty}^\infty
+ \TVM(u^s) \, ds$.}
Define the function
\begin{equation}
m(t) = \int_{\{ x \in \Omega : u(x) \leq t\}} \norm{\nabla u}_A
\, dx,
\end{equation}
- and note that $m(\infty) = J_A(u)$. Since $m(t)$ is non-decreasing
+ and note that $m(\infty) = \TVM(u)$. Since $m(t)$ is non-decreasing
with $t$, we can apply the existence theorem of Lebesgue \fixme{ref}
to conclude that $m'(t)$ exists a.e.\ (w.r.t.\ Lebesgue measure) and
that Lebesgue's inequality holds:
\begin{equation}
\int_{-\infty}^\infty m'(t)\, dt \leq m(\infty) - m(-\infty) =
- J_A(u).
+ \TVM(u).
\end{equation}
Next, fix an $s \in \mathbb{R}$ and define the function
\begin{equation}
\end{equation}
since $\eta_r(u) \to u^s$ when $r \to 0$. As this holds for any
$\norm{\xi}_A^* \leq 1$, we get from \fixme{ref def} that $m'(s)
- \geq J_A(u_s)$ almost everywhere and conclude
+ \geq \TVM(u_s)$ almost everywhere and conclude
\begin{equation}
- J_A(u) \geq \int_{-\infty}^\infty m'(t) \, dt \geq
- \int_{-\infty}^\infty J_A(u^s) \, ds.
+ \TVM(u) \geq \int_{-\infty}^\infty m'(t) \, dt \geq
+ \int_{-\infty}^\infty \TVM(u^s) \, ds.
\end{equation}
\fixme{something something set of measure zero, lebesgue integral
something something}
\end{proof}
This is where we magically introduce the perimeter, or more specifically
\begin{equation}
- \Per(U;\Omega) = J_A(\idfun_U).
+ \PerM(U;\Omega) = \TVM(\idfun_U).
\end{equation}
This is just a definition, the magic lies in the fact that we can
measure this perimeter as usual, just making sure to take the metric
tensor into account
\begin{equation}
- \Per(U;\Omega) = \int_{\partial^*U \cap \Omega} \sqrt{\nu_U(x)^T
+ \PerM(U;\Omega) = \int_{\partial^*U \cap \Omega} \sqrt{\nu_U(x)^T
A(x) \nu_U(x)} \, d \mathcal{H}^{n-1}.
\end{equation}
And then we say a little bit about this reduced boundary. It is a bit
where $c: [a,b] \to \Omega$ is a parametrized curve, such that $c'(t)$
is a tangent vector.
+We are now left with finding $u \in BV(\Omega)$ that minimizes
+\begin{equation}
+ E_v(u) =
+ \int_\Omega (u - v)^2 \, dx
+ + \int_{-\infty}^\infty \PerM( \{ u >
+ \lambda \}; \Omega) \, d\lambda.
+ \label{eq:continuous_energy}
+\end{equation}
+
\section{Cauchy--Crofton formulas}
If we parametrize straight lines as shown in Figure \fixme{ref}, we can
this set $d\mathcal{L} = \dpdr$ we are ready to introduce the
Cauchy--Crofton formula, which gives us a way to calculate the length of
a curve by looking at the measure of the set of lines that intersect the
-curve. Note that the measure $\dpdr$ is invariant under rigid motions,
-meaning combinations of translations and rotations.
+curve. Note that the measure $d\mathcal{L}$ is invariant under rigid
+motions, meaning combinations of translations and rotations.
\begin{theorem}[The Euclidean Cauchy--Crofton formula]
Given a 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} n_C(\phi, \rho) \, \dpdr = 2 \abs{C},
+ \int_\mathcal{L} n(L \cap C) \, d\mathcal{L}(L) = 2 \abs{C},
\end{equation}
- where $n_C(\phi, \rho)$ is the number of times the line $(\phi,
- \rho)$ intersects with the curve $C$.
+ where $n(L \cap C)$ is the number of times the line $L$ intersects
+ with the curve $C$.
\end{theorem}
This elegant formula is very useful when we later will discretize our
energy function. The set of lines $\mathcal{L}$ is then discretized
each point depends on a metric tensor $M(p)$ varying continuously over
our space, the Cauchy--Crofton formula becomes
\begin{equation}
- \abs{C}_R = \int_\mathcal{L} n_C(\phi, \rho) \, \frac{\det
- M(p)}{2 \left( u_L^T \cdot M(p) \cdot u_L \right)} \, \dpdr.
+ \abs{C}_R = \int_\mathcal{L} n(L \cap C) \, \frac{\det
+ M(p)}{2 \left( u_L^T \cdot M(p) \cdot u_L \right)} \,
+ d\mathcal{L}(L).
\end{equation}
- Here, $u_L$ is a unit vector along the line currently being
- integrated over.
+ Here, $u_L$ is a unit vector in the point $p$ along the line
+ currently being integrated over.
\end{theorem}
\begin{proof}
Assume first that our space is equipped with at constant metric
- tensor $A$. The length of our curve using this Riemannian metric can
+ tensor $M$. The length of our curve using this Riemannian metric can
then be calculated by transforming the curve back into euclidean
space and applying the Euclidean Cauchy--Crofton formula
\begin{align}
- \abs{C}_A &= \abs{A^{\sfrac{1}{2}}C} \\
- &= \int n(L \cap A^{\sfrac{1}{2}}C) \,
+ \abs{C}_M &= \abs{M^{\sfrac{1}{2}}C} \\
+ &= \int n(L \cap M^{\sfrac{1}{2}}C) \,
d\mathcal{L}(L) \\
- &= \int n(A^{-\sfrac{1}{2}}L \cap C) \,
+ &= \int n(M^{-\sfrac{1}{2}}L \cap C) \,
d\mathcal{L}(L) \\
- &= \int n(M \cap C) \abs{J(\Ahalf)}\,
- d\mathcal{L}(M)
+ &= \int n(P \cap C) \abs{J(\Mhalf)}\,
+ d\mathcal{L}(P)
\end{align}
- where $J(A^{\sfrac{1}{2}})$ is the Jacobian of our coordinate
- transformation $M = A^{-\sfrac{1}{2}} L$, which is what we need to
+ where $J(M^{\sfrac{1}{2}})$ is the Jacobian of our coordinate
+ transformation $P = M^{-\sfrac{1}{2}} L$, which is what we need to
find next.
- As $A\in \mathbb{R}^{2\times2}$ is symmetric, so is
- $A^{\sfrac{1}{2}}$, and it emits an eigendecomposition $\Ahalf
+ As $M\in \mathbb{R}^{2\times2}$ is symmetric, so is
+ $M^{\sfrac{1}{2}}$, and it emits an eigendecomposition $\Mhalf
= U\Lambda U^T$ where the components correspond to the following
coordinate transformations
\begin{align}
rotated by the $U^T$ operator, such that if $u$ is a unit vector
along the line being integrated over, then $v = U^Tu$ and
\begin{equation}
- \abs{J} = \frac{\det A}{\left(u^T A u\right)^{\sfrac{3}{2}}}
+ \abs{J} = \frac{\det M}{\left(u^T M u\right)^{\sfrac{3}{2}}}
\end{equation}
And now we argue that this also holds for a non-constant but
continuous metric tensor $M(x)$. By partitioning the domain into
\section{Discretization}
-It would probably be rewarding to look thorougly into the
-Cauchy--Crofton formula, especially in the Riemann metric case. We have
-to verify that all the approximations work, and that it indeed converges
-also in this case.
+We assume that our discrete images are given on a grid $\mathcal{G}$,
+where each discrete point is called a pixel and we further assume that
+each \emph{pixel} takes a value in the set of levels $\mathcal{L} = \{0,
+\hdots, L-1\}$. This is a reasonable assumption for grayscale images.
+(\fixme{repeated}). We now want to discretize the energy function in
+\eqref{eq:continuous_energy}. In addition, we want to decompose the
+energyfunction as a sum over all the levels $\mathcal{L}$.
-\subsection{Discrete Riemannian Cauchy--Crofton formula}
+\subsection{Fidelity term}
+
+For some pixel position $x \in \mathcal{G}$ and some level value $k \in
+\mathcal{L}$, we define the following function
+\begin{equation}
+ N_x(k) = \abs{k - v_x}^p
+\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 - v}^p\, dx \approx \sum_{x \in \mathcal{G}}
+ \abs{u_x - v_x}^p = \sum_{x \in \mathcal{G}} N_x(u_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{L}$:
+\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 - v_x}^p =
+ \sum_{x \in \mathcal{G}} N_x(u_x) = \sum_{\lambda = 0}^{L-2} \sum_x
+ \big( N_x(\lambda + 1) - N_x(\lambda) \big) \, u_x^\lambda + N_x(0).
+\end{equation}
+\fixme{Why stop at $L-2$}
+
+\subsection{Regularization term}
+
+Discretizing the regularization term is more challenging. We introduce
+the discrete levels to get
+\begin{equation}
+ \int_{-\infty}^\infty \PerM( \{ u > \lambda \}; \Omega) \, d\lambda
+ \approx \sum_{\lambda = 0}^{L-2} \PerM( \{ u > \lambda \}; \Omega).
+\end{equation}
+Note that we do not include any $\Delta \lambda$ difference, as we can
+just absorb this into the $\beta$ parameter. The perimeter is then
+calculated using a discretized version of the Cauchy--Crofton formula
+introduced in \fixme{ref}.
+\fixme{Why stop at $L-2$}
+
+\subsubsection{Discrete Riemannian Cauchy--Crofton formula}
\begin{figure}
\centering
where $\mathcal{L}_D$ is a discretization of the set of lines in the
plane $\mathcal{L}$ \fixme{figure}. Since the final goal is to work with
digital images, it makes sense to discretize our domain $\Omega$ as a
-regular (?) lattice? $\mathcal{S}$. The set of lines $\mathcal{L}_D$ can
+regular (?) lattice? $\mathcal{G}$. The set of lines $\mathcal{L}_D$ can
then be made up of lines through points on this lattice as shown in
Figure \fixme{some figure with lines}. Further working towards our graph
cut representation of the problem later, each line is made up of
to be a cut...}
\subsubsection{Consistency}
+
Consider a square centered around grid point $p$ with side lengths
$\sqrt{\delta}$. The size of this square will go to zero, while the
number of points along the outer edge $\lfloor 1 / \delta \rfloor$ of
obviously. \fixme{Some figure here aswell. And maybe some consideration
of $\sup \Delta \phi$ etc.}
+\subsubsection{Back to the energy...}
+
+\fixme{We talk about edges, but it's not really a graph yet...}
+The curve $C$ that we are measuring is the perimeter of a level set
+$X = \{ u > \lambda \}$. This means that the edges crossing the curve --
+$e \cap C$ -- are exactly those edges with one terminal inside $X$ and
+the other terminal outside $X$. \fixme{Terminal?}. This means we can
+write the discrete regularization energy term as
+\begin{equation}
+ \sum_{\lambda = 0}^{L-2} \PerM( \{ u > \lambda \}; \Omega)
+ \approx \sum_{\lambda = 0}^{L-2} \sum_{e} w_e
+ \abs{u_{\delta_1(e)}^\lambda - u_{\delta_2(e)}^\lambda}
+\end{equation}
+where the weight $w_e$ comes directly from the Riemannian
+Cauchy--Crofton formula
+\begin{equation}
+ w_e = \frac{\det M(e) \norm{e}^2 \delta^2 \Delta \phi}{2 (e^T \cdot
+ M(e) \cdot e)^{\sfrac{3}{2}}}.
+\end{equation}
+Now we just have to make sure that this fits well with the
+graph-representable function definition and all will be good.
+
+\subsection{Total energy}
+
+By combining the discretized fidelity and regularization terms and
+ignoring the constant term $N_x(0)$ we obtain an energy function which
+is decomposed into a sum over all the levels
+\begin{equation}
+ E_v(u) =
+ \sum_{\lambda=0}^{L-2} \sum_x E^x_\lambda(u^\lambda_x)
+ + \beta \sum_{\lambda = 0}^{L-2} \sum_{(x, y)} E^{x,y}(u^\lambda_x,
+ u^\lambda_y)
+ =: \sum_{\lambda=0}^{L-2} F_\lambda(u^\lambda)
+ \label{eq:total_energy}
+\end{equation}
+\fixme{oops label}
+where
+\begin{align}
+ E^x_\lambda(u^\lambda_x) &=
+ \big(
+ N_x(\lambda + 1) -
+ N_x(\lambda)
+ \big)
+ \, u^\lambda_x
+ \label{eq:fidelity_energy} \\
+ E^{x,y}(u^\lambda_x, u^\lambda_y) &=
+ w_{xy} \abs{u_x^\lambda - u_y^\lambda}.
+ \label{eq:neigh_energy}
+\end{align}
+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 (?) 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}
+In the following we will see two graph cut algorithms that find
+thresholded images minimizing each level, \emph{while guaranteeing that
+they meet this requirement.}
+
\section{Graph cut formulation}
Maybe this is more tightly connected with the previous section, but the
main point is that we describe how the edge weights are computed.
-\subsection{Networks}
-
-But first we have to say what a graph is. And also what a cut is!
+\subsection{Graphs}
\fixme{Maybe just call them graphs this time?}
-\subsection{Network representable energy functions}
+Using the notation of \cite{cormen2009introduction} we will denote a
+directed graph as $G = (V, E)$ where $V$ is a finite set of vertices, and
+$E$ is a binary relation on $V$. If $(u, v) \in E$ we say that there is
+an edge from $u$ to $v$ in the graph $G$.
+
+We introduce the non-negative capacity function $c : V \times V \to
+\left[0, \infty \right)$. Only edges $(u, v) \in E$ can have a positive
+capacity $c(u, v) = q > 0$ and it means that it is possible to send a
+\emph{flow} of maximum $q$ units from $u$ to $v$. For convenience we
+will let $c(u, v) = 0$ for any pair $(u, v) \not\in E$, and we do not
+allow self-loops in our graph. When a directed graph $G$ is equipped
+with capacity function $c$, one might call it a capacitated graph or a
+graph, but as all our graphs will be capacitated from this point, we
+will just call them graphs and we write $G = (V, E, c)$.
+
+There are two special vertices in the graph, the source $s$ and the
+sink $t$. Contrary to other vertices, which can neither produce nor
+receive excess flow, the source can produce and the sink can receive an
+unlimited amount of flow. The most basic problem in graph flow theory
+is the question of how much flow it is possible to send through the
+graph from the source to the sink.
+
+What we seek in our final graph is a minimum $s$-$t$-cut, a ``minimal''
+line through the graph that cuts a set of edges and
+divides the vertex set in two, separating the source from the sink.
+\begin{definition}[$s$-$t$-cut]
+ Given a graph $G = (V, E, c)$, an $s$-$t$-cut $(S, T)$ of $G$ is a
+ partition of $V$ into $S$ and $T = V - S$ such that $s \in S$ and $t
+ \in T$. The capacity of the cut is
+ \begin{equation}
+ c(S, T) = \sum_{u \in S} \sum_{v \in T} c(u, v),
+ \end{equation}
+ and a minimum $s$-$t$-cut is a cut whose capacity is minimum over
+ all $s$-$t$-cuts.
+ \label{def:s_t_cut}
+\end{definition}
+Note that there might exist several minimum $s$-$t$-cuts in a graph,
+resulting in different partitions of $V$. This is why we need to verify
+later that the cuts we obtain result in stackable thresholded images.
+
+\subsection{Graph representable energy functions}
+
+The next step is to find a way to construct a graph such that we can
+minimize the energy in \eqref{eq:total_energy} by finding a minimum
+$s$-$t$-cut. We will do this by creating small and simple graphs
+representing the separate summands of the energy. For these small
+graphs it will be easy to verify that the minimal cut also minimizes
+the corresponding part of the energy function, and they can then be
+merged giving a graph representing the complete energy function in
+\eqref{eq:total_energy}.
+
+First we need to establish the definition of a graph representable
+function, presented by Kolmogorov and Zabih in
+\cite{kolmogorov2002energy}.
+\begin{definition}[Graph representable functions]
+ A function $\mathcal{E}(x_1, \hdots, x_n)$ of $n$ binary
+ variables is graph-representable if
+ there exists a graph $G = (V, E, c)$ with terminals $s$ and $t$,
+ and a subset of vertices $V_0 = \{v_1, \hdots, v_n\} \subseteq V -
+ \{s, t\}$ such that, for any configuration $(x_1, \hdots, x_n) \in
+ \{0, 1\}^n$, the value of the energy $\mathcal{E}(x_1, \hdots, x_n)$
+ is equal to a constant plus the cost of the minimum $s$-$t$-cut
+ among all cuts $C = (S, T)$ where $x_i = 0 \Leftrightarrow v_i \in
+ S$ and $x_i = 1 \Leftrightarrow v_i \in T$, $\forall\, 1 \leq i \leq
+ n$.
+ \label{def:graph_representable}
+\end{definition}
+From this definition we see that if we have a graph-representable
+function $\mathcal{E}$ it is possible to find an exact global minimum
+of $\mathcal{E}$ by finding a minimal $s$-$t$-cut in a graph
+representing $\mathcal{E}$.
+
+Furthermore Kolmogorov and Zabih present an important result concerning
+what kinds of functions are graph-representable.
+\begin{theorem}[Identification of graph representable functions]
+ Given an energy function $\mathcal{E}$ of $n$ binary variables of
+ the form
+ \begin{equation}
+ \mathcal{E}(x_1, \hdots, x_n) = \sum_i \mathcal{E}^i(x_i)
+ + \sum_{i<j} \mathcal{E}^{i,j}(x_i, x_j),
+ \end{equation}
+ then $\mathcal{E}$ is graph representable if and only if each term
+ $\mathcal{E}^{i,j}$ satisfies the inequality
+ \begin{equation}
+ \mathcal{E}^{i,j}(0,0) + \mathcal{E}^{i,j}(1,1) \leq
+ \mathcal{E}^{i,j}(0,1) + \mathcal{E}^{i,j}(1,0).
+ \end{equation}
+ \label{thm:graph_rep_id}
+\end{theorem}
+This theorem will allow us to verify that our energy function actually
+is graph representable.
+
+Finally, the following theorem, proved by Kolmogorov and Zabih in
+\cite{kolmogorov2002energy}, will be crucial in our graph
+construction.
+\begin{theorem}[Additivity]
+ The sum of a finite number of graph-representable functions
+ \begin{equation}
+ \mathcal{E}(x_1, \hdots, x_n) = \sum_k \mathcal{E}^k(x_1,
+ \hdots, x_n),
+ \end{equation}
+ each represented by a graph $G^k = (V^k, E^k, c^k)$, is
+ graph-representable by $G = (V, E, c)$ where $V = \cup_k V^k$,
+ $E = \cup_k E^k$ and $c(u,v) = \sum_k c^k(u,v)$.
+\end{theorem}
+It allows us to construct small graphs representing the different
+summands of our energy function \eqref{eq:total_energy}, before adding
+them together to create a final graph representing the total energy.
+
+Note that when we apply this theorem later, we will assume that all the
+summands of \eqref{eq:total_energy} have the whole picture as their
+domain. It is unproblematic to extend $E_\lambda^x(u^\lambda_x)$ and
+$E^{x,y}(u^\lambda_x, u^\lambda_y)$ such that they take all the pixels
+their argument and then ignore all pixels except the ones they actually
+depend on.
+
+\subsection{Graph construction}
+
+\begin{figure}[t]
+ \input{fig/norm_subgraph}
+\end{figure}
+
+We will construct a graph in such a way that if a variable
+$u^\lambda_x$ ends up on the source side of the cut we set $u^\lambda_x
+= 0$, and if it ends up on the sink side we set $u^\lambda_x = 1$, as in
+Definition \ref{def:graph_representable}. This is an arbitrary choice,
+but still something we have to keep in mind through the rest of the
+section.
+
+We will now consider the two kinds of summands in the energy function in
+\eqref{eq:total_energy}. The fidelity term coming from our aim to keep
+the output image close to the original image, and the regularization
+term coming from our aim to minimize the total variation.
+
+\begin{table}[b]
+ \centering
+ \caption{Each row represents one of the two possible values of
+ $u^\lambda_x \in \{0,1\}$. The energy $E_\lambda^x(u^\lambda_x)$
+ and minimum cut obtaining this configuration is shown. The last
+ two columns show the capacities of the cut for each of the two
+ graph constructions in Figure \ref{fig:norm_subgraph}. We
+ verify that for each of the two graph constructions, the cut
+ capacities are equal to the energies, plus a constant.
+ }
+ \begin{tabular}{ccccc}
+ \hline
+ $u^\lambda_x$ & $E_\lambda^x(u^\lambda_x)$ & Min.\ cut $(S, T)$
+ & graph \subref{fig:norm_subgraph_pos} cut cap. & Graph
+ \subref{fig:norm_subgraph_neg} cut cap. \\
+ \hline
+ $0$ & $E_\lambda^x(0)$ & $(\{s, u^\lambda_x\}, \{t\})$ &
+ $E_\lambda^x(0)$ & $0$ \\
+ $1$ & $0$ & $(\{s\}, \{u^\lambda_x, t\})$ & $0$ & $-E_\lambda^x(0)$
+ \end{tabular}
+ \label{tab:fid_energy}
+\end{table}
+
+\subsubsection{Fidelity term}
+
+\begin{figure}[t]
+ \input{fig/neigh_subgraph}
+\end{figure}
+
+The fidelity term of our energy function in \eqref{eq:total_energy}
+simplifies to
+\begin{align}
+ E_\lambda^x(0) &=
+ N_x(\lambda + 1) -
+ N_x(\lambda)
+ \label{eq:fidelity_energy0} \\
+ E_\lambda^x(1) &= 0
+ \label{eq:fidelity_energy1}
+\end{align}
+where $E_\lambda^x(0)$ might be positive or negative depending on
+$\lambda$ and the pixel value $v_x$.
+
+Figure \ref{fig:norm_subgraph} shows how graphs can be constructed to
+represent this part of the total energy. The construction differs
+depending on whether $E_\lambda^x(0)$ is positive or negative. Table
+\ref{tab:fid_energy} shows how the cuts correspond to the values of
+$u^\lambda_x$ and we can easily verify that the constructed graph
+actually represents the fidelity term in the energy function.
+
+\subsubsection{Regularization term}
-Then, how we can represent an energy function as a graph, and how
-finding the minimum cut will give a minimum energy value.
+For our neighboring relation in \eqref{eq:total_energy} of the form
+\begin{equation}
+ E^{x,y}(u^\lambda_x, u^\lambda_y) =
+ w_{xy} \abs{u^\lambda_x - u^\lambda_y}
+\end{equation}
+we have
+\begin{equation}
+\begin{aligned}
+ E^{x,y}(0, 0) &= 0, \\
+ E^{x,y}(0, 1) &= w_{xy}, \\
+ E^{x,y}(1, 0) &= w_{xy}, \\
+ E^{x,y}(1, 1) &= 0,
+\end{aligned}
+\label{eq:neigh_energies}
+\end{equation}
+and by Theorem \ref{thm:graph_rep_id} our energy function is graph
+representable. In \cite{kolmogorov2002energy}, Kolmogorov and Zabih
+presents a way to construct a graph for any graph representable function
+on the form shown in Theorem \ref{thm:graph_rep_id}. Since the
+energies in \eqref{eq:neigh_energies} are especially simple, the
+construction and presentation is simplified.
+
+
+Figure \ref{fig:neigh_subgraph} shows two different ways of how a graph
+can be constructed to represent the regularization term. See Table
+\ref{tab:neigh_energy} for an overview of how the two values of
+$u^\lambda_x$ corresponds to cuts in the graph.
+
+\begin{table}[t]
+ \centering
+ \caption{An overview of the possible configurations of the variables
+ in the term $E^{x,y}(u^\lambda_x, u^\lambda_y)$. For each
+ configuration the corresponding energy and the cut yielding this
+ configuration is shown. The last two columns show the capacities
+ of the cut in the two alternative graph constructions shown in
+ Figure \ref{fig:neigh_subgraph}. We verify that for each of the
+ two graph constructions, the cut capacities are equal to the
+ energies, plus a constant.
+ }
+ \begin{tabular}{ccccc}
+ \hline
+ $(u^\lambda_x, u^\lambda_y)$
+ & $E^{x,y}(u^\lambda_x, u^\lambda_y)$
+ & Min.\ cut $(S, T)$
+ & Alt.\ \subref{fig:neigh_subgraph_alt1} cut cap. & Alt.\
+ \subref{fig:neigh_subgraph_alt2} cut cap. \\
+ \hline
+ $(0, 0)$ & $0$ & $(\{s, u^\lambda_x, u^\lambda_y\}, \{t\})$ & $w_{xy}$ & $0$ \\
+ $(0, 1)$ & $w_{xy}$ & $(\{s, u^\lambda_x\}, \{u^\lambda_y, t\})$ & $2w_{xy}$ & $w_{xy}$ \\
+ $(1, 0)$ & $w_{xy}$ & $(\{s, u^\lambda_y\}, \{u^\lambda_x, t\})$ & $2w_{xy}$ & $w_{xy}$ \\
+ $(1, 1)$ & $0$ & $(\{s\}, \{u^\lambda_x, u^\lambda_y, t\})$ & $w_{xy}$ & $0$
+ \end{tabular}
+ \label{tab:neigh_energy}
+\end{table}
-\subsection{Network construction}
+Figure \ref{fig:big_graph} shows a visualization of how the final
+graph might look with all its edges. The source will have a lot of
+outgoing edges, one for each pixel, while the sink has one incoming edge
+for each pixel. The vertices corresponding to the pixels are only
+connected to the source, the sink, and their neighboring pixels, so
+their edge degree is much smaller than for $s$ and $t$.
-Then we describe how we actually create the graph.
+\begin{figure}[bh]
+ \input{fig/big_graph}
+\end{figure}
\chapter{Maximum flow approach}
-Next up is the maximum flow. Here we can get a lot from the project.
+In the previous section we have seen how finding the minimum cut of
+carefully constructed graph can give us the thresholded image
+minimizing the energy function for one level value $\lambda$. We will in
+this and the next section see how such a minimum cut can be found by
+sending flow through the graph and trying to identify the
+``bottleneck''.
+
+\section{Flow graphs}
+We have already introduced capacities, and briefly mentioned the notion
+of flow as something limited by the capacity. In other words, flow is
+something we can send through our graph, but the capacity limits how
+much we can send along each edge. It is useful to imagine a water supply
+graph with pipes of different sizes and a water source and sink.
+
+Formally, we introduce the \emph{flow} as the function $f : V \times V
+\to [0, \infty)$. This function keeps count of how much flow we are
+sending through each edge of our graph and must satisfy the following
+two constraints
+\begin{description}
+ \item[Capacity constraint:] For all $u, v \in V, 0 \leq f(u, v) \leq c(u,
+ v)$, i.e., for every pair of vertices, the flow is less than or
+ equal to the capacity.
+ \item[Flow conservation:] For all $u \in V - \{s, t\}$
+ \begin{equation}
+ \sum_{v \in V} f(v, u) = \sum_{v \in V} f(u, v),
+ \end{equation}
+ i.e., for any vertex except the source and the sink, the
+ flow into the vertex must be equal to the flow out of the
+ vertex.
+\end{description}
+Note that we have defined $f$ with all pairs of vertices as its domain,
+even though it is only non-zero on edges $(u, v) \in E$. This makes it
+easier to write sums as in the flow conservation constraint. We say that
+an edge $(u, v)$ is \emph{saturated} if $f(u, v) = c(u, v)$.
+
+We define $\abs{f}$ as the net amount of flow from the source to the
+sink in the graph. Because of the flow conservation constraint, this
+can be calculated as the net amount of flow going out of the source
+\begin{equation}
+ \abs{f} = \sum_{v \in V} f(s, v) - \sum_{u \in V} f(u, s).
+ \label{eq:total_flow}
+\end{equation}
+Furthermore we denote the \emph{net flow} across an $s$-$t$-cut $C = (S,
+T)$ as
+\begin{equation}
+ f(S, T) = \sum_{u \in S} \sum_{v \in T} f(u, v) - \sum_{v \in T}
+ \sum_{u \in S} f(v, u).
+ \label{eq:net_flow_cut}
+\end{equation}
+Note how this definition differs from the capacity of a cut $c(S,T)$ in
+Definition \ref{def:s_t_cut}. While the capacity of a cut represents how
+much flow it is maximally possible to send from $S$ to $T$, the net flow
+across a cut represents the net amount of flow going across the cut,
+counting negatively the flow that goes back from $T$ to $S$.
+
+For any $s$-$t$-cut we have that
+\begin{equation}
+ \abs{f} = f(S, T).
+ \label{eq:total_flow_eq_cut}
+\end{equation}
+This is quite intuitive given the flow conservation constraint, and a
+full chain of arguments can be found in \cite{cormen2009introduction}.
+
+Following from \eqref{eq:net_flow_cut} and \eqref{eq:total_flow_eq_cut},
+we find that
+\begin{equation}
+ \begin{aligned}
+ \abs{f} &= \sum_{u \in S} \sum_{v \in T} f(u, v)
+ - \sum_{v \in T} \sum_{u \in S} f(v, u) \\
+ &\leq \sum_{u \in S} \sum_{v \in T} f(u, v) \\
+ &\leq \sum_{u \in S} \sum_{v \in T} c(u, v) \\
+ &\leq c(S, T)
+ \label{eq:flow_leq_cut}
+ \end{aligned}
+\end{equation}
+for any $s$-$t$-cut $C = (S, T)$. A very central result in graph flow
+theory called the max-flow min-cut theorem will be presented later. It
+states that the inequality of \eqref{eq:flow_leq_cut} becomes an equality
+when $f$ is a maximum flow for some cut $C = (S, T)$, and that all such
+cuts are minimum cuts.
+
+%Recall that we want to find the minimum $s$-$t$ cut in our graph. When
+%finding this minimum cut, we make use of an important duality theorem in
+%graph flow theory, stating that the capacity of a minimum $s$-$t$-cut
+%in a graph, is equal to the maximum flow from the source to the sink.
-\section{Flow networks}
+But how does this help us? We know that if we know the maximum flow
+value, and we have an $s$-$t$-cut with capacity equal to the maximum
+flow, we actually have a minimum cut. The question is then, how do we
+find a maximum flow, and how do we find a minimum cut?
-What is flow. Max--flow--min--cut theorem.
+\section{Augmenting path algorithms}
+The family of augmenting flow algorithms represent a popular approach to
+the maximum flow problem. The idea is simply to look for paths from the
+source to the sink that can carry additional flow, so-called augmenting
+paths, and then send the maximum possible amount of flow along such a
+path. When no such path exists anymore, no more flow can be sent from
+the source to the sink, and a maximum flow has been reached.
-\section{Augmenting flow algorithms}
+\subsection{Residual graph}
+When further discussing approaches to solving the maximum flow problem
+we will need the notion of a residual graph $G_f = (V_f, E_f, c_f)$,
+which is derived from the original graph $G$ and contains the edges
+along which it is possible to send additional flow. This means that
+$E_f$ contains the edges $(u, v)$ from $E$ where $f(u,v) < c(u,v)$. But
+that is not all; an important realisation is that it is also possible to
+push flow \emph{back} along an edge where the flow is already positive.
+In other words, sending flow from $v$ to $u$ by cancelling some or all
+of the flow that is already going from $u$ to $v$.
-\subsection{Ford--Fulkerson}
+Thus the capacity function $c_f$ of our residual graph becomes
+\begin{equation}
+ c_f(u, v) =
+ \begin{cases}
+ c(u, v) - f(u, v) & \text{if } (u, v) \in E, \\
+ f(v, u) & \text{if } (v, u) \in E, \\
+ 0 & \text{otherwise.}
+ \end{cases}
+\end{equation}
+The vertices $V_f$ of $G_f$ are the same as the original graph $G$,
+while the edges $E_f$ are taken to be all pairs of vertices $(u, v)$
+with $c_f(u, v) > 0$. Since we can in $E_f$ at most have all the
+original edges, and their reversals, we have $\abs{E_f} \leq 2 \abs{E}$.
+
+Note that there is ambiguity in the definition of the residual graph
+in the case where the original graph contains anti-parallel edges. One
+could avoid this by defining $c_f(u, v) = f(v, u) + c(u, v) - f(u, v)$
+instead, or as they do in \cite{cormen2009introduction}, disallow
+anti-parallel edges altogether. In any case it is not something we will
+have to think about in the implementation, since we will not actually
+construct the residual graph.
-\subsection{Edmonds--Karp}
+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 keeps two non-overlapping trees
+of non-saturated edges, one with the source and the other with the sink
+as its root. These trees are ``grown'' towards eachother and augmenting
+paths are found when their leaf nodes touch. In theory, the complexity
+of this algorithm is not great, but according to Boykov and Kolmogorov,
+it outperforms the other algorithms in experimental comparisons for this
+specific application.
\section{Push--relabel algorithm}
-\subsection{Network reuse}
+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(0)$
+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.
+
+\subsection{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.
\section{Boykov--Kolmogorov algorithm}
projects / master / commitdiff