publisher={Springer}\r
}\r
\r
+@article{rudin1992nonlinear,\r
+ title={Nonlinear total variation based noise removal algorithms},\r
+ author={Rudin, Leonid I and Osher, Stanley and Fatemi, Emad},\r
+ journal={Physica D: Nonlinear Phenomena},\r
+ volume={60},\r
+ number={1},\r
+ pages={259--268},\r
+ year={1992},\r
+ publisher={Elsevier}\r
+}\r
+\r
+@article{darbon2006image,\r
+ title={Image restoration with discrete constrained total variation part I: Fast and exact optimization},\r
+ author={Darbon, J{\'e}r{\^o}me and Sigelle, Marc},\r
+ journal={Journal of Mathematical Imaging and Vision},\r
+ volume={26},\r
+ number={3},\r
+ pages={261--276},\r
+ year={2006},\r
+ publisher={Springer}\r
+}\r
+\r
+@inproceedings{boykov2003computing,\r
+ title={Computing geodesics and minimal surfaces via graph cuts},\r
+ author={Boykov, Yuri and Kolmogorov, Vladimir},\r
+ booktitle={Computer Vision, 2003. Proceedings. Ninth IEEE International Conference on},\r
+ pages={26--33},\r
+ year={2003},\r
+ organization={IEEE}\r
+}\r
+\r
</bibtex>\r
--- /dev/null
+\newcommand{\abs}[1]{\lvert #1 \rvert}
+\newcommand{\norm}[1]{\lVert #1 \rVert}
+\newcommand{\TV}{\mathit{TV}}
+\newcommand{\BV}{\mathit{BV}}
+\newcommand{\idfun}{\mathbbm{1}}
+\definecolor{fixme-color}{rgb}{0.7,0.2,0}
+\newcommand{\fixme}[1]{\textcolor{fixme-color}{FIXME: #1}}
--- /dev/null
+\section{Graph cut formulation}
+Short introduction to graphs. Notation. What is a cut, what is a minimal
+s-t-cut. How do we get from the discrete energy to the graph and why
+does finding the graph cut result in minimizing the energy of the
+current label. Boykov and Kolmogorov.
+
+\subsection{Min-cut max-flow theorem}
+Want to find the minimal s-t-cut, so we introduce flow here maybe? Or
+earlier? Finding the max flow results in also finding the minimal cut,
+and there are «good» algorithms for finding the maximal flow.
+
+\subsection{Maximum flow algorithms}
+Two different ideas. One is to always maintain a valid flow in the
+network, Ford Foulkerson, always looking for new paths from s to t along
+which we can send additional flow. This leads to Edmonds-Karp, and
+Dinic/Dinitz. Second idea is push-relabel, allowing positive excess in
+the nodes. We lift the source, and let as much flow as possible flow
+towards the sink.
+
+\subsubsection{Dinic's algorithm}
+The «best» of the augmenting-path algorithms (citation needed). Makes a
+level graph with BFS, and finds a blocking flow through this graph which
+only goes from nodes with one label to nodes with greater labels. When
+this blocking flow is found we restart at the BFS.
+
+\subsubsection{The push-relabel algorithm}
+The nodes are allowed to have a positive excess but we still follow the
+capacity constraints. The nodes also have a labeling, which has to be
+valid. Push flow from active nodes, and relabel, preferrably in a
+specific order, until it is not possible anymore. The minimal cut is
+actually found before the flow becomes valid (has to be explained).
+
+Because of the way the graph evolves, we can reuse the graph between the
+different labels.
+
+Initially proposed in \cite{goldberg1988new} where one can also find
+proof that the algorithm maintains a valid labeling and that it
+terminates at at maximal flow. Also has the lemma that if you have a
+preflow, then there is no path from $s$ to $t$ in the residual graph.
+But what is a preflow? Also, the distance label never decreases. Also
+propose the global relabeling heuristic.
+
+Shiloach and Vishkin originally proposed the FIFO way of doing things
+(according to Golberg and Tarjan)?
+
+Cherkassky and Goldberg present some good implementation tips in
+\cite{cherkassky1997implementing}, for example the global and gap
+relabeling heuristics. Also has OK presentations of the different ways
+of choosing the order of discharge operations. In the results I should
+include some results with and without the heuristics.
+
+\cite{ahuja1997computational} has some good performance analysis.
+
+\cite{boykov2006graph} is a good source on how graph cuts are used in
+image analysis. It considers both one level graphs, and multilevel
+graphs. Which is nice. It has just a small section on the different
+algorithms that exist and that the special Boykov Kolmogorov algorithm
+is the best.
+
+\cite{derigs1989implementing} is also a good source for the different
+heuristics. And also on the neighbor list that should not be reset after
+each discharge (should try to implement this again, and reset at every
+relabel!). Also has some semi-usable notes on how the minimal cut is
+found after phase 1 of the algorithm is finished. Note that in this
+article the gap relabeling heuristic is also called RELABEL-GLOBAL.
+
+\begin{algorithm}
+\begin{algorithmic}
+ \Function{Push}{$u$}
+ \ForAll{$v$ neighbour of $u$}
+ \If{$res(u, v) > 0$}
+ \State $f \gets min(res(u, v), excess[u])$
+ \State $flow(u, v) \mathrel{+}= f$
+ \State $flow(v, u) \mathrel{-}= f$
+ \State $excess[u] \mathrel{-}= f$
+ \State $excess[v] \mathrel{+}= f$
+ \EndIf
+ \EndFor
+ \EndFunction
+\end{algorithmic}
+\caption{\sf The push procedure of the Push-Relabel algorithm}
+\end{algorithm}
+
+\begin{algorithm}
+\begin{algorithmic}
+ \Function{Relabel}{$u$}
+ \If{$u$ is only node at its height}
+ \Call{Gap}{$u$}
+ \Else
+ \State $height[u] \gets min(height[v] \, \forall v \in
+ neighbors[u] : res(u, v) > 0) + 1$
+ \EndIf
+ \EndFunction
+\end{algorithmic}
+\caption{\sf The relabel procedure of the Push-Relabel algorithm}
+\end{algorithm}
+
+\begin{algorithm}
+\begin{algorithmic}
+ \Function{Gap}{$u$}
+ \State $k \gets height[u]$
+ \ForAll{$v$ with height $\geq k$}
+ \State $height[v] \gets N$
+ \EndFor
+ \EndFunction
+\end{algorithmic}
+\caption{\sf The gap procedure of the Push-Relabel algorithm}
+\end{algorithm}
+
+\subsubsection{Other algorithms}
+Special made augmenting flow algorithm for this type of graphs. Used in
+software1, software2, etc. Is fast.
+
\usepackage{hyperref} % Lager hyperlinker i evt. pdf-dokument
% men har noen bugs, så den er kommentert
% bort her.
+
+\usepackage{tikz}
+\usetikzlibrary{matrix}
+\usetikzlibrary{arrows}
+\usetikzlibrary{positioning}
+\usetikzlibrary{calc}
+
+\usepackage{color}
+\usepackage{bbm} % Indicator function \mathbbm{1}
+\usepackage{caption}
+\usepackage{subcaption}
+
\usepackage{algorithm}
\usepackage{algpseudocode}
+
+\newtheorem{theorem}{Theorem}
+\newtheorem{definition}{Definition}
% Indeksgenerering er kommentert ut her. Ta bort prosenttegnene
% hvis du vil ha en indeks:
%\usepackage{makeidx}
%\makeindex
+\input{commands}
+
% Selve dokumentet begynner:
\begin{document}
% LaTeX om selv!
\input{introduction}
-\input{theory}
+\input{methods}
+\input{total}
+\input{cut}
+\input{noise}
\input{results}
% Bibliografi/referanseliste skal komme før appendiks
--- /dev/null
+\section{Methods in image restoration}
+
+\subsection{Gaussian filtering}
+Convolve the image with Gaussian function. Adds blur.
+
+\subsection{Anistropic diffusion}
+Some partial differential equation magic.
+
+\subsection{Non-local means}
+Somehow averaging over all the pixels in the image. Has a very good
+method noise.
+
+\subsection{Other filters}
+Median filter for example.
+
--- /dev/null
+\section{Noise theory}
+What kinds of noise occur in what kinds of settings. Which methods
+should be good for which
+
+\subsection{Gaussian noise}
+Common. Thermal noise? Central limit theorem. Read in Jähne.
+
+\subsection{Poisson noise}
+
+\subsection{Salt and pepper noise}
+
+\subsection{Shot noise}
+Same as Poisson noise? Quantum effects in the way light is emitted at
+the source. For regular photographs it will only be a determining factor
+in low light conditions.
+
+Dark current, or reverse bias leakage current, which is due to the
+random generation of electrons in the sensor, can also be modelled by a
+Poisson distribution.
+
% skriver i 'book'-dokumentklassen vil du kanskje bytte ut
% \section med \chapter, \subsection med \section, osv...
-\section{Theory}
-
-\subsection{Methods}
-
-\subsubsection{Anistrophic diffusion}
-Some partial differential equation magic.
-
-\subsubsection{Gaussian filtering}
-Convolve the image with Gaussian function. Adds blur.
-
-\subsubsection{Non-local means}
-Somehow averaging over all the pixels in the image. Has a very good
-method noise.
-
-\subsubsection{Non-linear filters}
-Median filter for example.
-
-\subsubsection{Total variation}
-Our method. WHAT IS SPLIT BREGMAN?
-
-\subsection{Continuous formulation}
-Formulate the energy function, describe the different terms. Total
-variation. \cite{caselles2011total} has a nice introduction to total
-variation. It describes the foundations, the spaces, and some methods,
-but mostly continuous.
+\section{Total variation}
+The method was initially introduced by L. Rudin, S. Osher, and E. Fatemi
+in \cite{rudin1992nonlinear} in the context of image restoration, but is
+now used in numerous other applications. One of its strong theoretical
+points is its ability to recover edges, compared to other restoration
+algorithms which might smooth over edges in the original image. The
+method can be formulated as minimizing the energy function
+\begin{equation}
+ E_v(u) =
+ \int_\Omega \abs{u - v}^p
+ + \beta \TV(u)
+ \label{eq:energy_function}
+\end{equation}
+where $v$ is the original image, and $\Omega$ is the image domain, in
+our case a subspace of $\mathbb{R}^2$. Normally, $p$ is chosen to be 1
+or 2. \fixme{Why? Theoretical niceness? Models the real world? There are
+also other norms describing other noise models?}
+
+The term $\TV(u)$ is the total variation of the image, and is defined as
+follows:
+\begin{definition}[Total variation]
+ Given a function $u \in L^1(\Omega)$, the total variation of $u$,
+ often written $\int_\Omega \abs{Du}$, where the $D$ is the gradient
+ taken in the distributional sense, is
+ \begin{equation}
+ \TV(u)
+ = \int_\Omega \abs{Du}
+ = \sup \left\{ \int_\Omega u \,\, \mathrm{div} \, \varphi :
+ \varphi \in C^\infty_c\left(\Omega, \mathbb{R}^N\right),
+ \norm{\varphi}_{L^\infty(\Omega)} \leq 1
+ \right\}.
+ \end{equation}
+ The test functions $\varphi$ are taken from
+ $C^\infty_c\left(\Omega, \mathbb{R}^N\right)$, the space of smooth
+ functions from $\Omega$ to $\mathbb{R}^N$ with compact support.
+\end{definition}
+
+Minimizing the total variationn will smooth out differences in the
+image. The term $\int_\Omega \abs{u - v}^p$ of
+\eqref{eq:energy_function} will constrain the solution image $u$ to be
+close to the input image $v$. Combining the two will give us an image
+close to the original, but more regular. The parameter $\beta$ controls
+how strongly we want to regularize the image.
+
+Following the notation used in \cite{darbon2006image}, we let $u_s$
+denote the value of the image $u$ at position $s \in \Omega$ and
+introduce the level sets $u^\lambda$ of the image
+\begin{equation}
+ u^\lambda = \idfun_{u \leq \lambda}
+ \label{eq:level_set}
+\end{equation}
+which will be important when dividing the problem into smaller
+subproblems.
+
+When discretizing the energy function \eqref{eq:energy_function} later
+on, we want to decompose it as a sum over the different levels, or pixel
+values, of the image. This is why we present the coarea formula here,
+and use it to write the total variation $\TV(u)$ as an integral over the
+range level values $\lambda$. The coarea formula and further references
+can be found in \cite{caselles2011total}.
+
+But first we need to introduce the space of functions with bounded
+variation and the perimeter of a set.
+\begin{definition}[Functions of bounded variation]
+ The space of functions with bounded variation $\BV(\Omega)$ is the
+ space of functions $u \in L^1(\Omega)$ such that the total variation
+ is finite, i.e.
+ \begin{equation}
+ \BV(\Omega) = \left\{ u \in L^1(\Omega) : \TV(u) < \infty
+ \right\}.
+ \end{equation}
+\end{definition}
+\begin{definition}[Set perimeter]
+ The perimeter of a set $E$ is defined as
+ \begin{equation}
+ P(E,\Omega) = \int_\Omega \abs{D \idfun_E}.
+ \end{equation}
+ A measurable set $E \subset \Omega$ is of finite perimeter in
+ $\Omega$ if $\idfun_E \in \BV(\Omega)$.
+\end{definition}
+These two definitions allows us to present the coarea formula
+\begin{theorem}[The coarea formula]
+ Let $u \in BV(\Omega)$. Then for almost every $\lambda$ the set
+ $u^\lambda$ is of finite perimeter, and one has the coarea formula
+ \begin{equation}
+ \TV(u)
+ = \int_\Omega \abs{D u}
+ = \int_\mathbb{R} P(u^\lambda) \, d\lambda
+ \label{eq:coarea_formula}
+ \end{equation}
+\end{theorem}
+This allows us to rewrite \eqref{eq:energy_function} to
+\begin{equation}
+ E_v(u) =
+ \int_\Omega \abs{u - v}^p
+ + \beta \int_\mathbb{R} P(u^\lambda) \, d\lambda,
+ \label{eq:energy_perimeter}
+\end{equation}
+which is a good starting point for the discretization. More analysis of
+the continuous problem is given in \cite{caselles2011total}, together
+with different numerical methods for solving it. Results on how the set
+of jumps in the resulting image $u$ is contained in the set of jumps in
+the original image $v$ are also presented.
+
+\fixme{\cite{caselles2011total} has a nice introduction to total
+ variation. It describes the foundations, the spaces, and some
+ methods, but mostly continuous. Could include some more theory here,
+ like the coarea stuff, and the number of discontinuities?}
+
+\fixme{Write something about the different methods that exist?}
\subsection{Discrete problem}
-How do we get to the discrete formulation. Why is it OK to do what we do
-with the gradient. How do we decompose the energy function into
-different levels and why is it OK to optimize each level separately. How
-do we choose the different neighbourhoods, and what could this mean for
-the result. Darbon and Sigelle. \cite{chan2011numerical} describes
-different methods for using total variation, and especially how one can
-use it in the discrete setting. It also goes into the min-cut method.
-Yay. Good book.
+Since digital images are given on a discrete grid, with values taken
+from a discrete and finite set of levels, we want to discretize
+\eqref{eq:energy_perimeter}, and solve the resulting discrete problem.
+We assume that the image is given on a discrete grid $S$,
+where the value of each pixel is taken from the set $\mathcal{L} = [0,
+\hdots, L-1]$. A reasonable assumption for grayscale images.
+
+\begin{figure}
+ \centering
+ \begin{tikzpicture}[
+ baseline=(M.center),
+ bg1/.style={fill=black!10},
+ bg2/.style={fill=black!25},
+ bg3/.style={fill=black!40},
+ bg4/.style={fill=black!55},
+ bg5/.style={fill=black!70},
+ ]
+ \matrix (M) [
+ matrix of nodes,
+ %column sep={6pt,between origins},
+ %row sep={6pt,between origins}
+ ]
+{
+ \node (A) {}; & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\
+ & 1 & 1 & 1 & 1 & 1 & 1 & 1 & & 1 & 1 & 1 & 1 \\
+ & 1 & 1 & & 1 & 1 & 1 & & & & 1 & 1 & 1 \\
+ & 1 & & & & 1 & 1 & & & & & 1 & 1 \\
+ \node (L) {$\mathcal{L}$}; & 1 & & & & & & & & & & & 1 \\
+ & |[bg5]| 5 & |[bg3]| 3 & |[bg2]| 2 & |[bg3]| 3 & |[bg4]| 4 & |[bg4]| 4 & |[bg2]| 2 & |[bg1]| 1 & |[bg2]| 2 & |[bg3]| 3 & |[bg4]| 4 & |[bg5]| 5 \\
+};
+
+\draw[->] (L) -- (A);
+\end{tikzpicture}
+\caption{The level sets of a one-dimensional image.}
+\label{fig:level_sets}
+\end{figure}
+After discretizing the image, we obtain a finite number of level sets
+$u^\lambda$. Figure \ref{fig:level_sets} shows an example of a
+one-dimensional image divided into its level sets.
+
+\begin{figure}[b]
+ \centering
+ \begin{subfigure}[b]{0.4\textwidth}
+ \centering
+ \begin{tikzpicture}[->,>=stealth',shorten >=1pt,auto,node distance=1.5cm,
+ thick,main node/.style={circle,draw}]
+
+ \node[main node] (1) {u};
+ \node[main node] (2) [below of=1] {};
+ \node[main node] (3) [right of=1] {v};
+ \node[main node] (4) [left of=1] {};
+ \node[main node] (5) [above of=1] {};
+
+ \path[every node/.style={font=\sffamily\small}]
+ (1) edge node [right] {} (2)
+ (1) edge node [right] {$w_{uv}$} (3)
+ (1) edge node [left] {} (4)
+ (1) edge node [left] {} (5);
+ \end{tikzpicture}
+ \caption{Size four neighborhood.}
+ \label{fig:c4_grid}
+ \end{subfigure}
+ ~
+ \begin{subfigure}[b]{0.4\textwidth}
+ \centering
+ \begin{tikzpicture}[->,>=stealth',shorten >=1pt,auto,node distance=1.5cm,
+ thick,main node/.style={circle,draw}]
+
+ \node[main node] (1) {u};
+ \node[main node] (2) [below of=1] {};
+ \node[main node] (3) [right of=1] {};
+ \node[main node] (4) [left of=1] {};
+ \node[main node] (5) [above of=1] {};
+ \node[main node] (6) [below right=1.5cm of 1] {};
+ \node[main node] (7) [above right=1.5cm of 1] {};
+ \node[main node] (8) [above left=1.5cm of 1] {};
+ \node[main node] (9) [below left=1.5cm of 1] {};
+
+ \path[every node/.style={font=\sffamily\small}]
+ (1) edge node [right] {} (2)
+ (1) edge node [right] {} (3)
+ (1) edge node [left] {} (4)
+ (1) edge node [left] {} (5)
+ (1) edge node [right] {} (6)
+ (1) edge node [right] {} (7)
+ (1) edge node [left] {} (8)
+ (1) edge node [left] {} (9);
+ \end{tikzpicture}
+ \caption{Size eight neighborhood.}
+ \label{fig:c8_grid}
+ \end{subfigure}
+ \caption{Two common neighborhood stencils.}
+ \label{fig:common_neighborhoods}
+\end{figure}
+We discretize the total variation in \eqref{eq:coarea_formula} as
+follows
+\begin{equation}
+ \TV(u)
+ = \sum_{\lambda = 0}^{L-2} P(u^\lambda)
+ = \sum_{\lambda = 0}^{L-2} \sum_{(u,v)} w_{uv}
+ \abs{u_u^\lambda - u_v^\lambda},
+ \label{eq:tv_discrete_int}
+\end{equation}
+where the second sum is over nodes in a neighborhood relation and the
+$w_{uv}$ is the weight of this relation. The sum only goes up to the
+level $L-2$ since $u^\lambda$ is constant equal to 1 in every pixel of
+the image. Figure \ref{fig:common_neighborhoods} shows the two most used
+neighbordhood stencils, which will be discussed more thorougly in
+\fixme{Section ??}. Boykov and Kolmogorov argue in
+\cite{boykov2003computing} that if the weight is chosen as
+\begin{equation}
+ w_{uv} = \frac{h^2 \cdot \Delta \phi_{uv}}{2 \cdot \abs{d_{uv}}},
+\end{equation}
+the discrete perimeter in \eqref{eq:tv_discrete_int} converges to the
+continuous perimeter in \eqref{eq:coarea_formula}. Here, $h$ is the grid
+size and $d_{uv}$ is the euclidean distance of the edge. The $\Delta
+\phi_{uv}$ is the difference between the angle of this edge and the
+next edge, if the edges are sorted by increasing angles. These
+parameters are also shown in Figure \ref{fig:common_neighborhoods}.
+Boykov and Kolmogorov prove we have convergence when all of $h$, $\Delta
+\phi_{uv}$, and $d_{uv}$ go to zero.
+\fixme{Maybe move (and expand) this somewhere else? Where the grids are
+discussed for example? Also, the section is a bit messy.}
+
+\fixme{BAD TRANSITION}
+
+For two binary variables $a$ and $b$ we can easily verify that $\abs{a -
+b} = a + b - 2 a b$. Using this we rewrite \eqref{eq:tv_discrete_int} to
+\begin{equation}
+ \TV(u) = \sum_{\lambda = 0}^{L-2} \sum_{(u,v)} w_{uv}
+ \left(
+ \left(1 - 2 u_t^\lambda \right) u_s^\lambda + u_t^\lambda
+ \right).
+ \label{eq:tv_discrete}
+\end{equation}
+
+We define the following function for some pixel value $x$ and some pixel
+position $u$ in the original image $v$
+\begin{equation}
+ N_u(x) = \abs{x - v_u}^p,
+\end{equation}
+such that the $\int_\Omega \abs{u - v}^p$ of \eqref{eq:energy_perimeter}
+can be discretized in the following way \fixme{Give this a name?}
+\begin{equation}
+ \sum_u \abs{u_u - v_u}^p = \sum_u N_u(u_u)
+ \label{eq:norm_discrete_int}.
+\end{equation}
+We want to write the energy as a sum over the different levels of the
+image, so we use the following decomposition formula which holds for any
+function $N(k)$ taking values $k \in \mathcal{L}$
+\begin{equation}
+\begin{aligned}
+ N(k) &= \sum_{\lambda=0}^{k-1} \big(
+ N(\lambda + 1) - N(\lambda)
+ \big) + N(0) \\
+ &= \sum_{\lambda=0}^{L-2} \big(
+ N(\lambda + 1) - N(\lambda)
+ \big) \mathbbm{1}_{\lambda < k} + N(0) \\
+ &= \sum_{\lambda=0}^{L-2} \big(
+ N(\lambda + 1) - N(\lambda)
+ \big) (1 - u^\lambda) + N(0).
+\end{aligned}
+\end{equation}
+This allows us to rewrite \eqref{eq:norm_discrete_int} and obtain
+\begin{equation}
+ \sum_{\lambda=0}^{L-2} \sum_u
+ \big(
+ N_u(\lambda + 1) -
+ N_u(\lambda)
+ \big)
+ (1 - u^\lambda) + N_u(0)
+\end{equation}
+We have now discretized the energy function and decomposed it into a sum
+over all the levels $\lambda$ of the image. If we minimize each level
+separately it is obvious that we also minimize the total energy. The
+question is if the obtained level sets $u^\lambda$ can be combined to
+make an output image. The level sets were defined as $u^\lambda =
+\mathbbm{1}_{u \leq \lambda}$, so we need them to be monotonely
+increasing in increasing level values, i.e.\
+\begin{equation}
+ u^\lambda_u \leq u^\mu_u \quad \forall \lambda \leq \mu,
+ \quad \forall u \in S.
+\end{equation}
+
+\fixme{Want to use $u$ and $v$ as pixels/nodes, what to use as image?
+Since $s$ and $t$ normally is the source and sink.}
+
+\fixme{Some like to call one term the fidelity term, and the other the
+regularization term?}
+
+\fixme{\cite{boykov2003computing} motivates the discrete estimation of
+ the perimeter of a set, and it makes sense. The size of the
+ perimeter should be proportional to the number of edges it crosses,
+ which is why we sum up $\abs{u^\lambda_s - u^\lambda_t}$ (which is
+ only different from zero at the cut).
+}
+
+\fixme{How do we get to the discrete formulation. Why is it OK to do
+ what we do with the gradient. How do we decompose the energy
+ function into different levels and why is it OK to optimize each
+ level separately. How do we choose the different neighbourhoods, and
+ what could this mean for the result. Darbon and Sigelle.
+ \cite{chan2011numerical} describes different methods for using total
+ variation, and especially how one can use it in the discrete
+ setting. It also goes into the min-cut method. Yay. Good book.
+}
\subsection{Graph cut formulation}
Short introduction to graphs. Notation. What is a cut, what is a minimal
Special made augmenting flow algorithm for this type of graphs. Used in
software1, software2, etc. Is fast.
-\subsection{Noise theory}
-What kinds of noise occur in what kinds of settings. Which methods
-should be good for which
-
-\subsubsection{Gaussian noise}
-Common. Thermal noise? Central limit theorem. Read in Jähne.
-
-\subsubsection{Poisson noise}
-
-\subsubsection{Salt and pepper noise}
-
-\subsubsection{Shot noise}
-Same as Poisson noise? Quantum effects in the way light is emitted at
-the source. For regular photographs it will only be a determining factor
-in low light conditions.
-
-Dark current, or reverse bias leakage current, which is due to the
-random generation of electrons in the sensor, can also be modelled by a
-Poisson distribution.
-
% NB: Bruken av \and mellom navn!
\titlepage
-\title{Image restoration using graph cuts}
+\title{Total variation based image restoration using graph cuts}
\author{Bjørn Rustad}
\date{\today}
\maketitle
--- /dev/null
+% Dette eksempelet er laget for article-dokumentklassen. Hvis
+% skriver i 'book'-dokumentklassen vil du kanskje bytte ut
+% \section med \chapter, \subsection med \section, osv...
+
+\section{Total variation}
+The method was initially introduced by L. Rudin, S. Osher, and E. Fatemi
+in \cite{rudin1992nonlinear} in the context of image restoration, but is
+now used in numerous other applications. One of its strong theoretical
+points is its ability to recover edges, compared to other restoration
+algorithms which might smooth over edges in the original image. The
+method can be formulated as minimizing the energy function
+\begin{equation}
+ E_v(u) =
+ \int_\Omega \abs{u - v}^p
+ + \beta \TV(u)
+ \label{eq:energy_function}
+\end{equation}
+where $v$ is the original image, and $\Omega$ is the image domain, in
+our case a subspace of $\mathbb{R}^2$. Normally, $p$ is chosen to be 1
+or 2. \fixme{Why? Theoretical niceness? Models the real world? There are
+also other norms describing other noise models?}
+
+The term $\TV(u)$ is the total variation of the image, and is defined as
+follows:
+\begin{definition}[Total variation]
+ Given a function $u \in L^1(\Omega)$, the total variation of $u$,
+ often written $\int_\Omega \abs{Du}$, where the $D$ is the gradient
+ taken in the distributional sense, is
+ \begin{equation}
+ \TV(u)
+ = \int_\Omega \abs{Du}
+ = \sup \left\{ \int_\Omega u \,\, \mathrm{div} \, \varphi :
+ \varphi \in C^\infty_c\left(\Omega, \mathbb{R}^N\right),
+ \norm{\varphi}_{L^\infty(\Omega)} \leq 1
+ \right\}.
+ \end{equation}
+ The test functions $\varphi$ are taken from
+ $C^\infty_c\left(\Omega, \mathbb{R}^N\right)$, the space of smooth
+ functions from $\Omega$ to $\mathbb{R}^N$ with compact support.
+\end{definition}
+
+Minimizing the total variationn will smooth out differences in the
+image. The term $\int_\Omega \abs{u - v}^p$ of
+\eqref{eq:energy_function} will constrain the solution image $u$ to be
+close to the input image $v$. Combining the two will give us an image
+close to the original, but more regular. The parameter $\beta$ controls
+how strongly we want to regularize the image.
+
+Following the notation used in \cite{darbon2006image}, we let $u_s$
+denote the value of the image $u$ at position $s \in \Omega$ and
+introduce the level sets $u^\lambda$ of the image
+\begin{equation}
+ u^\lambda = \idfun_{u \leq \lambda}
+ \label{eq:level_set}
+\end{equation}
+which will be important when dividing the problem into smaller
+subproblems.
+
+When discretizing the energy function \eqref{eq:energy_function} later
+on, we want to decompose it as a sum over the different levels, or pixel
+values, of the image. This is why we present the coarea formula here,
+and use it to write the total variation $\TV(u)$ as an integral over the
+range level values $\lambda$. The coarea formula and further references
+can be found in \cite{caselles2011total}.
+
+But first we need to introduce the space of functions with bounded
+variation and the perimeter of a set.
+\begin{definition}[Functions of bounded variation]
+ The space of functions with bounded variation $\BV(\Omega)$ is the
+ space of functions $u \in L^1(\Omega)$ such that the total variation
+ is finite, i.e.
+ \begin{equation}
+ \BV(\Omega) = \left\{ u \in L^1(\Omega) : \TV(u) < \infty
+ \right\}.
+ \end{equation}
+\end{definition}
+\begin{definition}[Set perimeter]
+ The perimeter of a set $E$ is defined as
+ \begin{equation}
+ P(E,\Omega) = \int_\Omega \abs{D \idfun_E}.
+ \end{equation}
+ A measurable set $E \subset \Omega$ is of finite perimeter in
+ $\Omega$ if $\idfun_E \in \BV(\Omega)$.
+\end{definition}
+These two definitions allows us to present the coarea formula
+\begin{theorem}[The coarea formula]
+ Let $u \in BV(\Omega)$. Then for almost every $\lambda$ the set
+ $u^\lambda$ is of finite perimeter, and one has the coarea formula
+ \begin{equation}
+ \TV(u)
+ = \int_\Omega \abs{D u}
+ = \int_\mathbb{R} P(u^\lambda) \, d\lambda
+ \label{eq:coarea_formula}
+ \end{equation}
+\end{theorem}
+This allows us to rewrite \eqref{eq:energy_function} to
+\begin{equation}
+ E_v(u) =
+ \int_\Omega \abs{u - v}^p
+ + \beta \int_\mathbb{R} P(u^\lambda) \, d\lambda,
+ \label{eq:energy_perimeter}
+\end{equation}
+which is a good starting point for the discretization. More analysis of
+the continuous problem is given in \cite{caselles2011total}, together
+with different numerical methods for solving it. Results on how the set
+of jumps in the resulting image $u$ is contained in the set of jumps in
+the original image $v$ are also presented.
+
+\fixme{\cite{caselles2011total} has a nice introduction to total
+ variation. It describes the foundations, the spaces, and some
+ methods, but mostly continuous. Could include some more theory here,
+ like the coarea stuff, and the number of discontinuities?}
+
+\fixme{Write something about the different methods that exist?}
+
+\subsection{Discrete problem}
+Since digital images are given on a discrete grid, with values taken
+from a discrete and finite set of levels, we want to discretize
+\eqref{eq:energy_perimeter}, and solve the resulting discrete problem.
+We assume that the image is given on a discrete grid $S$,
+where the value of each pixel is taken from the set $\mathcal{L} = [0,
+\hdots, L-1]$. A reasonable assumption for grayscale images.
+
+\begin{figure}
+ \centering
+ \begin{tikzpicture}[
+ baseline=(M.center),
+ bg1/.style={fill=black!10},
+ bg2/.style={fill=black!25},
+ bg3/.style={fill=black!40},
+ bg4/.style={fill=black!55},
+ bg5/.style={fill=black!70},
+ ]
+ \matrix (M) [
+ matrix of nodes,
+ %column sep={6pt,between origins},
+ %row sep={6pt,between origins}
+ ]
+{
+ \node (A) {}; & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\
+ & 1 & 1 & 1 & 1 & 1 & 1 & 1 & & 1 & 1 & 1 & 1 \\
+ & 1 & 1 & & 1 & 1 & 1 & & & & 1 & 1 & 1 \\
+ & 1 & & & & 1 & 1 & & & & & 1 & 1 \\
+ \node (L) {$\mathcal{L}$}; & 1 & & & & & & & & & & & 1 \\
+ & |[bg5]| 5 & |[bg3]| 3 & |[bg2]| 2 & |[bg3]| 3 & |[bg4]| 4 & |[bg4]| 4 & |[bg2]| 2 & |[bg1]| 1 & |[bg2]| 2 & |[bg3]| 3 & |[bg4]| 4 & |[bg5]| 5 \\
+};
+
+\draw[->] (L) -- (A);
+\end{tikzpicture}
+\caption{The level sets of a one-dimensional image.}
+\label{fig:level_sets}
+\end{figure}
+After discretizing the image, we obtain a finite number of level sets
+$u^\lambda$. Figure \ref{fig:level_sets} shows an example of a
+one-dimensional image divided into its level sets.
+
+\begin{figure}[b]
+ \centering
+ \begin{subfigure}[b]{0.4\textwidth}
+ \centering
+ \begin{tikzpicture}[->,>=stealth',shorten >=1pt,auto,node distance=1.5cm,
+ thick,main node/.style={circle,draw}]
+
+ \node[main node] (1) {u};
+ \node[main node] (2) [below of=1] {};
+ \node[main node] (3) [right of=1] {v};
+ \node[main node] (4) [left of=1] {};
+ \node[main node] (5) [above of=1] {};
+
+ \path[every node/.style={font=\sffamily\small}]
+ (1) edge node [right] {} (2)
+ (1) edge node [right] {$w_{uv}$} (3)
+ (1) edge node [left] {} (4)
+ (1) edge node [left] {} (5);
+ \end{tikzpicture}
+ \caption{Size four neighborhood.}
+ \label{fig:c4_grid}
+ \end{subfigure}
+ ~
+ \begin{subfigure}[b]{0.4\textwidth}
+ \centering
+ \begin{tikzpicture}[->,>=stealth',shorten >=1pt,auto,node distance=1.5cm,
+ thick,main node/.style={circle,draw}]
+
+ \node[main node] (1) {u};
+ \node[main node] (2) [below of=1] {};
+ \node[main node] (3) [right of=1] {};
+ \node[main node] (4) [left of=1] {};
+ \node[main node] (5) [above of=1] {};
+ \node[main node] (6) [below right=1.5cm of 1] {};
+ \node[main node] (7) [above right=1.5cm of 1] {};
+ \node[main node] (8) [above left=1.5cm of 1] {};
+ \node[main node] (9) [below left=1.5cm of 1] {};
+
+ \path[every node/.style={font=\sffamily\small}]
+ (1) edge node [right] {} (2)
+ (1) edge node [right] {} (3)
+ (1) edge node [left] {} (4)
+ (1) edge node [left] {} (5)
+ (1) edge node [right] {} (6)
+ (1) edge node [right] {} (7)
+ (1) edge node [left] {} (8)
+ (1) edge node [left] {} (9);
+ \end{tikzpicture}
+ \caption{Size eight neighborhood.}
+ \label{fig:c8_grid}
+ \end{subfigure}
+ \caption{Two common neighborhood stencils.}
+ \label{fig:common_neighborhoods}
+\end{figure}
+We discretize the total variation in \eqref{eq:coarea_formula} as
+follows
+\begin{equation}
+ \TV(u)
+ = \sum_{\lambda = 0}^{L-2} P(u^\lambda)
+ = \sum_{\lambda = 0}^{L-2} \sum_{(u,v)} w_{uv}
+ \abs{u_u^\lambda - u_v^\lambda},
+ \label{eq:tv_discrete_int}
+\end{equation}
+where the second sum is over nodes in a neighborhood relation and the
+$w_{uv}$ is the weight of this relation. The sum only goes up to the
+level $L-2$ since $u^\lambda$ is constant equal to 1 in every pixel of
+the image. Figure \ref{fig:common_neighborhoods} shows the two most used
+neighbordhood stencils, which will be discussed more thorougly in
+\fixme{Section ??}. Boykov and Kolmogorov argue in
+\cite{boykov2003computing} that if the weight is chosen as
+\begin{equation}
+ w_{uv} = \frac{h^2 \cdot \Delta \phi_{uv}}{2 \cdot \abs{d_{uv}}},
+\end{equation}
+the discrete perimeter in \eqref{eq:tv_discrete_int} converges to the
+continuous perimeter in \eqref{eq:coarea_formula}. Here, $h$ is the grid
+size and $d_{uv}$ is the euclidean distance of the edge. The $\Delta
+\phi_{uv}$ is the difference between the angle of this edge and the
+next edge, if the edges are sorted by increasing angles. These
+parameters are also shown in Figure \ref{fig:common_neighborhoods}.
+Boykov and Kolmogorov prove we have convergence when all of $h$, $\Delta
+\phi_{uv}$, and $d_{uv}$ go to zero.
+\fixme{Maybe move (and expand) this somewhere else? Where the grids are
+discussed for example? Also, the section is a bit messy.}
+
+\fixme{BAD TRANSITION}
+
+For two binary variables $a$ and $b$ we can easily verify that $\abs{a -
+b} = a + b - 2 a b$. Using this we rewrite \eqref{eq:tv_discrete_int} to
+\begin{equation}
+ \TV(u) = \sum_{\lambda = 0}^{L-2} \sum_{(u,v)} w_{uv}
+ \left(
+ \left(1 - 2 u_t^\lambda \right) u_s^\lambda + u_t^\lambda
+ \right).
+ \label{eq:tv_discrete}
+\end{equation}
+
+We define the following function for some pixel value $x$ and some pixel
+position $u$ in the original image $v$
+\begin{equation}
+ N_u(x) = \abs{x - v_u}^p,
+\end{equation}
+such that the $\int_\Omega \abs{u - v}^p$ of \eqref{eq:energy_perimeter}
+can be discretized in the following way \fixme{Give this a name?}
+\begin{equation}
+ \sum_u \abs{u_u - v_u}^p = \sum_u N_u(u_u)
+ \label{eq:norm_discrete_int}.
+\end{equation}
+We want to write the energy as a sum over the different levels of the
+image, so we use the following decomposition formula which holds for any
+function $N(k)$ taking values $k \in \mathcal{L}$
+\begin{equation}
+\begin{aligned}
+ N(k) &= \sum_{\lambda=0}^{k-1} \big(
+ N(\lambda + 1) - N(\lambda)
+ \big) + N(0) \\
+ &= \sum_{\lambda=0}^{L-2} \big(
+ N(\lambda + 1) - N(\lambda)
+ \big) \mathbbm{1}_{\lambda < k} + N(0) \\
+ &= \sum_{\lambda=0}^{L-2} \big(
+ N(\lambda + 1) - N(\lambda)
+ \big) (1 - u^\lambda) + N(0).
+\end{aligned}
+\end{equation}
+This allows us to rewrite \eqref{eq:norm_discrete_int} and obtain
+\begin{equation}
+ \sum_{\lambda=0}^{L-2} \sum_u
+ \big(
+ N_u(\lambda + 1) -
+ N_u(\lambda)
+ \big)
+ (1 - u^\lambda) + N_u(0)
+\end{equation}
+We have now discretized the energy function and decomposed it into a sum
+over all the levels $\lambda$ of the image. If we minimize each level
+separately it is obvious that we also minimize the total energy. The
+question is if the obtained level sets $u^\lambda$ can be combined to
+make an output image. The level sets were defined as $u^\lambda =
+\mathbbm{1}_{u \leq \lambda}$, so we need them to be monotonely
+increasing in increasing level values, i.e.\
+\begin{equation}
+ u^\lambda_u \leq u^\mu_u \quad \forall \lambda \leq \mu,
+ \quad \forall u \in S.
+\end{equation}
+
+\fixme{Want to use $u$ and $v$ as pixels/nodes, what to use as image?
+Since $s$ and $t$ normally is the source and sink.}
+
+\fixme{Some like to call one term the fidelity term, and the other the
+regularization term?}
+
+\fixme{\cite{boykov2003computing} motivates the discrete estimation of
+ the perimeter of a set, and it makes sense. The size of the
+ perimeter should be proportional to the number of edges it crosses,
+ which is why we sum up $\abs{u^\lambda_s - u^\lambda_t}$ (which is
+ only different from zero at the cut).
+}
+
+\fixme{How do we get to the discrete formulation. Why is it OK to do
+ what we do with the gradient. How do we decompose the energy
+ function into different levels and why is it OK to optimize each
+ level separately. How do we choose the different neighbourhoods, and
+ what could this mean for the result. Darbon and Sigelle.
+ \cite{chan2011numerical} describes different methods for using total
+ variation, and especially how one can use it in the discrete
+ setting. It also goes into the min-cut method. Yay. Good book.
+}
+
projects / prosjektoppgave / commitdiff