boost/graph/boykov_kolmogorov_max_flow.hpp
// Copyright (c) 2006, Stephan Diederich // // This code may be used under either of the following two licences: // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without // restriction, including without limitation the rights to use, // copy, modify, merge, publish, distribute, sublicense, and/or // sell copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: // // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. // // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND // NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT // HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, // WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR // OTHER DEALINGS IN THE SOFTWARE. OF SUCH DAMAGE. // // Or: // // Distributed under the Boost Software License, Version 1.0. // (See accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) #ifndef BOOST_BOYKOV_KOLMOGOROV_MAX_FLOW_HPP #define BOOST_BOYKOV_KOLMOGOROV_MAX_FLOW_HPP #include <boost/config.hpp> #include <boost/assert.hpp> #include <vector> #include <list> #include <utility> #include <iosfwd> #include <algorithm> // for std::min and std::max #include <boost/pending/queue.hpp> #include <boost/limits.hpp> #include <boost/property_map/property_map.hpp> #include <boost/none_t.hpp> #include <boost/graph/graph_concepts.hpp> #include <boost/graph/named_function_params.hpp> #include <boost/graph/lookup_edge.hpp> #include <boost/concept/assert.hpp> // The algorithm impelemented here is described in: // // Boykov, Y., Kolmogorov, V. "An Experimental Comparison of Min-Cut/Max-Flow // Algorithms for Energy Minimization in Vision", In IEEE Transactions on // Pattern Analysis and Machine Intelligence, vol. 26, no. 9, pp. 1124-1137, // Sep 2004. // // For further reading, also see: // // Kolmogorov, V. "Graph Based Algorithms for Scene Reconstruction from Two or // More Views". PhD thesis, Cornell University, Sep 2003. namespace boost { namespace detail { template < class Graph, class EdgeCapacityMap, class ResidualCapacityEdgeMap, class ReverseEdgeMap, class PredecessorMap, class ColorMap, class DistanceMap, class IndexMap > class bk_max_flow { typedef typename property_traits< EdgeCapacityMap >::value_type tEdgeVal; typedef graph_traits< Graph > tGraphTraits; typedef typename tGraphTraits::vertex_iterator vertex_iterator; typedef typename tGraphTraits::vertex_descriptor vertex_descriptor; typedef typename tGraphTraits::edge_descriptor edge_descriptor; typedef typename tGraphTraits::edge_iterator edge_iterator; typedef typename tGraphTraits::out_edge_iterator out_edge_iterator; typedef boost::queue< vertex_descriptor > tQueue; // queue of vertices, used in adoption-stage typedef typename property_traits< ColorMap >::value_type tColorValue; typedef color_traits< tColorValue > tColorTraits; typedef typename property_traits< DistanceMap >::value_type tDistanceVal; public: bk_max_flow(Graph& g, EdgeCapacityMap cap, ResidualCapacityEdgeMap res, ReverseEdgeMap rev, PredecessorMap pre, ColorMap color, DistanceMap dist, IndexMap idx, vertex_descriptor src, vertex_descriptor sink) : m_g(g) , m_index_map(idx) , m_cap_map(cap) , m_res_cap_map(res) , m_rev_edge_map(rev) , m_pre_map(pre) , m_tree_map(color) , m_dist_map(dist) , m_source(src) , m_sink(sink) , m_active_nodes() , m_in_active_list_vec(num_vertices(g), false) , m_in_active_list_map(make_iterator_property_map( m_in_active_list_vec.begin(), m_index_map)) , m_has_parent_vec(num_vertices(g), false) , m_has_parent_map( make_iterator_property_map(m_has_parent_vec.begin(), m_index_map)) , m_time_vec(num_vertices(g), 0) , m_time_map( make_iterator_property_map(m_time_vec.begin(), m_index_map)) , m_flow(0) , m_time(1) , m_last_grow_vertex(graph_traits< Graph >::null_vertex()) { // initialize the color-map with gray-values vertex_iterator vi, v_end; for (boost::tie(vi, v_end) = vertices(m_g); vi != v_end; ++vi) { set_tree(*vi, tColorTraits::gray()); } // Initialize flow to zero which means initializing // the residual capacity equal to the capacity edge_iterator ei, e_end; for (boost::tie(ei, e_end) = edges(m_g); ei != e_end; ++ei) { put(m_res_cap_map, *ei, get(m_cap_map, *ei)); BOOST_ASSERT(get(m_rev_edge_map, get(m_rev_edge_map, *ei)) == *ei); // check if the reverse edge map is build up // properly } // init the search trees with the two terminals set_tree(m_source, tColorTraits::black()); set_tree(m_sink, tColorTraits::white()); put(m_time_map, m_source, 1); put(m_time_map, m_sink, 1); } tEdgeVal max_flow() { // augment direct paths from SOURCE->SINK and SOURCE->VERTEX->SINK augment_direct_paths(); // start the main-loop while (true) { bool path_found; edge_descriptor connecting_edge; boost::tie(connecting_edge, path_found) = grow(); // find a path from source to sink if (!path_found) { // we're finished, no more paths were found break; } ++m_time; augment(connecting_edge); // augment that path adopt(); // rebuild search tree structure } return m_flow; } // the complete class is protected, as we want access to members in // derived test-class (see test/boykov_kolmogorov_max_flow_test.cpp) protected: void augment_direct_paths() { // in a first step, we augment all direct paths from // source->NODE->sink and additionally paths from source->sink. This // improves especially graphcuts for segmentation, as most of the // nodes have source/sink connects but shouldn't have an impact on // other maxflow problems (this is done in grow() anyway) out_edge_iterator ei, e_end; for (boost::tie(ei, e_end) = out_edges(m_source, m_g); ei != e_end; ++ei) { edge_descriptor from_source = *ei; vertex_descriptor current_node = target(from_source, m_g); if (current_node == m_sink) { tEdgeVal cap = get(m_res_cap_map, from_source); put(m_res_cap_map, from_source, 0); m_flow += cap; continue; } edge_descriptor to_sink; bool is_there; boost::tie(to_sink, is_there) = lookup_edge(current_node, m_sink, m_g); if (is_there) { tEdgeVal cap_from_source = get(m_res_cap_map, from_source); tEdgeVal cap_to_sink = get(m_res_cap_map, to_sink); if (cap_from_source > cap_to_sink) { set_tree(current_node, tColorTraits::black()); add_active_node(current_node); set_edge_to_parent(current_node, from_source); put(m_dist_map, current_node, 1); put(m_time_map, current_node, 1); // add stuff to flow and update residuals. we dont need // to update reverse_edges, as incoming/outgoing edges // to/from source/sink don't count for max-flow put(m_res_cap_map, from_source, get(m_res_cap_map, from_source) - cap_to_sink); put(m_res_cap_map, to_sink, 0); m_flow += cap_to_sink; } else if (cap_to_sink > 0) { set_tree(current_node, tColorTraits::white()); add_active_node(current_node); set_edge_to_parent(current_node, to_sink); put(m_dist_map, current_node, 1); put(m_time_map, current_node, 1); // add stuff to flow and update residuals. we dont need // to update reverse_edges, as incoming/outgoing edges // to/from source/sink don't count for max-flow put(m_res_cap_map, to_sink, get(m_res_cap_map, to_sink) - cap_from_source); put(m_res_cap_map, from_source, 0); m_flow += cap_from_source; } } else if (get(m_res_cap_map, from_source)) { // there is no sink connect, so we can't augment this path, // but to avoid adding m_source to the active nodes, we just // activate this node and set the approciate things set_tree(current_node, tColorTraits::black()); set_edge_to_parent(current_node, from_source); put(m_dist_map, current_node, 1); put(m_time_map, current_node, 1); add_active_node(current_node); } } for (boost::tie(ei, e_end) = out_edges(m_sink, m_g); ei != e_end; ++ei) { edge_descriptor to_sink = get(m_rev_edge_map, *ei); vertex_descriptor current_node = source(to_sink, m_g); if (get(m_res_cap_map, to_sink)) { set_tree(current_node, tColorTraits::white()); set_edge_to_parent(current_node, to_sink); put(m_dist_map, current_node, 1); put(m_time_map, current_node, 1); add_active_node(current_node); } } } /** * Returns a pair of an edge and a boolean. if the bool is true, the * edge is a connection of a found path from s->t , read "the link" and * source(returnVal, m_g) is the end of the path found in the * source-tree target(returnVal, m_g) is the beginning of the path found * in the sink-tree */ std::pair< edge_descriptor, bool > grow() { BOOST_ASSERT(m_orphans.empty()); vertex_descriptor current_node; while ((current_node = get_next_active_node()) != graph_traits< Graph >::null_vertex()) { // if there is one BOOST_ASSERT(get_tree(current_node) != tColorTraits::gray() && (has_parent(current_node) || current_node == m_source || current_node == m_sink)); if (get_tree(current_node) == tColorTraits::black()) { // source tree growing out_edge_iterator ei, e_end; if (current_node != m_last_grow_vertex) { m_last_grow_vertex = current_node; boost::tie(m_last_grow_edge_it, m_last_grow_edge_end) = out_edges(current_node, m_g); } for (; m_last_grow_edge_it != m_last_grow_edge_end; ++m_last_grow_edge_it) { edge_descriptor out_edge = *m_last_grow_edge_it; if (get(m_res_cap_map, out_edge) > 0) { // check if we have capacity left on this edge vertex_descriptor other_node = target(out_edge, m_g); if (get_tree(other_node) == tColorTraits::gray()) { // it's a free node set_tree(other_node, tColorTraits::black()); // aquire other node // to our search // tree set_edge_to_parent( other_node, out_edge); // set us as parent put(m_dist_map, other_node, get(m_dist_map, current_node) + 1); // and update the // distance-heuristic put(m_time_map, other_node, get(m_time_map, current_node)); add_active_node(other_node); } else if (get_tree(other_node) == tColorTraits::black()) { // we do this to get shorter paths. check if we // are nearer to the source as its parent is if (is_closer_to_terminal( current_node, other_node)) { set_edge_to_parent(other_node, out_edge); put(m_dist_map, other_node, get(m_dist_map, current_node) + 1); put(m_time_map, other_node, get(m_time_map, current_node)); } } else { BOOST_ASSERT(get_tree(other_node) == tColorTraits::white()); // kewl, found a path from one to the other // search tree, return // the connecting edge in src->sink dir return std::make_pair(out_edge, true); } } } // for all out-edges } // source-tree-growing else { BOOST_ASSERT( get_tree(current_node) == tColorTraits::white()); out_edge_iterator ei, e_end; if (current_node != m_last_grow_vertex) { m_last_grow_vertex = current_node; boost::tie(m_last_grow_edge_it, m_last_grow_edge_end) = out_edges(current_node, m_g); } for (; m_last_grow_edge_it != m_last_grow_edge_end; ++m_last_grow_edge_it) { edge_descriptor in_edge = get(m_rev_edge_map, *m_last_grow_edge_it); if (get(m_res_cap_map, in_edge) > 0) { // check if there is capacity left vertex_descriptor other_node = source(in_edge, m_g); if (get_tree(other_node) == tColorTraits::gray()) { // it's a free node set_tree(other_node, tColorTraits::white()); // aquire that node // to our search // tree set_edge_to_parent( other_node, in_edge); // set us as parent add_active_node( other_node); // activate that node put(m_dist_map, other_node, get(m_dist_map, current_node) + 1); // set its distance put(m_time_map, other_node, get(m_time_map, current_node)); // and time } else if (get_tree(other_node) == tColorTraits::white()) { if (is_closer_to_terminal( current_node, other_node)) { // we are closer to the sink than its parent // is, so we "adopt" him set_edge_to_parent(other_node, in_edge); put(m_dist_map, other_node, get(m_dist_map, current_node) + 1); put(m_time_map, other_node, get(m_time_map, current_node)); } } else { BOOST_ASSERT(get_tree(other_node) == tColorTraits::black()); // kewl, found a path from one to the other // search tree, // return the connecting edge in src->sink dir return std::make_pair(in_edge, true); } } } // for all out-edges } // sink-tree growing // all edges of that node are processed, and no more paths were // found. // remove if from the front of the active queue finish_node(current_node); } // while active_nodes not empty // no active nodes anymore and no path found, we're done return std::make_pair(edge_descriptor(), false); } /** * augments path from s->t and updates residual graph * source(e, m_g) is the end of the path found in the source-tree * target(e, m_g) is the beginning of the path found in the sink-tree * this phase generates orphans on satured edges, if the attached verts * are from different search-trees orphans are ordered in distance to * sink/source. first the farest from the source are front_inserted into * the orphans list, and after that the sink-tree-orphans are * front_inserted. when going to adoption stage the orphans are * popped_front, and so we process the nearest verts to the terminals * first */ void augment(edge_descriptor e) { BOOST_ASSERT(get_tree(target(e, m_g)) == tColorTraits::white()); BOOST_ASSERT(get_tree(source(e, m_g)) == tColorTraits::black()); BOOST_ASSERT(m_orphans.empty()); const tEdgeVal bottleneck = find_bottleneck(e); // now we push the found flow through the path // for each edge we saturate we have to look for the verts that // belong to that edge, one of them becomes an orphans now process // the connecting edge put(m_res_cap_map, e, get(m_res_cap_map, e) - bottleneck); BOOST_ASSERT(get(m_res_cap_map, e) >= 0); put(m_res_cap_map, get(m_rev_edge_map, e), get(m_res_cap_map, get(m_rev_edge_map, e)) + bottleneck); // now we follow the path back to the source vertex_descriptor current_node = source(e, m_g); while (current_node != m_source) { edge_descriptor pred = get_edge_to_parent(current_node); put(m_res_cap_map, pred, get(m_res_cap_map, pred) - bottleneck); BOOST_ASSERT(get(m_res_cap_map, pred) >= 0); put(m_res_cap_map, get(m_rev_edge_map, pred), get(m_res_cap_map, get(m_rev_edge_map, pred)) + bottleneck); if (get(m_res_cap_map, pred) == 0) { set_no_parent(current_node); m_orphans.push_front(current_node); } current_node = source(pred, m_g); } // then go forward in the sink-tree current_node = target(e, m_g); while (current_node != m_sink) { edge_descriptor pred = get_edge_to_parent(current_node); put(m_res_cap_map, pred, get(m_res_cap_map, pred) - bottleneck); BOOST_ASSERT(get(m_res_cap_map, pred) >= 0); put(m_res_cap_map, get(m_rev_edge_map, pred), get(m_res_cap_map, get(m_rev_edge_map, pred)) + bottleneck); if (get(m_res_cap_map, pred) == 0) { set_no_parent(current_node); m_orphans.push_front(current_node); } current_node = target(pred, m_g); } // and add it to the max-flow m_flow += bottleneck; } /** * returns the bottleneck of a s->t path (end_of_path is last vertex in * source-tree, begin_of_path is first vertex in sink-tree) */ inline tEdgeVal find_bottleneck(edge_descriptor e) { BOOST_USING_STD_MIN(); tEdgeVal minimum_cap = get(m_res_cap_map, e); vertex_descriptor current_node = source(e, m_g); // first go back in the source tree while (current_node != m_source) { edge_descriptor pred = get_edge_to_parent(current_node); minimum_cap = min BOOST_PREVENT_MACRO_SUBSTITUTION( minimum_cap, get(m_res_cap_map, pred)); current_node = source(pred, m_g); } // then go forward in the sink-tree current_node = target(e, m_g); while (current_node != m_sink) { edge_descriptor pred = get_edge_to_parent(current_node); minimum_cap = min BOOST_PREVENT_MACRO_SUBSTITUTION( minimum_cap, get(m_res_cap_map, pred)); current_node = target(pred, m_g); } return minimum_cap; } /** * rebuild search trees * empty the queue of orphans, and find new parents for them or just * drop them from the search trees */ void adopt() { while (!m_orphans.empty() || !m_child_orphans.empty()) { vertex_descriptor current_node; if (m_child_orphans.empty()) { // get the next orphan from the main-queue and remove it current_node = m_orphans.front(); m_orphans.pop_front(); } else { current_node = m_child_orphans.front(); m_child_orphans.pop(); } if (get_tree(current_node) == tColorTraits::black()) { // we're in the source-tree tDistanceVal min_distance = (std::numeric_limits< tDistanceVal >::max)(); edge_descriptor new_parent_edge; out_edge_iterator ei, e_end; for (boost::tie(ei, e_end) = out_edges(current_node, m_g); ei != e_end; ++ei) { const edge_descriptor in_edge = get(m_rev_edge_map, *ei); BOOST_ASSERT(target(in_edge, m_g) == current_node); // we should be the target of this // edge if (get(m_res_cap_map, in_edge) > 0) { vertex_descriptor other_node = source(in_edge, m_g); if (get_tree(other_node) == tColorTraits::black() && has_source_connect(other_node)) { if (get(m_dist_map, other_node) < min_distance) { min_distance = get(m_dist_map, other_node); new_parent_edge = in_edge; } } } } if (min_distance != (std::numeric_limits< tDistanceVal >::max)()) { set_edge_to_parent(current_node, new_parent_edge); put(m_dist_map, current_node, min_distance + 1); put(m_time_map, current_node, m_time); } else { put(m_time_map, current_node, 0); for (boost::tie(ei, e_end) = out_edges(current_node, m_g); ei != e_end; ++ei) { edge_descriptor in_edge = get(m_rev_edge_map, *ei); vertex_descriptor other_node = source(in_edge, m_g); if (get_tree(other_node) == tColorTraits::black() && other_node != m_source) { if (get(m_res_cap_map, in_edge) > 0) { add_active_node(other_node); } if (has_parent(other_node) && source( get_edge_to_parent(other_node), m_g) == current_node) { // we are the parent of that node // it has to find a new parent, too set_no_parent(other_node); m_child_orphans.push(other_node); } } } set_tree(current_node, tColorTraits::gray()); } // no parent found } // source-tree-adoption else { // now we should be in the sink-tree, check that... BOOST_ASSERT( get_tree(current_node) == tColorTraits::white()); out_edge_iterator ei, e_end; edge_descriptor new_parent_edge; tDistanceVal min_distance = (std::numeric_limits< tDistanceVal >::max)(); for (boost::tie(ei, e_end) = out_edges(current_node, m_g); ei != e_end; ++ei) { const edge_descriptor out_edge = *ei; if (get(m_res_cap_map, out_edge) > 0) { const vertex_descriptor other_node = target(out_edge, m_g); if (get_tree(other_node) == tColorTraits::white() && has_sink_connect(other_node)) if (get(m_dist_map, other_node) < min_distance) { min_distance = get(m_dist_map, other_node); new_parent_edge = out_edge; } } } if (min_distance != (std::numeric_limits< tDistanceVal >::max)()) { set_edge_to_parent(current_node, new_parent_edge); put(m_dist_map, current_node, min_distance + 1); put(m_time_map, current_node, m_time); } else { put(m_time_map, current_node, 0); for (boost::tie(ei, e_end) = out_edges(current_node, m_g); ei != e_end; ++ei) { const edge_descriptor out_edge = *ei; const vertex_descriptor other_node = target(out_edge, m_g); if (get_tree(other_node) == tColorTraits::white() && other_node != m_sink) { if (get(m_res_cap_map, out_edge) > 0) { add_active_node(other_node); } if (has_parent(other_node) && target( get_edge_to_parent(other_node), m_g) == current_node) { // we were it's parent, so it has to find a // new one, too set_no_parent(other_node); m_child_orphans.push(other_node); } } } set_tree(current_node, tColorTraits::gray()); } // no parent found } // sink-tree adoption } // while !orphans.empty() } // adopt /** * return next active vertex if there is one, otherwise a null_vertex */ inline vertex_descriptor get_next_active_node() { while (true) { if (m_active_nodes.empty()) return graph_traits< Graph >::null_vertex(); vertex_descriptor v = m_active_nodes.front(); // if it has no parent, this node can't be active (if its not // source or sink) if (!has_parent(v) && v != m_source && v != m_sink) { m_active_nodes.pop(); put(m_in_active_list_map, v, false); } else { BOOST_ASSERT(get_tree(v) == tColorTraits::black() || get_tree(v) == tColorTraits::white()); return v; } } } /** * adds v as an active vertex, but only if its not in the list already */ inline void add_active_node(vertex_descriptor v) { BOOST_ASSERT(get_tree(v) != tColorTraits::gray()); if (get(m_in_active_list_map, v)) { if (m_last_grow_vertex == v) { m_last_grow_vertex = graph_traits< Graph >::null_vertex(); } return; } else { put(m_in_active_list_map, v, true); m_active_nodes.push(v); } } /** * finish_node removes a node from the front of the active queue (its * called in grow phase, if no more paths can be found using this node) */ inline void finish_node(vertex_descriptor v) { BOOST_ASSERT(m_active_nodes.front() == v); m_active_nodes.pop(); put(m_in_active_list_map, v, false); m_last_grow_vertex = graph_traits< Graph >::null_vertex(); } /** * removes a vertex from the queue of active nodes (actually this does * nothing, but checks if this node has no parent edge, as this is the * criteria for being no more active) */ inline void remove_active_node(vertex_descriptor v) { BOOST_ASSERT(!has_parent(v)); } /** * returns the search tree of v; tColorValue::black() for source tree, * white() for sink tree, gray() for no tree */ inline tColorValue get_tree(vertex_descriptor v) const { return get(m_tree_map, v); } /** * sets search tree of v; tColorValue::black() for source tree, white() * for sink tree, gray() for no tree */ inline void set_tree(vertex_descriptor v, tColorValue t) { put(m_tree_map, v, t); } /** * returns edge to parent vertex of v; */ inline edge_descriptor get_edge_to_parent(vertex_descriptor v) const { return get(m_pre_map, v); } /** * returns true if the edge stored in m_pre_map[v] is a valid entry */ inline bool has_parent(vertex_descriptor v) const { return get(m_has_parent_map, v); } /** * sets edge to parent vertex of v; */ inline void set_edge_to_parent( vertex_descriptor v, edge_descriptor f_edge_to_parent) { BOOST_ASSERT(get(m_res_cap_map, f_edge_to_parent) > 0); put(m_pre_map, v, f_edge_to_parent); put(m_has_parent_map, v, true); } /** * removes the edge to parent of v (this is done by invalidating the * entry an additional map) */ inline void set_no_parent(vertex_descriptor v) { put(m_has_parent_map, v, false); } /** * checks if vertex v has a connect to the sink-vertex (@var m_sink) * @param v the vertex which is checked * @return true if a path to the sink was found, false if not */ inline bool has_sink_connect(vertex_descriptor v) { tDistanceVal current_distance = 0; vertex_descriptor current_vertex = v; while (true) { if (get(m_time_map, current_vertex) == m_time) { // we found a node which was already checked this round. use // it for distance calculations current_distance += get(m_dist_map, current_vertex); break; } if (current_vertex == m_sink) { put(m_time_map, m_sink, m_time); break; } if (has_parent(current_vertex)) { // it has a parent, so get it current_vertex = target(get_edge_to_parent(current_vertex), m_g); ++current_distance; } else { // no path found return false; } } current_vertex = v; while (get(m_time_map, current_vertex) != m_time) { put(m_dist_map, current_vertex, current_distance); --current_distance; put(m_time_map, current_vertex, m_time); current_vertex = target(get_edge_to_parent(current_vertex), m_g); } return true; } /** * checks if vertex v has a connect to the source-vertex (@var m_source) * @param v the vertex which is checked * @return true if a path to the source was found, false if not */ inline bool has_source_connect(vertex_descriptor v) { tDistanceVal current_distance = 0; vertex_descriptor current_vertex = v; while (true) { if (get(m_time_map, current_vertex) == m_time) { // we found a node which was already checked this round. use // it for distance calculations current_distance += get(m_dist_map, current_vertex); break; } if (current_vertex == m_source) { put(m_time_map, m_source, m_time); break; } if (has_parent(current_vertex)) { // it has a parent, so get it current_vertex = source(get_edge_to_parent(current_vertex), m_g); ++current_distance; } else { // no path found return false; } } current_vertex = v; while (get(m_time_map, current_vertex) != m_time) { put(m_dist_map, current_vertex, current_distance); --current_distance; put(m_time_map, current_vertex, m_time); current_vertex = source(get_edge_to_parent(current_vertex), m_g); } return true; } /** * returns true, if p is closer to a terminal than q */ inline bool is_closer_to_terminal( vertex_descriptor p, vertex_descriptor q) { // checks the timestamps first, to build no cycles, and after that // the real distance return (get(m_time_map, q) <= get(m_time_map, p) && get(m_dist_map, q) > get(m_dist_map, p) + 1); } //////// // member vars //////// Graph& m_g; IndexMap m_index_map; EdgeCapacityMap m_cap_map; ResidualCapacityEdgeMap m_res_cap_map; ReverseEdgeMap m_rev_edge_map; PredecessorMap m_pre_map; // stores paths found in the growth stage ColorMap m_tree_map; // maps each vertex into one of the two search tree // or none (gray()) DistanceMap m_dist_map; // stores distance to source/sink nodes vertex_descriptor m_source; vertex_descriptor m_sink; tQueue m_active_nodes; std::vector< bool > m_in_active_list_vec; iterator_property_map< std::vector< bool >::iterator, IndexMap > m_in_active_list_map; std::list< vertex_descriptor > m_orphans; tQueue m_child_orphans; // we use a second queuqe for child orphans, as // they are FIFO processed std::vector< bool > m_has_parent_vec; iterator_property_map< std::vector< bool >::iterator, IndexMap > m_has_parent_map; std::vector< long > m_time_vec; // timestamp of each node, used for // sink/source-path calculations iterator_property_map< std::vector< long >::iterator, IndexMap > m_time_map; tEdgeVal m_flow; long m_time; vertex_descriptor m_last_grow_vertex; out_edge_iterator m_last_grow_edge_it; out_edge_iterator m_last_grow_edge_end; }; } // namespace boost::detail /** * non-named-parameter version, given everything * this is the catch all version */ template < class Graph, class CapacityEdgeMap, class ResidualCapacityEdgeMap, class ReverseEdgeMap, class PredecessorMap, class ColorMap, class DistanceMap, class IndexMap > typename property_traits< CapacityEdgeMap >::value_type boykov_kolmogorov_max_flow(Graph& g, CapacityEdgeMap cap, ResidualCapacityEdgeMap res_cap, ReverseEdgeMap rev_map, PredecessorMap pre_map, ColorMap color, DistanceMap dist, IndexMap idx, typename graph_traits< Graph >::vertex_descriptor src, typename graph_traits< Graph >::vertex_descriptor sink) { typedef typename graph_traits< Graph >::vertex_descriptor vertex_descriptor; typedef typename graph_traits< Graph >::edge_descriptor edge_descriptor; // as this method is the last one before we instantiate the solver, we do // the concept checks here BOOST_CONCEPT_ASSERT( (VertexListGraphConcept< Graph >)); // to have vertices(), // num_vertices(), BOOST_CONCEPT_ASSERT((EdgeListGraphConcept< Graph >)); // to have edges() BOOST_CONCEPT_ASSERT( (IncidenceGraphConcept< Graph >)); // to have source(), target() and // out_edges() BOOST_CONCEPT_ASSERT((ReadablePropertyMapConcept< CapacityEdgeMap, edge_descriptor >)); // read flow-values from edges BOOST_CONCEPT_ASSERT((ReadWritePropertyMapConcept< ResidualCapacityEdgeMap, edge_descriptor >)); // write flow-values to residuals BOOST_CONCEPT_ASSERT((ReadablePropertyMapConcept< ReverseEdgeMap, edge_descriptor >)); // read out reverse edges BOOST_CONCEPT_ASSERT((ReadWritePropertyMapConcept< PredecessorMap, vertex_descriptor >)); // store predecessor there BOOST_CONCEPT_ASSERT((ReadWritePropertyMapConcept< ColorMap, vertex_descriptor >)); // write corresponding tree BOOST_CONCEPT_ASSERT((ReadWritePropertyMapConcept< DistanceMap, vertex_descriptor >)); // write distance to source/sink BOOST_CONCEPT_ASSERT((ReadablePropertyMapConcept< IndexMap, vertex_descriptor >)); // get index 0...|V|-1 BOOST_ASSERT(num_vertices(g) >= 2 && src != sink); detail::bk_max_flow< Graph, CapacityEdgeMap, ResidualCapacityEdgeMap, ReverseEdgeMap, PredecessorMap, ColorMap, DistanceMap, IndexMap > algo(g, cap, res_cap, rev_map, pre_map, color, dist, idx, src, sink); return algo.max_flow(); } /** * non-named-parameter version, given capacity, residual_capacity, * reverse_edges, and an index map. */ template < class Graph, class CapacityEdgeMap, class ResidualCapacityEdgeMap, class ReverseEdgeMap, class IndexMap > typename property_traits< CapacityEdgeMap >::value_type boykov_kolmogorov_max_flow(Graph& g, CapacityEdgeMap cap, ResidualCapacityEdgeMap res_cap, ReverseEdgeMap rev, IndexMap idx, typename graph_traits< Graph >::vertex_descriptor src, typename graph_traits< Graph >::vertex_descriptor sink) { typename graph_traits< Graph >::vertices_size_type n_verts = num_vertices(g); std::vector< typename graph_traits< Graph >::edge_descriptor > predecessor_vec(n_verts); std::vector< default_color_type > color_vec(n_verts); std::vector< typename graph_traits< Graph >::vertices_size_type > distance_vec(n_verts); return boykov_kolmogorov_max_flow(g, cap, res_cap, rev, make_iterator_property_map(predecessor_vec.begin(), idx), make_iterator_property_map(color_vec.begin(), idx), make_iterator_property_map(distance_vec.begin(), idx), idx, src, sink); } /** * non-named-parameter version, some given: capacity, residual_capacity, * reverse_edges, color_map and an index map. Use this if you are interested in * the minimum cut, as the color map provides that info. */ template < class Graph, class CapacityEdgeMap, class ResidualCapacityEdgeMap, class ReverseEdgeMap, class ColorMap, class IndexMap > typename property_traits< CapacityEdgeMap >::value_type boykov_kolmogorov_max_flow(Graph& g, CapacityEdgeMap cap, ResidualCapacityEdgeMap res_cap, ReverseEdgeMap rev, ColorMap color, IndexMap idx, typename graph_traits< Graph >::vertex_descriptor src, typename graph_traits< Graph >::vertex_descriptor sink) { typename graph_traits< Graph >::vertices_size_type n_verts = num_vertices(g); std::vector< typename graph_traits< Graph >::edge_descriptor > predecessor_vec(n_verts); std::vector< typename graph_traits< Graph >::vertices_size_type > distance_vec(n_verts); return boykov_kolmogorov_max_flow(g, cap, res_cap, rev, make_iterator_property_map(predecessor_vec.begin(), idx), color, make_iterator_property_map(distance_vec.begin(), idx), idx, src, sink); } /** * named-parameter version, some given */ template < class Graph, class P, class T, class R > typename detail::edge_capacity_value< Graph, P, T, R >::type boykov_kolmogorov_max_flow(Graph& g, typename graph_traits< Graph >::vertex_descriptor src, typename graph_traits< Graph >::vertex_descriptor sink, const bgl_named_params< P, T, R >& params) { return boykov_kolmogorov_max_flow(g, choose_const_pmap(get_param(params, edge_capacity), g, edge_capacity), choose_pmap(get_param(params, edge_residual_capacity), g, edge_residual_capacity), choose_const_pmap(get_param(params, edge_reverse), g, edge_reverse), choose_pmap( get_param(params, vertex_predecessor), g, vertex_predecessor), choose_pmap(get_param(params, vertex_color), g, vertex_color), choose_pmap(get_param(params, vertex_distance), g, vertex_distance), choose_const_pmap(get_param(params, vertex_index), g, vertex_index), src, sink); } /** * named-parameter version, none given */ template < class Graph > typename property_traits< typename property_map< Graph, edge_capacity_t >::const_type >::value_type boykov_kolmogorov_max_flow(Graph& g, typename graph_traits< Graph >::vertex_descriptor src, typename graph_traits< Graph >::vertex_descriptor sink) { bgl_named_params< int, buffer_param_t > params(0); // bogus empty param return boykov_kolmogorov_max_flow(g, src, sink, params); } } // namespace boost #endif // BOOST_BOYKOV_KOLMOGOROV_MAX_FLOW_HPP