git.sesse.net Git - casparcg/blob - dependencies/boost/boost/graph/boykov_kolmogorov_max_flow.hpp

   1 //  Copyright (c) 2006, Stephan Diederich
   2 //
   3 //  This code may be used under either of the following two licences:
   4 //
   5 //    Permission is hereby granted, free of charge, to any person
   6 //    obtaining a copy of this software and associated documentation
   7 //    files (the "Software"), to deal in the Software without
   8 //    restriction, including without limitation the rights to use,
   9 //    copy, modify, merge, publish, distribute, sublicense, and/or
  10 //    sell copies of the Software, and to permit persons to whom the
  11 //    Software is furnished to do so, subject to the following
  12 //    conditions:
  13 //
  14 //    The above copyright notice and this permission notice shall be
  15 //    included in all copies or substantial portions of the Software.
  16 //
  17 //    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  18 //    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
  19 //    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  20 //    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
  21 //    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
  22 //    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  23 //    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  24 //    OTHER DEALINGS IN THE SOFTWARE. OF SUCH DAMAGE.
  25 //
  26 //  Or:
  27 //
  28 //    Distributed under the Boost Software License, Version 1.0.
  29 //    (See accompanying file LICENSE_1_0.txt or copy at
  30 //    http://www.boost.org/LICENSE_1_0.txt)
  31
  32 #ifndef BOOST_BOYKOV_KOLMOGOROV_MAX_FLOW_HPP
  33 #define BOOST_BOYKOV_KOLMOGOROV_MAX_FLOW_HPP
  34
  35 #include <boost/config.hpp>
  36 #include <boost/assert.hpp>
  37 #include <vector>
  38 #include <list>
  39 #include <utility>
  40 #include <iosfwd>
  41 #include <algorithm> // for std::min and std::max
  42
  43 #include <boost/pending/queue.hpp>
  44 #include <boost/limits.hpp>
  45 #include <boost/property_map/property_map.hpp>
  46 #include <boost/none_t.hpp>
  47 #include <boost/graph/graph_concepts.hpp>
  48 #include <boost/graph/named_function_params.hpp>
  49 #include <boost/graph/lookup_edge.hpp>
  50
  51 // The algorithm impelemented here is described in:
  52 //
  53 // Boykov, Y., Kolmogorov, V. "An Experimental Comparison of Min-Cut/Max-Flow
  54 // Algorithms for Energy Minimization in Vision", In IEEE Transactions on
  55 // Pattern Analysis and Machine Intelligence, vol. 26, no. 9, pp. 1124-1137,
  56 // Sep 2004.
  57 //
  58 // For further reading, also see:
  59 //
  60 // Kolmogorov, V. "Graph Based Algorithms for Scene Reconstruction from Two or
  61 // More Views". PhD thesis, Cornell University, Sep 2003.
  62
  63 namespace boost {
  64
  65 namespace detail {
  66
  67 template <class Graph,
  68           class EdgeCapacityMap,
  69           class ResidualCapacityEdgeMap,
  70           class ReverseEdgeMap,
  71           class PredecessorMap,
  72           class ColorMap,
  73           class DistanceMap,
  74           class IndexMap>
  75 class bk_max_flow {
  76   typedef typename property_traits<EdgeCapacityMap>::value_type tEdgeVal;
  77   typedef graph_traits<Graph> tGraphTraits;
  78   typedef typename tGraphTraits::vertex_iterator vertex_iterator;
  79   typedef typename tGraphTraits::vertex_descriptor vertex_descriptor;
  80   typedef typename tGraphTraits::edge_descriptor edge_descriptor;
  81   typedef typename tGraphTraits::edge_iterator edge_iterator;
  82   typedef typename tGraphTraits::out_edge_iterator out_edge_iterator;
  83   typedef boost::queue<vertex_descriptor> tQueue;                               //queue of vertices, used in adoption-stage
  84   typedef typename property_traits<ColorMap>::value_type tColorValue;
  85   typedef color_traits<tColorValue> tColorTraits;
  86   typedef typename property_traits<DistanceMap>::value_type tDistanceVal;
  87
  88     public:
  89       bk_max_flow(Graph& g,
  90                   EdgeCapacityMap cap,
  91                   ResidualCapacityEdgeMap res,
  92                   ReverseEdgeMap rev,
  93                   PredecessorMap pre,
  94                   ColorMap color,
  95                   DistanceMap dist,
  96                   IndexMap idx,
  97                   vertex_descriptor src,
  98                   vertex_descriptor sink):
  99       m_g(g),
 100       m_index_map(idx),
 101       m_cap_map(cap),
 102       m_res_cap_map(res),
 103       m_rev_edge_map(rev),
 104       m_pre_map(pre),
 105       m_tree_map(color),
 106       m_dist_map(dist),
 107       m_source(src),
 108       m_sink(sink),
 109       m_active_nodes(),
 110       m_in_active_list_vec(num_vertices(g), false),
 111       m_in_active_list_map(make_iterator_property_map(m_in_active_list_vec.begin(), m_index_map)),
 112       m_has_parent_vec(num_vertices(g), false),
 113       m_has_parent_map(make_iterator_property_map(m_has_parent_vec.begin(), m_index_map)),
 114       m_time_vec(num_vertices(g), 0),
 115       m_time_map(make_iterator_property_map(m_time_vec.begin(), m_index_map)),
 116       m_flow(0),
 117       m_time(1),
 118       m_last_grow_vertex(graph_traits<Graph>::null_vertex()){
 119         // initialize the color-map with gray-values
 120         vertex_iterator vi, v_end;
 121         for(boost::tie(vi, v_end) = vertices(m_g); vi != v_end; ++vi){
 122           set_tree(*vi, tColorTraits::gray());
 123         }
 124         // Initialize flow to zero which means initializing
 125         // the residual capacity equal to the capacity
 126         edge_iterator ei, e_end;
 127         for(boost::tie(ei, e_end) = edges(m_g); ei != e_end; ++ei) {
 128           put(m_res_cap_map, *ei, get(m_cap_map, *ei));
 129           BOOST_ASSERT(get(m_rev_edge_map, get(m_rev_edge_map, *ei)) == *ei); //check if the reverse edge map is build up properly
 130         }
 131         //init the search trees with the two terminals
 132         set_tree(m_source, tColorTraits::black());
 133         set_tree(m_sink, tColorTraits::white());
 134         put(m_time_map, m_source, 1);
 135         put(m_time_map, m_sink, 1);
 136       }
 137
 138       tEdgeVal max_flow(){
 139         //augment direct paths from SOURCE->SINK and SOURCE->VERTEX->SINK
 140         augment_direct_paths();
 141         //start the main-loop
 142         while(true){
 143           bool path_found;
 144           edge_descriptor connecting_edge;
 145           boost::tie(connecting_edge, path_found) = grow(); //find a path from source to sink
 146           if(!path_found){
 147             //we're finished, no more paths were found
 148             break;
 149           }
 150           ++m_time;
 151           augment(connecting_edge); //augment that path
 152           adopt(); //rebuild search tree structure
 153         }
 154         return m_flow;
 155       }
 156
 157       // the complete class is protected, as we want access to members in
 158       // derived test-class (see test/boykov_kolmogorov_max_flow_test.cpp)
 159     protected:
 160       void augment_direct_paths(){
 161         // in a first step, we augment all direct paths from source->NODE->sink
 162         // and additionally paths from source->sink. This improves especially
 163         // graphcuts for segmentation, as most of the nodes have source/sink
 164         // connects but shouldn't have an impact on other maxflow problems
 165         // (this is done in grow() anyway)
 166         out_edge_iterator ei, e_end;
 167         for(boost::tie(ei, e_end) = out_edges(m_source, m_g); ei != e_end; ++ei){
 168           edge_descriptor from_source = *ei;
 169           vertex_descriptor current_node = target(from_source, m_g);
 170           if(current_node == m_sink){
 171             tEdgeVal cap = get(m_res_cap_map, from_source);
 172             put(m_res_cap_map, from_source, 0);
 173             m_flow += cap;
 174             continue;
 175           }
 176           edge_descriptor to_sink;
 177           bool is_there;
 178           boost::tie(to_sink, is_there) = lookup_edge(current_node, m_sink, m_g);
 179           if(is_there){
 180             tEdgeVal cap_from_source = get(m_res_cap_map, from_source);
 181             tEdgeVal cap_to_sink = get(m_res_cap_map, to_sink);
 182             if(cap_from_source > cap_to_sink){
 183               set_tree(current_node, tColorTraits::black());
 184               add_active_node(current_node);
 185               set_edge_to_parent(current_node, from_source);
 186               put(m_dist_map, current_node, 1);
 187               put(m_time_map, current_node, 1);
 188               // add stuff to flow and update residuals. we dont need to
 189               // update reverse_edges, as incoming/outgoing edges to/from
 190               // source/sink don't count for max-flow
 191               put(m_res_cap_map, from_source, get(m_res_cap_map, from_source) - cap_to_sink);
 192               put(m_res_cap_map, to_sink, 0);
 193               m_flow += cap_to_sink;
 194             } else if(cap_to_sink > 0){
 195               set_tree(current_node, tColorTraits::white());
 196               add_active_node(current_node);
 197               set_edge_to_parent(current_node, to_sink);
 198               put(m_dist_map, current_node, 1);
 199               put(m_time_map, current_node, 1);
 200               // add stuff to flow and update residuals. we dont need to update
 201               // reverse_edges, as incoming/outgoing edges to/from source/sink
 202               // don't count for max-flow
 203               put(m_res_cap_map, to_sink, get(m_res_cap_map, to_sink) - cap_from_source);
 204               put(m_res_cap_map, from_source, 0);
 205               m_flow += cap_from_source;
 206             }
 207           } else if(get(m_res_cap_map, from_source)){
 208             // there is no sink connect, so we can't augment this path, but to
 209             // avoid adding m_source to the active nodes, we just activate this
 210             // node and set the approciate things
 211             set_tree(current_node, tColorTraits::black());
 212             set_edge_to_parent(current_node, from_source);
 213             put(m_dist_map, current_node, 1);
 214             put(m_time_map, current_node, 1);
 215             add_active_node(current_node);
 216           }
 217         }
 218         for(boost::tie(ei, e_end) = out_edges(m_sink, m_g); ei != e_end; ++ei){
 219           edge_descriptor to_sink = get(m_rev_edge_map, *ei);
 220           vertex_descriptor current_node = source(to_sink, m_g);
 221           if(get(m_res_cap_map, to_sink)){
 222             set_tree(current_node, tColorTraits::white());
 223             set_edge_to_parent(current_node, to_sink);
 224             put(m_dist_map, current_node, 1);
 225             put(m_time_map, current_node, 1);
 226             add_active_node(current_node);
 227           }
 228         }
 229       }
 230
 231       /**
 232        * Returns a pair of an edge and a boolean. if the bool is true, the
 233        * edge is a connection of a found path from s->t , read "the link" and
 234        * source(returnVal, m_g) is the end of the path found in the source-tree
 235        * target(returnVal, m_g) is the beginning of the path found in the sink-tree
 236        */
 237       std::pair<edge_descriptor, bool> grow(){
 238         BOOST_ASSERT(m_orphans.empty());
 239         vertex_descriptor current_node;
 240         while((current_node = get_next_active_node()) != graph_traits<Graph>::null_vertex()){ //if there is one
 241           BOOST_ASSERT(get_tree(current_node) != tColorTraits::gray() &&
 242                        (has_parent(current_node) ||
 243                          current_node == m_source ||
 244                          current_node == m_sink));
 245
 246           if(get_tree(current_node) == tColorTraits::black()){
 247             //source tree growing
 248             out_edge_iterator ei, e_end;
 249             if(current_node != m_last_grow_vertex){
 250               m_last_grow_vertex = current_node;
 251               boost::tie(m_last_grow_edge_it, m_last_grow_edge_end) = out_edges(current_node, m_g);
 252             }
 253             for(; m_last_grow_edge_it != m_last_grow_edge_end; ++m_last_grow_edge_it) {
 254               edge_descriptor out_edge = *m_last_grow_edge_it;
 255               if(get(m_res_cap_map, out_edge) > 0){ //check if we have capacity left on this edge
 256                 vertex_descriptor other_node = target(out_edge, m_g);
 257                 if(get_tree(other_node) == tColorTraits::gray()){ //it's a free node
 258                   set_tree(other_node, tColorTraits::black()); //aquire other node to our search tree
 259                   set_edge_to_parent(other_node, out_edge);   //set us as parent
 260                   put(m_dist_map, other_node, get(m_dist_map, current_node) + 1);  //and update the distance-heuristic
 261                   put(m_time_map, other_node, get(m_time_map, current_node));
 262                   add_active_node(other_node);
 263                 } else if(get_tree(other_node) == tColorTraits::black()) {
 264                   // we do this to get shorter paths. check if we are nearer to
 265                   // the source as its parent is
 266                   if(is_closer_to_terminal(current_node, other_node)){
 267                     set_edge_to_parent(other_node, out_edge);
 268                     put(m_dist_map, other_node, get(m_dist_map, current_node) + 1);
 269                     put(m_time_map, other_node, get(m_time_map, current_node));
 270                   }
 271                 } else{
 272                   BOOST_ASSERT(get_tree(other_node)==tColorTraits::white());
 273                   //kewl, found a path from one to the other search tree, return
 274                   // the connecting edge in src->sink dir
 275                   return std::make_pair(out_edge, true);
 276                 }
 277               }
 278             } //for all out-edges
 279           } //source-tree-growing
 280           else{
 281             BOOST_ASSERT(get_tree(current_node) == tColorTraits::white());
 282             out_edge_iterator ei, e_end;
 283             if(current_node != m_last_grow_vertex){
 284               m_last_grow_vertex = current_node;
 285               boost::tie(m_last_grow_edge_it, m_last_grow_edge_end) = out_edges(current_node, m_g);
 286             }
 287             for(; m_last_grow_edge_it != m_last_grow_edge_end; ++m_last_grow_edge_it){
 288               edge_descriptor in_edge = get(m_rev_edge_map, *m_last_grow_edge_it);
 289               if(get(m_res_cap_map, in_edge) > 0){ //check if there is capacity left
 290                 vertex_descriptor other_node = source(in_edge, m_g);
 291                 if(get_tree(other_node) == tColorTraits::gray()){ //it's a free node
 292                   set_tree(other_node, tColorTraits::white());      //aquire that node to our search tree
 293                   set_edge_to_parent(other_node, in_edge);          //set us as parent
 294                   add_active_node(other_node);                      //activate that node
 295                   put(m_dist_map, other_node, get(m_dist_map, current_node) + 1); //set its distance
 296                   put(m_time_map, other_node, get(m_time_map, current_node));//and time
 297                 } else if(get_tree(other_node) == tColorTraits::white()){
 298                   if(is_closer_to_terminal(current_node, other_node)){
 299                     //we are closer to the sink than its parent is, so we "adopt" him
 300                     set_edge_to_parent(other_node, in_edge);
 301                     put(m_dist_map, other_node, get(m_dist_map, current_node) + 1);
 302                     put(m_time_map, other_node, get(m_time_map, current_node));
 303                   }
 304                 } else{
 305                   BOOST_ASSERT(get_tree(other_node)==tColorTraits::black());
 306                   //kewl, found a path from one to the other search tree,
 307                   // return the connecting edge in src->sink dir
 308                   return std::make_pair(in_edge, true);
 309                 }
 310               }
 311             } //for all out-edges
 312           } //sink-tree growing
 313
 314           //all edges of that node are processed, and no more paths were found.
 315           // remove if from the front of the active queue
 316           finish_node(current_node);
 317         } //while active_nodes not empty
 318
 319         //no active nodes anymore and no path found, we're done
 320         return std::make_pair(edge_descriptor(), false);
 321       }
 322
 323       /**
 324        * augments path from s->t and updates residual graph
 325        * source(e, m_g) is the end of the path found in the source-tree
 326        * target(e, m_g) is the beginning of the path found in the sink-tree
 327        * this phase generates orphans on satured edges, if the attached verts are
 328        * from different search-trees orphans are ordered in distance to
 329        * sink/source. first the farest from the source are front_inserted into
 330        * the orphans list, and after that the sink-tree-orphans are
 331        * front_inserted. when going to adoption stage the orphans are popped_front,
 332        * and so we process the nearest verts to the terminals first
 333        */
 334       void augment(edge_descriptor e) {
 335         BOOST_ASSERT(get_tree(target(e, m_g)) == tColorTraits::white());
 336         BOOST_ASSERT(get_tree(source(e, m_g)) == tColorTraits::black());
 337         BOOST_ASSERT(m_orphans.empty());
 338
 339         const tEdgeVal bottleneck = find_bottleneck(e);
 340         //now we push the found flow through the path
 341         //for each edge we saturate we have to look for the verts that belong to that edge, one of them becomes an orphans
 342         //now process the connecting edge
 343         put(m_res_cap_map, e, get(m_res_cap_map, e) - bottleneck);
 344         BOOST_ASSERT(get(m_res_cap_map, e) >= 0);
 345         put(m_res_cap_map, get(m_rev_edge_map, e), get(m_res_cap_map, get(m_rev_edge_map, e)) + bottleneck);
 346
 347         //now we follow the path back to the source
 348         vertex_descriptor current_node = source(e, m_g);
 349         while(current_node != m_source){
 350           edge_descriptor pred = get_edge_to_parent(current_node);
 351           put(m_res_cap_map, pred, get(m_res_cap_map, pred) - bottleneck);
 352           BOOST_ASSERT(get(m_res_cap_map, pred) >= 0);
 353           put(m_res_cap_map, get(m_rev_edge_map, pred), get(m_res_cap_map, get(m_rev_edge_map, pred)) + bottleneck);
 354           if(get(m_res_cap_map, pred) == 0){
 355             set_no_parent(current_node);
 356             m_orphans.push_front(current_node);
 357           }
 358           current_node = source(pred, m_g);
 359         }
 360         //then go forward in the sink-tree
 361         current_node = target(e, m_g);
 362         while(current_node != m_sink){
 363           edge_descriptor pred = get_edge_to_parent(current_node);
 364           put(m_res_cap_map, pred, get(m_res_cap_map, pred) - bottleneck);
 365           BOOST_ASSERT(get(m_res_cap_map, pred) >= 0);
 366           put(m_res_cap_map, get(m_rev_edge_map, pred), get(m_res_cap_map, get(m_rev_edge_map, pred)) + bottleneck);
 367           if(get(m_res_cap_map, pred) == 0){
 368             set_no_parent(current_node);
 369             m_orphans.push_front(current_node);
 370           }
 371           current_node = target(pred, m_g);
 372         }
 373         //and add it to the max-flow
 374         m_flow += bottleneck;
 375       }
 376
 377       /**
 378        * returns the bottleneck of a s->t path (end_of_path is last vertex in
 379        * source-tree, begin_of_path is first vertex in sink-tree)
 380        */
 381       inline tEdgeVal find_bottleneck(edge_descriptor e){
 382         BOOST_USING_STD_MIN();
 383         tEdgeVal minimum_cap = get(m_res_cap_map, e);
 384         vertex_descriptor current_node = source(e, m_g);
 385         //first go back in the source tree
 386         while(current_node != m_source){
 387           edge_descriptor pred = get_edge_to_parent(current_node);
 388           minimum_cap = min BOOST_PREVENT_MACRO_SUBSTITUTION(minimum_cap, get(m_res_cap_map, pred));
 389           current_node = source(pred, m_g);
 390         }
 391         //then go forward in the sink-tree
 392         current_node = target(e, m_g);
 393         while(current_node != m_sink){
 394           edge_descriptor pred = get_edge_to_parent(current_node);
 395           minimum_cap = min BOOST_PREVENT_MACRO_SUBSTITUTION(minimum_cap, get(m_res_cap_map, pred));
 396           current_node = target(pred, m_g);
 397         }
 398         return minimum_cap;
 399       }
 400
 401       /**
 402        * rebuild search trees
 403        * empty the queue of orphans, and find new parents for them or just drop
 404        * them from the search trees
 405        */
 406       void adopt(){
 407         while(!m_orphans.empty() || !m_child_orphans.empty()){
 408           vertex_descriptor current_node;
 409           if(m_child_orphans.empty()){
 410             //get the next orphan from the main-queue  and remove it
 411             current_node = m_orphans.front();
 412             m_orphans.pop_front();
 413           } else{
 414             current_node = m_child_orphans.front();
 415             m_child_orphans.pop();
 416           }
 417           if(get_tree(current_node) == tColorTraits::black()){
 418             //we're in the source-tree
 419             tDistanceVal min_distance = (std::numeric_limits<tDistanceVal>::max)();
 420             edge_descriptor new_parent_edge;
 421             out_edge_iterator ei, e_end;
 422             for(boost::tie(ei, e_end) = out_edges(current_node, m_g); ei != e_end; ++ei){
 423               const edge_descriptor in_edge = get(m_rev_edge_map, *ei);
 424               BOOST_ASSERT(target(in_edge, m_g) == current_node); //we should be the target of this edge
 425               if(get(m_res_cap_map, in_edge) > 0){
 426                 vertex_descriptor other_node = source(in_edge, m_g);
 427                 if(get_tree(other_node) == tColorTraits::black() && has_source_connect(other_node)){
 428                   if(get(m_dist_map, other_node) < min_distance){
 429                     min_distance = get(m_dist_map, other_node);
 430                     new_parent_edge = in_edge;
 431                   }
 432                 }
 433               }
 434             }
 435             if(min_distance != (std::numeric_limits<tDistanceVal>::max)()){
 436               set_edge_to_parent(current_node, new_parent_edge);
 437               put(m_dist_map, current_node, min_distance + 1);
 438               put(m_time_map, current_node, m_time);
 439             } else{
 440               put(m_time_map, current_node, 0);
 441               for(boost::tie(ei, e_end) = out_edges(current_node, m_g); ei != e_end; ++ei){
 442                 edge_descriptor in_edge = get(m_rev_edge_map, *ei);
 443                 vertex_descriptor other_node = source(in_edge, m_g);
 444                 if(get_tree(other_node) == tColorTraits::black() && has_parent(other_node)){
 445                   if(get(m_res_cap_map, in_edge) > 0){
 446                     add_active_node(other_node);
 447                   }
 448                   if(source(get_edge_to_parent(other_node), m_g) == current_node){
 449                     //we are the parent of that node
 450                     //it has to find a new parent, too
 451                     set_no_parent(other_node);
 452                     m_child_orphans.push(other_node);
 453                   }
 454                 }
 455               }
 456               set_tree(current_node, tColorTraits::gray());
 457             } //no parent found
 458           } //source-tree-adoption
 459           else{
 460             //now we should be in the sink-tree, check that...
 461             BOOST_ASSERT(get_tree(current_node) == tColorTraits::white());
 462             out_edge_iterator ei, e_end;
 463             edge_descriptor new_parent_edge;
 464             tDistanceVal min_distance = (std::numeric_limits<tDistanceVal>::max)();
 465             for(boost::tie(ei, e_end) = out_edges(current_node, m_g); ei != e_end; ++ei){
 466               const edge_descriptor out_edge = *ei;
 467               if(get(m_res_cap_map, out_edge) > 0){
 468                 const vertex_descriptor other_node = target(out_edge, m_g);
 469                 if(get_tree(other_node) == tColorTraits::white() && has_sink_connect(other_node))
 470                   if(get(m_dist_map, other_node) < min_distance){
 471                     min_distance = get(m_dist_map, other_node);
 472                     new_parent_edge = out_edge;
 473                   }
 474               }
 475             }
 476             if(min_distance != (std::numeric_limits<tDistanceVal>::max)()){
 477               set_edge_to_parent(current_node, new_parent_edge);
 478               put(m_dist_map, current_node, min_distance + 1);
 479               put(m_time_map, current_node, m_time);
 480             } else{
 481               put(m_time_map, current_node, 0);
 482               for(boost::tie(ei, e_end) = out_edges(current_node, m_g); ei != e_end; ++ei){
 483                 const edge_descriptor out_edge = *ei;
 484                 const vertex_descriptor other_node = target(out_edge, m_g);
 485                 if(get_tree(other_node) == tColorTraits::white() && has_parent(other_node)){
 486                   if(get(m_res_cap_map, out_edge) > 0){
 487                     add_active_node(other_node);
 488                   }
 489                   if(target(get_edge_to_parent(other_node), m_g) == current_node){
 490                     //we were it's parent, so it has to find a new one, too
 491                     set_no_parent(other_node);
 492                     m_child_orphans.push(other_node);
 493                   }
 494                 }
 495               }
 496               set_tree(current_node, tColorTraits::gray());
 497             } //no parent found
 498           } //sink-tree adoption
 499         } //while !orphans.empty()
 500       } //adopt
 501
 502       /**
 503        * return next active vertex if there is one, otherwise a null_vertex
 504        */
 505       inline vertex_descriptor get_next_active_node(){
 506         while(true){
 507           if(m_active_nodes.empty())
 508             return graph_traits<Graph>::null_vertex();
 509           vertex_descriptor v = m_active_nodes.front();
 510
 511       //if it has no parent, this node can't be active (if its not source or sink)
 512       if(!has_parent(v) && v != m_source && v != m_sink){
 513             m_active_nodes.pop();
 514             put(m_in_active_list_map, v, false);
 515           } else{
 516             BOOST_ASSERT(get_tree(v) == tColorTraits::black() || get_tree(v) == tColorTraits::white());
 517             return v;
 518           }
 519         }
 520       }
 521
 522       /**
 523        * adds v as an active vertex, but only if its not in the list already
 524        */
 525       inline void add_active_node(vertex_descriptor v){
 526         BOOST_ASSERT(get_tree(v) != tColorTraits::gray());
 527         if(get(m_in_active_list_map, v)){
 528           return;
 529         } else{
 530           put(m_in_active_list_map, v, true);
 531           m_active_nodes.push(v);
 532         }
 533       }
 534
 535       /**
 536        * 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)
 537        */
 538       inline void finish_node(vertex_descriptor v){
 539         BOOST_ASSERT(m_active_nodes.front() == v);
 540         m_active_nodes.pop();
 541         put(m_in_active_list_map, v, false);
 542         m_last_grow_vertex = graph_traits<Graph>::null_vertex();
 543       }
 544
 545       /**
 546        * removes a vertex from the queue of active nodes (actually this does nothing,
 547        * but checks if this node has no parent edge, as this is the criteria for
 548        * being no more active)
 549        */
 550       inline void remove_active_node(vertex_descriptor v){
 551         BOOST_ASSERT(!has_parent(v));
 552       }
 553
 554       /**
 555        * returns the search tree of v; tColorValue::black() for source tree,
 556        * white() for sink tree, gray() for no tree
 557        */
 558       inline tColorValue get_tree(vertex_descriptor v) const {
 559         return get(m_tree_map, v);
 560       }
 561
 562       /**
 563        * sets search tree of v; tColorValue::black() for source tree, white()
 564        * for sink tree, gray() for no tree
 565        */
 566       inline void set_tree(vertex_descriptor v, tColorValue t){
 567         put(m_tree_map, v, t);
 568       }
 569
 570       /**
 571        * returns edge to parent vertex of v;
 572        */
 573       inline edge_descriptor get_edge_to_parent(vertex_descriptor v) const{
 574         return get(m_pre_map, v);
 575       }
 576
 577       /**
 578        * returns true if the edge stored in m_pre_map[v] is a valid entry
 579        */
 580       inline bool has_parent(vertex_descriptor v) const{
 581         return get(m_has_parent_map, v);
 582       }
 583
 584       /**
 585        * sets edge to parent vertex of v;
 586        */
 587       inline void set_edge_to_parent(vertex_descriptor v, edge_descriptor f_edge_to_parent){
 588         BOOST_ASSERT(get(m_res_cap_map, f_edge_to_parent) > 0);
 589         put(m_pre_map, v, f_edge_to_parent);
 590         put(m_has_parent_map, v, true);
 591       }
 592
 593       /**
 594        * removes the edge to parent of v (this is done by invalidating the
 595        * entry an additional map)
 596        */
 597       inline void set_no_parent(vertex_descriptor v){
 598         put(m_has_parent_map, v, false);
 599       }
 600
 601       /**
 602        * checks if vertex v has a connect to the sink-vertex (@var m_sink)
 603        * @param v the vertex which is checked
 604        * @return true if a path to the sink was found, false if not
 605        */
 606       inline bool has_sink_connect(vertex_descriptor v){
 607         tDistanceVal current_distance = 0;
 608         vertex_descriptor current_vertex = v;
 609         while(true){
 610           if(get(m_time_map, current_vertex) == m_time){
 611             //we found a node which was already checked this round. use it for distance calculations
 612             current_distance += get(m_dist_map, current_vertex);
 613             break;
 614           }
 615           if(current_vertex == m_sink){
 616             put(m_time_map, m_sink, m_time);
 617             break;
 618           }
 619           if(has_parent(current_vertex)){
 620             //it has a parent, so get it
 621             current_vertex = target(get_edge_to_parent(current_vertex), m_g);
 622             ++current_distance;
 623           } else{
 624             //no path found
 625             return false;
 626           }
 627         }
 628         current_vertex=v;
 629         while(get(m_time_map, current_vertex) != m_time){
 630           put(m_dist_map, current_vertex, current_distance);
 631           --current_distance;
 632           put(m_time_map, current_vertex, m_time);
 633           current_vertex = target(get_edge_to_parent(current_vertex), m_g);
 634         }
 635         return true;
 636       }
 637
 638       /**
 639        * checks if vertex v has a connect to the source-vertex (@var m_source)
 640        * @param v the vertex which is checked
 641        * @return true if a path to the source was found, false if not
 642        */
 643       inline bool has_source_connect(vertex_descriptor v){
 644         tDistanceVal current_distance = 0;
 645         vertex_descriptor current_vertex = v;
 646         while(true){
 647           if(get(m_time_map, current_vertex) == m_time){
 648             //we found a node which was already checked this round. use it for distance calculations
 649             current_distance += get(m_dist_map, current_vertex);
 650             break;
 651           }
 652           if(current_vertex == m_source){
 653             put(m_time_map, m_source, m_time);
 654             break;
 655           }
 656           if(has_parent(current_vertex)){
 657             //it has a parent, so get it
 658             current_vertex = source(get_edge_to_parent(current_vertex), m_g);
 659             ++current_distance;
 660           } else{
 661             //no path found
 662             return false;
 663           }
 664         }
 665         current_vertex=v;
 666         while(get(m_time_map, current_vertex) != m_time){
 667             put(m_dist_map, current_vertex, current_distance);
 668             --current_distance;
 669             put(m_time_map, current_vertex, m_time);
 670             current_vertex = source(get_edge_to_parent(current_vertex), m_g);
 671         }
 672         return true;
 673       }
 674
 675       /**
 676        * returns true, if p is closer to a terminal than q
 677        */
 678       inline bool is_closer_to_terminal(vertex_descriptor p, vertex_descriptor q){
 679         //checks the timestamps first, to build no cycles, and after that the real distance
 680         return (get(m_time_map, q) <= get(m_time_map, p) &&
 681                 get(m_dist_map, q) > get(m_dist_map, p)+1);
 682       }
 683
 684       ////////
 685       // member vars
 686       ////////
 687       Graph& m_g;
 688       IndexMap m_index_map;
 689       EdgeCapacityMap m_cap_map;
 690       ResidualCapacityEdgeMap m_res_cap_map;
 691       ReverseEdgeMap m_rev_edge_map;
 692       PredecessorMap m_pre_map; //stores paths found in the growth stage
 693       ColorMap m_tree_map; //maps each vertex into one of the two search tree or none (gray())
 694       DistanceMap m_dist_map; //stores distance to source/sink nodes
 695       vertex_descriptor m_source;
 696       vertex_descriptor m_sink;
 697
 698       tQueue m_active_nodes;
 699       std::vector<bool> m_in_active_list_vec;
 700       iterator_property_map<std::vector<bool>::iterator, IndexMap> m_in_active_list_map;
 701
 702       std::list<vertex_descriptor> m_orphans;
 703       tQueue m_child_orphans; // we use a second queuqe for child orphans, as they are FIFO processed
 704
 705       std::vector<bool> m_has_parent_vec;
 706       iterator_property_map<std::vector<bool>::iterator, IndexMap> m_has_parent_map;
 707
 708       std::vector<long> m_time_vec; //timestamp of each node, used for sink/source-path calculations
 709       iterator_property_map<std::vector<long>::iterator, IndexMap> m_time_map;
 710       tEdgeVal m_flow;
 711       long m_time;
 712       vertex_descriptor m_last_grow_vertex;
 713       out_edge_iterator m_last_grow_edge_it;
 714       out_edge_iterator m_last_grow_edge_end;
 715 };
 716
 717 } //namespace boost::detail
 718
 719 /**
 720   * non-named-parameter version, given everything
 721   * this is the catch all version
 722   */
 723 template<class Graph,
 724          class CapacityEdgeMap,
 725          class ResidualCapacityEdgeMap,
 726          class ReverseEdgeMap, class PredecessorMap,
 727          class ColorMap,
 728          class DistanceMap,
 729          class IndexMap>
 730 typename property_traits<CapacityEdgeMap>::value_type
 731 boykov_kolmogorov_max_flow(Graph& g,
 732                            CapacityEdgeMap cap,
 733                            ResidualCapacityEdgeMap res_cap,
 734                            ReverseEdgeMap rev_map,
 735                            PredecessorMap pre_map,
 736                            ColorMap color,
 737                            DistanceMap dist,
 738                            IndexMap idx,
 739                            typename graph_traits<Graph>::vertex_descriptor src,
 740                            typename graph_traits<Graph>::vertex_descriptor sink)
 741 {
 742   typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
 743   typedef typename graph_traits<Graph>::edge_descriptor edge_descriptor;
 744
 745   //as this method is the last one before we instantiate the solver, we do the concept checks here
 746   function_requires<VertexListGraphConcept<Graph> >(); //to have vertices(), num_vertices(),
 747   function_requires<EdgeListGraphConcept<Graph> >(); //to have edges()
 748   function_requires<IncidenceGraphConcept<Graph> >(); //to have source(), target() and out_edges()
 749   function_requires<ReadablePropertyMapConcept<CapacityEdgeMap, edge_descriptor> >(); //read flow-values from edges
 750   function_requires<ReadWritePropertyMapConcept<ResidualCapacityEdgeMap, edge_descriptor> >(); //write flow-values to residuals
 751   function_requires<ReadablePropertyMapConcept<ReverseEdgeMap, edge_descriptor> >(); //read out reverse edges
 752   function_requires<ReadWritePropertyMapConcept<PredecessorMap, vertex_descriptor> >(); //store predecessor there
 753   function_requires<ReadWritePropertyMapConcept<ColorMap, vertex_descriptor> >(); //write corresponding tree
 754   function_requires<ReadWritePropertyMapConcept<DistanceMap, vertex_descriptor> >(); //write distance to source/sink
 755   function_requires<ReadablePropertyMapConcept<IndexMap, vertex_descriptor> >(); //get index 0...|V|-1
 756   BOOST_ASSERT(num_vertices(g) >= 2 && src != sink);
 757
 758   detail::bk_max_flow<
 759     Graph, CapacityEdgeMap, ResidualCapacityEdgeMap, ReverseEdgeMap,
 760     PredecessorMap, ColorMap, DistanceMap, IndexMap
 761   > algo(g, cap, res_cap, rev_map, pre_map, color, dist, idx, src, sink);
 762
 763   return algo.max_flow();
 764 }
 765
 766 /**
 767  * non-named-parameter version, given capacity, residucal_capacity,
 768  * reverse_edges, and an index map.
 769  */
 770 template<class Graph,
 771          class CapacityEdgeMap,
 772          class ResidualCapacityEdgeMap,
 773          class ReverseEdgeMap,
 774          class IndexMap>
 775 typename property_traits<CapacityEdgeMap>::value_type
 776 boykov_kolmogorov_max_flow(Graph& g,
 777                            CapacityEdgeMap cap,
 778                            ResidualCapacityEdgeMap res_cap,
 779                            ReverseEdgeMap rev,
 780                            IndexMap idx,
 781                            typename graph_traits<Graph>::vertex_descriptor src,
 782                            typename graph_traits<Graph>::vertex_descriptor sink)
 783 {
 784   typename graph_traits<Graph>::vertices_size_type n_verts = num_vertices(g);
 785   std::vector<typename graph_traits<Graph>::edge_descriptor> predecessor_vec(n_verts);
 786   std::vector<default_color_type> color_vec(n_verts);
 787   std::vector<typename graph_traits<Graph>::vertices_size_type> distance_vec(n_verts);
 788   return
 789     boykov_kolmogorov_max_flow(
 790       g, cap, res_cap, rev,
 791       make_iterator_property_map(predecessor_vec.begin(), idx),
 792       make_iterator_property_map(color_vec.begin(), idx),
 793       make_iterator_property_map(distance_vec.begin(), idx),
 794       idx, src, sink);
 795 }
 796
 797 /**
 798  * non-named-parameter version, some given: capacity, residual_capacity,
 799  * reverse_edges, color_map and an index map. Use this if you are interested in
 800  * the minimum cut, as the color map provides that info.
 801  */
 802 template<class Graph,
 803          class CapacityEdgeMap,
 804          class ResidualCapacityEdgeMap,
 805          class ReverseEdgeMap,
 806          class ColorMap,
 807          class IndexMap>
 808 typename property_traits<CapacityEdgeMap>::value_type
 809 boykov_kolmogorov_max_flow(Graph& g,
 810                            CapacityEdgeMap cap,
 811                            ResidualCapacityEdgeMap res_cap,
 812                            ReverseEdgeMap rev,
 813                            ColorMap color,
 814                            IndexMap idx,
 815                            typename graph_traits<Graph>::vertex_descriptor src,
 816                            typename graph_traits<Graph>::vertex_descriptor sink)
 817 {
 818   typename graph_traits<Graph>::vertices_size_type n_verts = num_vertices(g);
 819   std::vector<typename graph_traits<Graph>::edge_descriptor> predecessor_vec(n_verts);
 820   std::vector<typename graph_traits<Graph>::vertices_size_type> distance_vec(n_verts);
 821   return
 822     boykov_kolmogorov_max_flow(
 823       g, cap, res_cap, rev,
 824       make_iterator_property_map(predecessor_vec.begin(), idx),
 825       color,
 826       make_iterator_property_map(distance_vec.begin(), idx),
 827       idx, src, sink);
 828 }
 829
 830 /**
 831  * named-parameter version, some given
 832  */
 833 template<class Graph, class P, class T, class R>
 834 typename property_traits<typename property_map<Graph, edge_capacity_t>::const_type>::value_type
 835 boykov_kolmogorov_max_flow(Graph& g,
 836                            typename graph_traits<Graph>::vertex_descriptor src,
 837                            typename graph_traits<Graph>::vertex_descriptor sink,
 838                            const bgl_named_params<P, T, R>& params)
 839 {
 840   return
 841   boykov_kolmogorov_max_flow(
 842     g,
 843     choose_const_pmap(get_param(params, edge_capacity), g, edge_capacity),
 844     choose_pmap(get_param(params, edge_residual_capacity), g, edge_residual_capacity),
 845     choose_const_pmap(get_param(params, edge_reverse), g, edge_reverse),
 846     choose_pmap(get_param(params, vertex_predecessor), g, vertex_predecessor),
 847     choose_pmap(get_param(params, vertex_color), g, vertex_color),
 848     choose_pmap(get_param(params, vertex_distance), g, vertex_distance),
 849     choose_const_pmap(get_param(params, vertex_index), g, vertex_index),
 850     src, sink);
 851 }
 852
 853 /**
 854  * named-parameter version, none given
 855  */
 856 template<class Graph>
 857 typename property_traits<typename property_map<Graph, edge_capacity_t>::const_type>::value_type
 858 boykov_kolmogorov_max_flow(Graph& g,
 859                            typename graph_traits<Graph>::vertex_descriptor src,
 860                            typename graph_traits<Graph>::vertex_descriptor sink)
 861 {
 862   bgl_named_params<int, buffer_param_t> params(0); // bogus empty param
 863   return boykov_kolmogorov_max_flow(g, src, sink, params);
 864 }
 865
 866 } // namespace boost
 867
 868 #endif // BOOST_BOYKOV_KOLMOGOROV_MAX_FLOW_HPP
 869