From 24b84b2759fccef4131d4dccced3da328101fb70 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Bj=C3=B8rn=20Rustad?= Date: Tue, 10 Jun 2014 10:15:55 +0200 Subject: [PATCH] Fixes here and there --- bib.bib | 14 ++-- cut.tex | 121 ++++++++++++------------------- fig/big_network.tex | 41 +++++++++++ fig/norm_evolution.tex | 6 +- fig/norm_subgraph.tex | 12 ++-- flow.tex | 67 +++++++++--------- introduction.tex | 2 - main.tex | 16 +++-- methods.tex | 37 +++++----- pushrelabel.tex | 157 +++++++++++++++++++++++------------------ results.tex | 24 ++++--- titlepage.tex | 37 +++++++++- total.tex | 121 ++++++++++++++++--------------- 13 files changed, 364 insertions(+), 291 deletions(-) create mode 100644 fig/big_network.tex diff --git a/bib.bib b/bib.bib index 6d3cf18..711861b 100644 --- a/bib.bib +++ b/bib.bib @@ -139,12 +139,16 @@ article{derigs1989implementing, } @incollection{caselles2011total, - title={Total variation in imaging}, - author={Caselles, Vicent and Chambolle, Antonin and Novaga, Matteo}, - journal={Handbook of Mathematical Methods in Imaging}, - pages={1016--1057}, year={2011}, - publisher={Springer} + isbn={978-0-387-92919-4}, + booktitle={Handbook of Mathematical Methods in Imaging}, + editor={Scherzer, Otmar}, + doi={10.1007/978-0-387-92920-0_23}, + title={Total Variation in Imaging}, + url={http://dx.doi.org/10.1007/978-0-387-92920-0_23}, + publisher={Springer New York}, + author={Caselles, V. and Chambolle, A. and Novaga, M.}, + pages={1016-1057} } @incollection{chan2011numerical, diff --git a/cut.tex b/cut.tex index 147851d..6f76169 100644 --- a/cut.tex +++ b/cut.tex @@ -8,15 +8,15 @@ together is the fact that we should be able to stack the thresholded images in the end. We will in this section present a graph cut approach to the level -subproblem. A graph is constructed containing all the pixels as -vertices, with carefully chosen edge weights such that a minimum cut of -the graph will also minimize the energy function of +subproblem. A graph is constructed containing one vertex per pixel, with +carefully chosen edge weights such that a minimum cut of the graph will +also minimize the energy function of \eqref{eq:total_energy}. \subsection{Networks} 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$. We say that if $(u, v) \in E$ there is +$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 @@ -30,7 +30,7 @@ c)$. There are two special nodes in the network, the source $s$ and the sink $t$. Contrary to other nodes, which can neither produce nor receive -excess flow, the source can produce and the sink can receive, an +excess flow, the source can produce and the sink can receive an unlimited amount of flow. The most basic problem in flow network theory is the question of how much flow it is possible to send through the network from the source to the sink. @@ -49,24 +49,20 @@ divides the vertex set in two, separating the source from the sink. all $s$-$t$-cuts. \label{def:s_t_cut} \end{definition} -We will see later that if a vertex is in the set $S$ or $T$ is what -indicates whether we should set the binary variable $u^\lambda_x$ to 0 -or 1. Note that there might exist several minimum $s$-$t$-cuts in a network, -resulting in different partitions of $V$. This is why we need to be -careful later that we are sure the cuts we obtain will result in -stackable thresholded images. +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{Network representable energy functions} The next step is to find a way to construct a network 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 networks -representing the separate summands of the energy function. For these -small networks it will be easy to verify that the minimal cut also -minimizes the corresponding part of the energy function. These small -networks can then be merged giving a network representing the complete -energy function \eqref{eq:total_energy}. +representing the separate summands of the energy. For these small +networks 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 network representing the complete energy function in +\eqref{eq:total_energy}. First we need to establish the definition of a network representable function, presented by Kolmogorov and Zabih in @@ -79,7 +75,7 @@ function, presented by Kolmogorov and Zabih in \{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 + 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:network_representable} @@ -92,7 +88,7 @@ representing $\mathcal{E}$. Furthermore Kolmogorov and Zabih present an important result concerning what kinds of functions are network-representable. \begin{theorem}[Identification of network representable functions] - Given an energy function $\mathcal{E}$ of $n$ binary variables on + 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) @@ -110,7 +106,8 @@ This theorem will allow us to verify that our energy function actually is network representable. Finally, the following theorem, proved by Kolmogorov and Zabih in -\cite{kolmogorov2002energy}, will be crucial in our graph construction. +\cite{kolmogorov2002energy}, will be crucial in our network +construction. \begin{theorem}[Additivity] The sum of a finite number of network-representable functions \begin{equation} @@ -135,8 +132,8 @@ depend on. \subsection{Network construction} \label{sec:network_construction} We will construct a network in such a way that if a variable -$u^\lambda_x$ ends up on the source side of the cut, we set $u^\lambda_x -= 0$, and if it ends up on the sink side, we set $u^\lambda_x = 1$ as in +$u^\lambda_x$ ends up on the source side of the cut we set $u^\lambda_x += 0$, and if it ends up on the sink side we set $u^\lambda_x = 1$, as in Definition \ref{def:network_representable}. This is an arbitrary choice, but still something we have to keep in mind through the rest of the section. @@ -150,15 +147,15 @@ term coming from our aim to minimize the total variation. The fidelity term of our energy function in \eqref{eq:total_energy} simplifies to \begin{align} - E^x(0) &= + E_\lambda^x(0) &= N_x(\lambda + 1) - N_x(\lambda) \label{eq:fidelity_energy0} \\ - E^x(1) &= 0 + E_\lambda^x(1) &= 0 \label{eq:fidelity_energy1} \end{align} -where $E^x(0)$ might be positive or negative depending on $\lambda$ and -the pixel value $v_x$. +where $E_\lambda^x(0)$ might be positive or negative depending on +$\lambda$ and the pixel value $v_x$. \begin{figure} \input{fig/norm_subgraph} @@ -166,7 +163,7 @@ the pixel value $v_x$. Figure \ref{fig:norm_subgraph} shows how networks can be constructed to represent this part of the total energy. The construction differs -depending on whether $E^x(0)$ is positive or negative. Table +depending on whether $E_\lambda^x(0)$ is positive or negative. Table \ref{tab:fid_energy} shows how the cuts correspond to the values of $u^\lambda_x$ and we can easily verify that the constructed network actually represents the fidelity term in the energy function. @@ -174,22 +171,22 @@ actually represents the fidelity term in the energy function. \begin{table} \centering \caption{Each row represents one of the two possible values of - $u^\lambda_x \in \{0,1\}$. The energy $E^x(u^\lambda_x)$ and - minimum cut obtaining this configuration is shown. The last two - columns show the capacities of the cut for each of the two - network constructions in Figure - \ref{fig:norm_subgraph}. We verify that, for each of the two - network constructions, the cut capacities are equal to the - energies plus a constant. + $u^\lambda_x \in \{0,1\}$. The energy $E_\lambda^x(u^\lambda_x)$ + and minimum cut obtaining this configuration is shown. The last + two columns show the capacities of the cut for each of the two + network constructions in Figure \ref{fig:norm_subgraph}. We + verify that for each of the two network constructions, the cut + capacities are equal to the energies, plus a constant. } \begin{tabular}{ccccc} \hline - $u^\lambda_x$ & $E^x(u^\lambda_x)$ & Min.\ cut $(S, T)$ & - Network \subref{fig:norm_subgraph_pos} cut cap. & Network + $u^\lambda_x$ & $E_\lambda^x(u^\lambda_x)$ & Min.\ cut $(S, T)$ + & Network \subref{fig:norm_subgraph_pos} cut cap. & Network \subref{fig:norm_subgraph_neg} cut cap. \\ \hline - $0$ & $E^x(0)$ & $(\{s, u^\lambda_x\}, \{t\})$ & $E^x(0)$ & $0$ \\ - $1$ & $0$ & $(\{s\}, \{u^\lambda_x, t\})$ & $0$ & $-E^x(0)$ + $0$ & $E_\lambda^x(0)$ & $(\{s, u^\lambda_x\}, \{t\})$ & + $E_\lambda^x(0)$ & $0$ \\ + $1$ & $0$ & $(\{s\}, \{u^\lambda_x, t\})$ & $0$ & $-E_\lambda^x(0)$ \end{tabular} \label{tab:fid_energy} \end{table} @@ -198,10 +195,7 @@ actually represents the fidelity term in the energy function. For our neighboring relation in \eqref{eq:total_energy} of the form \begin{equation} E^{x,y}(u^\lambda_x, u^\lambda_y) = - w_{xy} - \left( - \left(1 - 2 u_y^\lambda \right) u_x^\lambda + u_y^\lambda - \right). + w_{xy} \abs{u^\lambda_x - u^\lambda_y} \end{equation} we have \begin{equation} @@ -236,7 +230,7 @@ $u^\lambda_x$ corresponds to cuts in the network. configuration the corresponding energy and the cut yielding this configuration is shown. The last two columns show the capacities of the cut in the two alternative network constructions shown in - Figure \ref{fig:neigh_subgraph}. We verify that, for each of the + Figure \ref{fig:neigh_subgraph}. We verify that for each of the two network constructions, the cut capacities are equal to the energies, plus a constant. } @@ -257,40 +251,13 @@ $u^\lambda_x$ corresponds to cuts in the network. \end{table} \begin{figure} - \centering - \begin{tikzpicture}[scale=1.3] - - \node[vertex] (s) at (1, 0) {s}; - \node[vertex] (t) at (1, -5) {t}; - - \begin{scope}[ - yshift=-100,every node/.append style={ - yscale=.3,yslant=.3,xslant=-0.6},yscale=.3,yslant=.4,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] (s) to (\x\y); - \begin{pgfonlayer}{background} - \path[edge] (\x\y) to (t); - \end{pgfonlayer} - } - } - \end{tikzpicture} - \caption{HELLO} - \label{fig:image_flow} + \input{fig/big_network} \end{figure} +Figure \ref{fig:big_network} shows a visualization of how the final +network might look with all its edges. The source will have a lot of +outgoing edges, one for each pixel, while the sink has one incoming edge +for each pixel. The vertices corresponding to the pixels are only +connected to the source, the sink, and their neighboring pixels, so +their edge degree is much smaller than for $s$ and $t$. + diff --git a/fig/big_network.tex b/fig/big_network.tex new file mode 100644 index 0000000..bada059 --- /dev/null +++ b/fig/big_network.tex @@ -0,0 +1,41 @@ +\centering +\begin{tikzpicture}[scale=1.3] + + \node[vertex] (s) at (.6, 0) {s}; + \node[vertex] (t) at (.6, -6) {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 network, with edges between neighboring + pixels, and edges connecting the source and sink to the rest of the + network. The edge capacities are not visualized, and of course some + of them might be zero while others might be very high. +} +\label{fig:big_network} diff --git a/fig/norm_evolution.tex b/fig/norm_evolution.tex index 21f65a6..ca85d77 100644 --- a/fig/norm_evolution.tex +++ b/fig/norm_evolution.tex @@ -3,7 +3,7 @@ \centering \begin{tikzpicture} \draw[->] (0,0) -- (4,0) node[right] {$\lambda$}; - \draw[->] (0,-2) -- (0,2) node[above] {$E^x(0)$}; + \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}); @@ -29,7 +29,7 @@ \centering \begin{tikzpicture} \draw[->] (0,0) -- (4,0) node[right] {$\lambda$}; - \draw[->] (0,-2) -- (0,2) node[above] {$E^x(0)$}; + \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[ @@ -48,7 +48,7 @@ \caption{$L^2$ fidelity term.} \label{fig:l2_norm_evolution} \end{subfigure} -\caption{Two figures showing how the fidelity energy term $E^x(0)$ +\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 index 43f44d7..9cc18f4 100644 --- a/fig/norm_subgraph.tex +++ b/fig/norm_subgraph.tex @@ -6,9 +6,9 @@ \node[vertex] (u) at (0, -1) {$u^\lambda_x$}; \node[vertex] (t) at (0, -2) {t}; - \path[edge] (u) -- node[weight noslope] {$E^x(0)$} (t); + \path[edge] (u) -- node[weight noslope] {$E_\lambda^x(0)$} (t); \end{tikzpicture} - \caption{The network when $E^x(0) > 0$, with constant equal to + \caption{The network when $E_\lambda^x(0) > 0$, with constant equal to 0. } \label{fig:norm_subgraph_pos} @@ -21,15 +21,15 @@ \node[vertex] (u) at (0, -1) {$u^\lambda_x$}; \node[vertex] (t) at (0, -2) {t}; - \path[edge] (s) -- node[weight noslope] {$-E^x(0)$} (u); + \path[edge] (s) -- node[weight noslope] {$-E_\lambda^x(0)$} (u); \end{tikzpicture} - \caption{Network when $E^x(0) < 0$, with constant equal to - $-E^x(0)$. + \caption{Network when $E_\lambda^x(0) < 0$, with constant equal to + $-E_\lambda^x(0)$. } \label{fig:norm_subgraph_neg} \end{subfigure} \caption{ - The network construction for the fidelity term $E^x(u^\lambda_x)$. + The network 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. } diff --git a/flow.tex b/flow.tex index fffc835..567a868 100644 --- a/flow.tex +++ b/flow.tex @@ -14,8 +14,7 @@ We have already introduced capacities, and briefly mentioned the notion of flow as something limited by the capacity. In other words, flow is something we can send through our network, but the capacity limits how much we can send along each edge. It is useful to imagine a water supply -network with pipes of different sizes. The source and sink would then be -the water source, and water sink in the network. +network with pipes of different sizes and a water source and sink. Formally, we introduce the \emph{flow} as the function $f : V \times V \to [0, \infty)$. This function keeps count of how much flow we are @@ -78,10 +77,11 @@ we find that \label{eq:flow_leq_cut} \end{aligned} \end{equation} -for any $s$-$t$-cut $C = (S, T)$. The max-flow min-cut theorem, which -will be presented later, states that this inequality becomes an equality -when $f$ is a maximum flow, for some cut $C = (S, T)$, and this will -then be a minimum cut. +for any $s$-$t$-cut $C = (S, T)$. A very central result in network flow +theory called the max-flow min-cut theorem will be presented later. It +states that the inequality of \eqref{eq:flow_leq_cut} becomes an equality +when $f$ is a maximum flow for some cut $C = (S, T)$, and that all such +cuts are minimum cuts. %Recall that we want to find the minimum $s$-$t$ cut in our network. When %finding this minimum cut, we make use of an important duality theorem in @@ -96,18 +96,18 @@ find a maximum flow, and how do we find a minimum cut? \subsection{Augmenting path algorithms} The family of augmenting flow algorithms represent a popular approach to the maximum flow problem. The idea is simply to look for paths from the -source to the sink that can carry additional flow, a so-called -augmenting path, and then send the maximum possible amount of flow along -this path. When no such path exists anymore, no more flow can be sent -from the source to the sink, and a maximum flow has been reached. +source to the sink that can carry additional flow, so-called augmenting +paths, and then send the maximum possible amount of flow along such a +path. When no such path exists anymore, no more flow can be sent from +the source to the sink, and a maximum flow has been reached. \subsubsection{Residual network} -When further discussing approaches to solving the max-flow problem we -will need the notion of a residual network $G_f = (V_f, E_f, c_f)$, +When further discussing approaches to solving the maximum flow problem +we will need the notion of a residual network $G_f = (V_f, E_f, c_f)$, which is derived from the original network $G$ and contains the edges along which it is possible to send additional flow. This means that $E_f$ contains the edges $(u, v)$ from $E$ where $f(u,v) < c(u,v)$. But -that is not all, an important realisation is that it is also possible to +that is not all; an important realisation is that it is also possible to push flow \emph{back} along an edge where the flow is already positive. In other words, sending flow from $v$ to $u$ by cancelling some or all of the flow that is already going from $u$ to $v$. @@ -138,7 +138,7 @@ With the residual network 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 network $G = (V, E, c)$ with source $s$ and - sink $t$, then the following is equivalent: + sink $t$, then the following statements are equivalent: \begin{enumerate} \item $f$ is a maximum flow in $G$. \item The residual network $G_f$ contains no augmenting paths. @@ -147,8 +147,11 @@ max-flow min-cut theorem. \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. +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 +network and $T = V - S$. Figure \ref{fig:aug_flow} shows a simple network which already has five units flowing from $s$ to $t$. The marked path is a possible augmenting @@ -187,14 +190,13 @@ possible. \label{alg:ford_fulkerson} \end{algorithm} -Of course, a choice has to be made on how to search for a path in the -residual network. A common choice is to do a breadth-first search from -the source node until the sink node 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 -algorithm. +There are different ways to find augmenting paths, and a common choice +is to do a breadth-first search from the source node until the sink node +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. \subsection{Other algorithms} There are many different maximum flow algorithms that fall into the @@ -207,12 +209,12 @@ a variant of the augmenting path algorithm. It maintains a distance labeling $d(u)$ of the vertices $u \in V$ in the network, where $d(u)$ is the shortest path from the $s$ to $u$ in the residual network. This can be computed with a simple breadth-first search. The next step is to -construct a 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 has been found, the -distance labels are recalculated, and the label of the sink will be -increased by at least one. +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 @@ -220,6 +222,7 @@ 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 the Boykov and Kolmogorov, -it outperforms the other algorithms in experimental comparisons. +of this algorithm is not great, but according to Boykov and Kolmogorov, +it outperforms the other algorithms in experimental comparisons for this +specific application. diff --git a/introduction.tex b/introduction.tex index 275cb84..212d743 100644 --- a/introduction.tex +++ b/introduction.tex @@ -44,5 +44,3 @@ this project, but we will also see some image restoration results in the end, showing how the method performs for different kinds of noise, and different input parameters. -\newpage - diff --git a/main.tex b/main.tex index ba9cbdb..3915944 100644 --- a/main.tex +++ b/main.tex @@ -1,16 +1,16 @@ % Først spesifiserer vi hvilken dokumentklasse vi vil ha og noen % globale opsjoner. Bytt ut 'article' med 'book' hvis du vil ha % med kapitler. -\documentclass[a4paper, twoside, titlepage, 11pt]{article} +\documentclass[a4paper, titlepage, twoside, 11pt]{article} % Så sier vi fra om hvilke tilleggspakker vi trenger % til dokumentet vårt. De som du ikke trenger (se kommentaren) % kan det være en fordel å kommentere ut (sett prosenttegn foran), % da vil kompilering gå raskere. -%\usepackage[T1]{fontenc} % norsk tegnsett (æøå) +%\usepackage[T1]{fontenc} % norsk tegnsett (æøå) %\usepackage[utf8]{inputenc} % norsk tegnsett -\usepackage{geometry} % anbefalt pakke for å styre marger. +\usepackage{geometry} % anbefalt pakke for å styre marger. %\usepackage{amsfonts,amssymb} % matematikksymboler \usepackage{amsmath} @@ -27,6 +27,8 @@ % men har noen bugs, så den er kommentert % bort her. +\usepackage[binary-units=true]{siunitx} + \usepackage{tikz} \usetikzlibrary{matrix} \usetikzlibrary{arrows} @@ -41,7 +43,8 @@ \usepackage{listings} \lstset{ - basicstyle=\footnotesize, + %basicstyle=\footnotesize\ttfamily, + basicstyle=\footnotesize\ttfamily, keepspaces=true, tabsize=4, %inputencoding=utf8/latin1, @@ -57,7 +60,7 @@ \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,->,fill] +\tikzstyle{edge} = [draw,thick,->] \tikzstyle{dash} = [draw,dashed,-] \tikzstyle{bent edge} = [draw,bend left=80,thick,->] \tikzstyle{double edge} = [draw,thick,<->] @@ -94,6 +97,8 @@ % Inkluder forsida: \input{titlepage} +\cleardoublepage + % Romerske tall på alt før selve rapporten starter er pent. \pagenumbering{roman} @@ -119,6 +124,7 @@ % LaTeX om selv! \input{introduction} +\clearpage \input{methods} \input{total} \input{cut} diff --git a/methods.tex b/methods.tex index 5d959e4..62a7d87 100644 --- a/methods.tex +++ b/methods.tex @@ -17,7 +17,7 @@ blur amounts to setting each pixel value as the weighted average of its neighbors in the original image. This will, in addition to smoothing out possible noise, add blur and remove details from the image. -The Gaussian happens to be the fundamental solution of the heat +The Gaussian function happens to be the fundamental solution of the heat equation $\partial_t u = \Delta u$. Convolving it with the original image $v$ is therefore equivalent to solving the heat equation with $v$ as initial value, for some time $t > 0$ depending on $\sigma$. Care must @@ -30,18 +30,18 @@ example \cite{weickert1998anisotropic} for more information. \subsection{Anisotropic diffusion} Since the Gaussian filter will blur out both noise and details of the -image, it would be nice to somehow reduce the amount of blur that is -added in areas of the image with edges, or what we believe to be edges. -This can be done by controlling the thermal diffusivity $\alpha(u)$ of -the image in the heat equation +image, it would be nice to somehow reduce the amount of blurring done +around edges, or what we believe to be edges. This can be done by +controlling the thermal diffusivity $\alpha(u)$ of the image in the heat +equation \begin{equation} \begin{cases} - \partial_t u &= \mathrm{div} \big( a(u) \nabla u\big) \\ - u|_{t=0} &= v. + \partial_t u = \mathrm{div} \big( \alpha(u) \nabla u\big) \\ + u|_{t=0} = v. \end{cases} \end{equation} The problem is then, how to detect the edges such that we can reduce -$\alpha(u)$ in those areas. +$\alpha(u)$ in those parts of the image. Perona and Malik proposed a modified heat equation \begin{equation} @@ -75,20 +75,19 @@ edges, and at the same time a large diffusivity \emph{along} the edges. See for example \cite{weickert1998anisotropic} for more information on these diffusion tensor methods. -\subsection{Non-local means} +\subsection{Other filters} The non-local means method of image restoration takes into account the -fact that two separate parts of an image might be very similar, for -example in the case of a regular texture. One defines some neighborhood -for each pixel, and when restoring a pixel, one considers \emph{all} -other pixels whose neighborhoods are similar to the neighborhood of the -current pixel, for some notion of similarity. One then averages all -these pixels to obtain the new pixel value. +fact that two spatially separated parts of an image might be very +similar, for example in the case of a regular texture. As in many other +methods, one defines some neighborhood for each pixel, but when +restoring a pixel, one considers all other pixels whose neighborhoods +are similar to the neighborhood of the current pixel, for some notion of +similarity. The new pixel value is then obtained by taking the average +of all these similar pixels. -\subsection{Other filters} A different filter much used in real-world image processing is the median filter. It replaces the value of each pixel with the median of -the values of the surrounding pixels. As with other methods, one has to -choose what the neighborhood of each pixel should be. It can perform -very well, especially with salt and pepper noise or images where some +the values of the neighboring pixels. It can perform very well in +practice, especially with salt and pepper noise or images where some pixels are missing. diff --git a/pushrelabel.tex b/pushrelabel.tex index 537b4cd..1f89abf 100644 --- a/pushrelabel.tex +++ b/pushrelabel.tex @@ -4,11 +4,11 @@ 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 network at all times, but still obtains a valid maximum flow -when the algorithm is finished. +when the algorithm terminates. \subsection{Preflow} Instead of maintaining a valid flow, we introduce the concept of a -\emph{preflow} as we relax the flow conservation constraint from +\emph{preflow} by relaxing the flow conservation constraint from earlier. For a preflow $f$, the capacity constraint still holds so the flow must always be less than the capacity, but we allow positive excess in the vertices. The flow conservation constraint from before then @@ -35,17 +35,20 @@ $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 of the network where connected vertices can not have a large +vertices in the network 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 network 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. -The height map $d : V \to \mathbb{N}$ is also often called a distance -labeling and it satisfies $d(t) = 0$, and for every edge $(u, v)$ in the -residual network, i.e.\ every edge with $c_f(u, v) > 0$, we require that -$d(u) \leq d(v) + 1$. For every vertex $u$ for which there is a path +The height map $d : V \to \mathbb{N}$ satisfies $d(t) = 0$, and for +every edge $(u, v)$ in the residual network, 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. + +For every vertex $u$ for which there is a path through the residual network to the sink, $d(u)$ will be a lower bound on the length of such a path. @@ -62,7 +65,7 @@ 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$ such that $d(u) = d(v) + 1$, i.e.\ to a vertex with a +$(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. @@ -77,16 +80,13 @@ admissible edge $(u, v)$. %increases for $v$, remains non-negative for $u$ and remains the same for %all other vertices. -The distance labels are not changed during the push procedure, however, -the residual network can change. The edge $(u, v)$ might disappear, and -an edge $(v, u)$ will surely appear if it does not already exist. Assume -that $d$ is a valid labeling. When starting the push procedure we have -$d(u) = d(v) + 1$, so we also have the following -\begin{align} - d(u) &\leq d(v) + 1 \\ - d(v) &\leq d(u) + 1, -\end{align} -and $d$ remains a valid labeling, even if the edge $(v, u)$ appears. +The residual network 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} @@ -112,7 +112,7 @@ and $d$ remains a valid labeling, even if the edge $(v, u)$ appears. \subsubsection{The relabel procedure} The relabel procedure is our tool for changing the distance labels of -our vertices. It changes the label of a vertex to the greatest possible +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 rather mathematical pseudocode implementation. @@ -149,8 +149,8 @@ network. 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 run the push and -relabel procedures when applicable until there are no more active +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 @@ -178,18 +178,18 @@ is relabeled. \label{alg:discharge} \end{algorithm} -In the second phase of the algorithm this preflow is turned into a -maximum flow by sending excess flow which did not reach the sink, from +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 network 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 network $G = (V, E, c)$, assume that the first phase of - the push-relabel algorithm has terminated and no more active - vertices remain. Then there exists a $k \in \mathbb{N} \cap \left(0, - N \right)$ such that there is no vertex with label $k$, and 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$. + 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} @@ -198,19 +198,19 @@ first phase is finished, and the following theorem allows us to do that. $N-1$ labels in $\{1, \ldots, N-1\}$, so there must exist an $k$ as described. - These can be no edge $(u, v) \in E_f$ such that $u \in S$ and $v \in - T$, as this would imply $d(v) \leq d(u) + 2$ which contradicts the - labeling constraint. 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$ have no flow. This means that the capacity of the cut is - equal to the flow through the cut, i.e.\ $c(S, T) = f(S, T)$. + 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 from $s$ to $t$. + equal to the flow. \end{proof} We will see later that with the gap relabeling heuristic, there will @@ -230,10 +230,10 @@ 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 (and then in the -opposite direction). Since the number of relabelings is bounded by -$\abs{V}$, and we have $\abs{E}$ edges, this gives the stated number of -saturating pushes. +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 @@ -243,48 +243,50 @@ idea is to define \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 the $u$ -inactive (but might activate $v$). The total increase in $\phi$ due 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 + c \abs{V} \underbrace{\abs{V} + \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 just stated that the discharge procedure is run on -active nodes until there are no more active nodes left. The choice of -the order in which to discharge these active nodes remain, and multiple +Until now we have stated that the discharge procedure is run on active +nodes until there are no more active nodes left. The choice of the order +in which to discharge these active vertices remain, and multiple possibilities exist. -The 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. +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 +network. -The highest level selection rule always discharges the vertex with the -largest distance label by always keeping track of which active vertices -has what 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}. +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 how the input network looks. +it, varies a lot with the structure of the input network. -\fixme{No reference to ON IMPLEMENTING THE PUSH-RELABEL?} +Cherkassky and Goldberg describe the algorithm along with different +selection rules, heuristics and their implementation in +\cite{cherkassky1997implementing}. \subsection{Heuristics} -Different heuristics exists that can speed up the algorithm +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 networks. The most used heuristics are the gap and global relabeling heuristics, both aiming to @@ -355,12 +357,12 @@ procedure will not change the validity of these two things. When running the push-relabel algorithm with the gap heuristic, we can be sure that there will never be a node $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$. Using the same reasoning as in Theorem -\ref{thm:cut_identification}, if there was a vertex with label $N-1$, -there would only be $N-3$ possibly having labels in $[1, 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$. +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 @@ -437,12 +439,12 @@ we reset the flow and labels. \subsubsection{Output image construction} -We mentioned already in section \ref{sec:total_energy} that in order to +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}. We will now show that through reuse of the -distance labels from the last iteration of the push-relabel algorithm, -we can guarantee that it is possible to stack the thresholded images. +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. \fixme{More precise statement than ``stack on top of each other''} @@ -496,6 +498,21 @@ The idea is now that we have halved the possible $\lambda$ interval for 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 network 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 +network 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. diff --git a/results.tex b/results.tex index af86901..27db990 100644 --- a/results.tex +++ b/results.tex @@ -159,14 +159,14 @@ the noisy images from Figure \ref{fig:noisy_lena} restored for different parameters of $\beta$ and the exponent $p$ in the fidelity term. First we note that the method works quite well on smooth surfaces, like Lena's shoulder and the background of the image, while her hair, and the -feather hanging from her hat has lost details, especially for higher +feather boa hanging from her hat has lost details, especially for higher values of $\beta$. We see in the method noise, especially in Figure -\ref{fig:laplace_restored_1}, that the feather in Lena's hat has lost -details. Being especially bright, the method noise here tells us that -the algorithm has corrected this area more than the rest, all the while -the noise was spread evenly over the whole image. +\ref{fig:laplace_restored_1}, that the feather boa in Lena's hat has +lost details. Being especially bright, the method noise here tells us +that the algorithm has corrected this area more than the rest, all the +while the noise was spread evenly over the whole image. As mentioned earlier, one of the theoretically nice properties of the total variation method is that it preserves edges. Looking at the @@ -181,10 +181,11 @@ again. \caption{ Performance comparison of the implemented algorithms, and the algorithm of Darbon and Sigelle \cite{darbon2006image} using - sections of the Lena test image of increasing size. The - Push-Relabel algorithm is run both using the highest level - selection rule (HL) and the first-in first-out selection rule - (FIFO). All numbers are in seconds. + sections of the Lena test image of increasing size, with the + $L^2$ norm in the fidelity term, $\beta = 10$ and a size four + neighborhood. The Push-Relabel algorithm is run both using the + highest level selection rule (HL) and the first-in first-out + selection rule (FIFO). All numbers are in seconds. } \begin{tabular}{lrrrr} Method & $64\times 64$ & $128\times 128$ & $256\times 256$ & $512\times @@ -204,5 +205,8 @@ Sigelle at \cite{darbonsigellebinary}, implemented following the description in \cite{darbon2006image}, has also been timed for comparison. They use the specialized algorithm of Boykov and Kolmogorov \cite{boykov2004experimental}, and additionally a divide and conquer -approach which could explain the superior performance. +approach which could explain the superior performance. All performance +experiments were performed on a Lenovo Thinkpad X230 Laptop with an +Intel\textregistered\ Core\texttrademark\ i5-3320M CPU at +\SI{2.60}{\giga\hertz} and \SI{8}{\giga\byte} of memory. diff --git a/titlepage.tex b/titlepage.tex index 8de7ea6..2801c39 100644 --- a/titlepage.tex +++ b/titlepage.tex @@ -11,7 +11,42 @@ \today \\[2cm] -\includegraphics[width=0.5\textwidth]{lena512.png} +%\includegraphics[width=0.5\textwidth]{lena512.png} +\centering +\begin{tikzpicture}[scale=1.3] + + \node[vertex] (s) at (.6, 0) {s}; + \node[vertex] (t) at (.6, -6) {t}; + + \begin{scope}[ + yshift=-100,every node/.append style={ + yscale=.3,yslant=.3,xslant=-0.6},yscale=.3,yslant=.3,xslant=-0.6 + ] + \node[scale=.3] at (1.5,1.5) {\includegraphics{lena512.png}}; + \foreach \x in {0,...,3} { + \foreach \y in {0,...,3} { + \node[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} \end{titlepage} diff --git a/total.tex b/total.tex index 75c8f94..acda826 100644 --- a/total.tex +++ b/total.tex @@ -13,10 +13,10 @@ $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. -How we can integrate the gradient of an image which can contain -discontinuities will be clarified later, but first we will look at how -this model can be grounded in spatial statistics; a field which forms -an important part of image analysis. +How we compute and integrate the gradient the image will be clarified +later, but first we will look at how this model can be grounded in +spatial statistics; a field which forms an important part of image +analysis. \subsection{Probabilistic background} \label{sec:prob} @@ -24,9 +24,8 @@ an important part of image analysis. Image processing is pushing forward in many fields at once, analysis, spatial statistics and discrete optimization to name a few. Methods can often be related across fields, and in this short section we will see -how the total variation method can be related to one instance of the -so-called maximum a posteriori estimation (MAP for short) of spatial -statistics. +how the total variation method can be related to the maximum a +posteriori estimation (MAP) approach in spatial statistics. There are many ways to model noise in an image, and different models take different physical phenomena into account. We will in this brief @@ -38,8 +37,8 @@ noise}, which given an image $u$ results in a captured image where $\epsilon$ represents the noise that can be modeled by different probability distributions. -For a captured noisy image $v$ and we want to recover $u$, or something -as close to $u$ as possible, so we introduce the likelihood function +Given a captured noisy image $v$ we want to recover $u$, so we introduce +the likelihood \begin{equation} p(u \mid v) = \frac{p(v \mid u) \cdot p(u)}{p(v)} \end{equation} @@ -55,12 +54,12 @@ increasing function, the maximization problem is equivalent to \min_u \big( - \ln p(v \mid u) - \ln p(u) \big). \label{eq:minimize_neg_log} \end{equation} -For discrete images consisting of pixels, we use $x$ and $y$ to denote -pixel positions such that $u_x$ and $u_y$ are the values of those -pixels, \emph{not to be confused with} partial derivatives $\partial_x -u$ and $\partial_y u$ which will not be used. If the noise term $p(v -\mid u)$ is modeled to be independently and identically distributed -Gaussian noise on each pixel, we have +For discrete images consisting of pixels, we use $x, y \in \mathbb{R}^2$ +to denote pixel positions such that $u_x$ and $u_y$ are the values of +those pixels, \emph{not to be confused with} partial derivatives +$\partial_x u$ and $\partial_y u$ which will not be used. If the noise +term $p(v \mid u)$ is modeled to be independently and identically +distributed Gaussian noise in each pixel, we have \begin{equation} p(v \mid u) \propto \prod_x \exp \left(-\frac{(u_x - v_x)^2}{2\sigma^2} \right). @@ -70,8 +69,8 @@ Gaussian noise on each pixel, we have The prior probability $p(u)$ depends on what kinds of images are considered, and how they are obtained. In the theory of Markov random fields, one lets the probability of a pixel value depend only on some -neighborhood $\mathcal{N}$. A common choice for the relation between -neighboring pixels is the Laplace distribution +neighborhood $\mathcal{N}$ of the pixel. A common choice for the +relation between neighboring pixels is the Laplace distribution \begin{equation} p(u) \propto \prod_x \prod_{y \in \mathcal{N}(x)} \exp \left( @@ -109,9 +108,9 @@ We see that the two discrete problems \eqref{eq:minimizing_probability_l2} and \eqref{eq:minimizing_probability_l1} resemble the original formulation in \eqref{eq:first_min_presentation}. Before presenting the -discretization, we will look at a more rigorous presentation of -\eqref{eq:first_min_presentation}, and the definitions, spaces and -theorems needed to treat it. +discretization of \eqref{eq:first_min_presentation}, we will look at a +more rigorous presentation of \eqref{eq:first_min_presentation}, and the +definitions, spaces and theorems needed to treat it. \subsection{Continuous formulation} The method of total variation was initially introduced to image @@ -127,7 +126,8 @@ minimizing the energy function + \beta\; \TV(u) \label{eq:energy_function} \end{equation} -where $v$ is the original image, and $\Omega$ is the image domain, in +over $L^p(\Omega)$. +Here $v$ is the original image, and $\Omega$ is the image domain, in our case an open subset of $\mathbb{R}^2$. Normally, $p$ is chosen to be 1 or 2, the two choices that correspond to the two different noise models in \eqref{eq:minimizing_probability_l1} and @@ -140,9 +140,9 @@ problem. The term $\TV(u)$ in \eqref{eq:energy_function} is the total variation of the image, introduced earlier as $\int_\Omega \abs{\nabla u}$. Since -it is problematic to calculate the traditional gradient, the following -definition is introduced, using Green's identity to transfer the -derivative to a smooth test function. +it can be problematic to calculate the traditional gradient, the +following definition is introduced, using Green's identity to transfer +the derivative to a smooth test function. \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 @@ -179,10 +179,10 @@ space of functions with bounded variation. %range of the level values $\lambda$. The coarea formula and further %references can be found in \cite{caselles2011total}. -Before stating the coarea formula, we need the formal definition of a -set perimeter. Here, $\idfun_E$ signifies the characteristic function of -the set $E$, the function which is equal to one in every point in $E$, -and zero elsewhere. +Before introducing the coarea formula, we need the formal definition of +a set perimeter. Here, $\idfun_E$ signifies the characteristic function +of the set $E$, the function which is equal to one in every point in +$E$, and zero elsewhere. \begin{definition}[Set perimeter] The perimeter of a set $E$ is defined as \begin{equation} @@ -191,14 +191,14 @@ and zero elsewhere. A measurable set $E \subset \Omega$ is of finite perimeter in $\Omega$ if $\idfun_E \in \BV(\Omega)$. \end{definition} -The set perimeter behaves mostly as one might expect, just note that if -the boundary of $E$ overlaps with the boundary of $\Omega$, then -$P(E,\Omega)$ will not include the overlapping part of the boundary. - -\fixme{Meh. Level sets.} - -For an image $u$ we denote by $u \leq \lambda$ the set $\{ x \in \Omega -: u_x \leq \lambda\}$ and the thresholded images +The set perimeter behaves mostly as one might expect for ``nice'' sets, +just note that if the boundary of $E$ overlaps with the boundary of +$\Omega$, then $P(E,\Omega)$ will not include the overlapping parts of +the boundary. + +For an image $u$ and a level $\lambda$ we denote the \emph{level set} by +$\{u \leq \lambda\}$ defined as the set $\{ x \in \Omega : u_x \leq +\lambda\}$. The thresholded images are defined as \begin{equation} u^\lambda = \idfun_{u \leq \lambda}. \label{eq:level_set} @@ -208,21 +208,22 @@ coarea formula, which will let us write the total variation as an integral over the range of the image levels $\lambda$. \begin{theorem}[The coarea formula] Let $u \in BV(\Omega)$. Then for almost every $\lambda$ the set - $u \leq \lambda$ is of finite perimeter, and one has the coarea + $\{u \leq \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 \leq \lambda,\Omega) \, d\lambda + = \int_\mathbb{R} P(\{u \leq \lambda\},\Omega) \, d\lambda \label{eq:coarea_formula} \end{equation} \end{theorem} -Assuming that our image $u$ is in $\BV(\Omega)$, this allows us to +Note that the equality in \eqref{eq:coarea_formula} holds for all $u \in +L^1(\Omega)$, not only when the total variation is finite. We can now rewrite \eqref{eq:energy_function} as \begin{equation} E_v(u) = \int_\Omega \abs{u - v}^p - + \beta \int_\mathbb{R} P(u \leq \lambda,\Omega) \, d\lambda. + + \beta \int_\mathbb{R} P(\{u \leq \lambda\},\Omega) \, d\lambda. \label{eq:energy_perimeter} \end{equation} This will be the starting point for our discretization, and it will @@ -357,7 +358,7 @@ setting, such that the total variation can be written as a sum over the level values \begin{equation} \TV(u) - = \sum_{\lambda = 0}^{L-2} P(u \leq \lambda, \Omega) + = \sum_{\lambda = 0}^{L-2} P(\{u \leq \lambda\}, \Omega) \approx \sum_{\lambda = 0}^{L-2} \sum_{(x,y)} w_{xy} \abs{u_x^\lambda - u_y^\lambda}. \label{eq:tv_discrete_int} @@ -368,11 +369,11 @@ neighborhood relation. The $w_{xy}$ is a weight parameter, and the sum over $\lambda$ ends at $L-2$ since $u^{L-1}$ is equal to 1 in every pixel of the image. -Intuitively, the perimeter of the level set $u^\lambda$ is proportional -to the number of pixels at the boundary of the set. This is again -proportional to the number of neighborhood relations crossing the +Intuitively, the perimeter of the level set $\{u \leq \lambda\}$ is +proportional to the number of pixels at the boundary of the set. This is +again proportional to the number of neighborhood relations crossing the boundary, and $\abs{u^\lambda_x - u^\lambda_y}$ only contributes to the -sum when one node is in $u^\lambda$ and the other is not. +sum when one pixel is in $\{u \leq \lambda\}$ and the other is not. Boykov and Kolmogorov argue in \cite{boykov2003computing} that if the weight is chosen as @@ -387,7 +388,7 @@ 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_{xy}$, and $d_{xy}$ go to -zero. +zero. \fixme{BLOOP} \subsubsection{Fidelity term} Now that the total variation is discretized, we need to take care of the @@ -405,7 +406,7 @@ This allows us to discretize $\int_\Omega \abs{u - v}^p$ of \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}$ +function $N(k)$ taking values $k \in \mathcal{L}$: \begin{equation} \begin{aligned} N(k) &= \sum_{\lambda=0}^{k-1} \big( @@ -413,7 +414,7 @@ function $N(k)$ taking values $k \in \mathcal{L}$ \big) + N(0) \\ &= \sum_{\lambda=0}^{L-2} \big( N(\lambda + 1) - N(\lambda) - \big) \mathbbm{1}_{\lambda < k} + N(0) + \big) \mathbbm{1}_{\lambda < k} + N(0). \end{aligned} \end{equation} Since $\mathbbm{1}_{\lambda < u_x} = (1 - u^\lambda_x)$ we rewrite @@ -425,8 +426,9 @@ Since $\mathbbm{1}_{\lambda < u_x} = (1 - u^\lambda_x)$ we rewrite N_x(\lambda + 1) - N_x(\lambda) \big) - (1 - u^\lambda_x) + N_x(0) + (1 - u^\lambda_x) + N_x(0). \end{equation} +\fixme{BLOOP} \subsubsection{Total energy} \label{sec:total_energy} @@ -439,7 +441,7 @@ are left with \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) + =: \sum_{\lambda=0}^{L-2} F_\lambda(u^\lambda) \label{eq:total_energy} \end{equation} where @@ -452,21 +454,18 @@ where (1 - u^\lambda_x) \label{eq:fidelity_energy} \\ E^{x,y}(u^\lambda_x, u^\lambda_y) &= - w_{xy} - \left( - \left(1 - 2 u_y^\lambda \right) u_x^\lambda + u_y^\lambda - \right). + w_{xy} \abs{u_x^\lambda - u_y^\lambda}. \label{eq:neigh_energy} \end{align} If we minimize each level separately it is obvious that we also minimize -the total energy. In other words if we for every $\lambda$ find a level -set $u^\lambda$ that minimizes $F_\lambda(u^\lambda)$, the sum in -\eqref{eq:total_energy} will be minimized. But the question remains if -these obtained level sets $u^\lambda$ can be combined to produce an -output image. The level sets were defined as $u^\lambda = \mathbbm{1}_{u -\leq \lambda}$, so we need them to be monotonically increasing in -increasing level values, i.e.\, +the total energy. In other words if we for every $\lambda$ find a +thresholded image $u^\lambda$ that minimizes $F_\lambda(u^\lambda)$, the +sum in \eqref{eq:total_energy} will be minimized. But the question +remains if these obtained $u^\lambda$ can be combined to produce an +output image. They were defined as $u^\lambda = \mathbbm{1}_{u \leq +\lambda}$, so we need them to be monotonically increasing in increasing +level values, i.e.\, \begin{equation} u^\lambda_x \leq u^\mu_x \quad \forall \lambda \leq \mu, \quad \forall x \in S. -- 2.47.3