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