diff --git a/dlib/graph_cuts.h b/dlib/graph_cuts.h new file mode 100644 index 000000000..10c677a04 --- /dev/null +++ b/dlib/graph_cuts.h @@ -0,0 +1,13 @@ +// Copyright (C) 2012 Davis E. King (davis@dlib.net) +// License: Boost Software License See LICENSE.txt for the full license. +#ifndef DLIB_GRAPH_CUTs_HEADER_ +#define DLIB_GRAPH_CUTs_HEADER_ + +#include "graph_cuts/min_cut.h" +#include "graph_cuts/general_flow_graph.h" +#include "graph_cuts/find_max_factor_graph_potts.h" + +#endif // DLIB_GRAPH_CUTs_HEADER_ + + + diff --git a/dlib/graph_cuts/find_max_factor_graph_potts.h b/dlib/graph_cuts/find_max_factor_graph_potts.h new file mode 100644 index 000000000..a2edbcaa9 --- /dev/null +++ b/dlib/graph_cuts/find_max_factor_graph_potts.h @@ -0,0 +1,477 @@ +// Copyright (C) 2012 Davis E. King (davis@dlib.net) +// License: Boost Software License See LICENSE.txt for the full license. +#ifndef DLIB_FIND_MAX_FACTOR_GRAPH_PoTTS_H__ +#define DLIB_FIND_MAX_FACTOR_GRAPH_PoTTS_H__ + +#include "find_max_factor_graph_potts_abstract.h" +#include "../matrix.h" +#include "min_cut.h" + +namespace dlib +{ + +// ---------------------------------------------------------------------------------------- + + namespace impl + { + + template < + typename potts_problem, + typename T = void + > + class flows_container + { + /* + This object notionally represents a matrix of flow values. It's + overloaded to represent this matrix efficiently though. In this case + it represents the matrix using a sparse representation. + */ + + typedef typename potts_problem::value_type edge_type; + std::vector > flows; + public: + + void setup( + const potts_problem& p + ) + { + flows.resize(p.number_of_nodes()); + for (unsigned long i = 0; i < flows.size(); ++i) + { + flows[i].resize(p.number_of_neighbors(i)); + } + } + + edge_type& operator() ( + const long r, + const long c + ) { return flows[r][c]; } + + const edge_type& operator() ( + const long r, + const long c + ) const { return flows[r][c]; } + }; + +// ---------------------------------------------------------------------------------------- + + template < + typename potts_problem + > + class flows_container::type> + { + /* + This object notionally represents a matrix of flow values. It's + overloaded to represent this matrix efficiently though. In this case + it represents the matrix using a dense representation. + + */ + typedef typename potts_problem::value_type edge_type; + const static unsigned long max_number_of_neighbors = potts_problem::max_number_of_neighbors; + matrix flows; + public: + + void setup( + const potts_problem& p + ) + { + flows.set_size(p.number_of_nodes(), max_number_of_neighbors); + } + + edge_type& operator() ( + const long r, + const long c + ) { return flows(r,c); } + + const edge_type& operator() ( + const long r, + const long c + ) const { return flows(r,c); } + }; + +// ---------------------------------------------------------------------------------------- + + template < + typename potts_problem + > + class potts_flow_graph + { + public: + typedef typename potts_problem::value_type edge_type; + private: + /*! + This is a utility class used by dlib::min_cut to convert a potts_problem + into the kind of flow graph expected by the min_cut object's main block + of code. + + Within this object, we will use the convention that one past + potts_problem::number_of_nodes() is the source node and two past is + the sink node. + !*/ + + potts_problem& g; + + // flows(i,j) == the flow from node id i to it's jth neighbor + flows_container flows; + // source_flows(i,0) == flow from source to node i, + // source_flows(i,1) == flow from node i to source + matrix source_flows; + + // sink_flows(i,0) == flow from sink to node i, + // sink_flows(i,1) == flow from node i to sink + matrix sink_flows; + + node_label source_label, sink_label; + public: + + potts_flow_graph( + potts_problem& g_ + ) : g(g_) + { + flows.setup(g); + + source_flows.set_size(g.number_of_nodes(), 2); + sink_flows.set_size(g.number_of_nodes(), 2); + source_flows = 0; + sink_flows = 0; + + source_label = FREE_NODE; + sink_label = FREE_NODE; + + // setup flows based on factor potentials + edge_type min_value = std::numeric_limits::max(); + for (unsigned long i = 0; i < g.number_of_nodes(); ++i) + { + source_flows(i,0) = -g.factor_value(i, true); + sink_flows(i,1) = -g.factor_value(i, false); + + source_flows(i,1) = 0; + sink_flows(i,0) = 0; + + if (source_flows(i,0) < min_value) + min_value = source_flows(i,0); + if (sink_flows(i,1) < min_value) + min_value = sink_flows(i,1); + + for (unsigned long j = 0; j < g.number_of_neighbors(i); ++j) + { + flows(i,j) = g.factor_value_disagreement(i, g.get_neighbor(i,j)); + } + } + + // The flows can't be negative, so if they are then just adjust them + // so they are positive. Note that doing this doesn't change the optimal + // labeling assignment. + if (min_value < 0) + { + source_flows -= min_value; + sink_flows -= min_value; + } + } + + class out_edge_iterator + { + friend class potts_flow_graph; + unsigned long idx; // base node idx + unsigned long cnt; // count over the neighbors of idx + public: + + out_edge_iterator( + ):idx(0),cnt(0){} + + out_edge_iterator( + unsigned long idx_, + unsigned long cnt_ + ):idx(idx_),cnt(cnt_) + {} + + bool operator!= ( + const out_edge_iterator& item + ) const { return cnt != item.cnt; } + + out_edge_iterator& operator++( + ) + { + ++cnt; + return *this; + } + }; + + class in_edge_iterator + { + friend class potts_flow_graph; + unsigned long idx; // base node idx + unsigned long cnt; // count over the neighbors of idx + public: + + in_edge_iterator( + ):idx(0),cnt(0) + {} + + + in_edge_iterator( + unsigned long idx_, + unsigned long cnt_ + ):idx(idx_),cnt(cnt_) + {} + + bool operator!= ( + const in_edge_iterator& item + ) const { return cnt != item.cnt; } + + in_edge_iterator& operator++( + ) + { + ++cnt; + return *this; + } + }; + + unsigned long number_of_nodes ( + ) const { return g.number_of_nodes() + 2; } + + out_edge_iterator out_begin( + const unsigned long& it + ) const { return out_edge_iterator(it, 0); } + + in_edge_iterator in_begin( + const unsigned long& it + ) const { return in_edge_iterator(it, 0); } + + out_edge_iterator out_end( + const unsigned long& it + ) const + { + if (it >= g.number_of_nodes()) + return out_edge_iterator(it, g.number_of_nodes()); + else + return out_edge_iterator(it, g.number_of_neighbors(it)+2); + } + + in_edge_iterator in_end( + const unsigned long& it + ) const + { + if (it >= g.number_of_nodes()) + return in_edge_iterator(it, g.number_of_nodes()); + else + return in_edge_iterator(it, g.number_of_neighbors(it)+2); + } + + + template + unsigned long node_id ( + const iterator_type& it + ) const + { + // if this isn't an iterator over the source or sink nodes + if (it.idx < g.number_of_nodes()) + { + const unsigned long num = g.number_of_neighbors(it.idx); + if (it.cnt < num) + return g.get_neighbor(it.idx, it.cnt); + else if (it.cnt == num) + return g.number_of_nodes(); + else + return g.number_of_nodes()+1; + } + else + { + return it.cnt; + } + } + + + edge_type get_flow ( + const unsigned long& it1, + const unsigned long& it2 + ) const + { + if (it1 >= g.number_of_nodes()) + { + // if it1 is the source + if (it1 == g.number_of_nodes()) + { + return source_flows(it2,0); + } + else // if it1 is the sink + { + return sink_flows(it2,0); + } + } + else if (it2 >= g.number_of_nodes()) + { + // if it2 is the source + if (it2 == g.number_of_nodes()) + { + return source_flows(it1,1); + } + else // if it2 is the sink + { + return sink_flows(it1,1); + } + } + else + { + return flows(it1, g.get_neighbor_idx(it1, it2)); + } + + } + + edge_type get_flow ( + const out_edge_iterator& it + ) const + { + if (it.idx < g.number_of_nodes()) + { + const unsigned long num = g.number_of_neighbors(it.idx); + if (it.cnt < num) + return flows(it.idx, it.cnt); + else if (it.cnt == num) + return source_flows(it.idx,1); + else + return sink_flows(it.idx,1); + } + else + { + // if it.idx is the source + if (it.idx == g.number_of_nodes()) + { + return source_flows(it.cnt,0); + } + else // if it.idx is the sink + { + return sink_flows(it.cnt,0); + } + } + } + + edge_type get_flow ( + const in_edge_iterator& it + ) const + { + return get_flow(node_id(it), it.idx); + } + + void adjust_flow ( + const unsigned long& it1, + const unsigned long& it2, + const edge_type& value + ) + { + if (it1 >= g.number_of_nodes()) + { + // if it1 is the source + if (it1 == g.number_of_nodes()) + { + source_flows(it2,0) += value; + source_flows(it2,1) -= value; + } + else // if it1 is the sink + { + sink_flows(it2,0) += value; + sink_flows(it2,1) -= value; + } + } + else if (it2 >= g.number_of_nodes()) + { + // if it2 is the source + if (it2 == g.number_of_nodes()) + { + source_flows(it1,1) += value; + source_flows(it1,0) -= value; + } + else // if it2 is the sink + { + sink_flows(it1,1) += value; + sink_flows(it1,0) -= value; + } + } + else + { + flows(it1, g.get_neighbor_idx(it1, it2)) += value; + flows(it2, g.get_neighbor_idx(it2, it1)) -= value; + } + + } + + void set_label ( + const unsigned long& it, + node_label value + ) + { + if (it < g.number_of_nodes()) + g.set_label(it, value); + else if (it == g.number_of_nodes()) + source_label = value; + else + sink_label = value; + } + + node_label get_label ( + const unsigned long& it + ) const + { + if (it < g.number_of_nodes()) + return g.get_label(it); + if (it == g.number_of_nodes()) + return source_label; + else + return sink_label; + } + + }; + } + +// ---------------------------------------------------------------------------------------- +// ---------------------------------------------------------------------------------------- +// ---------------------------------------------------------------------------------------- + + template < + typename potts_model + > + typename potts_model::value_type potts_model_score ( + const potts_model& g + ) + { + typename potts_model::value_type score = 0; + for (unsigned long i = 0; i < g.number_of_nodes(); ++i) + { + const bool label = (g.get_label(i)!=0); + score += g.factor_value(i, label); + } + + for (unsigned long i = 0; i < g.number_of_nodes(); ++i) + { + for (unsigned long n = 0; n < g.number_of_neighbors(i); ++n) + { + const unsigned long idx2 = g.get_neighbor(i,n); + if (g.get_label(i) != g.get_label(idx2) && i < idx2) + score -= g.factor_value_disagreement(i, idx2); + } + } + + return score; + } + +// ---------------------------------------------------------------------------------------- + + template < + typename potts_model + > + void find_max_factor_graph_potts ( + potts_model& m + ) + { + min_cut mc; + dlib::impl::potts_flow_graph pfg(m); + mc(pfg, m.number_of_nodes(), m.number_of_nodes()+1); + } + +// ---------------------------------------------------------------------------------------- + +} + +#endif // DLIB_FIND_MAX_FACTOR_GRAPH_PoTTS_H__ + diff --git a/dlib/graph_cuts/find_max_factor_graph_potts_abstract.h b/dlib/graph_cuts/find_max_factor_graph_potts_abstract.h new file mode 100644 index 000000000..daf5c7ee6 --- /dev/null +++ b/dlib/graph_cuts/find_max_factor_graph_potts_abstract.h @@ -0,0 +1,223 @@ +// Copyright (C) 2012 Davis E. King (davis@dlib.net) +// License: Boost Software License See LICENSE.txt for the full license. +#undef DLIB_FIND_MAX_FACTOR_GRAPH_PoTTS_ABSTRACT_H__ +#ifdef DLIB_FIND_MAX_FACTOR_GRAPH_PoTTS_ABSTRACT_H__ + +#include "../matrix.h" +#include "min_cut_abstract.h" + +namespace dlib +{ + +// ---------------------------------------------------------------------------------------- + + class potts_problem + { + /*! + WHAT THIS OBJECT REPRESENTS + This object represents a boolean valued factor graph or graphical model + that can be efficiently operated on using graph cuts. In particular, this + object defines the interface a MAP problem on a factor graph must + implement if it is to be solved using the find_max_factor_graph_potts() + routine defined at the bottom of this file. + + Note that there is no dlib::potts_problem object. What you are looking + at here is simply the interface definition for a Potts problem. You must + implement your own version of this object for the problem you wish to + solve and then pass it to the find_max_factor_graph_potts() routine. + + Note also that a factor graph should not have any nodes which are + neighbors with themselves. Additionally, the graph is undirected. This + mean that if A is a neighbor of B then B must be a neighbor of A for + the MAP problem to be valid. + !*/ + + public: + + unsigned long number_of_nodes ( + ) const; + /*! + ensures + - returns the number of nodes in the factor graph. Or in other words, + returns the number of variables in the MAP problem/Potts model. + !*/ + + unsigned long number_of_neighbors ( + unsigned long idx + ) const; + /*! + requires + - idx < number_of_nodes() + ensures + - returns the number of neighbors of node idx. + !*/ + + // This is an optional variable which specifies a number that is always + // greater than or equal to number_of_neighbors(idx). If you don't know + // the value at compile time then either don't include max_number_of_neighbors + // in your potts_problem object or set it to 0. + const static unsigned long max_number_of_neighbors = 0; + + unsigned long get_neighbor ( + unsigned long idx, + unsigned long n + ) const; + /*! + requires + - idx < number_of_nodes() + - n < number_of_neighbors(idx) + ensures + - returns the node index value of the n-th neighbor of + the node with index value idx. + - The neighbor relationship is reciprocal. That is, if + get_neighbor(A,i)==B then there is a value of j such + that get_neighbor(B,j)==A. + - A node is never its own neighbor. That is, there is + no i such that get_neighbor(idx,i)==idx. + !*/ + + unsigned long get_neighbor_idx ( + unsigned long idx1, + unsigned long idx2 + ) const; + /*! + requires + - idx1 < number_of_nodes() + - idx2 < number_of_nodes() + ensures + - This function is basically the inverse of get_neighbor(). + - returns a number IDX such that: + - get_neighbor(idx1,IDX) == idx2 + - IDX < number_of_neighbors(node_idx1) + !*/ + + void set_label ( + const unsigned long& idx, + node_label value + ); + /*! + requires + - idx < number_of_nodes() + ensures + - #get_label(idx) == value + !*/ + + node_label get_label ( + const unsigned long& idx + ) const; + /*! + requires + - idx < number_of_nodes() + ensures + - returns the current label for the idx-th node. This is a value which is + 0 if the node's label is false and is any other value if it is true. + + Note that this value is not used by factor_value() or factor_value_disagreement(). + It is simply here to provide a mechanism for find_max_factor_graph_potts() + to return its labeled result. Additionally, the reason it returns a + node_label rather than a bool is because doing it this way facilitates + use of a graph cut algorithm for the solution of the MAP problem. For + more of an explanation you should read the paper referenced by the min_cut + object. + !*/ + + // This typedef should be for a type like int or double. It + // must also be capable of representing signed values. + typedef an_integer_or_real_type value_type; + + value_type factor_value ( + unsigned long idx, + bool value + ) const; + /*! + requires + - idx < number_of_nodes() + ensures + - returns a value which indicates how "good" it is to assign the idx-node + a label equal to value. The larger the value, the more desirable + the label contained by value. + !*/ + + value_type factor_value_disagreement ( + unsigned long idx1, + unsigned long idx2 + ) const; + /*! + requires + - idx1 < number_of_nodes() + - idx2 < number_of_nodes() + - idx1 != idx2 + - the idx1-th node and idx2-th node are neighbors in the graph. That is, + get_neighbor(idx1,i)==idx2 for some value of i. + ensures + - returns a number >= 0. This is the penalty for giving node idx1 and idx2 + different labels. Larger values indicate a larger penalty. + - this function is symmetric. That is, it is true that: + factor_value_disagreement(i,j) == factor_value_disagreement(j,i) + !*/ + + }; + +// ---------------------------------------------------------------------------------------- + + template < + typename potts_problem + > + typename potts_problem::value_type potts_model_score ( + const potts_problem& prob + ); + /*! + requires + - potts_problem == an object with an interface compatible with the potts_problem + object defined at the top of this file. + - for all valid i and j: + - prob.factor_value_disagreement(i,j) == prob.factor_value_disagreement(j,i) + ensures + - computes the model score for the given potts_problem. To define this + precisely: + - let L(i) == the boolean label of the ith variable in prob. Or in other + words, L(i) == (prob.get_label(i) != 0). + - let F == the sum over valid i of prob.factor_value(i, L(i)). + - Let D == the sum of all values of prob.factor_value_disagreement(i,j) + whenever the following conditions are true about i and j: + - i and j are neighbors in the graph defined by prob, that is, + it is valid to call prob.factor_value_disagreement(i,j). + - L(i) != L(j) + - i < j + (i.e. We want to make sure to only count the edge between i and j once) + + - Then this function returns F - D + !*/ + +// ---------------------------------------------------------------------------------------- + + template < + typename potts_problem + > + void find_max_factor_graph_potts ( + potts_problem& prob + ) + /*! + requires + - potts_problem == an object with an interface compatible with the potts_problem + object defined at the top of this file. + - for all valid i and j: + - prob.factor_value_disagreement(i,j) >= 0 + - prob.factor_value_disagreement(i,j) == prob.factor_value_disagreement(j,i) + ensures + - This function is a tool for exactly solving the MAP problem in a Potts + model. In particular, this means that this function finds the assignments + to all the labels in prob which maximizes potts_model_score(#prob). + - Note that this routine is a little bit faster if all the values + returned by prob.factor_value() are negative. So if you can arrange for that + to be true without spending any extra CPU cycles then it might be a good idea + to do so. + !*/ + +// ---------------------------------------------------------------------------------------- + +} + +#endif // DLIB_FIND_MAX_FACTOR_GRAPH_PoTTS_ABSTRACT_H__ + + diff --git a/dlib/graph_cuts/general_flow_graph.h b/dlib/graph_cuts/general_flow_graph.h new file mode 100644 index 000000000..cf4184b7f --- /dev/null +++ b/dlib/graph_cuts/general_flow_graph.h @@ -0,0 +1,172 @@ +// Copyright (C) 2012 Davis E. King (davis@dlib.net) +// License: Boost Software License See LICENSE.txt for the full license. +#ifndef DLIB_GENERAL_FLOW_GRaPH_H__ +#define DLIB_GENERAL_FLOW_GRaPH_H__ + +#include "../graph_utils.h" + +namespace dlib +{ + +// ---------------------------------------------------------------------------------------- + + namespace impl + { + + template < + typename directed_graph_type + > + class general_flow_graph + { + /*! + this is a utility class used by dlib::min_cut to convert a directed_graph + into the kind of flow graph expected by the min_cut object's main block + of code. + !*/ + + directed_graph_type& g; + + typedef typename directed_graph_type::node_type node_type; + typedef typename directed_graph_type::type node_label; + + public: + + general_flow_graph( + directed_graph_type& g_ + ) : g(g_) + { + } + + class out_edge_iterator + { + friend class general_flow_graph; + unsigned long idx; // base node idx + unsigned long cnt; // count over the neighbors of idx + public: + + out_edge_iterator( + ):idx(0),cnt(0) {} + + out_edge_iterator( + unsigned long idx_, + unsigned long cnt_ + ):idx(idx_),cnt(cnt_) + {} + + bool operator!= ( + const out_edge_iterator& item + ) const { return cnt != item.cnt; } + + out_edge_iterator& operator++( + ) + { + ++cnt; + return *this; + } + }; + + class in_edge_iterator + { + + friend class general_flow_graph; + unsigned long idx; // base node idx + unsigned long cnt; // count over the neighbors of idx + public: + + in_edge_iterator( + ):idx(0),cnt(0) {} + + in_edge_iterator( + unsigned long idx_, + unsigned long cnt_ + ):idx(idx_),cnt(cnt_) + {} + + bool operator!= ( + const in_edge_iterator& item + ) const { return cnt != item.cnt; } + + in_edge_iterator& operator++( + ) + { + ++cnt; + return *this; + } + }; + + unsigned long number_of_nodes ( + ) const { return g.number_of_nodes(); } + + out_edge_iterator out_begin( + const unsigned long& it + ) const { return out_edge_iterator(it, 0); } + + in_edge_iterator in_begin( + const unsigned long& it + ) const { return in_edge_iterator(it, 0); } + + out_edge_iterator out_end( + const unsigned long& it + ) const { return out_edge_iterator(it, g.node(it).number_of_children()); } + + in_edge_iterator in_end( + const unsigned long& it + ) const { return in_edge_iterator(it, g.node(it).number_of_parents()); } + + unsigned long node_id ( + const out_edge_iterator& it + ) const { return g.node(it.idx).child(it.cnt).index(); } + unsigned long node_id ( + const in_edge_iterator& it + ) const { return g.node(it.idx).parent(it.cnt).index(); } + + typedef typename directed_graph_type::edge_type edge_type; + + edge_type get_flow (const unsigned long& it1, const unsigned long& it2) const + { + return edge(g, it1, it2); + } + edge_type get_flow (const out_edge_iterator& it) const + { + return g.node(it.idx).child_edge(it.cnt); + } + edge_type get_flow (const in_edge_iterator& it) const + { + return g.node(it.idx).parent_edge(it.cnt); + } + + void adjust_flow ( + const unsigned long& it1, + const unsigned long& it2, + const edge_type& value + ) + { + edge(g, it1, it2) += value; + edge(g, it2, it1) -= value; + } + + void set_label ( + const unsigned long& it, + node_label value + ) + { + g.node(it).data = value; + } + + node_label get_label ( + const unsigned long& it + ) const + { + return g.node(it).data; + } + + }; + + } + +// ---------------------------------------------------------------------------------------- + +} + +#endif // DLIB_GENERAL_FLOW_GRaPH_H__ + diff --git a/dlib/graph_cuts/min_cut.h b/dlib/graph_cuts/min_cut.h new file mode 100644 index 000000000..bd6a3dfbd --- /dev/null +++ b/dlib/graph_cuts/min_cut.h @@ -0,0 +1,572 @@ +// Copyright (C) 2012 Davis E. King (davis@dlib.net) +// License: Boost Software License See LICENSE.txt for the full license. +#ifndef DLIB_MIN_CuT_H__ +#define DLIB_MIN_CuT_H__ + +#include "min_cut_abstract.h" +#include "../matrix.h" +#include "general_flow_graph.h" +#include "../is_kind.h" + +#include +#include +#include + + +// ---------------------------------------------------------------------------------------- + + +namespace dlib +{ + + typedef unsigned char node_label; + +// ---------------------------------------------------------------------------------------- + + const node_label SOURCE_CUT = 0; + const node_label SINK_CUT = 1; + const node_label FREE_NODE = 2; + +// ---------------------------------------------------------------------------------------- + + template + typename disable_if,typename flow_graph::edge_type>::type + graph_cut_score ( + const flow_graph& g + ) + { + typedef typename flow_graph::edge_type edge_weight_type; + edge_weight_type score = 0; + typedef typename flow_graph::out_edge_iterator out_edge_iterator; + for (unsigned long i = 0; i < g.number_of_nodes(); ++i) + { + if (g.get_label(i) != SOURCE_CUT) + continue; + + for (out_edge_iterator n = g.out_begin(i); n != g.out_end(i); ++n) + { + if (g.get_label(g.node_id(n)) != SOURCE_CUT) + { + score += g.get_flow(n); + } + } + } + + return score; + } + + template + typename enable_if,typename directed_graph::edge_type>::type + graph_cut_score ( + const directed_graph& g + ) + { + return graph_cut_score(dlib::impl::general_flow_graph(g)); + } + +// ---------------------------------------------------------------------------------------- + + class min_cut + { + + public: + + min_cut() + { + } + + min_cut( const min_cut& ) + { + // Intentionally left empty since all the member variables + // don't logically contribute to the state of this object. + // This copy constructor is here to explicitly void the overhead + // of copying these transient variables. + } + + template < + typename directed_graph + > + typename enable_if >::type operator() ( + directed_graph& g, + const unsigned long source_node, + const unsigned long sink_node + ) const + { + DLIB_ASSERT(graph_contains_length_one_cycle(g) == false, + "\t void min_cut::operator()" + << "\n\t Invalid arguments were given to this function." + ); + DLIB_ASSERT(graph_has_symmetric_edges(g) == true, + "\t void min_cut::operator()" + << "\n\t Invalid arguments were given to this function." + ); + + dlib::impl::general_flow_graph temp(g); + (*this)(temp, source_node, sink_node); + } + + template < + typename flow_graph + > + typename disable_if >::type operator() ( + flow_graph& g, + const unsigned long source_node, + const unsigned long sink_node + ) const + { +#ifdef ENABLE_ASSERTS + DLIB_ASSERT(source_node != sink_node && + source_node < g.number_of_nodes() && + sink_node < g.number_of_nodes(), + "\t void min_cut::operator()" + << "\n\t Invalid arguments were given to this function." + << "\n\t g.number_of_nodes(): " << g.number_of_nodes() + << "\n\t source_node: " << source_node + << "\n\t sink_node: " << sink_node + << "\n\t this: " << this + ); + + for (unsigned long i = 0; i < g.number_of_nodes(); ++i) + { + typename flow_graph::out_edge_iterator j, end = g.out_end(i); + for (j = g.out_begin(i); j != end; ++j) + { + const unsigned long jj = g.node_id(j); + DLIB_ASSERT(g.get_flow(i,jj) >= 0, + "\t void min_cut::operator()" + << "\n\t i: "<< i + << "\n\t jj: "<< jj + << "\n\t g.get_flow(i,jj): "<< g.get_flow(i,jj) + << "\n\t this: "<< this + ); + + } + } +#endif + parent.clear(); + active.clear(); + orphans.clear(); + + typedef typename flow_graph::edge_type edge_type; + typedef typename flow_graph::out_edge_iterator out_edge_iterator; + typedef typename flow_graph::in_edge_iterator in_edge_iterator; + + // initialize labels + for (unsigned long i = 0; i < g.number_of_nodes(); ++i) + g.set_label(i, FREE_NODE); + + g.set_label(source_node, SOURCE_CUT); + g.set_label(sink_node, SINK_CUT); + + // used to indicate "no parent" + const unsigned long nil = g.number_of_nodes(); + + parent.assign(g.number_of_nodes(), nil); + + time = 1; + dist.assign(g.number_of_nodes(), 0); + ts.assign(g.number_of_nodes(), time); + + active.push_back(source_node); + active.push_back(sink_node); + + + in_edge_iterator in_begin = g.in_begin(active[0]); + out_edge_iterator out_begin = g.out_begin(active[0]); + + unsigned long source_side, sink_side; + while (grow(g,source_side,sink_side, in_begin, out_begin)) + { + ++time; + ts[source_node] = time; + ts[sink_node] = time; + + augment(g, source_node, sink_node, source_side, sink_side); + adopt(g, source_node, sink_node); + } + + } + + + private: + + unsigned long distance_to_origin ( + const unsigned long nil, + unsigned long p, + unsigned long + ) const + { + unsigned long start = p; + unsigned long count = 0; + while (p != nil) + { + if (ts[p] == time) + { + count += dist[p]; + + unsigned long count_down = count; + // adjust the dist and ts for the nodes on this path. + while (start != p) + { + ts[start] = time; + dist[start] = count_down; + --count_down; + start = parent[start]; + } + + return count; + } + p = parent[p]; + ++count; + } + + return std::numeric_limits::max(); + } + + template + void adopt ( + flow_graph& g, + const unsigned long source, + const unsigned long sink + ) const + { + typedef typename flow_graph::edge_type edge_type; + typedef typename flow_graph::out_edge_iterator out_edge_iterator; + typedef typename flow_graph::in_edge_iterator in_edge_iterator; + + // used to indicate "no parent" + const unsigned long nil = g.number_of_nodes(); + + while (orphans.size() > 0) + { + const unsigned long p = orphans.front(); + orphans.pop_front(); + + const unsigned char label_p = g.get_label(p); + + // Try to find a valid parent for p. + if (label_p == SOURCE_CUT) + { + const in_edge_iterator begin(g.in_begin(p)); + const in_edge_iterator end(g.in_end(p)); + unsigned long best_dist = std::numeric_limits::max(); + unsigned long best_node = 0; + for(in_edge_iterator q = begin; q != end; ++q) + { + const unsigned long id = g.node_id(q); + + if (g.get_label(id) != label_p || g.get_flow(q) <= 0 ) + continue; + + unsigned long temp = distance_to_origin(nil, id,source); + if (temp < best_dist) + { + best_dist = temp; + best_node = id; + } + + } + if (best_dist != std::numeric_limits::max()) + { + parent[p] = best_node; + dist[p] = dist[best_node] + 1; + ts[p] = time; + } + + // if we didn't find a parent for p + if (parent[p] == nil) + { + for(in_edge_iterator q = begin; q != end; ++q) + { + const unsigned long id = g.node_id(q); + + if (g.get_label(id) != SOURCE_CUT) + continue; + + if (g.get_flow(q) > 0) + active.push_back(id); + + if (parent[id] == p) + { + parent[id] = nil; + orphans.push_back(id); + } + } + g.set_label(p, FREE_NODE); + } + } + else + { + unsigned long best_node = 0; + unsigned long best_dist = std::numeric_limits::max(); + const out_edge_iterator begin(g.out_begin(p)); + const out_edge_iterator end(g.out_end(p)); + for(out_edge_iterator q = begin; q != end; ++q) + { + const unsigned long id = g.node_id(q); + if (g.get_label(id) != label_p || g.get_flow(q) <= 0) + continue; + + unsigned long temp = distance_to_origin(nil, id,sink); + + if (temp < best_dist) + { + best_dist = temp; + best_node = id; + } + } + + if (best_dist != std::numeric_limits::max()) + { + parent[p] = best_node; + dist[p] = dist[best_node] + 1; + ts[p] = time; + } + + // if we didn't find a parent for p + if (parent[p] == nil) + { + for(out_edge_iterator q = begin; q != end; ++q) + { + const unsigned long id = g.node_id(q); + + if (g.get_label(id) != SINK_CUT) + continue; + + if (g.get_flow(q) > 0) + active.push_back(id); + + if (parent[id] == p) + { + parent[id] = nil; + orphans.push_back(id); + } + } + + g.set_label(p, FREE_NODE); + } + } + + + } + + } + + template + void augment ( + flow_graph& g, + const unsigned long& source, + const unsigned long& sink, + const unsigned long& source_side, + const unsigned long& sink_side + ) const + { + typedef typename flow_graph::edge_type edge_type; + typedef typename flow_graph::out_edge_iterator out_edge_iterator; + typedef typename flow_graph::in_edge_iterator in_edge_iterator; + + // used to indicate "no parent" + const unsigned long nil = g.number_of_nodes(); + + unsigned long s = source_side; + unsigned long t = sink_side; + edge_type min_cap = g.get_flow(s,t); + + // find the bottleneck capacity on the current path. + + // check from source_side back to the source for the min capacity link. + t = s; + while (t != source) + { + s = parent[t]; + const edge_type temp = g.get_flow(s, t); + if (temp < min_cap) + { + min_cap = temp; + } + t = s; + } + + // check from sink_side back to the sink for the min capacity link + s = sink_side; + while (s != sink) + { + t = parent[s]; + const edge_type temp = g.get_flow(s, t); + if (temp < min_cap) + { + min_cap = temp; + } + s = t; + } + + + // now push the max possible amount of flow though the path + s = source_side; + t = sink_side; + g.adjust_flow(t,s, min_cap); + + // trace back towards the source + t = s; + while (t != source) + { + s = parent[t]; + g.adjust_flow(t,s, min_cap); + if (g.get_flow(s,t) <= 0) + { + parent[t] = nil; + orphans.push_back(t); + } + + t = s; + } + + // trace back towards the sink + s = sink_side; + while (s != sink) + { + t = parent[s]; + g.adjust_flow(t,s, min_cap); + if (g.get_flow(s,t) <= 0) + { + parent[s] = nil; + orphans.push_back(s); + } + s = t; + } + } + + template + bool grow ( + flow_graph& g, + unsigned long& source_side, + unsigned long& sink_side, + typename flow_graph::in_edge_iterator& in_begin, + typename flow_graph::out_edge_iterator& out_begin + ) const + /*! + ensures + - if (an augmenting path was found) then + - returns true + - (#source_side, #sink_side) == the point where the two trees meet. + #source_side is part of the source tree and #sink_side is part of + the sink tree. + - else + - returns false + !*/ + { + typedef typename flow_graph::edge_type edge_type; + typedef typename flow_graph::out_edge_iterator out_edge_iterator; + typedef typename flow_graph::in_edge_iterator in_edge_iterator; + + + while (active.size() != 0) + { + // pick an active node + const unsigned long A = active[0]; + + const unsigned char label_A = g.get_label(A); + + // process its neighbors + if (label_A == SOURCE_CUT) + { + const out_edge_iterator out_end = g.out_end(A); + for(out_edge_iterator& i = out_begin; i != out_end; ++i) + { + if (g.get_flow(i) > 0) + { + const unsigned long id = g.node_id(i); + const unsigned char label_i = g.get_label(id); + if (label_i == FREE_NODE) + { + active.push_back(id); + g.set_label(id,SOURCE_CUT); + parent[id] = A; + ts[id] = ts[A]; + dist[id] = dist[A] + 1; + } + else if (label_A != label_i) + { + source_side = A; + sink_side = id; + return true; + } + else if (is_closer(A, id)) + { + parent[id] = A; + ts[id] = ts[A]; + dist[id] = dist[A] + 1; + } + } + } + } + else if (label_A == SINK_CUT) + { + const in_edge_iterator in_end = g.in_end(A); + for(in_edge_iterator& i = in_begin; i != in_end; ++i) + { + if (g.get_flow(i) > 0) + { + const unsigned long id = g.node_id(i); + const unsigned char label_i = g.get_label(id); + if (label_i == FREE_NODE) + { + active.push_back(id); + g.set_label(id,SINK_CUT); + parent[id] = A; + ts[id] = ts[A]; + dist[id] = dist[A] + 1; + } + else if (label_A != label_i) + { + sink_side = A; + source_side = id; + return true; + } + else if (is_closer(A, id)) + { + parent[id] = A; + ts[id] = ts[A]; + dist[id] = dist[A] + 1; + } + } + } + } + + active.pop_front(); + if (active.size() != 0) + { + in_begin = g.in_begin(active[0]); + out_begin = g.out_begin(active[0]); + } + } + + return false; + } + + inline bool is_closer ( + unsigned long p, + unsigned long q + ) const + { + // return true if p is closer to a terminal than q + return ts[q] <= ts[p] && dist[q] > dist[p]; + } + + mutable std::vector dist; + mutable std::vector ts; + mutable uint32 time; + mutable std::vector parent; + + mutable std::deque active; + mutable std::deque orphans; + }; + +// ---------------------------------------------------------------------------------------- + +} + +// ---------------------------------------------------------------------------------------- + +#endif // DLIB_MIN_CuT_H__ + diff --git a/dlib/graph_cuts/min_cut_abstract.h b/dlib/graph_cuts/min_cut_abstract.h new file mode 100644 index 000000000..3b7869f1d --- /dev/null +++ b/dlib/graph_cuts/min_cut_abstract.h @@ -0,0 +1,472 @@ +// Copyright (C) 2012 Davis E. King (davis@dlib.net) +// License: Boost Software License See LICENSE.txt for the full license. +#undef DLIB_MIN_CuT_ABSTRACT_H__ +#ifdef DLIB_MIN_CuT_ABSTRACT_H__ + +#include "../graph_utils.h" + +// ---------------------------------------------------------------------------------------- + +namespace dlib +{ + + typedef unsigned char node_label; + +// ---------------------------------------------------------------------------------------- + + template < + typename directed_graph_type + > + class flow_graph + { + /*! + WHAT THIS OBJECT REPRESENTS + This object represents a flow capacity graph for use with the + min_cut algorithm defined below. In particular, this object + is a kind of directed graph where the edge weights specify the + flow capacities. + + Note that there is no dlib::flow_graph object. What you are + looking at here is simply the interface definition for a graph + which can be used with the min_cut algorithm. You must implement + your own version of this object for the graph you wish to work with + and then pass it to the min_cut::operator() routine. + + It's also worth pointing out that this graph has symmetric edge + connections. That is, if there is an edge from node A to node B + then there must also be an edge from node B to node A. + !*/ + + public: + + class out_edge_iterator + { + /*! + WHAT THIS OBJECT REPRESENTS + This is a simple forward iterator for iterating over the neighbors + of a node in the graph. It also represents the fact that the neighbors + are on the end of an outgoing edge. That is, the edge represents + the amount of flow which can flow towards the neighbor. + !*/ + + public: + out_edge_iterator( + ); + /*! + ensures + - constructs an iterator in an undefined state. It can't + be used until assigned with a valid iterator. + !*/ + + out_edge_iterator( + const out_edge_iterator& item + ); + /*! + ensures + - #*this is a copy of item + !*/ + + out_edge_iterator& operator=( + const out_edge_iterator& item + ); + /*! + ensures + - #*this is a copy of item + - returns #*this + !*/ + + bool operator!= ( + const out_edge_iterator& item + ) const; + /*! + requires + - *this and item are iterators over the neighbors for the + same node. + ensures + - returns false if *this and item both reference the same + node in the graph and true otherwise. + !*/ + + out_edge_iterator& operator++( + ); + /*! + ensures + - advances *this to the next neighbor node. + - returns a reference to the updated *this + (i.e. this is the ++object form of the increment operator) + !*/ + }; + + class in_edge_iterator + { + /*! + WHAT THIS OBJECT REPRESENTS + This is a simple forward iterator for iterating over the neighbors + of a node in the graph. It also represents the fact that the neighbors + are on the end of an incoming edge. That is, the edge represents + the amount of flow which can flow out of the neighbor node. + !*/ + + public: + + in_edge_iterator( + ); + /*! + ensures + - constructs an iterator in an undefined state. It can't + be used until assigned with a valid iterator. + !*/ + + in_edge_iterator( + const in_edge_iterator& item + ); + /*! + ensures + - #*this is a copy of item + !*/ + + in_edge_iterator& operator=( + const in_edge_iterator& item + ); + /*! + ensures + - #*this is a copy of item + - returns #*this + !*/ + + bool operator!= ( + const in_edge_iterator& item + ) const; + /*! + requires + - *this and item are iterators over the neighbors for the + same node. + ensures + - returns false if *this and item both reference the same + node in the graph and true otherwise. + !*/ + + in_edge_iterator& operator++( + ); + /*! + ensures + - advances *this to the next neighbor node. + - returns a reference to the updated *this + (i.e. this is the ++object form of the increment operator) + !*/ + }; + + + + unsigned long number_of_nodes ( + ) const; + /*! + ensures + - returns the number of nodes in the graph. + !*/ + + out_edge_iterator out_begin( + const unsigned long& idx + ) const; + /*! + requires + - idx < number_of_nodes() + ensures + - returns an iterator pointing to the first neighboring node of + the idx-th node. If no such node exists then returns out_end(idx). + - The returned iterator also represents the directed edge going from + node idx to the neighbor. + !*/ + + in_edge_iterator in_begin( + const unsigned long& idx + ) const; + /*! + requires + - idx < number_of_nodes() + ensures + - returns an iterator pointing to the first neighboring node of + the idx-th node. If no such node exists then returns in_end(idx). + - The returned iterator also represents the directed edge going from + the neighbor to node idx. + !*/ + + out_edge_iterator out_end( + const unsigned long& idx + ) const; + /*! + requires + - idx < number_of_nodes() + ensures + - returns an iterator to one past the last neighboring node of + the idx-th node. + !*/ + + in_edge_iterator in_end( + const unsigned long& idx + ) const; + /*! + requires + - idx < number_of_nodes() + ensures + - returns an iterator to one past the last neighboring node of + the idx-th node. + !*/ + + + unsigned long node_id ( + const out_edge_iterator& it + ) const; + /*! + requires + - it == a valid iterator (i.e. it must be in the range [out_begin(idx), out_end(idx)) + for some valid idx) + ensures + - returns a number IDX such that: + - 0 <= IDX < number_of_nodes() + - IDX == The index which uniquely identifies the node pointed to by the + iterator it. This number can be used with any member function in this + object which expect a node index. (e.g. get_label(IDX) == the label for the + node pointed to by it) + !*/ + + unsigned long node_id ( + const in_edge_iterator& it + ) const; + /*! + requires + - it == a valid iterator (i.e. it must be in the range [in_begin(idx), in_end(idx)) + for some valid idx) + ensures + - returns a number IDX such that: + - 0 <= IDX < number_of_nodes() + - IDX == The index which uniquely identifies the node pointed to by the + iterator it. This number can be used with any member function in this + object which expect a node index. (e.g. get_label(IDX) == the label for the + node pointed to by it) + !*/ + + // This typedef should be for a type like int or double. It + // must also be capable of representing signed values. + typedef an_integer_or_real_type edge_type; + + edge_type get_flow ( + const unsigned long& idx1, + const unsigned long& idx2 + ) const; + /*! + requires + - idx1 < number_of_nodes() + - idx2 < number_of_nodes() + - idx1 and idx2 are neighbors in the graph + ensures + - returns the residual flow capacity from the idx1-th + node to the idx2-th node. + !*/ + + edge_type get_flow ( + const out_edge_iterator& it + ) const; + /*! + requires + - it == a valid iterator (i.e. it must be in the range [out_begin(idx), out_end(idx)) + for some valid idx) + ensures + - let IDX = node_id(it) + - it represents the directed edge from a node, call it H, to the node IDX. Therefore, + this function returns get_flow(H,IDX) + !*/ + + edge_type get_flow ( + const in_edge_iterator& it + ) const; + /*! + requires + - it == a valid iterator (i.e. it must be in the range [in_begin(idx), in_end(idx)) + for some valid idx) + ensures + - let IDX = node_id(it) + - it represents the directed edge from node IDX to another node, call it H. Therefore, + this function returns get_flow(IDX,H) + !*/ + + void adjust_flow ( + const unsigned long& idx1, + const unsigned long& idx2, + const edge_type& value + ); + /*! + requires + - idx1 < number_of_nodes() + - idx2 < number_of_nodes() + - idx1 and idx2 are neighbors in the graph + ensures + - #get_flow(idx1,idx2) == get_flow(idx1,idx2) + value + - #get_flow(idx2,idx1) == get_flow(idx2,idx1) - value + !*/ + + void set_label ( + const unsigned long& idx, + node_label value + ); + /*! + requires + - idx < number_of_nodes() + ensures + - #get_label(idx) == value + !*/ + + node_label get_label ( + const unsigned long& idx + ) const; + /*! + requires + - idx < number_of_nodes() + ensures + - returns the label for the idx-th node in the graph. + !*/ + + }; + +// ---------------------------------------------------------------------------------------- + + const node_label SOURCE_CUT = 0; + const node_label SINK_CUT = 1; + const node_label FREE_NODE = 2; + +// ---------------------------------------------------------------------------------------- + + template < + typename flow_graph + > + typename flow_graph::edge_type graph_cut_score ( + const flow_graph& g + ); + /*! + requires + - flow_graph == an object with an interface compatible with the flow_graph + object defined at the top of this file, or, an implementation of + dlib/directed_graph/directed_graph_kernel_abstract.h. + ensures + - returns the sum of the outgoing edges from nodes with a label of SOURCE_CUT + to nodes with a label != SOURCE_CUT. Note that for a directed_graph object, + the labels are stored in the node's data field. + !*/ + +// ---------------------------------------------------------------------------------------- + + class min_cut + { + /*! + WHAT THIS OBJECT REPRESENTS + This is a function object which can be used to find the min cut + on a graph. + + The implementation is based on the method described in the following + paper: + An Experimental Comparison of Min-Cut/Max-Flow Algorithms for + Energy Minimization in Vision, by Yuri Boykov and Vladimir Kolmogorov, + in PAMI 2004. + + !*/ + + public: + + min_cut( + ); + /*! + ensures + - this object is properly initialized + !*/ + + template < + typename flow_graph + > + void operator() ( + flow_graph& g, + const unsigned long source_node, + const unsigned long sink_node + ) const; + /*! + requires + - flow_graph == an object with an interface compatible with the flow_graph + object defined at the top of this file. + - source_node != sink_node + - source_node < g.number_of_nodes() + - sink_node < g.number_of_nodes() + - for all valid i and j: + - g.get_flow(i,j) >= 0 + (i.e. all the flow capacities/edge weights are positive) + - g does not contain any self loops. That is, no nodes are neighbors with + themselves. + ensures + - Finds the minimum cut on the given graph. That is, this function finds + a labeling of nodes in g such that graph_cut_score(g) would be minimized. Note + that the flow values in #g are modified by this algorithm so if you want + to obtain the min cut score you must call min_cut::operator(), then copy + the flow values back into #g, and then call graph_cut_score(#g). But in most + cases you don't care about the value of the min cut score, rather, you + just want the labels in #g. + - #g.get_label(source_node) == SOURCE_CUT + - #g.get_label(sink_node) == SINK_CUT + - for all valid i: + - #g.get_label(i) == SOURCE_CUT, SINK_CUT, or FREE_NODE + - if (#g.get_label(i) == SOURCE_CUT) then + - The minimum cut of g places node i into the source side of the cut. + - if (#g.get_label(i) == SINK_CUT) then + - The minimum cut of g places node i into the sink side of the cut. + - if (#g.get_label(i) == FREE_NODE) then + - Node i can be labeled SOURCE_CUT or SINK_CUT. Both labelings + result in the same cut score. + - When interpreting g as a graph of flow capacities from the source_node + to the sink_node we can say that the min cut problem is equivalent to + the max flow problem. This equivalent problem is to find out how to push + as much "flow" from the source node to the sink node as possible. + Upon termination, #g will contain the final flow residuals in addition to + the graph cut labels. That is, for all valid i and j: + - #g.get_flow(i,j) == the residual flow capacity left after the max + possible amount of flow is passing from the source node to the sink + node. For example, this means that #g.get_flow(i,j) == 0 whenever + node i and j are in different cuts. + - #g.get_flow(i,j) >= 0 + !*/ + + template < + typename directed_graph + > + void operator() ( + directed_graph& g, + const unsigned long source_node, + const unsigned long sink_node + ) const; + /*! + requires + - directed_graph == an implementation of dlib/directed_graph/directed_graph_kernel_abstract.h + - directed_graph::type == node_label + - directed_graph::edge_type == and integer or double type + - source_node != sink_node + - source_node < g.number_of_nodes() + - sink_node < g.number_of_nodes() + - for all valid i and j: + - edge(g,i,j) >= 0 + (i.e. all the flow capacities/edge weights are positive) + - graph_contains_length_one_cycle(g) == false + - graph_has_symmetric_edges(g) == true + ensures + - This routine simply converts g into a flow graph and calls the version + of operator() defined above. Note that the conversion is done in O(1) + time, it's just an interface adaptor. + - edge weights in g correspond to network flows while the .data field of + each node in g corresponds to the graph node labels. + - upon termination, the flows and labels in g will have been modified + as described in the above operator() routine. + !*/ + }; + +// ---------------------------------------------------------------------------------------- + +} + +#endif // DLIB_MIN_CuT_ABSTRACT_H__ + +