\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}
--- /dev/null
+#include <iostream>
+#include <queue>
+#include <stack>
+#include <vector>
+#include <algorithm>
+#include <cassert>
+#include <cstdlib>
+
+#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;
+}
+
--- /dev/null
+#pragma once
+
+#include <iostream>
+#include <queue>
+#include <vector>
+#include <list>
+#include <set>
+#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<Edge> 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<Vertex > G;
+ std::vector<int> count;
+ SelectionRule& rule;
+
+ /* BK stuff */
+ std::queue<int> bkq;
+ std::queue<int> orphans;
+
+ int lastGrowVertex;
+ size_t lastIndex;
+public:
+ std::vector<char> 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();
+};
+
``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}
made possible by the construction of the residual graph just
described.
-\begin{figure}
+\begin{figure}[b]
\input{fig/aug_flow}
\end{figure}
\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
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.
-
projects / master / commitdiff