From fe2422332603977755e953cbd2a1d0ce082da58e Mon Sep 17 00:00:00 2001 From: =?utf8?q?Bj=C3=B8rn=20Rustad?= Date: Sun, 8 Feb 2015 22:20:40 +0100 Subject: [PATCH] FINAL --- appendix.tex | 71 ++++++- discrete.tex | 2 + graph.cpp | 529 +++++++++++++++++++++++++++++++++++++++++++++++++++ graph.hpp | 104 ++++++++++ main.tex | 2 +- maxflow.tex | 37 +--- 6 files changed, 709 insertions(+), 36 deletions(-) create mode 100644 graph.cpp create mode 100644 graph.hpp diff --git a/appendix.tex b/appendix.tex index a7381cc..654afd1 100644 --- a/appendix.tex +++ b/appendix.tex @@ -1,6 +1,73 @@ \chapter{\cpp{} implementation} +\label{app:impl} -\fixme{THIS SHOULD NOT BE INCLUDED ANYMORE!} +A \cpp{} implementation of the method described is included here and can +also be found online at \cite{github}. It uses +the open computer vision library OpenCV \cite{opencv_library} to load +and save image files and contains compilation and usage instructions. +The implementation has been tested on an installation of the Ubuntu +Linux distribution, but it should in theory be portable to other +platforms supported by OpenCV. A rudimentary graphics interface has also +been made, to make it easier to play with the parameters of the +algorithm. -\inputminted{c++}{../image-restoration/graph.cpp} +Both the push-relabel and the Boykov--Kolmogorov algorithms have been +implemented. Although an effort has been made to improve the performance +of both implementations, they are not ment to beat the fastest. The +focus has rather been on clarity and understanding. + +Note that when implementing maximum flow algorithms it is not a good +idea, memory- and performance-wise, to actually construct the residual +graph $G_f$. Instead, every time we update the flow $f(u,v)$ we set +the flow in the opposite direction to its negative value $f(v,u) = +-f(u,v)$. Then we can at any time, consider the value $c(u,v) - f(u,v)$ +in the place of the residual capacity $c_f(u,v)$. + +For the gap relabeling heuristic of the push-relabel algorithm, we need +to have a easy way of finding when a gap occurs. This is done by keeping +track of how many vertices exist with each label. + +When the capacities have been updated in the Boykov--Kolmogorov +algorithm, flow is sent along all two-edge paths such that they do not +have to be considered by the main loop of the algorithm. + +\usemintedstyle{borland} + +\section{main.cpp} +Main function of the command line executable. Here we read and parse command +line parameters. +\inputminted[linenos,fontsize=\fontsize{9}{9}]{cpp}{/home/burk/dev/image-restoration/main.cpp} + +\section{image.hpp} +Calculates and sets up graph weights based on the input image and +parameters. +\inputminted[linenos,fontsize=\fontsize{9}{9}]{cpp}{/home/burk/dev/image-restoration/image.hpp} + +\section{image.cpp} +\inputminted[linenos,fontsize=\fontsize{9}{9}]{cpp}{/home/burk/dev/image-restoration/image.cpp} + +\section{anisotropy.hpp} +Construction of anisotropy tensor. +\inputminted[linenos,fontsize=\fontsize{9}{9}]{cpp}{/home/burk/dev/image-restoration/anisotropy.hpp} + +\section{anisotropy.cpp} +\inputminted[linenos,fontsize=\fontsize{9}{9}]{cpp}{/home/burk/dev/image-restoration/anisotropy.cpp} + +\section{graph.hpp} +Graph class with the maximum flow algorithms. +\inputminted[linenos,fontsize=\fontsize{9}{9}]{cpp}{/home/burk/dev/master/graph.hpp} + +\section{graph.cpp} +\inputminted[linenos,fontsize=\fontsize{9}{9}]{cpp}{/home/burk/dev/master/graph.cpp} + +\section{selectionrule.hpp} +Highest level and FIFO selection rules for the push-relabel algorithm. +\inputminted[linenos,fontsize=\fontsize{9}{9}]{cpp}{/home/burk/dev/image-restoration/selectionrule.hpp} + +\section{selectionrule.cpp} +\inputminted[linenos,fontsize=\fontsize{9}{9}]{cpp}{/home/burk/dev/image-restoration/selectionrule.cpp} + +\section{neighborhood.hpp} +Neighborhood class that also calculates angular differences. +\inputminted[linenos,fontsize=\fontsize{9}{9}]{cpp}{/home/burk/dev/image-restoration/neighborhood.hpp} diff --git a/discrete.tex b/discrete.tex index 67f0a82..e69d7e3 100644 --- a/discrete.tex +++ b/discrete.tex @@ -455,6 +455,8 @@ In this section we will look at how these graphs are constructed such that their minimum cuts correspond to the minimizers of the functional $F^\lambda$. The description is taken with some small adjustments from my project work \cite{project}, and is included here for completeness. +An implementation of the described approach can be found in +Appendix~\ref{app:impl}. \subsection{Graphs} diff --git a/graph.cpp b/graph.cpp new file mode 100644 index 0000000..e4443bc --- /dev/null +++ b/graph.cpp @@ -0,0 +1,529 @@ +#include +#include +#include +#include +#include +#include +#include + +#include "graph.hpp" + +using namespace std; + +/* Add an edge from one vertex to another. */ +void FlowGraph::addEdge(int from, int to, int cap) { + G[from].e.push_back(Edge(from, to, cap, G[to].e.size())); + if (from == to) G[from].e.back().index++; + int index = G[from].e.size() - 1; + G[to].e.push_back(Edge(to, from, 0, index)); + + if (from == source) + G[to].si = index; + + if (to == sink) + G[from].ti = index; +} + +/* + * Add an edge and at the same time an antiparallel edge + * with the same capacity. + */ +void FlowGraph::addDoubleEdge(int from, int to, int cap) { + G[from].e.push_back(Edge(from, to, cap, G[to].e.size())); + G[to].e.push_back(Edge(to, from, cap, G[from].e.size() - 1)); +} + +#ifdef PUSH_RELABEL + +/* + * Change the capacity of an edge. Need the from-vertex and + * the index of the edge in its edge list (returned from addEdge. + */ +void FlowGraph::changeCapacity(int from, int to, int cap) { + int index; + if (from == source) { + index = G[to].si; + } else if (to == sink) { + index = G[from].ti; + } else { + exit(1); + } + + int diff = G[from].e[index].flow - cap; + + G[from].e[index].cap = cap; + + /* Check if we need to reduce the flow. */ + if (diff > 0) { + G[from].excess += diff; + G[to].excess -= diff; + G[from].e[index].flow = cap; + G[to].e[G[from].e[index].index].flow = -cap; + rule.add(from, G[from].height, G[from].excess); + } +} + +#else + +/* + * Change the capacity of an edge. Need the from-vertex and + * the index of the edge in its edge list (returned from addEdge) + */ +void FlowGraph::changeCapacity(int from, int to, int cap) { + int index; + if (from == source) { + index = G[to].si; + } else if (to == sink) { + index = G[from].ti; + } else { + exit(1); + } + + /* Nodes in the S set can not send any more flow anyways. */ + if (from == source && G[to].c == SOURCE) + return; + + G[from].e[index].cap = cap; + + if (from != source) + return; + + /* Push flow along two-edged path. */ + int si = G[to].si; + int ti = G[to].ti; + + Edge *sv, *vt; + sv = &G[source].e[si]; + vt = &G[to].e[ti]; + + int rs = sv->cap - sv->flow; + int rt = vt->cap - vt->flow; + + if (rs > 0 && rt > 0) { + int m = min(rs, rt); + push(*sv, m); + push(*vt, m); + + if (m == rs && G[to].p == sv) { + G[to].p = NULL; + orphans.push(to); + } + + if (m == rt && G[to].p == vt) { + G[to].p = NULL; + orphans.push(to); + } + } + + /* Activate vertices if the new edge was not saturated. */ + if (G[from].e[si].flow != G[from].e[si].cap) { + if (!G[from].active) { + bkq.push(from); + G[from].active = true; + } + if (!G[to].active) { + bkq.push(to); + G[to].active = true; + } + } +} + +#endif + +/* Reset all flow and excess. */ +void FlowGraph::resetFlow() { + for (size_t i = 0; i < G.size(); ++i) { + for (size_t j = 0; j < G[i].e.size(); ++j) { + G[i].e[j].flow = 0; + } + G[i].excess = 0; + } +} + +/* Reset all distance labels. */ +void FlowGraph::resetHeights() { + for (size_t i = 0; i < G.size(); ++i) { + G[i].height = 0; + } + fill(count.begin(), count.end(), 0); +} + +/* Push along an edge. */ +void FlowGraph::push(Edge &e) { + int flow = min(e.cap - e.flow, G[e.from].excess); + G[e.from].excess -= flow; + G[e.to].excess += flow; + e.flow += flow; + G[e.to].e[e.index].flow -= flow; + + rule.add(e.to, G[e.to].height, G[e.to].excess); +} + +/* Push given flow along an edge. */ +void FlowGraph::push(Edge &e, int f) { + e.flow += f; + G[e.to].e[e.index].flow -= f; +} + +/* Relabel a vertex. */ +void FlowGraph::relabel(int u) { + count[G[u].height]--; + G[u].height = 2*N; + + for (size_t i = 0; i < G[u].e.size(); ++i) { + if (G[u].e[i].cap > G[u].e[i].flow) { + G[u].height = min(G[u].height, G[G[u].e[i].to].height + 1); + } + } + + if (G[u].height >= N) { + G[u].height = N; + } + else { + count[G[u].height]++; + rule.add(u, G[u].height, G[u].excess); + } +} + +/* Relabel all vertices over the gap h to label N. */ +void FlowGraph::gap(int h) { + for (size_t i = 0; i < G.size(); ++i) { + if (G[i].height < h) continue; + if (G[i].height >= N) continue; + + rule.deactivate(i); + + count[G[i].height]--; + G[i].height = N; + } + + rule.gap(h); +} + +/* Discharge a vertex. */ +void FlowGraph::discharge(int u) { + size_t i; + for (i = 0; i < G[u].e.size() && G[u].excess > 0; ++i) { + if (G[u].e[i].cap > G[u].e[i].flow + && G[u].height == G[G[u].e[i].to].height + 1) { + push(G[u].e[i]); + } + } + + if (G[u].excess > 0) { + /* Check if a gap will appear. */ + if (count[G[u].height] == 1) + gap(G[u].height); + else + relabel(u); + } +} + +/* Run the push-relabel algorithm to find the min-cut. */ +void FlowGraph::minCutPushRelabel(int source, int sink) { + G[source].height = N; + + rule.activate(source); + rule.activate(sink); + + for (size_t i = 0; i < G[source].e.size(); ++i) { + G[source].excess = G[source].e[i].cap; + push(G[source].e[i]); + } + G[source].excess = 0; + + int c = 0; + /* Loop over active nodes using selection rule. */ + while (!rule.empty()) { + c++; + int u = rule.next(); + discharge(u); + } + + /* Output the cut based on vertex heights. */ + for (size_t i = 0; i < cut.size(); ++i) { + cut[i] = G[i].height >= N; + } +} + +/* The capacity of the edge, in the direction given by the tree. */ +int FlowGraph::treeCap(const Edge& e, Color col) const { + if (col == SOURCE) + return e.cap - e.flow; + else if (col == SINK) + return G[e.to].e[e.index].cap - G[e.to].e[e.index].flow; + else + return 0; +} + +/* Try to grow the trees S and T from the active vertices. */ +Edge *FlowGraph::grow() { + while (!bkq.empty()) { + int p = bkq.front(); + if (!G[p].active) { + bkq.pop(); + continue; + } + + size_t i = 0; + + /* + * If we're growing from the same vertex as before, + * continue with the same index, if not restart at 0. + */ + if (lastGrowVertex == p) + i = lastIndex; + + lastGrowVertex = p; + + for (; i < G[p].e.size(); ++i) { + lastIndex = i; + Edge *e = &G[p].e[i]; + int q = e->to; + + if (G[p].c == G[q].c) + continue; + + if (treeCap(*e, G[p].c) <= 0) + continue; + + if (G[q].c == FREE) { + /* Found a free vertex, add it. */ + G[q].c = G[p].c; + + int len; + if (G[p].c == SOURCE) { + G[q].p = e; + assert(treeOrigin(p, len) == source); + } else if (G[p].c == SINK) { + G[q].p = &G[q].e[e->index]; + assert(treeOrigin(p, len) == sink); + } else { + cout << G[p].c << endl; + exit(1); + } + (void)len; + + G[q].active = 1; + bkq.push(q); + } + else if (G[q].c != G[p].c) { + /* The trees meet! */ + if (G[p].c == SOURCE) { + return e; + } + else if (G[p].c == SINK) { + return &G[q].e[e->index]; + } + else { + exit(1); + } + + return NULL; + } + } + + bkq.pop(); + G[p].active = 0; + } + + /* Path is empty */ + return NULL; +} + +/* Augment along path given by the edge e, and implicitly by the trees. */ +int FlowGraph::augment(Edge* e) { + int m = e->cap - e->flow; + + /* Find maximum flow we can send. */ + Edge *cur = e; + while (cur != NULL) { + m = min(m, cur->cap - cur->flow); + cur = G[cur->from].p; + } + + cur = e; + while (cur != NULL) { + m = min(m, cur->cap - cur->flow); + cur = G[cur->to].p; + } + + cur = e; + bool back = true; + int len = 0; + /* Loop through path and update flow. */ + while (cur != NULL) { + /* If saturated, we must orphanize. */ + if (cur->cap - cur->flow == m) { + int u = cur->from; + int v = cur->to; + + if (G[u].c == SOURCE && G[v].c == SOURCE) { + if (v != source && v != sink) { + orphans.push(v); + G[v].p = NULL; + } + } + if (G[u].c == SINK && G[v].c == SINK) { + if (u != source && u != sink) { + orphans.push(u); + G[u].p = NULL; + } + } + } + len++; + push(*cur, m); + + /* + * If we reach the source, we must start again + * in e, and go towards the sink. + */ + if (back) { + cur = G[cur->from].p; + if (cur == NULL) { + back = false; + cur = G[e->to].p; + } + } else { + cur = G[cur->to].p; + } + } + return len; +} + +/* Find the origin of vertex u. */ +int FlowGraph::treeOrigin(int u, int &len) const { + int cur = u; + len = 0; + + if (G[cur].c == SOURCE) { + while (G[cur].p != NULL) { + cur = G[cur].p->from; + len++; + } + } else if (G[cur].c == SINK) { + while (G[cur].p != NULL) { + cur = G[cur].p->to; + len++; + } + } else { + exit(1); + } + + return cur; +} + +/* Adopt orphans. */ +void FlowGraph::adopt() { + while (orphans.size() > 0) { + int u = orphans.front(); + orphans.pop(); + + assert(G[u].c != FREE); + + int minlen = 1000000000; + int minidx = -1; + /* Aim to find parent close to the root of the tree. */ + for (size_t i = 0; i < G[u].e.size(); ++i) { + int v = G[u].e[i].to; + + if (G[u].c != G[v].c) + continue; + + if (treeCap(G[v].e[G[u].e[i].index], G[u].c) <= 0) + continue; + + int len; + int origin = treeOrigin(v, len); + if (origin != source && origin != sink) + continue; + + if (len < minlen) { + minlen = len; + minidx = i; + } + if (minlen <= 2) break; + /* Found a possible parent */ + } + + bool found = false; + if (minidx != -1) { + int i = minidx; + int v = G[u].e[i].to; + int len; + int origin = treeOrigin(v, len); + + if (origin == source) { + G[u].p = &G[v].e[G[u].e[i].index]; + found = true; + } else if (origin == sink) { + G[u].p = &G[u].e[i]; + found = true; + } else { + exit(1); + } + } + + /* If not found, free vertex, and orphanize possible children. */ + if (!found) { + for (size_t i = 0; i < G[u].e.size(); ++i) { + int v = G[u].e[i].to; + + if (G[u].c != G[v].c) + continue; + + if (treeCap(G[v].e[G[u].e[i].index], G[v].c) > 0) { + G[v].active = true; + bkq.push(v); + } + + if (v == source || v == sink) + continue; + + if (G[v].p + && (G[v].p->to == u + || G[v].p->from == u)) { + orphans.push(v); + G[v].p = NULL; + } + } + + G[u].c = FREE; + + G[u].active = false; + /* We might still have u in the queue */ + } + } +} + +void FlowGraph::minCutBK(int source, int sink) { + lastGrowVertex = -1; + adopt(); + + int numpaths = 0; + double totlen = 0; + while (true) { + Edge *e; + e = grow(); + + if (e == NULL) { + /* Empty path. */ + break; + } + + totlen += augment(e); + numpaths++; + adopt(); + } + + cout << "Avg length: " << double(totlen) / double(numpaths) << endl; + + int size1 = 0, size2 = 0; + for (size_t i = 0; i < cut.size(); ++i) { + if (G[i].c == SOURCE) size1++; + else if (G[i].c == SINK) size2++; + cut[i] = G[i].c == SOURCE; + } + cout << "Inbetweeners: " << cut.size() - size1 - size2 << endl; +} + diff --git a/graph.hpp b/graph.hpp new file mode 100644 index 0000000..ea183c9 --- /dev/null +++ b/graph.hpp @@ -0,0 +1,104 @@ +#pragma once + +#include +#include +#include +#include +#include +#include "selectionrule.hpp" + +enum Color { FREE = 0, SOURCE, SINK }; + +class Edge { +private: + +public: + int from, to; + int cap; + int flow; + int index; + + Edge(int f, int t, int c, int i) : + from(f), to(t), cap(c), flow(0), index(i) {} +}; + +class Vertex { +private: + +public: + std::vector e; + Color c; + bool active; + Edge *p; + int height; + int excess; + int si; + int ti; + + Vertex(Color c, bool a) : + c(c), active(a), p(NULL), height(0), excess(0), si(0), ti(0) {} + + Vertex() : + c(FREE), active(false), p(NULL), height(0), excess(0), si(0), ti(0) {} +}; + +class FlowGraph { +private: + int N; + int source, sink; + std::vector G; + std::vector count; + SelectionRule& rule; + + /* BK stuff */ + std::queue bkq; + std::queue orphans; + + int lastGrowVertex; + size_t lastIndex; +public: + std::vector cut; + + FlowGraph(int N, int source, int sink, SelectionRule& rule) : + N(N), + source(source), + sink(sink), + G(N), + count(N+1), + rule(rule), + lastGrowVertex(-1), + cut(N) { + +#ifdef BOYKOV_KOLMOGOROV + bkq.push(source); + bkq.push(sink); + G[source].c = SOURCE; + G[sink].c = SINK; + G[source].active = true; + G[sink].active = true; +#endif + } + + int getSource() const { return source; } + int getSink() const { return sink; } + void addEdge(int from, int to, int cap); + void addDoubleEdge(int from, int to, int cap); + void changeCapacity(int from, int index, int cap); + void resetFlow(); + void resetHeights(); + + void push(Edge &e); + void push(Edge &e, int f); + void relabel(int u); + void gap(int h); + void discharge(int u); + void minCutPushRelabel(int source, int sink); + + void minCutBK(int source, int sink); + int augment(Edge *e); + int treeCap(const Edge& e, Color col) const; + int treeOrigin(int u, int &len) const; + void adopt(); + Edge *grow(); +}; + diff --git a/main.tex b/main.tex index a7c4662..1ef24ef 100644 --- a/main.tex +++ b/main.tex @@ -21,7 +21,7 @@ %\mathtoolsset{showonlyrefs=true} -\newminted{cpp}{fontsize=\tiny} +%\newminted{cpp}{fontsize=\tiny} \usepackage{polyglossia} \setmainlanguage[variant=american]{english} diff --git a/maxflow.tex b/maxflow.tex index 0fbc607..e775362 100644 --- a/maxflow.tex +++ b/maxflow.tex @@ -9,6 +9,8 @@ sending flow through the graph and trying to identify the ``bottleneck''. This chapter, except for the description of the Boykov--Kolmogorov algorithm is taken with some adjustments from my project work \cite{project} and is included here for completeness. +Implementations of the push-relabel and Boykov--Kolmogorov algorithms +can be found in Appendix~\ref{app:impl}. \section{Flow graphs} @@ -163,7 +165,7 @@ path, and note that it follows an edge in $E$ in the reverse direction, made possible by the construction of the residual graph just described. -\begin{figure} +\begin{figure}[b] \input{fig/aug_flow} \end{figure} @@ -517,6 +519,7 @@ selection rules, heuristics and their implementation in \cite{cherkassky1997implementing}. \subsection{Heuristics} + Different heuristics exist that can speed up the algorithm considerably. Being heuristics, they are not guaranteed to work, and might perform differently on different kinds of graphs. The most used @@ -1066,35 +1069,3 @@ paths from $s$ to $t$ as there are pixels in the image. When increasing the capacity of edges $(v, t)$, a quick sweep over these two-edged paths to send any possible flow may speed up the algorithm. -\section{Implementation} - -A \cpp{} implementation was submitted together with this thesis and can -also be found online at \cite{github}. It uses -the open computer vision library OpenCV \cite{opencv_library} to load -and save image files and contains compilation and usage instructions. -The implementation has been tested on an installation of the Ubuntu -Linux distribution, but it should in theory be portable to other -platforms supported by OpenCV. A rudimentary graphics interface has also -been made, to make it easier to play with the parameters of the -algorithm. - -Both the push-relabel and the Boykov--Kolmogorov algorithms have been -implemented. Although an effort has been made to improve the performance -of both implementations, they are not ment to beat the fastest. The -focus has rather been on clarity and understanding. - -Note that when implementing maximum flow algorithms it is not a good -idea, memory- and performance-wise, to actually construct the residual -graph $G_f$. Instead, every time we update the flow $f(u,v)$ we set -the flow in the opposite direction to its negative value $f(v,u) = --f(u,v)$. Then we can at any time, consider the value $c(u,v) - f(u,v)$ -in the place of the residual capacity $c_f(u,v)$. - -For the gap relabeling heuristic of the push-relabel algorithm, we need -to have a easy way of finding when a gap occurs. This is done by keeping -track of how many vertices exist with each label. - -When the capacities have been updated in the Boykov--Kolmogorov -algorithm, flow is sent along all two-edge paths such that they do not -have to be considered by the main loop of the algorithm. - -- 2.47.3