}\r
\r
@incollection{caselles2011total,\r
- title={Total variation in imaging},\r
- author={Caselles, Vicent and Chambolle, Antonin and Novaga, Matteo},\r
- journal={Handbook of Mathematical Methods in Imaging},\r
- pages={1016--1057},\r
year={2011},\r
- publisher={Springer}\r
+ isbn={978-0-387-92919-4},\r
+ booktitle={Handbook of Mathematical Methods in Imaging},\r
+ editor={Scherzer, Otmar},\r
+ doi={10.1007/978-0-387-92920-0_23},\r
+ title={Total Variation in Imaging},\r
+ url={http://dx.doi.org/10.1007/978-0-387-92920-0_23},\r
+ publisher={Springer New York},\r
+ author={Caselles, V. and Chambolle, A. and Novaga, M.},\r
+ pages={1016-1057}\r
}\r
\r
@incollection{chan2011numerical,\r
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
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.
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
\{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}
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)
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}
\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.
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}
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.
\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}
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}
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.
}
\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$.
+
--- /dev/null
+\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}
\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});
\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[
\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}
\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}
\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.
}
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
\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
\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$.
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.
\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
\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
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
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.
end, showing how the method performs for different kinds of noise, and
different input parameters.
-\newpage
-
% 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}
% men har noen bugs, så den er kommentert
% bort her.
+\usepackage[binary-units=true]{siunitx}
+
\usepackage{tikz}
\usetikzlibrary{matrix}
\usetikzlibrary{arrows}
\usepackage{listings}
\lstset{
- basicstyle=\footnotesize,
+ %basicstyle=\footnotesize\ttfamily,
+ basicstyle=\footnotesize\ttfamily,
keepspaces=true,
tabsize=4,
%inputencoding=utf8/latin1,
\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,<->]
% Inkluder forsida:
\input{titlepage}
+\cleardoublepage
+
% Romerske tall på alt før selve rapporten starter er pent.
\pagenumbering{roman}
% LaTeX om selv!
\input{introduction}
+\clearpage
\input{methods}
\input{total}
\input{cut}
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
\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}
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.
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
$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.
\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.
%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}
\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.
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
\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}
$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
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
\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
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
\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''}
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.
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
\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
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.
\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}
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}
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
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}
\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).
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(
\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
+ \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
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
%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}
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}
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
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}
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
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
\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(
\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
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}
\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
(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.
projects / prosjektoppgave / commitdiff