From f2c6f401dcbc6bcda950541aa961d5c8d8c6c90a Mon Sep 17 00:00:00 2001 From: =?utf8?q?Bj=C3=B8rn=20Rustad?= Date: Wed, 12 Mar 2014 18:44:02 +0100 Subject: [PATCH] Wrote a lot, and some more --- bib.bib | 31 +++++ commands.tex | 7 + cut.tex | 113 ++++++++++++++++ main.tex | 22 ++- methods.tex | 15 ++ noise.tex | 20 +++ theory.tex | 369 ++++++++++++++++++++++++++++++++++++++++++-------- titlepage.tex | 2 +- total.tex | 323 +++++++++++++++++++++++++++++++++++++++++++ 9 files changed, 847 insertions(+), 55 deletions(-) create mode 100644 commands.tex create mode 100644 cut.tex create mode 100644 methods.tex create mode 100644 noise.tex create mode 100644 total.tex diff --git a/bib.bib b/bib.bib index a22d23b..fb93e21 100644 --- a/bib.bib +++ b/bib.bib @@ -90,4 +90,35 @@ publisher={Springer} } +@article{rudin1992nonlinear, + title={Nonlinear total variation based noise removal algorithms}, + author={Rudin, Leonid I and Osher, Stanley and Fatemi, Emad}, + journal={Physica D: Nonlinear Phenomena}, + volume={60}, + number={1}, + pages={259--268}, + year={1992}, + publisher={Elsevier} +} + +@article{darbon2006image, + title={Image restoration with discrete constrained total variation part I: Fast and exact optimization}, + author={Darbon, J{\'e}r{\^o}me and Sigelle, Marc}, + journal={Journal of Mathematical Imaging and Vision}, + volume={26}, + number={3}, + pages={261--276}, + year={2006}, + publisher={Springer} +} + +@inproceedings{boykov2003computing, + title={Computing geodesics and minimal surfaces via graph cuts}, + author={Boykov, Yuri and Kolmogorov, Vladimir}, + booktitle={Computer Vision, 2003. Proceedings. Ninth IEEE International Conference on}, + pages={26--33}, + year={2003}, + organization={IEEE} +} + diff --git a/commands.tex b/commands.tex new file mode 100644 index 0000000..8150234 --- /dev/null +++ b/commands.tex @@ -0,0 +1,7 @@ +\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}} diff --git a/cut.tex b/cut.tex new file mode 100644 index 0000000..da8f62b --- /dev/null +++ b/cut.tex @@ -0,0 +1,113 @@ +\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. + diff --git a/main.tex b/main.tex index dc75594..46ab832 100644 --- a/main.tex +++ b/main.tex @@ -19,14 +19,31 @@ \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} @@ -64,7 +81,10 @@ % LaTeX om selv! \input{introduction} -\input{theory} +\input{methods} +\input{total} +\input{cut} +\input{noise} \input{results} % Bibliografi/referanseliste skal komme før appendiks diff --git a/methods.tex b/methods.tex new file mode 100644 index 0000000..8ac6b82 --- /dev/null +++ b/methods.tex @@ -0,0 +1,15 @@ +\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. + diff --git a/noise.tex b/noise.tex new file mode 100644 index 0000000..e545274 --- /dev/null +++ b/noise.tex @@ -0,0 +1,20 @@ +\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. + diff --git a/theory.tex b/theory.tex index 57516e0..4a438da 100644 --- a/theory.tex +++ b/theory.tex @@ -2,41 +2,324 @@ % 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 @@ -151,23 +434,3 @@ article the gap relabeling heuristic is also called RELABEL-GLOBAL. 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. - diff --git a/titlepage.tex b/titlepage.tex index 82d2256..5a4f093 100644 --- a/titlepage.tex +++ b/titlepage.tex @@ -2,7 +2,7 @@ % 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 diff --git a/total.tex b/total.tex new file mode 100644 index 0000000..64e7e0b --- /dev/null +++ b/total.tex @@ -0,0 +1,323 @@ +% 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. +} + -- 2.47.3