ANIMA  4.0
animaGraph.h
Go to the documentation of this file.
1 /* graph.h */
2 /*
3 ###################################################################
4 # #
5 # MAXFLOW - software for computing mincut/maxflow in a graph #
6 # Version 3.0 #
7 # http://www.cs.adastral.ucl.ac.uk/~vnk/software.html #
8 # #
9 # Yuri Boykov (yuri@csd.uwo.ca) #
10 # Vladimir Kolmogorov (vnk@adastral.ucl.ac.uk) #
11 # 2001-2006 #
12 # #
13 ###################################################################
14 
15 1. Introduction.
16 
17 This software library implements the maxflow algorithm described in
18 
19  "An Experimental Comparison of Min-Cut/Max-Flow Algorithms for Energy Minimization in Vision."
20  Yuri Boykov and Vladimir Kolmogorov.
21  In IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI),
22  September 2004
23 
24 This algorithm was developed by Yuri Boykov and Vladimir Kolmogorov
25 at Siemens Corporate Research. To make it available for public use,
26 it was later reimplemented by Vladimir Kolmogorov based on open publications.
27 
28 If you use this software for research purposes, you should cite
29 the aforementioned paper in any resulting publication.
30 
31 ----------------------------------------------------------------------
32 
33 REUSING TREES:
34 
35 Starting with version 3.0, there is a also an option of reusing search
36 trees from one maxflow computation to the next, as described in
37 
38  "Efficiently Solving Dynamic Markov Random Fields Using Graph Cuts."
39  Pushmeet Kohli and Philip H.S. Torr
40  International Conference on Computer Vision (ICCV), 2005
41 
42 If you use this option, you should cite
43 the aforementioned paper in any resulting publication.
44 
45 Tested under windows, Visual C++ 6.0 compiler and unix (SunOS 5.8
46 and RedHat Linux 7.0, GNU c++ compiler).
47 
48 ##################################################################
49 
50 2. License & disclaimer.
51 
52  Copyright 2001 Vladimir Kolmogorov (vnk@adastral.ucl.ac.uk), Yuri Boykov (yuri@csd.uwo.ca).
53 
54  This software can be used for research purposes only.
55 
56  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
57  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
58  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
59  A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
60  OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
61  SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
62  LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
63  DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
64  THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
65  (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
66  OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
67 
68 ##################################################################
69 */
70 
71 
72 #pragma once
73 
74 #include <string.h>
75 #include "animaBlock.h"
76 
77 #include <assert.h>
78 // NOTE: in UNIX you need to use -DNDEBUG preprocessor option to supress assert's!!!
79 
80 
81 
82 // captype: type of edge capacities (excluding t-links)
83 // tcaptype: type of t-links (edges between nodes and terminals)
84 // flowtype: type of total flow
85 //
86 // Current instantiations are in instances.inc
87 
88 
89 namespace anima {
90 
91 template <typename captype, typename tcaptype, typename flowtype> class Graph
92 {
93 public:
94  typedef enum
95  {
96  SOURCE = 0,
97  SINK = 1
98  } termtype; // terminals
99  typedef int node_id;
100 
102  // BASIC INTERFACE FUNCTIONS //
103  // (should be enough for most applications) //
105 
106  // Constructor.
107  // The first argument gives an estimate of the maximum number of nodes that can be added
108  // to the graph, and the second argument is an estimate of the maximum number of edges.
109  // The last (optional) argument is the pointer to the function which will be called
110  // if an error occurs; an error message is passed to this function.
111  // If this argument is omitted, exit(1) will be called.
112  //
113  // IMPORTANT: It is possible to add more nodes to the graph than node_num_max
114  // (and node_num_max can be zero). However, if the count is exceeded, then
115  // the internal memory is reallocated (increased by 50%) which is expensive.
116  // Also, temporarily the amount of allocated memory would be more than twice than needed.
117  // Similarly for edges.
118  // If you wish to avoid this overhead, you can download version 2.2, where nodes and edges are stored in blocks.
119  Graph(int node_num_max, int edge_num_max, void (*err_function)(char *) = NULL);
120 
121  // Destructor
122  ~Graph();
123 
124  // Adds node(s) to the graph. By default, one node is added (num=1); then first call returns 0, second call returns 1, and so on.
125  // If num>1, then several nodes are added, and node_id of the first one is returned.
126  // IMPORTANT: see note about the constructor
127  node_id add_node(int num = 1);
128 
129  // Adds a bidirectional edge between 'i' and 'j' with the weights 'cap' and 'rev_cap'.
130  // IMPORTANT: see note about the constructor
131  void add_edge(node_id i, node_id j, captype cap, captype rev_cap);
132 
133  // Adds new edges 'SOURCE->i' and 'i->SINK' with corresponding weights.
134  // Can be called multiple times for each node.
135  // Weights can be negative.
136  // NOTE: the number of such edges is not counted in edge_num_max.
137  // No internal memory is allocated by this call.
138  void add_tweights(node_id i, tcaptype cap_source, tcaptype cap_sink);
139 
140 
141  // Computes the maxflow. Can be called several times.
142  // FOR DESCRIPTION OF reuse_trees, SEE mark_node().
143  // FOR DESCRIPTION OF changed_list, SEE remove_from_changed_list().
144  flowtype maxflow(bool reuse_trees = false, Block<node_id>* changed_list = NULL);
145 
146  // After the maxflow is computed, this function returns to which
147  // segment the node 'i' belongs (Graph<captype,tcaptype,flowtype>::SOURCE or Graph<captype,tcaptype,flowtype>::SINK).
148  //
149  // Occasionally there may be several minimum cuts. If a node can be assigned
150  // to both the source and the sink, then default_segm is returned.
151  termtype what_segment(node_id i, termtype default_segm = SOURCE);
152 
153 
154 
156  // ADVANCED INTERFACE FUNCTIONS //
157  // (provide access to the graph) //
159 
160  //private:
161 public:
162  struct node;
163  struct arc;
164 
165 public:
166 
168  // 1. Reallocating graph. //
170 
171  // Removes all nodes and edges.
172  // After that functions add_node() and add_edge() must be called again.
173  //
174  // Advantage compared to deleting Graph and allocating it again:
175  // no calls to delete/new (which could be quite slow).
176  //
177  // If the graph structure stays the same, then an alternative
178  // is to go through all nodes/edges and set new residual capacities
179  // (see functions below).
180  void reset();
181 
183  // 2. Functions for getting pointers to arcs and for reading graph structure. //
184  // NOTE: adding new arcs may invalidate these pointers (if reallocation //
185  // happens). So it's best not to add arcs while reading graph structure. //
187 
188  // The following two functions return arcs in the same order that they
189  // were added to the graph. NOTE: for each call add_edge(i,j,cap,cap_rev)
190  // the first arc returned will be i->j, and the second j->i.
191  // If there are no more arcs, then the function can still be called, but
192  // the returned arc_id is undetermined.
193  typedef arc* arc_id;
194  arc_id get_first_arc();
195  arc_id get_next_arc(arc_id a);
196 
197  // other functions for reading graph structure
198  int get_node_num() { return node_num; }
199  int get_arc_num() { return (int)(arc_last - arcs); }
200  void get_arc_ends(arc_id a, node_id& i, node_id& j); // returns i,j to that a = i->j
201 
203  // 3. Functions for reading residual capacities. //
205 
206  // returns residual capacity of SOURCE->i minus residual capacity of i->SINK
207  tcaptype get_trcap(node_id i);
208  // returns residual capacity of arc a
209  captype get_rcap(arc* a);
210 
212  // 4. Functions for setting residual capacities. //
213  // NOTE: If these functions are used, the value of the flow //
214  // returned by maxflow() will not be valid! //
216 
217  void set_trcap(node_id i, tcaptype trcap);
218  void set_rcap(arc* a, captype rcap);
219 
221  // 5. Functions related to reusing trees & list of changed nodes. //
223 
224  // If flag reuse_trees is true while calling maxflow(), then search trees
225  // are reused from previous maxflow computation.
226  // In this case before calling maxflow() the user must
227  // specify which parts of the graph have changed by calling mark_node():
228  // add_tweights(i),set_trcap(i) => call mark_node(i)
229  // add_edge(i,j),set_rcap(a) => call mark_node(i); mark_node(j)
230  //
231  // This option makes sense only if a small part of the graph is changed.
232  // The initialization procedure goes only through marked nodes then.
233  //
234  // mark_node(i) can either be called before or after graph modification.
235  // Can be called more than once per node, but calls after the first one
236  // do not have any effect.
237  //
238  // NOTE:
239  // - This option cannot be used in the first call to maxflow().
240  // - It is not necessary to call mark_node() if the change is ``not essential'',
241  // i.e. sign(trcap) is preserved for a node and zero/nonzero status is preserved for an arc.
242  // - To check that you marked all necessary nodes, you can call maxflow(false) after calling maxflow(true).
243  // If everything is correct, the two calls must return the same value of flow. (Useful for debugging).
244  void mark_node(node_id i);
245 
246  // If changed_list is not NULL while calling maxflow(), then the algorithm
247  // keeps a list of nodes which could potentially have changed their segmentation label.
248  // Nodes which are not in the list are guaranteed to keep their old segmentation label (SOURCE or SINK).
249  // Example usage:
250  //
251  // typedef Graph<int,int,int> G;
252  // G* g = new Graph(nodeNum, edgeNum);
253  // Block<G::node_id>* changed_list = new Block<G::node_id>(128);
254  //
255  // ... // add nodes and edges
256  //
257  // g->maxflow(); // first call should be without arguments
258  // for (int iter=0; iter<10; iter++)
259  // {
260  // ... // change graph, call mark_node() accordingly
261  //
262  // g->maxflow(true, changed_list);
263  // G::node_id* ptr;
264  // for (ptr=changed_list->ScanFirst(); ptr; ptr=changed_list->ScanNext())
265  // {
266  // G::node_id i = *ptr; assert(i>=0 && i<nodeNum);
267  // g->remove_from_changed_list(i);
268  // // do something with node i...
269  // if (g->what_segment(i) == G::SOURCE) { ... }
270  // }
271  // changed_list->Reset();
272  // }
273  // delete changed_list;
274  //
275  // NOTE:
276  // - If changed_list option is used, then reuse_trees must be used as well.
277  // - In the example above, the user may omit calls g->remove_from_changed_list(i) and changed_list->Reset() in a given iteration.
278  // Then during the next call to maxflow(true, &changed_list) new nodes will be added to changed_list.
279  // - If the next call to maxflow() does not use option reuse_trees, then calling remove_from_changed_list()
280  // is not necessary. ("changed_list->Reset()" or "delete changed_list" should still be called, though).
281  void remove_from_changed_list(node_id i)
282  {
283  assert(i>=0 && i<node_num && nodes[i].is_in_changed_list);
284  nodes[i].is_in_changed_list = 0;
285  }
286 
287 
288 
289 
290 
291 
295 
296  //private:
297 public:
298  // internal variables and functions
299 
300  struct node
301  {
302  arc *first; // first outcoming arc
303 
304  arc *parent; // node's parent
305  node *next; // pointer to the next active node
306  // (or to itself if it is the last node in the list)
307  int TS; // timestamp showing when DIST was computed
308  int DIST; // distance to the terminal
309  int is_sink : 1; // flag showing whether the node is in the source or in the sink tree (if parent!=NULL)
310  int is_marked : 1; // set by mark_node()
311  int is_in_changed_list : 1; // set by maxflow if
312 
313  tcaptype tr_cap; // if tr_cap > 0 then tr_cap is residual capacity of the arc SOURCE->node
314  // otherwise -tr_cap is residual capacity of the arc node->SINK
315 
316  };
317 
318  struct arc
319  {
320  node *head; // node the arc points to
321  arc *next; // next arc with the same originating node
322  arc *sister; // reverse arc
323 
324  captype r_cap; // residual capacity
325  };
326 
327  struct nodeptr
328  {
331  };
332  static const int NODEPTR_BLOCK_SIZE = 128;
333 
334  node *nodes, *node_last, *node_max; // node_last = nodes+node_num, node_max = nodes+node_num_max;
335  arc *arcs, *arc_last, *arc_max; // arc_last = arcs+2*edge_num, arc_max = arcs+2*edge_num_max;
336 
337  int node_num;
338 
340 
341  void (*error_function)(char *); // this function is called if a error occurs,
342  // with a corresponding error message
343  // (or exit(1) is called if it's NULL)
344 
345  flowtype flow; // total flow
346 
347  // reusing trees & list of changed pixels
348  int maxflow_iteration; // counter
350 
352 
353  node *queue_first[2], *queue_last[2]; // list of active nodes
354  nodeptr *orphan_first, *orphan_last; // list of pointers to orphans
355  int TIME; // monotonically increasing global counter
356 
358 
359  void reallocate_nodes(int num); // num is the number of new nodes
360  void reallocate_arcs();
361 
362  // functions for processing active list
363  void set_active(node *i);
364  node *next_active();
365 
366  // functions for processing orphans list
367  void set_orphan_front(node* i); // add to the beginning of the list
368  void set_orphan_rear(node* i); // add to the end of the list
369 
370  void add_to_changed_list(node* i);
371 
372  void maxflow_init(); // called if reuse_trees == false
373  void maxflow_reuse_trees_init(); // called if reuse_trees == true
374  void augment(arc *middle_arc);
375  void process_source_orphan(node *i);
376  void process_sink_orphan(node *i);
377 
378  void test_consistency(node* current_node=NULL); // debug function
379 };
380 
381 
382 
383 
384 
385 
386 
387 
388 
389 
390 
392 // Implementation - inline functions //
394 
395 template <typename captype, typename tcaptype, typename flowtype>
396  inline Graph<captype, tcaptype, flowtype>::Graph(int node_num_max, int edge_num_max, void (*err_function)(char *))
397  : node_num(0),
398  nodeptr_block(NULL),
399  error_function(err_function)
400 {
401  if (node_num_max < 16) node_num_max = 16;
402  if (edge_num_max < 16) edge_num_max = 16;
403  arcs = (arc*) malloc(2*edge_num_max*sizeof(arc));
404  nodes = (node*) malloc(node_num_max*sizeof(node));
405  if (!nodes || !arcs) {
406  free(arcs);
407  throw std::bad_alloc();
408  }
409  node_last = nodes;
410  node_max = nodes + node_num_max;
411  arc_last = arcs;
412  arc_max = arcs + 2*edge_num_max;
413 
414  maxflow_iteration = 0;
415  flow = 0;
416 }
417 
418 template <typename captype, typename tcaptype, typename flowtype>
420 {
421  if (nodeptr_block)
422  {
423  delete nodeptr_block;
424  nodeptr_block = NULL;
425  }
426  free(nodes);
427  free(arcs);
428 }
429 
430 template <typename captype, typename tcaptype, typename flowtype>
432 {
433  node_last = nodes;
434  arc_last = arcs;
435  node_num = 0;
436 
437  if (nodeptr_block)
438  {
439  delete nodeptr_block;
440  nodeptr_block = NULL;
441  }
442 
443  maxflow_iteration = 0;
444  flow = 0;
445 }
446 
447 template <typename captype, typename tcaptype, typename flowtype>
449 {
450  int node_num_max = (int)(node_max - nodes);
451  node* nodes_old = nodes;
452 
453  node_num_max += node_num_max / 2;
454  if (node_num_max < node_num + num) node_num_max = node_num + num;
455  nodes = (node*) realloc(nodes_old, node_num_max*sizeof(node));
456  if (!nodes) { /*if (error_function) (*error_function)("Not enough memory!");*/ exit(1); }
457 
459  node_max = nodes + node_num_max;
460 
461  if (nodes != nodes_old)
462  {
463  arc* a;
464  for (a=arcs; a<arc_last; a++)
465  {
466  a->head = (node*) ((char*)a->head + (((char*) nodes) - ((char*) nodes_old)));
467  }
468  }
469 }
470 
471 template <typename captype, typename tcaptype, typename flowtype>
473 {
474 
475  int arc_num_max = (int)(arc_max - arcs);
476  int arc_num = (int)(arc_last - arcs);
477  arc* arcs_old = arcs;
478 
479  arc_num_max += arc_num_max / 2; if (arc_num_max & 1) arc_num_max ++;
480  arcs = (arc*) realloc(arcs_old, arc_num_max*sizeof(arc));
481  if (!arcs) { /*if (error_function) (*error_function)("Not enough memory!"); */ exit(1); }
482 
483  arc_last = arcs + arc_num;
484  arc_max = arcs + arc_num_max;
485 
486  if (arcs != arcs_old)
487  {
488  node* i;
489  arc* a;
490  for (i=nodes; i<node_last; i++)
491  {
492  if (i->first) i->first = (arc*) ((char*)i->first + (((char*) arcs) - ((char*) arcs_old)));
493  }
494  for (a=arcs; a<arc_last; a++)
495  {
496  if (a->next) a->next = (arc*) ((char*)a->next + (((char*) arcs) - ((char*) arcs_old)));
497  a->sister = (arc*) ((char*)a->sister + (((char*) arcs) - ((char*) arcs_old)));
498  }
499  }
500 }
501 
502 template <typename captype, typename tcaptype, typename flowtype>
504 {
505  assert(num > 0);
506 
507  if (node_last + num > node_max) reallocate_nodes(num);
508 
509  if (num == 1)
510  {
511  node_last -> first = NULL;
512  node_last -> tr_cap = 0;
513  node_last -> is_marked = 0;
514  node_last -> is_in_changed_list = 0;
515 
516  node_last ++;
517  return node_num ++;
518  }
519  else
520  {
521  memset(node_last, 0, num*sizeof(node));
522 
523  node_id i = node_num;
524  node_num += num;
525  node_last += num;
526  return i;
527  }
528 }
529 
530 template <typename captype, typename tcaptype, typename flowtype>
531  inline void Graph<captype,tcaptype,flowtype>::add_tweights(node_id i, tcaptype cap_source, tcaptype cap_sink)
532 {
533  assert(i >= 0 && i < node_num);
534 
535  tcaptype delta = nodes[i].tr_cap;
536  if (delta > 0) cap_source += delta;
537  else cap_sink -= delta;
538  flow += (cap_source < cap_sink) ? cap_source : cap_sink;
539  nodes[i].tr_cap = cap_source - cap_sink;
540 }
541 
542 template <typename captype, typename tcaptype, typename flowtype>
543  inline void Graph<captype,tcaptype,flowtype>::add_edge(node_id _i, node_id _j, captype cap, captype rev_cap)
544 {
545  assert(_i >= 0 && _i < node_num);
546  assert(_j >= 0 && _j < node_num);
547  assert(_i != _j);
548  assert(cap >= 0);
549  assert(rev_cap >= 0);
550 
551  if (arc_last == arc_max) reallocate_arcs();
552 
553  arc *a = arc_last ++;
554  arc *a_rev = arc_last ++;
555 
556  node* i = nodes + _i;
557  node* j = nodes + _j;
558 
559  a -> sister = a_rev;
560  a_rev -> sister = a;
561  a -> next = i -> first;
562  i -> first = a;
563  a_rev -> next = j -> first;
564  j -> first = a_rev;
565  a -> head = j;
566  a_rev -> head = i;
567  a -> r_cap = cap;
568  a_rev -> r_cap = rev_cap;
569 }
570 
571 template <typename captype, typename tcaptype, typename flowtype>
573 {
574  return arcs;
575 }
576 
577 template <typename captype, typename tcaptype, typename flowtype>
579 {
580  return a + 1;
581 }
582 
583 template <typename captype, typename tcaptype, typename flowtype>
584  inline void Graph<captype,tcaptype,flowtype>::get_arc_ends(arc* a, node_id& i, node_id& j)
585 {
586  assert(a >= arcs && a < arc_last);
587  i = (node_id) (a->sister->head - nodes);
588  j = (node_id) (a->head - nodes);
589 }
590 
591 template <typename captype, typename tcaptype, typename flowtype>
593 {
594  assert(i>=0 && i<node_num);
595  return nodes[i].tr_cap;
596 }
597 
598 template <typename captype, typename tcaptype, typename flowtype>
600 {
601  assert(a >= arcs && a < arc_last);
602  return a->r_cap;
603 }
604 
605 template <typename captype, typename tcaptype, typename flowtype>
606  inline void Graph<captype,tcaptype,flowtype>::set_trcap(node_id i, tcaptype trcap)
607 {
608  assert(i>=0 && i<node_num);
609  nodes[i].tr_cap = trcap;
610 }
611 
612 template <typename captype, typename tcaptype, typename flowtype>
613  inline void Graph<captype,tcaptype,flowtype>::set_rcap(arc* a, captype rcap)
614 {
615  assert(a >= arcs && a < arc_last);
616  a->r_cap = rcap;
617 }
618 
619 
620 template <typename captype, typename tcaptype, typename flowtype>
622 {
623  if (nodes[i].parent)
624  {
625  return (nodes[i].is_sink) ? SINK : SOURCE;
626  }
627  else
628  {
629  return default_segm;
630  }
631 }
632 
633 template <typename captype, typename tcaptype, typename flowtype>
635 {
636  node* i = nodes + _i;
637  if (!i->next)
638  {
639  /* it's not in the list yet */
640  if (queue_last[1]) queue_last[1] -> next = i;
641  else queue_first[1] = i;
642  queue_last[1] = i;
643  i -> next = i;
644  }
645  i->is_marked = 1;
646 }
647 
648 /*
649  special constants for node->parent
650 */
651 #define TERMINAL ( (arc *) 1 ) /* to terminal */
652 #define ORPHAN ( (arc *) 2 ) /* orphan */
653 
654 
655 #define INFINITE_D ((int)(((unsigned)-1)/2)) /* infinite distance to the terminal */
656 
657 /***********************************************************************/
658 
659 /*
660  Functions for processing active list.
661  i->next points to the next node in the list
662  (or to i, if i is the last node in the list).
663  If i->next is NULL iff i is not in the list.
664 
665  There are two queues. Active nodes are added
666  to the end of the second queue and read from
667  the front of the first queue. If the first queue
668  is empty, it is replaced by the second queue
669  (and the second queue becomes empty).
670 */
671 
672 
673 template <typename captype, typename tcaptype, typename flowtype>
675 {
676  if (!i->next)
677  {
678  /* it's not in the list yet */
679  if (queue_last[1]) queue_last[1] -> next = i;
680  else queue_first[1] = i;
681  queue_last[1] = i;
682  i -> next = i;
683  }
684 }
685 
686 /*
687  Returns the next active node.
688  If it is connected to the sink, it stays in the list,
689  otherwise it is removed from the list
690 */
691 template <typename captype, typename tcaptype, typename flowtype>
693 {
694  node *i;
695 
696  while ( 1 )
697  {
698  if (!(i=queue_first[0]))
699  {
700  queue_first[0] = i = queue_first[1];
701  queue_last[0] = queue_last[1];
702  queue_first[1] = NULL;
703  queue_last[1] = NULL;
704  if (!i) return NULL;
705  }
706 
707  /* remove it from the active list */
708  if (i->next == i) queue_first[0] = queue_last[0] = NULL;
709  else queue_first[0] = i -> next;
710  i -> next = NULL;
711 
712  /* a node in the list is active iff it has a parent */
713  if (i->parent) return i;
714  }
715 }
716 
717 /***********************************************************************/
718 
719 template <typename captype, typename tcaptype, typename flowtype>
721 {
722  nodeptr *np;
723  i -> parent = ORPHAN;
724  np = nodeptr_block -> New();
725  np -> ptr = i;
726  np -> next = orphan_first;
727  orphan_first = np;
728 }
729 
730 template <typename captype, typename tcaptype, typename flowtype>
732 {
733  nodeptr *np;
734  i -> parent = ORPHAN;
735  np = nodeptr_block -> New();
736  np -> ptr = i;
737  if (orphan_last) orphan_last -> next = np;
738  else orphan_first = np;
739  orphan_last = np;
740  np -> next = NULL;
741 }
742 
743 /***********************************************************************/
744 
745 template <typename captype, typename tcaptype, typename flowtype>
747 {
748  if (changed_list && !i->is_in_changed_list)
749  {
750  node_id* ptr = changed_list->New();
751  *ptr = (node_id)(i - nodes);
752  i->is_in_changed_list = true;
753  }
754 }
755 
756 /***********************************************************************/
757 
758 template <typename captype, typename tcaptype, typename flowtype>
760 {
761  node *i;
762 
763  queue_first[0] = queue_last[0] = NULL;
764  queue_first[1] = queue_last[1] = NULL;
765  orphan_first = NULL;
766 
767  TIME = 0;
768 
769  for (i=nodes; i<node_last; i++)
770  {
771  i -> next = NULL;
772  i -> is_marked = 0;
773  i -> is_in_changed_list = 0;
774  i -> TS = TIME;
775  if (i->tr_cap > 0)
776  {
777  /* i is connected to the source */
778  i -> is_sink = 0;
779  i -> parent = TERMINAL;
780  set_active(i);
781  i -> DIST = 1;
782  }
783  else if (i->tr_cap < 0)
784  {
785  /* i is connected to the sink */
786  i -> is_sink = 1;
787  i -> parent = TERMINAL;
788  set_active(i);
789  i -> DIST = 1;
790  }
791  else
792  {
793  i -> parent = NULL;
794  }
795  }
796 }
797 
798 template <typename captype, typename tcaptype, typename flowtype>
800 {
801  node* i;
802  node* j;
803  node* queue = queue_first[1];
804  arc* a;
805  nodeptr* np;
806 
807  queue_first[0] = queue_last[0] = NULL;
808  queue_first[1] = queue_last[1] = NULL;
809  orphan_first = orphan_last = NULL;
810 
811  TIME ++;
812 
813  while ((i=queue))
814  {
815  queue = i->next;
816  if (queue == i) queue = NULL;
817  i->next = NULL;
818  i->is_marked = 0;
819  set_active(i);
820 
821  if (i->tr_cap == 0)
822  {
823  if (i->parent) set_orphan_rear(i);
824  continue;
825  }
826 
827  if (i->tr_cap > 0)
828  {
829  if (!i->parent || i->is_sink)
830  {
831  i->is_sink = 0;
832  for (a=i->first; a; a=a->next)
833  {
834  j = a->head;
835  if (!j->is_marked)
836  {
837  if (j->parent == a->sister) set_orphan_rear(j);
838  if (j->parent && j->is_sink && a->r_cap > 0) set_active(j);
839  }
840  }
842  }
843  }
844  else
845  {
846  if (!i->parent || !i->is_sink)
847  {
848  i->is_sink = 1;
849  for (a=i->first; a; a=a->next)
850  {
851  j = a->head;
852  if (!j->is_marked)
853  {
854  if (j->parent == a->sister) set_orphan_rear(j);
855  if (j->parent && !j->is_sink && a->sister->r_cap > 0) set_active(j);
856  }
857  }
859  }
860  }
861  i->parent = TERMINAL;
862  i -> TS = TIME;
863  i -> DIST = 1;
864  }
865 
866  //test_consistency();
867 
868  /* adoption */
869  while ((np=orphan_first))
870  {
871  orphan_first = np -> next;
872  i = np -> ptr;
873  nodeptr_block -> Delete(np);
874  if (!orphan_first) orphan_last = NULL;
875  if (i->is_sink) process_sink_orphan(i);
876  else process_source_orphan(i);
877  }
878  /* adoption end */
879 
880  //test_consistency();
881 }
882 
883 template <typename captype, typename tcaptype, typename flowtype>
885 {
886  node *i;
887  arc *a;
888  tcaptype bottleneck;
889 
890 
891  /* 1. Finding bottleneck capacity */
892  /* 1a - the source tree */
893  bottleneck = middle_arc -> r_cap;
894  for (i=middle_arc->sister->head; ; i=a->head)
895  {
896  a = i -> parent;
897  if (a == TERMINAL) break;
898  if (bottleneck > a->sister->r_cap) bottleneck = a -> sister -> r_cap;
899  }
900  if (bottleneck > i->tr_cap) bottleneck = i -> tr_cap;
901  /* 1b - the sink tree */
902  for (i=middle_arc->head; ; i=a->head)
903  {
904  a = i -> parent;
905  if (a == TERMINAL) break;
906  if (bottleneck > a->r_cap) bottleneck = a -> r_cap;
907  }
908  if (bottleneck > - i->tr_cap) bottleneck = - i -> tr_cap;
909 
910 
911  /* 2. Augmenting */
912  /* 2a - the source tree */
913  middle_arc -> sister -> r_cap += bottleneck;
914  middle_arc -> r_cap -= bottleneck;
915  for (i=middle_arc->sister->head; ; i=a->head)
916  {
917  a = i -> parent;
918  if (a == TERMINAL) break;
919  a -> r_cap += bottleneck;
920  a -> sister -> r_cap -= bottleneck;
921  if (!a->sister->r_cap)
922  {
923  set_orphan_front(i); // add i to the beginning of the adoption list
924  }
925  }
926  i -> tr_cap -= bottleneck;
927  if (!i->tr_cap)
928  {
929  set_orphan_front(i); // add i to the beginning of the adoption list
930  }
931  /* 2b - the sink tree */
932  for (i=middle_arc->head; ; i=a->head)
933  {
934  a = i -> parent;
935  if (a == TERMINAL) break;
936  a -> sister -> r_cap += bottleneck;
937  a -> r_cap -= bottleneck;
938  if (!a->r_cap)
939  {
940  set_orphan_front(i); // add i to the beginning of the adoption list
941  }
942  }
943  i -> tr_cap += bottleneck;
944  if (!i->tr_cap)
945  {
946  set_orphan_front(i); // add i to the beginning of the adoption list
947  }
948 
949 
950  flow += bottleneck;
951 }
952 
953 /***********************************************************************/
954 
955 template <typename captype, typename tcaptype, typename flowtype>
957 {
958  node *j;
959  arc *a0, *a0_min = NULL, *a;
960  int d, d_min = INFINITE_D;
961 
962  /* trying to find a new parent */
963  for (a0=i->first; a0; a0=a0->next)
964  if (a0->sister->r_cap)
965  {
966  j = a0 -> head;
967  if (!j->is_sink && (a=j->parent))
968  {
969  /* checking the origin of j */
970  d = 0;
971  while ( 1 )
972  {
973  if (j->TS == TIME)
974  {
975  d += j -> DIST;
976  break;
977  }
978  a = j -> parent;
979  d ++;
980  if (a==TERMINAL)
981  {
982  j -> TS = TIME;
983  j -> DIST = 1;
984  break;
985  }
986  if (a==ORPHAN) { d = INFINITE_D; break; }
987  j = a -> head;
988  }
989  if (d<INFINITE_D) /* j originates from the source - done */
990  {
991  if (d<d_min)
992  {
993  a0_min = a0;
994  d_min = d;
995  }
996  /* set marks along the path */
997  for (j=a0->head; j->TS!=TIME; j=j->parent->head)
998  {
999  j -> TS = TIME;
1000  j -> DIST = d --;
1001  }
1002  }
1003  }
1004  }
1005 
1006  if ((i->parent = a0_min))
1007  {
1008  i -> TS = TIME;
1009  i -> DIST = d_min + 1;
1010  }
1011  else
1012  {
1013  /* no parent is found */
1015 
1016  /* process neighbors */
1017  for (a0=i->first; a0; a0=a0->next)
1018  {
1019  j = a0 -> head;
1020  if (!j->is_sink && (a=j->parent))
1021  {
1022  if (a0->sister->r_cap) set_active(j);
1023  if (a!=TERMINAL && a!=ORPHAN && a->head==i)
1024  {
1025  set_orphan_rear(j); // add j to the end of the adoption list
1026  }
1027  }
1028  }
1029  }
1030 }
1031 
1032 template <typename captype, typename tcaptype, typename flowtype>
1034 {
1035  node *j;
1036  arc *a0, *a0_min = NULL, *a;
1037  int d, d_min = INFINITE_D;
1038 
1039  /* trying to find a new parent */
1040  for (a0=i->first; a0; a0=a0->next)
1041  if (a0->r_cap)
1042  {
1043  j = a0 -> head;
1044  if (j->is_sink && (a=j->parent))
1045  {
1046  /* checking the origin of j */
1047  d = 0;
1048  while ( 1 )
1049  {
1050  if (j->TS == TIME)
1051  {
1052  d += j -> DIST;
1053  break;
1054  }
1055  a = j -> parent;
1056  d ++;
1057  if (a==TERMINAL)
1058  {
1059  j -> TS = TIME;
1060  j -> DIST = 1;
1061  break;
1062  }
1063  if (a==ORPHAN) { d = INFINITE_D; break; }
1064  j = a -> head;
1065  }
1066  if (d<INFINITE_D) /* j originates from the sink - done */
1067  {
1068  if (d<d_min)
1069  {
1070  a0_min = a0;
1071  d_min = d;
1072  }
1073  /* set marks along the path */
1074  for (j=a0->head; j->TS!=TIME; j=j->parent->head)
1075  {
1076  j -> TS = TIME;
1077  j -> DIST = d --;
1078  }
1079  }
1080  }
1081  }
1082 
1083  if ((i->parent = a0_min))
1084  {
1085  i -> TS = TIME;
1086  i -> DIST = d_min + 1;
1087  }
1088  else
1089  {
1090  /* no parent is found */
1092 
1093  /* process neighbors */
1094  for (a0=i->first; a0; a0=a0->next)
1095  {
1096  j = a0 -> head;
1097  if (j->is_sink && (a=j->parent))
1098  {
1099  if (a0->r_cap) set_active(j);
1100  if (a!=TERMINAL && a!=ORPHAN && a->head==i)
1101  {
1102  set_orphan_rear(j); // add j to the end of the adoption list
1103  }
1104  }
1105  }
1106  }
1107 }
1108 
1109 /***********************************************************************/
1110 
1111 template <typename captype, typename tcaptype, typename flowtype>
1112  inline flowtype Graph<captype,tcaptype,flowtype>::maxflow(bool reuse_trees, Block<node_id>* _changed_list)
1113 {
1114  node *i, *j, *current_node = NULL;
1115  arc *a;
1116  nodeptr *np, *np_next;
1117 
1118  if (!nodeptr_block)
1119  {
1121  }
1122 
1123  changed_list = _changed_list;
1124  if (maxflow_iteration == 0 && reuse_trees) { /*if (error_function) (*error_function)("reuse_trees cannot be used in the first call to maxflow()!");*/ exit(1); }
1125  if (changed_list && !reuse_trees) { /*if (error_function) (*error_function)("changed_list cannot be used without reuse_trees!");*/ exit(1); }
1126 
1127  if (reuse_trees) maxflow_reuse_trees_init();
1128  else maxflow_init();
1129 
1130  // main loop
1131  while ( 1 )
1132  {
1133  // test_consistency(current_node);
1134 
1135  if ((i=current_node))
1136  {
1137  i -> next = NULL; /* remove active flag */
1138  if (!i->parent) i = NULL;
1139  }
1140  if (!i)
1141  {
1142  if (!(i = next_active())) break;
1143  }
1144 
1145  /* growth */
1146  if (!i->is_sink)
1147  {
1148  /* grow source tree */
1149  for (a=i->first; a; a=a->next)
1150  if (a->r_cap)
1151  {
1152  j = a -> head;
1153  if (!j->parent)
1154  {
1155  j -> is_sink = 0;
1156  j -> parent = a -> sister;
1157  j -> TS = i -> TS;
1158  j -> DIST = i -> DIST + 1;
1159  set_active(j);
1161  }
1162  else if (j->is_sink) break;
1163  else if (j->TS <= i->TS &&
1164  j->DIST > i->DIST)
1165  {
1166  /* heuristic - trying to make the distance from j to the source shorter */
1167  j -> parent = a -> sister;
1168  j -> TS = i -> TS;
1169  j -> DIST = i -> DIST + 1;
1170  }
1171  }
1172  }
1173  else
1174  {
1175  /* grow sink tree */
1176  for (a=i->first; a; a=a->next)
1177  if (a->sister->r_cap)
1178  {
1179  j = a -> head;
1180  if (!j->parent)
1181  {
1182  j -> is_sink = 1;
1183  j -> parent = a -> sister;
1184  j -> TS = i -> TS;
1185  j -> DIST = i -> DIST + 1;
1186  set_active(j);
1188  }
1189  else if (!j->is_sink) { a = a -> sister; break; }
1190  else if (j->TS <= i->TS &&
1191  j->DIST > i->DIST)
1192  {
1193  /* heuristic - trying to make the distance from j to the sink shorter */
1194  j -> parent = a -> sister;
1195  j -> TS = i -> TS;
1196  j -> DIST = i -> DIST + 1;
1197  }
1198  }
1199  }
1200 
1201  TIME ++;
1202 
1203  if (a)
1204  {
1205  i -> next = i; /* set active flag */
1206  current_node = i;
1207 
1208  /* augmentation */
1209  augment(a);
1210  /* augmentation end */
1211 
1212  /* adoption */
1213  while ((np=orphan_first))
1214  {
1215  np_next = np -> next;
1216  np -> next = NULL;
1217 
1218  while ((np=orphan_first))
1219  {
1220  orphan_first = np -> next;
1221  i = np -> ptr;
1222  nodeptr_block -> Delete(np);
1223  if (!orphan_first) orphan_last = NULL;
1224  if (i->is_sink) process_sink_orphan(i);
1225  else process_source_orphan(i);
1226  }
1227 
1228  orphan_first = np_next;
1229  }
1230  /* adoption end */
1231  }
1232  else current_node = NULL;
1233  }
1234  // test_consistency();
1235 
1236  if (!reuse_trees || (maxflow_iteration % 64) == 0)
1237  {
1238  delete nodeptr_block;
1239  nodeptr_block = NULL;
1240  }
1241 
1242  maxflow_iteration ++;
1243  return flow;
1244 }
1245 
1246 /***********************************************************************/
1247 
1248 
1249 template <typename captype, typename tcaptype, typename flowtype>
1251 {
1252  node *i;
1253  arc *a;
1254  int r;
1255  int num1 = 0, num2 = 0;
1256 
1257  // test whether all nodes i with i->next!=NULL are indeed in the queue
1258  for (i=nodes; i<node_last; i++)
1259  {
1260  if (i->next || i==current_node) num1 ++;
1261  }
1262  for (r=0; r<3; r++)
1263  {
1264  i = (r == 2) ? current_node : queue_first[r];
1265  if (i)
1266  for ( ; ; i=i->next)
1267  {
1268  num2 ++;
1269  if (i->next == i)
1270  {
1271  if (r<2) assert(i == queue_last[r]);
1272  else assert(i == current_node);
1273  break;
1274  }
1275  }
1276  }
1277  assert(num1 == num2);
1278 
1279  for (i=nodes; i<node_last; i++)
1280  {
1281  // test whether all edges in seach trees are non-saturated
1282  if (i->parent == NULL) {}
1283  else if (i->parent == ORPHAN) {}
1284  else if (i->parent == TERMINAL)
1285  {
1286  if (!i->is_sink) assert(i->tr_cap > 0);
1287  else assert(i->tr_cap < 0);
1288  }
1289  else
1290  {
1291  if (!i->is_sink) assert (i->parent->sister->r_cap > 0);
1292  else assert (i->parent->r_cap > 0);
1293  }
1294  // test whether passive nodes in search trees have neighbors in
1295  // a different tree through non-saturated edges
1296  if (i->parent && !i->next)
1297  {
1298  if (!i->is_sink)
1299  {
1300  assert(i->tr_cap >= 0);
1301  for (a=i->first; a; a=a->next)
1302  {
1303  if (a->r_cap > 0) assert(a->head->parent && !a->head->is_sink);
1304  }
1305  }
1306  else
1307  {
1308  assert(i->tr_cap <= 0);
1309  for (a=i->first; a; a=a->next)
1310  {
1311  if (a->sister->r_cap > 0) assert(a->head->parent && a->head->is_sink);
1312  }
1313  }
1314  }
1315  // test marking invariants
1316  if (i->parent && i->parent!=ORPHAN && i->parent!=TERMINAL)
1317  {
1318  assert(i->TS <= i->parent->head->TS);
1319  if (i->TS == i->parent->head->TS) assert(i->DIST > i->parent->head->DIST);
1320  }
1321  }
1322 }
1323 
1324 }
void process_sink_orphan(node *i)
Definition: animaGraph.h:1033
int maxflow_iteration
Definition: animaGraph.h:348
arc * arc_last
Definition: animaGraph.h:335
node * node_max
Definition: animaGraph.h:334
Block< node_id > * changed_list
Definition: animaGraph.h:349
void set_orphan_front(node *i)
Definition: animaGraph.h:720
node * next_active()
Definition: animaGraph.h:692
arc * arc_max
Definition: animaGraph.h:335
nodeptr * orphan_first
Definition: animaGraph.h:354
void set_active(node *i)
Definition: animaGraph.h:674
node * node_last
Definition: animaGraph.h:334
void maxflow_init()
Definition: animaGraph.h:759
int get_arc_num()
Definition: animaGraph.h:199
DBlock< nodeptr > * nodeptr_block
Definition: animaGraph.h:339
node * queue_first[2]
Definition: animaGraph.h:353
flowtype flow
Definition: animaGraph.h:345
void set_trcap(node_id i, tcaptype trcap)
Definition: animaGraph.h:606
void mark_node(node_id i)
Definition: animaGraph.h:634
#define INFINITE_D
Definition: animaGraph.h:655
void augment(arc *middle_arc)
Definition: animaGraph.h:884
void reallocate_nodes(int num)
Definition: animaGraph.h:448
void process_source_orphan(node *i)
Definition: animaGraph.h:956
static const int NODEPTR_BLOCK_SIZE
Definition: animaGraph.h:332
void get_arc_ends(arc_id a, node_id &i, node_id &j)
Definition: animaGraph.h:584
void set_orphan_rear(node *i)
Definition: animaGraph.h:731
nodeptr * orphan_last
Definition: animaGraph.h:354
void reallocate_arcs()
Definition: animaGraph.h:472
node * queue_last[2]
Definition: animaGraph.h:353
void add_edge(node_id i, node_id j, captype cap, captype rev_cap)
Definition: animaGraph.h:543
#define TERMINAL
Definition: animaGraph.h:651
void(* error_function)(char *)
Definition: animaGraph.h:341
captype get_rcap(arc *a)
Definition: animaGraph.h:599
flowtype maxflow(bool reuse_trees=false, Block< node_id > *changed_list=NULL)
Definition: animaGraph.h:1112
node_id add_node(int num=1)
Definition: animaGraph.h:503
void reset()
Definition: animaGraph.h:431
void remove_from_changed_list(node_id i)
Definition: animaGraph.h:281
arc_id get_next_arc(arc_id a)
Definition: animaGraph.h:578
void maxflow_reuse_trees_init()
Definition: animaGraph.h:799
void test_consistency(node *current_node=NULL)
Definition: animaGraph.h:1250
Type * New(int num=1)
Definition: animaBlock.h:116
void add_to_changed_list(node *i)
Definition: animaGraph.h:746
Graph(int node_num_max, int edge_num_max, void(*err_function)(char *)=NULL)
Definition: animaGraph.h:396
void add_tweights(node_id i, tcaptype cap_source, tcaptype cap_sink)
Definition: animaGraph.h:531
arc_id get_first_arc()
Definition: animaGraph.h:572
int get_node_num()
Definition: animaGraph.h:198
#define ORPHAN
Definition: animaGraph.h:652
arc * arc_id
Definition: animaGraph.h:193
void set_rcap(arc *a, captype rcap)
Definition: animaGraph.h:613
tcaptype get_trcap(node_id i)
Definition: animaGraph.h:592
termtype what_segment(node_id i, termtype default_segm=SOURCE)
Definition: animaGraph.h:621
node * nodes
Definition: animaGraph.h:334