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 #include <boost/concept/assert.hpp>
51
52 // The algorithm impelemented here is described in:
53 //
54 // Boykov, Y., Kolmogorov, V. "An Experimental Comparison of Min-Cut/Max-Flow
55 // Algorithms for Energy Minimization in Vision", In IEEE Transactions on
56 // Pattern Analysis and Machine Intelligence, vol. 26, no. 9, pp. 1124-1137,
57 // Sep 2004.
58 //
59 // For further reading, also see:
60 //
61 // Kolmogorov, V. "Graph Based Algorithms for Scene Reconstruction from Two or
62 // More Views". PhD thesis, Cornell University, Sep 2003.
63
64 namespace boost
65 {
66
67 namespace detail
68 {
69
70 template < class Graph, class EdgeCapacityMap,
71 class ResidualCapacityEdgeMap, class ReverseEdgeMap,
72 class PredecessorMap, class ColorMap, class DistanceMap,
73 class IndexMap >
74 class bk_max_flow
75 {
76 typedef
77 typename property_traits< EdgeCapacityMap >::value_type tEdgeVal;
78 typedef graph_traits< Graph > tGraphTraits;
79 typedef typename tGraphTraits::vertex_iterator vertex_iterator;
80 typedef typename tGraphTraits::vertex_descriptor vertex_descriptor;
81 typedef typename tGraphTraits::edge_descriptor edge_descriptor;
82 typedef typename tGraphTraits::edge_iterator edge_iterator;
83 typedef typename tGraphTraits::out_edge_iterator out_edge_iterator;
84 typedef boost::queue< vertex_descriptor >
85 tQueue; // queue of vertices, used in adoption-stage
86 typedef typename property_traits< ColorMap >::value_type tColorValue;
87 typedef color_traits< tColorValue > tColorTraits;
88 typedef
89 typename property_traits< DistanceMap >::value_type tDistanceVal;
90
91 public:
bk_max_flow(Graph & g,EdgeCapacityMap cap,ResidualCapacityEdgeMap res,ReverseEdgeMap rev,PredecessorMap pre,ColorMap color,DistanceMap dist,IndexMap idx,vertex_descriptor src,vertex_descriptor sink)92 bk_max_flow(Graph& g, EdgeCapacityMap cap, ResidualCapacityEdgeMap res,
93 ReverseEdgeMap rev, PredecessorMap pre, ColorMap color,
94 DistanceMap dist, IndexMap idx, vertex_descriptor src,
95 vertex_descriptor sink)
96 : m_g(g)
97 , m_index_map(idx)
98 , m_cap_map(cap)
99 , m_res_cap_map(res)
100 , m_rev_edge_map(rev)
101 , m_pre_map(pre)
102 , m_tree_map(color)
103 , m_dist_map(dist)
104 , m_source(src)
105 , m_sink(sink)
106 , m_active_nodes()
107 , m_in_active_list_vec(num_vertices(g), false)
108 , m_in_active_list_map(make_iterator_property_map(
109 m_in_active_list_vec.begin(), m_index_map))
110 , m_has_parent_vec(num_vertices(g), false)
111 , m_has_parent_map(
112 make_iterator_property_map(m_has_parent_vec.begin(), m_index_map))
113 , m_time_vec(num_vertices(g), 0)
114 , m_time_map(
115 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 {
120 // initialize the color-map with gray-values
121 vertex_iterator vi, v_end;
122 for (boost::tie(vi, v_end) = vertices(m_g); vi != v_end; ++vi)
123 {
124 set_tree(*vi, tColorTraits::gray());
125 }
126 // Initialize flow to zero which means initializing
127 // the residual capacity equal to the capacity
128 edge_iterator ei, e_end;
129 for (boost::tie(ei, e_end) = edges(m_g); ei != e_end; ++ei)
130 {
131 put(m_res_cap_map, *ei, get(m_cap_map, *ei));
132 BOOST_ASSERT(get(m_rev_edge_map, get(m_rev_edge_map, *ei))
133 == *ei); // check if the reverse edge map is build up
134 // properly
135 }
136 // init the search trees with the two terminals
137 set_tree(m_source, tColorTraits::black());
138 set_tree(m_sink, tColorTraits::white());
139 put(m_time_map, m_source, 1);
140 put(m_time_map, m_sink, 1);
141 }
142
max_flow()143 tEdgeVal max_flow()
144 {
145 // augment direct paths from SOURCE->SINK and SOURCE->VERTEX->SINK
146 augment_direct_paths();
147 // start the main-loop
148 while (true)
149 {
150 bool path_found;
151 edge_descriptor connecting_edge;
152 boost::tie(connecting_edge, path_found)
153 = grow(); // find a path from source to sink
154 if (!path_found)
155 {
156 // we're finished, no more paths were found
157 break;
158 }
159 ++m_time;
160 augment(connecting_edge); // augment that path
161 adopt(); // rebuild search tree structure
162 }
163 return m_flow;
164 }
165
166 // the complete class is protected, as we want access to members in
167 // derived test-class (see test/boykov_kolmogorov_max_flow_test.cpp)
168 protected:
augment_direct_paths()169 void augment_direct_paths()
170 {
171 // in a first step, we augment all direct paths from
172 // source->NODE->sink and additionally paths from source->sink. This
173 // improves especially graphcuts for segmentation, as most of the
174 // nodes have source/sink connects but shouldn't have an impact on
175 // other maxflow problems (this is done in grow() anyway)
176 out_edge_iterator ei, e_end;
177 for (boost::tie(ei, e_end) = out_edges(m_source, m_g); ei != e_end;
178 ++ei)
179 {
180 edge_descriptor from_source = *ei;
181 vertex_descriptor current_node = target(from_source, m_g);
182 if (current_node == m_sink)
183 {
184 tEdgeVal cap = get(m_res_cap_map, from_source);
185 put(m_res_cap_map, from_source, 0);
186 m_flow += cap;
187 continue;
188 }
189 edge_descriptor to_sink;
190 bool is_there;
191 boost::tie(to_sink, is_there)
192 = lookup_edge(current_node, m_sink, m_g);
193 if (is_there)
194 {
195 tEdgeVal cap_from_source = get(m_res_cap_map, from_source);
196 tEdgeVal cap_to_sink = get(m_res_cap_map, to_sink);
197 if (cap_from_source > cap_to_sink)
198 {
199 set_tree(current_node, tColorTraits::black());
200 add_active_node(current_node);
201 set_edge_to_parent(current_node, from_source);
202 put(m_dist_map, current_node, 1);
203 put(m_time_map, current_node, 1);
204 // add stuff to flow and update residuals. we dont need
205 // to update reverse_edges, as incoming/outgoing edges
206 // to/from source/sink don't count for max-flow
207 put(m_res_cap_map, from_source,
208 get(m_res_cap_map, from_source) - cap_to_sink);
209 put(m_res_cap_map, to_sink, 0);
210 m_flow += cap_to_sink;
211 }
212 else if (cap_to_sink > 0)
213 {
214 set_tree(current_node, tColorTraits::white());
215 add_active_node(current_node);
216 set_edge_to_parent(current_node, to_sink);
217 put(m_dist_map, current_node, 1);
218 put(m_time_map, current_node, 1);
219 // add stuff to flow and update residuals. we dont need
220 // to update reverse_edges, as incoming/outgoing edges
221 // to/from source/sink don't count for max-flow
222 put(m_res_cap_map, to_sink,
223 get(m_res_cap_map, to_sink) - cap_from_source);
224 put(m_res_cap_map, from_source, 0);
225 m_flow += cap_from_source;
226 }
227 }
228 else if (get(m_res_cap_map, from_source))
229 {
230 // there is no sink connect, so we can't augment this path,
231 // but to avoid adding m_source to the active nodes, we just
232 // activate this node and set the approciate things
233 set_tree(current_node, tColorTraits::black());
234 set_edge_to_parent(current_node, from_source);
235 put(m_dist_map, current_node, 1);
236 put(m_time_map, current_node, 1);
237 add_active_node(current_node);
238 }
239 }
240 for (boost::tie(ei, e_end) = out_edges(m_sink, m_g); ei != e_end;
241 ++ei)
242 {
243 edge_descriptor to_sink = get(m_rev_edge_map, *ei);
244 vertex_descriptor current_node = source(to_sink, m_g);
245 if (get(m_res_cap_map, to_sink))
246 {
247 set_tree(current_node, tColorTraits::white());
248 set_edge_to_parent(current_node, to_sink);
249 put(m_dist_map, current_node, 1);
250 put(m_time_map, current_node, 1);
251 add_active_node(current_node);
252 }
253 }
254 }
255
256 /**
257 * Returns a pair of an edge and a boolean. if the bool is true, the
258 * edge is a connection of a found path from s->t , read "the link" and
259 * source(returnVal, m_g) is the end of the path found in the
260 * source-tree target(returnVal, m_g) is the beginning of the path found
261 * in the sink-tree
262 */
grow()263 std::pair< edge_descriptor, bool > grow()
264 {
265 BOOST_ASSERT(m_orphans.empty());
266 vertex_descriptor current_node;
267 while ((current_node = get_next_active_node())
268 != graph_traits< Graph >::null_vertex())
269 { // if there is one
270 BOOST_ASSERT(get_tree(current_node) != tColorTraits::gray()
271 && (has_parent(current_node) || current_node == m_source
272 || current_node == m_sink));
273
274 if (get_tree(current_node) == tColorTraits::black())
275 {
276 // source tree growing
277 out_edge_iterator ei, e_end;
278 if (current_node != m_last_grow_vertex)
279 {
280 m_last_grow_vertex = current_node;
281 boost::tie(m_last_grow_edge_it, m_last_grow_edge_end)
282 = out_edges(current_node, m_g);
283 }
284 for (; m_last_grow_edge_it != m_last_grow_edge_end;
285 ++m_last_grow_edge_it)
286 {
287 edge_descriptor out_edge = *m_last_grow_edge_it;
288 if (get(m_res_cap_map, out_edge) > 0)
289 { // check if we have capacity left on this edge
290 vertex_descriptor other_node
291 = target(out_edge, m_g);
292 if (get_tree(other_node) == tColorTraits::gray())
293 { // it's a free node
294 set_tree(other_node,
295 tColorTraits::black()); // aquire other node
296 // to our search
297 // tree
298 set_edge_to_parent(
299 other_node, out_edge); // set us as parent
300 put(m_dist_map, other_node,
301 get(m_dist_map, current_node)
302 + 1); // and update the
303 // distance-heuristic
304 put(m_time_map, other_node,
305 get(m_time_map, current_node));
306 add_active_node(other_node);
307 }
308 else if (get_tree(other_node)
309 == tColorTraits::black())
310 {
311 // we do this to get shorter paths. check if we
312 // are nearer to the source as its parent is
313 if (is_closer_to_terminal(
314 current_node, other_node))
315 {
316 set_edge_to_parent(other_node, out_edge);
317 put(m_dist_map, other_node,
318 get(m_dist_map, current_node) + 1);
319 put(m_time_map, other_node,
320 get(m_time_map, current_node));
321 }
322 }
323 else
324 {
325 BOOST_ASSERT(get_tree(other_node)
326 == tColorTraits::white());
327 // kewl, found a path from one to the other
328 // search tree, return
329 // the connecting edge in src->sink dir
330 return std::make_pair(out_edge, true);
331 }
332 }
333 } // for all out-edges
334 } // source-tree-growing
335 else
336 {
337 BOOST_ASSERT(
338 get_tree(current_node) == tColorTraits::white());
339 out_edge_iterator ei, e_end;
340 if (current_node != m_last_grow_vertex)
341 {
342 m_last_grow_vertex = current_node;
343 boost::tie(m_last_grow_edge_it, m_last_grow_edge_end)
344 = out_edges(current_node, m_g);
345 }
346 for (; m_last_grow_edge_it != m_last_grow_edge_end;
347 ++m_last_grow_edge_it)
348 {
349 edge_descriptor in_edge
350 = get(m_rev_edge_map, *m_last_grow_edge_it);
351 if (get(m_res_cap_map, in_edge) > 0)
352 { // check if there is capacity left
353 vertex_descriptor other_node = source(in_edge, m_g);
354 if (get_tree(other_node) == tColorTraits::gray())
355 { // it's a free node
356 set_tree(other_node,
357 tColorTraits::white()); // aquire that node
358 // to our search
359 // tree
360 set_edge_to_parent(
361 other_node, in_edge); // set us as parent
362 add_active_node(
363 other_node); // activate that node
364 put(m_dist_map, other_node,
365 get(m_dist_map, current_node)
366 + 1); // set its distance
367 put(m_time_map, other_node,
368 get(m_time_map, current_node)); // and time
369 }
370 else if (get_tree(other_node)
371 == tColorTraits::white())
372 {
373 if (is_closer_to_terminal(
374 current_node, other_node))
375 {
376 // we are closer to the sink than its parent
377 // is, so we "adopt" him
378 set_edge_to_parent(other_node, in_edge);
379 put(m_dist_map, other_node,
380 get(m_dist_map, current_node) + 1);
381 put(m_time_map, other_node,
382 get(m_time_map, current_node));
383 }
384 }
385 else
386 {
387 BOOST_ASSERT(get_tree(other_node)
388 == tColorTraits::black());
389 // kewl, found a path from one to the other
390 // search tree,
391 // return the connecting edge in src->sink dir
392 return std::make_pair(in_edge, true);
393 }
394 }
395 } // for all out-edges
396 } // sink-tree growing
397
398 // all edges of that node are processed, and no more paths were
399 // found.
400 // remove if from the front of the active queue
401 finish_node(current_node);
402 } // while active_nodes not empty
403
404 // no active nodes anymore and no path found, we're done
405 return std::make_pair(edge_descriptor(), false);
406 }
407
408 /**
409 * augments path from s->t and updates residual graph
410 * source(e, m_g) is the end of the path found in the source-tree
411 * target(e, m_g) is the beginning of the path found in the sink-tree
412 * this phase generates orphans on satured edges, if the attached verts
413 * are from different search-trees orphans are ordered in distance to
414 * sink/source. first the farest from the source are front_inserted into
415 * the orphans list, and after that the sink-tree-orphans are
416 * front_inserted. when going to adoption stage the orphans are
417 * popped_front, and so we process the nearest verts to the terminals
418 * first
419 */
augment(edge_descriptor e)420 void augment(edge_descriptor e)
421 {
422 BOOST_ASSERT(get_tree(target(e, m_g)) == tColorTraits::white());
423 BOOST_ASSERT(get_tree(source(e, m_g)) == tColorTraits::black());
424 BOOST_ASSERT(m_orphans.empty());
425
426 const tEdgeVal bottleneck = find_bottleneck(e);
427 // now we push the found flow through the path
428 // for each edge we saturate we have to look for the verts that
429 // belong to that edge, one of them becomes an orphans now process
430 // the connecting edge
431 put(m_res_cap_map, e, get(m_res_cap_map, e) - bottleneck);
432 BOOST_ASSERT(get(m_res_cap_map, e) >= 0);
433 put(m_res_cap_map, get(m_rev_edge_map, e),
434 get(m_res_cap_map, get(m_rev_edge_map, e)) + bottleneck);
435
436 // now we follow the path back to the source
437 vertex_descriptor current_node = source(e, m_g);
438 while (current_node != m_source)
439 {
440 edge_descriptor pred = get_edge_to_parent(current_node);
441 put(m_res_cap_map, pred, get(m_res_cap_map, pred) - bottleneck);
442 BOOST_ASSERT(get(m_res_cap_map, pred) >= 0);
443 put(m_res_cap_map, get(m_rev_edge_map, pred),
444 get(m_res_cap_map, get(m_rev_edge_map, pred)) + bottleneck);
445 if (get(m_res_cap_map, pred) == 0)
446 {
447 set_no_parent(current_node);
448 m_orphans.push_front(current_node);
449 }
450 current_node = source(pred, m_g);
451 }
452 // then go forward in the sink-tree
453 current_node = target(e, m_g);
454 while (current_node != m_sink)
455 {
456 edge_descriptor pred = get_edge_to_parent(current_node);
457 put(m_res_cap_map, pred, get(m_res_cap_map, pred) - bottleneck);
458 BOOST_ASSERT(get(m_res_cap_map, pred) >= 0);
459 put(m_res_cap_map, get(m_rev_edge_map, pred),
460 get(m_res_cap_map, get(m_rev_edge_map, pred)) + bottleneck);
461 if (get(m_res_cap_map, pred) == 0)
462 {
463 set_no_parent(current_node);
464 m_orphans.push_front(current_node);
465 }
466 current_node = target(pred, m_g);
467 }
468 // and add it to the max-flow
469 m_flow += bottleneck;
470 }
471
472 /**
473 * returns the bottleneck of a s->t path (end_of_path is last vertex in
474 * source-tree, begin_of_path is first vertex in sink-tree)
475 */
find_bottleneck(edge_descriptor e)476 inline tEdgeVal find_bottleneck(edge_descriptor e)
477 {
478 BOOST_USING_STD_MIN();
479 tEdgeVal minimum_cap = get(m_res_cap_map, e);
480 vertex_descriptor current_node = source(e, m_g);
481 // first go back in the source tree
482 while (current_node != m_source)
483 {
484 edge_descriptor pred = get_edge_to_parent(current_node);
485 minimum_cap = min BOOST_PREVENT_MACRO_SUBSTITUTION(
486 minimum_cap, get(m_res_cap_map, pred));
487 current_node = source(pred, m_g);
488 }
489 // then go forward in the sink-tree
490 current_node = target(e, m_g);
491 while (current_node != m_sink)
492 {
493 edge_descriptor pred = get_edge_to_parent(current_node);
494 minimum_cap = min BOOST_PREVENT_MACRO_SUBSTITUTION(
495 minimum_cap, get(m_res_cap_map, pred));
496 current_node = target(pred, m_g);
497 }
498 return minimum_cap;
499 }
500
501 /**
502 * rebuild search trees
503 * empty the queue of orphans, and find new parents for them or just
504 * drop them from the search trees
505 */
adopt()506 void adopt()
507 {
508 while (!m_orphans.empty() || !m_child_orphans.empty())
509 {
510 vertex_descriptor current_node;
511 if (m_child_orphans.empty())
512 {
513 // get the next orphan from the main-queue and remove it
514 current_node = m_orphans.front();
515 m_orphans.pop_front();
516 }
517 else
518 {
519 current_node = m_child_orphans.front();
520 m_child_orphans.pop();
521 }
522 if (get_tree(current_node) == tColorTraits::black())
523 {
524 // we're in the source-tree
525 tDistanceVal min_distance
526 = (std::numeric_limits< tDistanceVal >::max)();
527 edge_descriptor new_parent_edge;
528 out_edge_iterator ei, e_end;
529 for (boost::tie(ei, e_end) = out_edges(current_node, m_g);
530 ei != e_end; ++ei)
531 {
532 const edge_descriptor in_edge
533 = get(m_rev_edge_map, *ei);
534 BOOST_ASSERT(target(in_edge, m_g)
535 == current_node); // we should be the target of this
536 // edge
537 if (get(m_res_cap_map, in_edge) > 0)
538 {
539 vertex_descriptor other_node = source(in_edge, m_g);
540 if (get_tree(other_node) == tColorTraits::black()
541 && has_source_connect(other_node))
542 {
543 if (get(m_dist_map, other_node) < min_distance)
544 {
545 min_distance = get(m_dist_map, other_node);
546 new_parent_edge = in_edge;
547 }
548 }
549 }
550 }
551 if (min_distance
552 != (std::numeric_limits< tDistanceVal >::max)())
553 {
554 set_edge_to_parent(current_node, new_parent_edge);
555 put(m_dist_map, current_node, min_distance + 1);
556 put(m_time_map, current_node, m_time);
557 }
558 else
559 {
560 put(m_time_map, current_node, 0);
561 for (boost::tie(ei, e_end)
562 = out_edges(current_node, m_g);
563 ei != e_end; ++ei)
564 {
565 edge_descriptor in_edge = get(m_rev_edge_map, *ei);
566 vertex_descriptor other_node = source(in_edge, m_g);
567 if (get_tree(other_node) == tColorTraits::black()
568 && other_node != m_source)
569 {
570 if (get(m_res_cap_map, in_edge) > 0)
571 {
572 add_active_node(other_node);
573 }
574 if (has_parent(other_node)
575 && source(
576 get_edge_to_parent(other_node), m_g)
577 == current_node)
578 {
579 // we are the parent of that node
580 // it has to find a new parent, too
581 set_no_parent(other_node);
582 m_child_orphans.push(other_node);
583 }
584 }
585 }
586 set_tree(current_node, tColorTraits::gray());
587 } // no parent found
588 } // source-tree-adoption
589 else
590 {
591 // now we should be in the sink-tree, check that...
592 BOOST_ASSERT(
593 get_tree(current_node) == tColorTraits::white());
594 out_edge_iterator ei, e_end;
595 edge_descriptor new_parent_edge;
596 tDistanceVal min_distance
597 = (std::numeric_limits< tDistanceVal >::max)();
598 for (boost::tie(ei, e_end) = out_edges(current_node, m_g);
599 ei != e_end; ++ei)
600 {
601 const edge_descriptor out_edge = *ei;
602 if (get(m_res_cap_map, out_edge) > 0)
603 {
604 const vertex_descriptor other_node
605 = target(out_edge, m_g);
606 if (get_tree(other_node) == tColorTraits::white()
607 && has_sink_connect(other_node))
608 if (get(m_dist_map, other_node) < min_distance)
609 {
610 min_distance = get(m_dist_map, other_node);
611 new_parent_edge = out_edge;
612 }
613 }
614 }
615 if (min_distance
616 != (std::numeric_limits< tDistanceVal >::max)())
617 {
618 set_edge_to_parent(current_node, new_parent_edge);
619 put(m_dist_map, current_node, min_distance + 1);
620 put(m_time_map, current_node, m_time);
621 }
622 else
623 {
624 put(m_time_map, current_node, 0);
625 for (boost::tie(ei, e_end)
626 = out_edges(current_node, m_g);
627 ei != e_end; ++ei)
628 {
629 const edge_descriptor out_edge = *ei;
630 const vertex_descriptor other_node
631 = target(out_edge, m_g);
632 if (get_tree(other_node) == tColorTraits::white()
633 && other_node != m_sink)
634 {
635 if (get(m_res_cap_map, out_edge) > 0)
636 {
637 add_active_node(other_node);
638 }
639 if (has_parent(other_node)
640 && target(
641 get_edge_to_parent(other_node), m_g)
642 == current_node)
643 {
644 // we were it's parent, so it has to find a
645 // new one, too
646 set_no_parent(other_node);
647 m_child_orphans.push(other_node);
648 }
649 }
650 }
651 set_tree(current_node, tColorTraits::gray());
652 } // no parent found
653 } // sink-tree adoption
654 } // while !orphans.empty()
655 } // adopt
656
657 /**
658 * return next active vertex if there is one, otherwise a null_vertex
659 */
get_next_active_node()660 inline vertex_descriptor get_next_active_node()
661 {
662 while (true)
663 {
664 if (m_active_nodes.empty())
665 return graph_traits< Graph >::null_vertex();
666 vertex_descriptor v = m_active_nodes.front();
667
668 // if it has no parent, this node can't be active (if its not
669 // source or sink)
670 if (!has_parent(v) && v != m_source && v != m_sink)
671 {
672 m_active_nodes.pop();
673 put(m_in_active_list_map, v, false);
674 }
675 else
676 {
677 BOOST_ASSERT(get_tree(v) == tColorTraits::black()
678 || get_tree(v) == tColorTraits::white());
679 return v;
680 }
681 }
682 }
683
684 /**
685 * adds v as an active vertex, but only if its not in the list already
686 */
add_active_node(vertex_descriptor v)687 inline void add_active_node(vertex_descriptor v)
688 {
689 BOOST_ASSERT(get_tree(v) != tColorTraits::gray());
690 if (get(m_in_active_list_map, v))
691 {
692 if (m_last_grow_vertex == v)
693 {
694 m_last_grow_vertex = graph_traits< Graph >::null_vertex();
695 }
696 return;
697 }
698 else
699 {
700 put(m_in_active_list_map, v, true);
701 m_active_nodes.push(v);
702 }
703 }
704
705 /**
706 * finish_node removes a node from the front of the active queue (its
707 * called in grow phase, if no more paths can be found using this node)
708 */
finish_node(vertex_descriptor v)709 inline void finish_node(vertex_descriptor v)
710 {
711 BOOST_ASSERT(m_active_nodes.front() == v);
712 m_active_nodes.pop();
713 put(m_in_active_list_map, v, false);
714 m_last_grow_vertex = graph_traits< Graph >::null_vertex();
715 }
716
717 /**
718 * removes a vertex from the queue of active nodes (actually this does
719 * nothing, but checks if this node has no parent edge, as this is the
720 * criteria for being no more active)
721 */
remove_active_node(vertex_descriptor v)722 inline void remove_active_node(vertex_descriptor v)
723 {
724 BOOST_ASSERT(!has_parent(v));
725 }
726
727 /**
728 * returns the search tree of v; tColorValue::black() for source tree,
729 * white() for sink tree, gray() for no tree
730 */
get_tree(vertex_descriptor v) const731 inline tColorValue get_tree(vertex_descriptor v) const
732 {
733 return get(m_tree_map, v);
734 }
735
736 /**
737 * sets search tree of v; tColorValue::black() for source tree, white()
738 * for sink tree, gray() for no tree
739 */
set_tree(vertex_descriptor v,tColorValue t)740 inline void set_tree(vertex_descriptor v, tColorValue t)
741 {
742 put(m_tree_map, v, t);
743 }
744
745 /**
746 * returns edge to parent vertex of v;
747 */
get_edge_to_parent(vertex_descriptor v) const748 inline edge_descriptor get_edge_to_parent(vertex_descriptor v) const
749 {
750 return get(m_pre_map, v);
751 }
752
753 /**
754 * returns true if the edge stored in m_pre_map[v] is a valid entry
755 */
has_parent(vertex_descriptor v) const756 inline bool has_parent(vertex_descriptor v) const
757 {
758 return get(m_has_parent_map, v);
759 }
760
761 /**
762 * sets edge to parent vertex of v;
763 */
set_edge_to_parent(vertex_descriptor v,edge_descriptor f_edge_to_parent)764 inline void set_edge_to_parent(
765 vertex_descriptor v, edge_descriptor f_edge_to_parent)
766 {
767 BOOST_ASSERT(get(m_res_cap_map, f_edge_to_parent) > 0);
768 put(m_pre_map, v, f_edge_to_parent);
769 put(m_has_parent_map, v, true);
770 }
771
772 /**
773 * removes the edge to parent of v (this is done by invalidating the
774 * entry an additional map)
775 */
set_no_parent(vertex_descriptor v)776 inline void set_no_parent(vertex_descriptor v)
777 {
778 put(m_has_parent_map, v, false);
779 }
780
781 /**
782 * checks if vertex v has a connect to the sink-vertex (@var m_sink)
783 * @param v the vertex which is checked
784 * @return true if a path to the sink was found, false if not
785 */
has_sink_connect(vertex_descriptor v)786 inline bool has_sink_connect(vertex_descriptor v)
787 {
788 tDistanceVal current_distance = 0;
789 vertex_descriptor current_vertex = v;
790 while (true)
791 {
792 if (get(m_time_map, current_vertex) == m_time)
793 {
794 // we found a node which was already checked this round. use
795 // it for distance calculations
796 current_distance += get(m_dist_map, current_vertex);
797 break;
798 }
799 if (current_vertex == m_sink)
800 {
801 put(m_time_map, m_sink, m_time);
802 break;
803 }
804 if (has_parent(current_vertex))
805 {
806 // it has a parent, so get it
807 current_vertex
808 = target(get_edge_to_parent(current_vertex), m_g);
809 ++current_distance;
810 }
811 else
812 {
813 // no path found
814 return false;
815 }
816 }
817 current_vertex = v;
818 while (get(m_time_map, current_vertex) != m_time)
819 {
820 put(m_dist_map, current_vertex, current_distance);
821 --current_distance;
822 put(m_time_map, current_vertex, m_time);
823 current_vertex
824 = target(get_edge_to_parent(current_vertex), m_g);
825 }
826 return true;
827 }
828
829 /**
830 * checks if vertex v has a connect to the source-vertex (@var m_source)
831 * @param v the vertex which is checked
832 * @return true if a path to the source was found, false if not
833 */
has_source_connect(vertex_descriptor v)834 inline bool has_source_connect(vertex_descriptor v)
835 {
836 tDistanceVal current_distance = 0;
837 vertex_descriptor current_vertex = v;
838 while (true)
839 {
840 if (get(m_time_map, current_vertex) == m_time)
841 {
842 // we found a node which was already checked this round. use
843 // it for distance calculations
844 current_distance += get(m_dist_map, current_vertex);
845 break;
846 }
847 if (current_vertex == m_source)
848 {
849 put(m_time_map, m_source, m_time);
850 break;
851 }
852 if (has_parent(current_vertex))
853 {
854 // it has a parent, so get it
855 current_vertex
856 = source(get_edge_to_parent(current_vertex), m_g);
857 ++current_distance;
858 }
859 else
860 {
861 // no path found
862 return false;
863 }
864 }
865 current_vertex = v;
866 while (get(m_time_map, current_vertex) != m_time)
867 {
868 put(m_dist_map, current_vertex, current_distance);
869 --current_distance;
870 put(m_time_map, current_vertex, m_time);
871 current_vertex
872 = source(get_edge_to_parent(current_vertex), m_g);
873 }
874 return true;
875 }
876
877 /**
878 * returns true, if p is closer to a terminal than q
879 */
is_closer_to_terminal(vertex_descriptor p,vertex_descriptor q)880 inline bool is_closer_to_terminal(
881 vertex_descriptor p, vertex_descriptor q)
882 {
883 // checks the timestamps first, to build no cycles, and after that
884 // the real distance
885 return (get(m_time_map, q) <= get(m_time_map, p)
886 && get(m_dist_map, q) > get(m_dist_map, p) + 1);
887 }
888
889 ////////
890 // member vars
891 ////////
892 Graph& m_g;
893 IndexMap m_index_map;
894 EdgeCapacityMap m_cap_map;
895 ResidualCapacityEdgeMap m_res_cap_map;
896 ReverseEdgeMap m_rev_edge_map;
897 PredecessorMap m_pre_map; // stores paths found in the growth stage
898 ColorMap m_tree_map; // maps each vertex into one of the two search tree
899 // or none (gray())
900 DistanceMap m_dist_map; // stores distance to source/sink nodes
901 vertex_descriptor m_source;
902 vertex_descriptor m_sink;
903
904 tQueue m_active_nodes;
905 std::vector< bool > m_in_active_list_vec;
906 iterator_property_map< std::vector< bool >::iterator, IndexMap >
907 m_in_active_list_map;
908
909 std::list< vertex_descriptor > m_orphans;
910 tQueue m_child_orphans; // we use a second queuqe for child orphans, as
911 // they are FIFO processed
912
913 std::vector< bool > m_has_parent_vec;
914 iterator_property_map< std::vector< bool >::iterator, IndexMap >
915 m_has_parent_map;
916
917 std::vector< long > m_time_vec; // timestamp of each node, used for
918 // sink/source-path calculations
919 iterator_property_map< std::vector< long >::iterator, IndexMap >
920 m_time_map;
921 tEdgeVal m_flow;
922 long m_time;
923 vertex_descriptor m_last_grow_vertex;
924 out_edge_iterator m_last_grow_edge_it;
925 out_edge_iterator m_last_grow_edge_end;
926 };
927
928 } // namespace boost::detail
929
930 /**
931 * non-named-parameter version, given everything
932 * this is the catch all version
933 */
934 template < class Graph, class CapacityEdgeMap, class ResidualCapacityEdgeMap,
935 class ReverseEdgeMap, class PredecessorMap, class ColorMap,
936 class DistanceMap, class IndexMap >
937 typename property_traits< CapacityEdgeMap >::value_type
boykov_kolmogorov_max_flow(Graph & g,CapacityEdgeMap cap,ResidualCapacityEdgeMap res_cap,ReverseEdgeMap rev_map,PredecessorMap pre_map,ColorMap color,DistanceMap dist,IndexMap idx,typename graph_traits<Graph>::vertex_descriptor src,typename graph_traits<Graph>::vertex_descriptor sink)938 boykov_kolmogorov_max_flow(Graph& g, CapacityEdgeMap cap,
939 ResidualCapacityEdgeMap res_cap, ReverseEdgeMap rev_map,
940 PredecessorMap pre_map, ColorMap color, DistanceMap dist, IndexMap idx,
941 typename graph_traits< Graph >::vertex_descriptor src,
942 typename graph_traits< Graph >::vertex_descriptor sink)
943 {
944 typedef typename graph_traits< Graph >::vertex_descriptor vertex_descriptor;
945 typedef typename graph_traits< Graph >::edge_descriptor edge_descriptor;
946
947 // as this method is the last one before we instantiate the solver, we do
948 // the concept checks here
949 BOOST_CONCEPT_ASSERT(
950 (VertexListGraphConcept< Graph >)); // to have vertices(),
951 // num_vertices(),
952 BOOST_CONCEPT_ASSERT((EdgeListGraphConcept< Graph >)); // to have edges()
953 BOOST_CONCEPT_ASSERT(
954 (IncidenceGraphConcept< Graph >)); // to have source(), target() and
955 // out_edges()
956 BOOST_CONCEPT_ASSERT((ReadablePropertyMapConcept< CapacityEdgeMap,
957 edge_descriptor >)); // read flow-values from edges
958 BOOST_CONCEPT_ASSERT((ReadWritePropertyMapConcept< ResidualCapacityEdgeMap,
959 edge_descriptor >)); // write flow-values to residuals
960 BOOST_CONCEPT_ASSERT((ReadablePropertyMapConcept< ReverseEdgeMap,
961 edge_descriptor >)); // read out reverse edges
962 BOOST_CONCEPT_ASSERT((ReadWritePropertyMapConcept< PredecessorMap,
963 vertex_descriptor >)); // store predecessor there
964 BOOST_CONCEPT_ASSERT((ReadWritePropertyMapConcept< ColorMap,
965 vertex_descriptor >)); // write corresponding tree
966 BOOST_CONCEPT_ASSERT((ReadWritePropertyMapConcept< DistanceMap,
967 vertex_descriptor >)); // write distance to source/sink
968 BOOST_CONCEPT_ASSERT((ReadablePropertyMapConcept< IndexMap,
969 vertex_descriptor >)); // get index 0...|V|-1
970 BOOST_ASSERT(num_vertices(g) >= 2 && src != sink);
971
972 detail::bk_max_flow< Graph, CapacityEdgeMap, ResidualCapacityEdgeMap,
973 ReverseEdgeMap, PredecessorMap, ColorMap, DistanceMap, IndexMap >
974 algo(g, cap, res_cap, rev_map, pre_map, color, dist, idx, src, sink);
975
976 return algo.max_flow();
977 }
978
979 /**
980 * non-named-parameter version, given capacity, residucal_capacity,
981 * reverse_edges, and an index map.
982 */
983 template < class Graph, class CapacityEdgeMap, class ResidualCapacityEdgeMap,
984 class ReverseEdgeMap, class IndexMap >
985 typename property_traits< CapacityEdgeMap >::value_type
boykov_kolmogorov_max_flow(Graph & g,CapacityEdgeMap cap,ResidualCapacityEdgeMap res_cap,ReverseEdgeMap rev,IndexMap idx,typename graph_traits<Graph>::vertex_descriptor src,typename graph_traits<Graph>::vertex_descriptor sink)986 boykov_kolmogorov_max_flow(Graph& g, CapacityEdgeMap cap,
987 ResidualCapacityEdgeMap res_cap, ReverseEdgeMap rev, IndexMap idx,
988 typename graph_traits< Graph >::vertex_descriptor src,
989 typename graph_traits< Graph >::vertex_descriptor sink)
990 {
991 typename graph_traits< Graph >::vertices_size_type n_verts
992 = num_vertices(g);
993 std::vector< typename graph_traits< Graph >::edge_descriptor >
994 predecessor_vec(n_verts);
995 std::vector< default_color_type > color_vec(n_verts);
996 std::vector< typename graph_traits< Graph >::vertices_size_type >
997 distance_vec(n_verts);
998 return boykov_kolmogorov_max_flow(g, cap, res_cap, rev,
999 make_iterator_property_map(predecessor_vec.begin(), idx),
1000 make_iterator_property_map(color_vec.begin(), idx),
1001 make_iterator_property_map(distance_vec.begin(), idx), idx, src, sink);
1002 }
1003
1004 /**
1005 * non-named-parameter version, some given: capacity, residual_capacity,
1006 * reverse_edges, color_map and an index map. Use this if you are interested in
1007 * the minimum cut, as the color map provides that info.
1008 */
1009 template < class Graph, class CapacityEdgeMap, class ResidualCapacityEdgeMap,
1010 class ReverseEdgeMap, class ColorMap, class IndexMap >
1011 typename property_traits< CapacityEdgeMap >::value_type
boykov_kolmogorov_max_flow(Graph & g,CapacityEdgeMap cap,ResidualCapacityEdgeMap res_cap,ReverseEdgeMap rev,ColorMap color,IndexMap idx,typename graph_traits<Graph>::vertex_descriptor src,typename graph_traits<Graph>::vertex_descriptor sink)1012 boykov_kolmogorov_max_flow(Graph& g, CapacityEdgeMap cap,
1013 ResidualCapacityEdgeMap res_cap, ReverseEdgeMap rev, ColorMap color,
1014 IndexMap idx, typename graph_traits< Graph >::vertex_descriptor src,
1015 typename graph_traits< Graph >::vertex_descriptor sink)
1016 {
1017 typename graph_traits< Graph >::vertices_size_type n_verts
1018 = num_vertices(g);
1019 std::vector< typename graph_traits< Graph >::edge_descriptor >
1020 predecessor_vec(n_verts);
1021 std::vector< typename graph_traits< Graph >::vertices_size_type >
1022 distance_vec(n_verts);
1023 return boykov_kolmogorov_max_flow(g, cap, res_cap, rev,
1024 make_iterator_property_map(predecessor_vec.begin(), idx), color,
1025 make_iterator_property_map(distance_vec.begin(), idx), idx, src, sink);
1026 }
1027
1028 /**
1029 * named-parameter version, some given
1030 */
1031 template < class Graph, class P, class T, class R >
1032 typename property_traits<
1033 typename property_map< Graph, edge_capacity_t >::const_type >::value_type
boykov_kolmogorov_max_flow(Graph & g,typename graph_traits<Graph>::vertex_descriptor src,typename graph_traits<Graph>::vertex_descriptor sink,const bgl_named_params<P,T,R> & params)1034 boykov_kolmogorov_max_flow(Graph& g,
1035 typename graph_traits< Graph >::vertex_descriptor src,
1036 typename graph_traits< Graph >::vertex_descriptor sink,
1037 const bgl_named_params< P, T, R >& params)
1038 {
1039 return boykov_kolmogorov_max_flow(g,
1040 choose_const_pmap(get_param(params, edge_capacity), g, edge_capacity),
1041 choose_pmap(get_param(params, edge_residual_capacity), g,
1042 edge_residual_capacity),
1043 choose_const_pmap(get_param(params, edge_reverse), g, edge_reverse),
1044 choose_pmap(
1045 get_param(params, vertex_predecessor), g, vertex_predecessor),
1046 choose_pmap(get_param(params, vertex_color), g, vertex_color),
1047 choose_pmap(get_param(params, vertex_distance), g, vertex_distance),
1048 choose_const_pmap(get_param(params, vertex_index), g, vertex_index),
1049 src, sink);
1050 }
1051
1052 /**
1053 * named-parameter version, none given
1054 */
1055 template < class Graph >
1056 typename property_traits<
1057 typename property_map< Graph, edge_capacity_t >::const_type >::value_type
boykov_kolmogorov_max_flow(Graph & g,typename graph_traits<Graph>::vertex_descriptor src,typename graph_traits<Graph>::vertex_descriptor sink)1058 boykov_kolmogorov_max_flow(Graph& g,
1059 typename graph_traits< Graph >::vertex_descriptor src,
1060 typename graph_traits< Graph >::vertex_descriptor sink)
1061 {
1062 bgl_named_params< int, buffer_param_t > params(0); // bogus empty param
1063 return boykov_kolmogorov_max_flow(g, src, sink, params);
1064 }
1065
1066 } // namespace boost
1067
1068 #endif // BOOST_BOYKOV_KOLMOGOROV_MAX_FLOW_HPP
1069