]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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(*vi, 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(m_source, tColorTraits::black()); | |
134 | set_tree(m_sink, 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(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(current_node, tColorTraits::black()); | |
185 | add_active_node(current_node); | |
186 | set_edge_to_parent(current_node, 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(current_node, tColorTraits::white()); | |
197 | add_active_node(current_node); | |
198 | set_edge_to_parent(current_node, 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(current_node, tColorTraits::black()); | |
213 | set_edge_to_parent(current_node, from_source); | |
214 | put(m_dist_map, current_node, 1); | |
215 | put(m_time_map, current_node, 1); | |
216 | add_active_node(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(current_node, tColorTraits::white()); | |
224 | set_edge_to_parent(current_node, to_sink); | |
225 | put(m_dist_map, current_node, 1); | |
226 | put(m_time_map, current_node, 1); | |
227 | add_active_node(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(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(other_node) == tColorTraits::gray()){ //it's a free node | |
259 | set_tree(other_node, tColorTraits::black()); //aquire other node to our search tree | |
260 | set_edge_to_parent(other_node, 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(other_node); | |
264 | } else if(get_tree(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(current_node, other_node)){ | |
268 | set_edge_to_parent(other_node, 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(other_node) == tColorTraits::gray()){ //it's a free node | |
293 | set_tree(other_node, tColorTraits::white()); //aquire that node to our search tree | |
294 | set_edge_to_parent(other_node, in_edge); //set us as parent | |
295 | add_active_node(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(other_node) == tColorTraits::white()){ | |
299 | if(is_closer_to_terminal(current_node, other_node)){ | |
300 | //we are closer to the sink than its parent is, so we "adopt" him | |
301 | set_edge_to_parent(other_node, 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(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(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(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(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(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(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(other_node) == tColorTraits::black() && has_source_connect(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(current_node, 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(other_node) == tColorTraits::black() && other_node != m_source){ | |
446 | if(get(m_res_cap_map, in_edge) > 0){ | |
447 | add_active_node(other_node); | |
448 | } | |
449 | if(has_parent(other_node) && source(get_edge_to_parent(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(current_node, 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(other_node) == tColorTraits::white() && has_sink_connect(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(current_node, 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(other_node) == tColorTraits::white() && other_node != m_sink){ | |
487 | if(get(m_res_cap_map, out_edge) > 0){ | |
488 | add_active_node(other_node); | |
489 | } | |
490 | if(has_parent(other_node) && target(get_edge_to_parent(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(current_node, 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(current_vertex)){ | |
624 | //it has a parent, so get it | |
625 | current_vertex = target(get_edge_to_parent(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(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(current_vertex)){ | |
661 | //it has a parent, so get it | |
662 | current_vertex = source(get_edge_to_parent(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(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 |