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

source code of boost/boost/graph/boykov_kolmogorov_max_flow.hpp