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 | namespace detail { |
67 | |
68 | template <class Graph, |
69 | class EdgeCapacityMap, |
70 | class ResidualCapacityEdgeMap, |
71 | class ReverseEdgeMap, |
72 | class PredecessorMap, |
73 | class ColorMap, |
74 | class DistanceMap, |
75 | class IndexMap> |
76 | class 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 | */ |
727 | template<class Graph, |
728 | class CapacityEdgeMap, |
729 | class ResidualCapacityEdgeMap, |
730 | class ReverseEdgeMap, class PredecessorMap, |
731 | class ColorMap, |
732 | class DistanceMap, |
733 | class IndexMap> |
734 | typename property_traits<CapacityEdgeMap>::value_type |
735 | boykov_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 | */ |
774 | template<class Graph, |
775 | class CapacityEdgeMap, |
776 | class ResidualCapacityEdgeMap, |
777 | class ReverseEdgeMap, |
778 | class IndexMap> |
779 | typename property_traits<CapacityEdgeMap>::value_type |
780 | boykov_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 | */ |
806 | template<class Graph, |
807 | class CapacityEdgeMap, |
808 | class ResidualCapacityEdgeMap, |
809 | class ReverseEdgeMap, |
810 | class ColorMap, |
811 | class IndexMap> |
812 | typename property_traits<CapacityEdgeMap>::value_type |
813 | boykov_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 | */ |
837 | template<class Graph, class P, class T, class R> |
838 | typename property_traits<typename property_map<Graph, edge_capacity_t>::const_type>::value_type |
839 | boykov_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 | */ |
860 | template<class Graph> |
861 | typename property_traits<typename property_map<Graph, edge_capacity_t>::const_type>::value_type |
862 | boykov_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 | |