From 02ce8fa0f07d85d5ebcd62a7266a2dc037b13a9a Mon Sep 17 00:00:00 2001 From: =?utf8?q?Bj=C3=B8rn=20Rustad?= Date: Wed, 22 Oct 2014 19:18:57 +0200 Subject: [PATCH] Copy stuff from project --- .gitignore | 1 + commands.tex | 1 + fig/aug_flow.tex | 22 + fig/big_graph.tex | 41 ++ fig/neigh_subgraph.tex | 44 ++ fig/norm_evolution.tex | 54 ++ fig/norm_subgraph.tex | 36 + main.tex | 20 + theory.tex | 1427 ++++++++++++++++++++++++++++++++++++---- 9 files changed, 1510 insertions(+), 136 deletions(-) create mode 100644 fig/aug_flow.tex create mode 100644 fig/big_graph.tex create mode 100644 fig/neigh_subgraph.tex create mode 100644 fig/norm_evolution.tex create mode 100644 fig/norm_subgraph.tex diff --git a/.gitignore b/.gitignore index 080e206..859c7b3 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,4 @@ *.pdf *.toc *.swp +*.xdv diff --git a/commands.tex b/commands.tex index 0d9b8e4..6110fdb 100644 --- a/commands.tex +++ b/commands.tex @@ -9,6 +9,7 @@ \newcommand{\idfun}{\chi} \newcommand{\dpdr}{d\phi \, d\rho} \newcommand{\Ahalf}{A^{\sfrac{1}{2}}} +\newcommand{\Mhalf}{M^{\sfrac{1}{2}}} \definecolor{fixme-color}{rgb}{0.7,0.2,0} \newcommand{\fixme}[1]{\textcolor{fixme-color}{FIXME: #1}} \def\cpp{C{}\texttt{++}} diff --git a/fig/aug_flow.tex b/fig/aug_flow.tex new file mode 100644 index 0000000..458818a --- /dev/null +++ b/fig/aug_flow.tex @@ -0,0 +1,22 @@ +\centering +\begin{tikzpicture}[scale=1.5] + \node[vertex] (s) at (0, 0) {s}; + \node[vertex] (u) at (2, 1) {u}; + \node[vertex] (v) at (2, -1) {v}; + \node[vertex] (t) at (4, 0) {t}; + + \path[edge] (s) -- node[weight] {5/5} (u); + \path[selected edge] (s) -- node[weight] {} (v); + \path[edge] (s) -- node[weight] {0/5} (v); + \path[selected edge] (u) -- node[weight] {} (v); + \path[edge] (u) -- node[weight] {5/5} (v); + \path[selected edge] (u) -- node[weight] {} (t); + \path[edge] (u) -- node[weight] {0/5} (t); + \path[edge] (v) -- node[weight] {5/5} (t); +\end{tikzpicture} +\caption{A graph with the flow and capacity of each edge shown as + flow/capacity. The marked path is a valid augmenting path from + $s$ to $t$, and by increasing the flow along it, we will push + flow back from $v$ to $u$. +} +\label{fig:aug_flow} diff --git a/fig/big_graph.tex b/fig/big_graph.tex new file mode 100644 index 0000000..d905378 --- /dev/null +++ b/fig/big_graph.tex @@ -0,0 +1,41 @@ +\centering +\begin{tikzpicture}[scale=1.3] + + \node[vertex] (s) at (.6, -1) {s}; + \node[vertex] (t) at (.6, -5) {t}; + + \begin{scope}[ + yshift=-100,every node/.append style={ + yscale=.3,yslant=.3,xslant=-0.6},yscale=.3,yslant=.3,xslant=-0.6 + ] + \foreach \x in {0,...,3} { + \foreach \y in {0,...,3} { + \node[fill=white,vertex,draw=black] (\x\y) at (\x, \y) {}; + } + } + \foreach \x in {0,...,3} { + \foreach \y [count=\yi] in {0,...,2} { + \path[double edge] (\x\y)--(\x\yi); + \path[double edge] (\y\x)--(\yi\x); + } + } + \end{scope} + + + \foreach \x in {0,...,3} { + \foreach \y in {0,...,3} { + \path[edge,out=-90,in=90] (s) to + (\x\y.center); + \begin{pgfonlayer}{background} + \path[edge,out=-90,in=90] (\x\y) to (t); + \end{pgfonlayer} + } + } +\end{tikzpicture} +\caption{ + A visualization of the final graph, with edges between neighboring + pixels, and edges connecting the source and sink to the rest of the + graph. The edge capacities are not visualized, and of course some + of them might be zero while others might be very high. +} +\label{fig:big_graph} diff --git a/fig/neigh_subgraph.tex b/fig/neigh_subgraph.tex new file mode 100644 index 0000000..54935b4 --- /dev/null +++ b/fig/neigh_subgraph.tex @@ -0,0 +1,44 @@ +\centering +\begin{subfigure}[t]{0.4\textwidth} + \centering + \begin{tikzpicture}[scale=1.5] + \node[vertex] (s) at (0, 0) {s}; + \node[vertex] (u) at (-1,-2) {$u^\lambda_x$}; + \node[vertex] (v) at (1, -2) {$u^\lambda_y$}; + \node[vertex] (t) at (0, -4) {t}; + + \path[edge] (s) -- node[weight] {$w_{xy}$} (u); + \path[edge] (u) -- node[weight] {$2w_{xy}$} (v); + \path[edge] (v) -- node[weight] {$w_{xy}$} (t); + \end{tikzpicture} + \caption{Representing $E^{x,y}(u^\lambda_x, u^\lambda_y)$ with a + constant term of $w_{xy}$.} + \label{fig:neigh_subgraph_alt1} +\end{subfigure} +~ +\begin{subfigure}[t]{0.4\textwidth} + \centering + \begin{tikzpicture}[scale=1.5] + \node[vertex] (s) at (0, 0) {s}; + \node[vertex] (u) at (-1,-2) {$u^\lambda_x$}; + \node[vertex] (v) at (1, -2) {$u^\lambda_y$}; + \node[vertex] (t) at (0, -4) {t}; + + \path[edge] ([yshift=1.5pt]u.east) -- node[weight around] + {$w_{xy}$} + ([yshift=1.5pt]v.west); + \path[edge] ([yshift=-1.5pt]v.west) -- node[weight around] + {$w_{xy}$} + ([yshift=-1.5pt]u.east); + \end{tikzpicture} + \caption{Representing $E^{x,y}(u^\lambda_x, u^\lambda_y)$ with a + constant term of 0.} + \label{fig:neigh_subgraph_alt2} +\end{subfigure} +\caption{Two alternative ways of constructing a graph representing + the energy term $E^{x,y}(u^\lambda_x, u^\lambda_y)$. See Table + \ref{tab:neigh_energy} for an overview of the different possible + configurations of $(u^\lambda_x, u^\lambda_y)$ and how they + correspond to cuts through the graph. +} +\label{fig:neigh_subgraph} diff --git a/fig/norm_evolution.tex b/fig/norm_evolution.tex new file mode 100644 index 0000000..ca85d77 --- /dev/null +++ b/fig/norm_evolution.tex @@ -0,0 +1,54 @@ +\centering +\begin{subfigure}[t]{0.4\textwidth} + \centering + \begin{tikzpicture} + \draw[->] (0,0) -- (4,0) node[right] {$\lambda$}; + \draw[->] (0,-2) -- (0,2) node[above] {$E_\lambda^x(0)$}; + \draw[line,domain=0:1] plot ({\x},{-1}); + \draw[line,domain=1:2,dashed] plot ({\x},{2*\x-3}); + \draw[line,domain=2:4] plot ({\x},{1}); + + \draw (1,2pt) -- (1,-2pt) node[ + font=\small, + anchor=south + ] + {$v_x-1$}; + + \draw (2,2pt) -- (2,-2pt) node[ + font=\small, + anchor=north + ] + {$v_x$}; + + \end{tikzpicture} + \caption{$L^1$ fidelity term.} + \label{fig:l1_norm_evolution} +\end{subfigure} +~ +\begin{subfigure}[t]{0.4\textwidth} + \centering + \begin{tikzpicture} + \draw[->] (0,0) -- (4,0) node[right] {$\lambda$}; + \draw[->] (0,-2) -- (0,2) node[above] {$E_\lambda^x(0)$}; + \draw[line,domain=0:4] plot ({\x},{.8*\x - 1.6}); + + \draw (2.5,2pt) -- (2.5,-2pt) node[ + font=\small, + anchor=north + ] + {$v_x$}; + + \draw (1.5,2pt) -- (1.5,-2pt) node[ + font=\small, + anchor=south + ] + {$v_x-1$}; + + \end{tikzpicture} + \caption{$L^2$ fidelity term.} + \label{fig:l2_norm_evolution} +\end{subfigure} +\caption{Two figures showing how the fidelity energy term $E_\lambda^x(0)$ + in \eqref{eq:total_energy} increases monotonically with $\lambda$. +} +\label{fig:norm_evolution} diff --git a/fig/norm_subgraph.tex b/fig/norm_subgraph.tex new file mode 100644 index 0000000..3e95ca1 --- /dev/null +++ b/fig/norm_subgraph.tex @@ -0,0 +1,36 @@ +\centering +\begin{subfigure}[t]{0.45\textwidth} + \centering + \begin{tikzpicture}[scale=2.0] + \node[vertex] (s) at (0, 0) {s}; + \node[vertex] (u) at (0, -1) {$u^\lambda_x$}; + \node[vertex] (t) at (0, -2) {t}; + + \path[edge] (u) -- node[weight noslope] {$E_\lambda^x(0)$} (t); + \end{tikzpicture} + \caption{The graph when $E_\lambda^x(0) > 0$, with constant equal to + 0. + } + \label{fig:norm_subgraph_pos} +\end{subfigure} +~ +\begin{subfigure}[t]{0.45\textwidth} + \centering + \begin{tikzpicture}[scale=2.0] + \node[vertex] (s) at (0, 0) {s}; + \node[vertex] (u) at (0, -1) {$u^\lambda_x$}; + \node[vertex] (t) at (0, -2) {t}; + + \path[edge] (s) -- node[weight noslope] {$-E_\lambda^x(0)$} (u); + \end{tikzpicture} + \caption{Graph when $E_\lambda^x(0) < 0$, with constant equal to + $-E_\lambda^x(0)$. + } + \label{fig:norm_subgraph_neg} +\end{subfigure} +\caption{ + The graph construction for the fidelity term $E_\lambda^x(u^\lambda_x)$. + See Table \ref{tab:fid_energy} for an overview of the different + possible cuts, and on why this construction works. +} +\label{fig:norm_subgraph} diff --git a/main.tex b/main.tex index aec555a..21139dd 100644 --- a/main.tex +++ b/main.tex @@ -22,6 +22,26 @@ \usetikzlibrary{positioning} \usetikzlibrary{calc} +\tikzstyle{vertex} = [circle,minimum size=20pt,inner sep=0pt,draw] +\tikzstyle{small vertex} = [ + font=\footnotesize,circle, + minimum size=9pt,inner sep=0pt,fill=gray!40 +] +\tikzstyle{selected vertex} = [vertex, fill=red!24] +\tikzstyle{edge} = [draw,thick,->] +\tikzstyle{dash} = [draw,dashed,-] +\tikzstyle{bent edge} = [draw,bend left=80,thick,->] +\tikzstyle{double edge} = [draw,thick,<->] +\tikzstyle{weight} = [font=\small,sloped,fill=white] +\tikzstyle{weight noslope} = [font=\small,fill=white] +\tikzstyle{weight around} = [font=\small,sloped,auto] +\tikzstyle{selected edge} = [draw,line width=5pt,-,red!50] +\tikzstyle{ignored edge} = [draw,line width=5pt,-,black!20] +\tikzstyle{line} = [ultra thick,scale=1,smooth,variable=\x,blue] +\pgfdeclarelayer{background} +\pgfdeclarelayer{foreground} +\pgfsetlayers{background,main,foreground} + \usepackage{color} %\usepackage[margin=0.8cm]{caption} \usepackage{caption} diff --git a/theory.tex b/theory.tex index aa26c15..fbd2aa7 100644 --- a/theory.tex +++ b/theory.tex @@ -16,10 +16,33 @@ introduce the anisotropic total variation \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 @@ -51,69 +74,104 @@ Our wish for the metric tensor is that it down-weighs the variation 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} @@ -146,89 +204,62 @@ over all the layers 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} @@ -262,23 +293,23 @@ If we let $\norm{\xi}_A = \sqrt{\xi^T A \xi}$ and $\norm{\eta}_A^* = \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 @@ -291,6 +322,15 @@ the tangent vector 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 @@ -299,16 +339,16 @@ define the set of all staight lines as $\mathcal{L} = \{ (\phi, \rho) : 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 @@ -320,32 +360,33 @@ can be approximated by a sum. 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} @@ -420,7 +461,7 @@ can be approximated by a sum. 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 @@ -491,12 +532,64 @@ tested later). Then the stencil shape can be adjusted as well. \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 @@ -525,7 +618,7 @@ the approximation 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 @@ -575,6 +668,7 @@ length is approximated by 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 @@ -589,43 +683,1104 @@ cross $e$ is maximum $1$ so that's great. Also, $M(e) \to M(p)$ 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 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} -- 2.47.3