1 // Copyright (c) 2006, Stephan Diederich
3 // This code may be used under either of the following two licences:
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
14 // The above copyright notice and this permission notice shall be
15 // included in all copies or substantial portions of the Software.
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.
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)
32 #ifndef BOOST_BOYKOV_KOLMOGOROV_MAX_FLOW_HPP
33 #define BOOST_BOYKOV_KOLMOGOROV_MAX_FLOW_HPP
35 #include <boost/config.hpp>
36 #include <boost/assert.hpp>
41 #include <algorithm> // for std::min and std::max
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>
51 // The algorithm impelemented here is described in:
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,
58 // For further reading, also see:
60 // Kolmogorov, V. "Graph Based Algorithms for Scene Reconstruction from Two or
61 // More Views". PhD thesis, Cornell University, Sep 2003.
67 template <class Graph,
68 class EdgeCapacityMap,
69 class ResidualCapacityEdgeMap,
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;
91 ResidualCapacityEdgeMap res,
97 vertex_descriptor src,
98 vertex_descriptor sink):
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)),
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());
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
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);
139 //augment direct paths from SOURCE->SINK and SOURCE->VERTEX->SINK
140 augment_direct_paths();
141 //start the main-loop
144 edge_descriptor connecting_edge;
145 boost::tie(connecting_edge, path_found) = grow(); //find a path from source to sink
147 //we're finished, no more paths were found
151 augment(connecting_edge); //augment that path
152 adopt(); //rebuild search tree structure
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)
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);
176 edge_descriptor to_sink;
178 boost::tie(to_sink, is_there) = lookup_edge(current_node, m_sink, m_g);
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;
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);
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);
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
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));
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);
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));
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);
278 } //for all out-edges
279 } //source-tree-growing
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);
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));
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);
311 } //for all out-edges
312 } //sink-tree growing
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
319 //no active nodes anymore and no path found, we're done
320 return std::make_pair(edge_descriptor(), false);
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
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());
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);
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);
358 current_node = source(pred, m_g);
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);
371 current_node = target(pred, m_g);
373 //and add it to the max-flow
374 m_flow += bottleneck;
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)
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);
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);
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
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();
414 current_node = m_child_orphans.front();
415 m_child_orphans.pop();
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;
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);
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);
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);
456 set_tree(current_node, tColorTraits::gray());
458 } //source-tree-adoption
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;
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);
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);
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);
496 set_tree(current_node, tColorTraits::gray());
498 } //sink-tree adoption
499 } //while !orphans.empty()
503 * return next active vertex if there is one, otherwise a null_vertex
505 inline vertex_descriptor get_next_active_node(){
507 if(m_active_nodes.empty())
508 return graph_traits<Graph>::null_vertex();
509 vertex_descriptor v = m_active_nodes.front();
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);
516 BOOST_ASSERT(get_tree(v) == tColorTraits::black() || get_tree(v) == tColorTraits::white());
523 * adds v as an active vertex, but only if its not in the list already
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)){
530 put(m_in_active_list_map, v, true);
531 m_active_nodes.push(v);
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)
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();
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)
550 inline void remove_active_node(vertex_descriptor v){
551 BOOST_ASSERT(!has_parent(v));
555 * returns the search tree of v; tColorValue::black() for source tree,
556 * white() for sink tree, gray() for no tree
558 inline tColorValue get_tree(vertex_descriptor v) const {
559 return get(m_tree_map, v);
563 * sets search tree of v; tColorValue::black() for source tree, white()
564 * for sink tree, gray() for no tree
566 inline void set_tree(vertex_descriptor v, tColorValue t){
567 put(m_tree_map, v, t);
571 * returns edge to parent vertex of v;
573 inline edge_descriptor get_edge_to_parent(vertex_descriptor v) const{
574 return get(m_pre_map, v);
578 * returns true if the edge stored in m_pre_map[v] is a valid entry
580 inline bool has_parent(vertex_descriptor v) const{
581 return get(m_has_parent_map, v);
585 * sets edge to parent vertex of v;
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);
594 * removes the edge to parent of v (this is done by invalidating the
595 * entry an additional map)
597 inline void set_no_parent(vertex_descriptor v){
598 put(m_has_parent_map, v, false);
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
606 inline bool has_sink_connect(vertex_descriptor v){
607 tDistanceVal current_distance = 0;
608 vertex_descriptor current_vertex = v;
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);
615 if(current_vertex == m_sink){
616 put(m_time_map, m_sink, m_time);
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);
629 while(get(m_time_map, current_vertex) != m_time){
630 put(m_dist_map, current_vertex, current_distance);
632 put(m_time_map, current_vertex, m_time);
633 current_vertex = target(get_edge_to_parent(current_vertex), m_g);
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
643 inline bool has_source_connect(vertex_descriptor v){
644 tDistanceVal current_distance = 0;
645 vertex_descriptor current_vertex = v;
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);
652 if(current_vertex == m_source){
653 put(m_time_map, m_source, m_time);
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);
666 while(get(m_time_map, current_vertex) != m_time){
667 put(m_dist_map, current_vertex, current_distance);
669 put(m_time_map, current_vertex, m_time);
670 current_vertex = source(get_edge_to_parent(current_vertex), m_g);
676 * returns true, if p is closer to a terminal than q
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);
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;
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;
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
705 std::vector<bool> m_has_parent_vec;
706 iterator_property_map<std::vector<bool>::iterator, IndexMap> m_has_parent_map;
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;
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;
717 } //namespace boost::detail
720 * non-named-parameter version, given everything
721 * this is the catch all version
723 template<class Graph,
724 class CapacityEdgeMap,
725 class ResidualCapacityEdgeMap,
726 class ReverseEdgeMap, class PredecessorMap,
730 typename property_traits<CapacityEdgeMap>::value_type
731 boykov_kolmogorov_max_flow(Graph& g,
733 ResidualCapacityEdgeMap res_cap,
734 ReverseEdgeMap rev_map,
735 PredecessorMap pre_map,
739 typename graph_traits<Graph>::vertex_descriptor src,
740 typename graph_traits<Graph>::vertex_descriptor sink)
742 typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
743 typedef typename graph_traits<Graph>::edge_descriptor edge_descriptor;
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);
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);
763 return algo.max_flow();
767 * non-named-parameter version, given capacity, residucal_capacity,
768 * reverse_edges, and an index map.
770 template<class Graph,
771 class CapacityEdgeMap,
772 class ResidualCapacityEdgeMap,
773 class ReverseEdgeMap,
775 typename property_traits<CapacityEdgeMap>::value_type
776 boykov_kolmogorov_max_flow(Graph& g,
778 ResidualCapacityEdgeMap res_cap,
781 typename graph_traits<Graph>::vertex_descriptor src,
782 typename graph_traits<Graph>::vertex_descriptor sink)
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);
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),
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.
802 template<class Graph,
803 class CapacityEdgeMap,
804 class ResidualCapacityEdgeMap,
805 class ReverseEdgeMap,
808 typename property_traits<CapacityEdgeMap>::value_type
809 boykov_kolmogorov_max_flow(Graph& g,
811 ResidualCapacityEdgeMap res_cap,
815 typename graph_traits<Graph>::vertex_descriptor src,
816 typename graph_traits<Graph>::vertex_descriptor sink)
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);
822 boykov_kolmogorov_max_flow(
823 g, cap, res_cap, rev,
824 make_iterator_property_map(predecessor_vec.begin(), idx),
826 make_iterator_property_map(distance_vec.begin(), idx),
831 * named-parameter version, some given
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)
841 boykov_kolmogorov_max_flow(
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),
854 * named-parameter version, none given
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)
862 bgl_named_params<int, buffer_param_t> params(0); // bogus empty param
863 return boykov_kolmogorov_max_flow(g, src, sink, params);
868 #endif // BOOST_BOYKOV_KOLMOGOROV_MAX_FLOW_HPP