Surface Mesh Segmentation
 All Classes Namespaces Files Functions Typedefs Pages
graph.h
1 /*
2 ###################################################################
3 # #
4 # MAXFLOW - software for computing mincut/maxflow in a graph #
5 # Version 2.21 #
6 # http://www.cs.ucl.ac.uk/staff/V.Kolmogorov/software.html #
7 # #
8 # Yuri Boykov (yuri@csd.uwo.ca) #
9 # Vladimir Kolmogorov (v.kolmogorov@cs.ucl.ac.uk) #
10 # 2001 #
11 # #
12 ###################################################################
13 
14 1. Introduction.
15 
16 This software library implements the maxflow algorithm
17 described in
18 
19  An Experimental Comparison of Min-Cut/Max-Flow Algorithms
20  for Energy Minimization in Vision.
21  Yuri Boykov and Vladimir Kolmogorov.
22  In IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI),
23  September 2004
24 
25 This algorithm was developed by Yuri Boykov and Vladimir Kolmogorov
26 at Siemens Corporate Research. To make it available for public use,
27 it was later reimplemented by Vladimir Kolmogorov based on open publications.
28 
29 If you use this software for research purposes, you should cite
30 the aforementioned paper in any resulting publication.
31 
32 Tested under windows, Visual C++ 6.0 compiler and unix (SunOS 5.8
33 and RedHat Linux 7.0, GNU c++ compiler).
34 
35 ##################################################################
36 
37 2. Licence.
38 
39 Copyright UCL Business PLC
40 
41 This program is available under dual licence:
42 
43 1) Under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.
44 Note that any program that incorporates the code under this licence must, under the terms of the GNU GPL, be released under a licence compatible with the GPL. GNU GPL does not permit incorporating this program into proprietary programs. If you wish to do this, please see the alternative licence available below.
45 GNU General Public License can be found at http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
46 
47 2) Proprietary Licence from UCL Business PLC.
48 To enable programers to include the MaxFlow software in a proprietary system (which is not allowed by the GNU GPL), this licence gives you the right to incorporate the software in your program and distribute under any licence of your choosing. The full terms of the licence and applicable fee, are available from the Licensors at: http://www.uclb-elicensing.com/optimisation_software/maxflow_computervision.html
49 
50 ##################################################################
51 
52 3. Graph representation.
53 
54 There are two versions of the algorithm using different
55 graph representations (adjacency list and forward star).
56 The former one uses more than twice as much memory as the
57 latter one but is 10-20% faster.
58 
59 Memory allocation (assuming that all capacities are 'short' - 2 bytes):
60 
61  | Nodes | Arcs
62 ------------------------------------------
63 Adjacency list | *24 bytes | *14 bytes
64 Forward star | *28 bytes | 6 bytes
65 
66 (* means that often it should be rounded up to be a multiple of 4
67 - some compilers (e.g. Visual C++) seem to round up elements
68 of arrays unless the are structures containing only char[].)
69 
70 Note that arcs are always added in pairs - in forward and reverse directions.
71 Arcs between nodes and terminals (the source and the sink) are
72 not stored as arcs, but rather as a part of nodes.
73 
74 The assumption for the forward star representation is that
75 the maximum number of arcs per node (except the source
76 and the sink) is much less than ARC_BLOCK_SIZE (1024 by default).
77 
78 Both versions have the same interface.
79 
80 ##################################################################
81 
82 4. Example usage.
83 
84 This section shows how to use the library to compute
85 a minimum cut on the following graph:
86 
87  SOURCE
88  / \
89  1/ \2
90  / 3 \
91  node0 -----> node1
92  | <----- |
93  | 4 |
94  \ /
95  5\ /6
96  \ /
97  SINK
98 
100 
101 #include <stdio.h>
102 #include "graph.h"
103 
104 void main()
105 {
106  Graph::node_id nodes[2];
107  Graph *g = new Graph();
108 
109  nodes[0] = g -> add_node();
110  nodes[1] = g -> add_node();
111  g -> set_tweights(nodes[0], 1, 5);
112  g -> set_tweights(nodes[1], 2, 6);
113  g -> add_edge(nodes[0], nodes[1], 3, 4);
114 
115  Graph::flowtype flow = g -> maxflow();
116 
117  printf("Flow = %d\n", flow);
118  printf("Minimum cut:\n");
119  if (g->what_segment(nodes[0]) == Graph::SOURCE)
120  printf("node0 is in the SOURCE set\n");
121  else
122  printf("node0 is in the SINK set\n");
123  if (g->what_segment(nodes[1]) == Graph::SOURCE)
124  printf("node1 is in the SOURCE set\n");
125  else
126  printf("node1 is in the SINK set\n");
127 
128  delete g;
129 }
130 
132 */
133 
134 /* block.h */
135 /*
136  Template classes Block and DBlock
137  Implement adding and deleting items of the same type in blocks.
138 
139  If there there are many items then using Block or DBlock
140  is more efficient than using 'new' and 'delete' both in terms
141  of memory and time since
142  (1) On some systems there is some minimum amount of memory
143  that 'new' can allocate (e.g., 64), so if items are
144  small that a lot of memory is wasted.
145  (2) 'new' and 'delete' are designed for items of varying size.
146  If all items has the same size, then an algorithm for
147  adding and deleting can be made more efficient.
148  (3) All Block and DBlock functions are inline, so there are
149  no extra function calls.
150 
151  Differences between Block and DBlock:
152  (1) DBlock allows both adding and deleting items,
153  whereas Block allows only adding items.
154  (2) Block has an additional operation of scanning
155  items added so far (in the order in which they were added).
156  (3) Block allows to allocate several consecutive
157  items at a time, whereas DBlock can add only a single item.
158 
159  Note that no constructors or destructors are called for items.
160 
161  Example usage for items of type 'MyType':
162 
164  #include "block.h"
165  #define BLOCK_SIZE 1024
166  typedef struct { int a, b; } MyType;
167  MyType *ptr, *array[10000];
168 
169  ...
170 
171  Block<MyType> *block = new Block<MyType>(BLOCK_SIZE);
172 
173  // adding items
174  for (int i=0; i<sizeof(array); i++)
175  {
176  ptr = block -> New();
177  ptr -> a = ptr -> b = rand();
178  }
179 
180  // reading items
181  for (ptr=block->ScanFirst(); ptr; ptr=block->ScanNext())
182  {
183  printf("%d %d\n", ptr->a, ptr->b);
184  }
185 
186  delete block;
187 
188  ...
189 
190  DBlock<MyType> *dblock = new DBlock<MyType>(BLOCK_SIZE);
191 
192  // adding items
193  for (int i=0; i<sizeof(array); i++)
194  {
195  array[i] = dblock -> New();
196  }
197 
198  // deleting items
199  for (int i=0; i<sizeof(array); i+=2)
200  {
201  dblock -> Delete(array[i]);
202  }
203 
204  // adding items
205  for (int i=0; i<sizeof(array); i++)
206  {
207  array[i] = dblock -> New();
208  }
209 
210  delete dblock;
211 
213 
214  Note that DBlock deletes items by marking them as
215  empty (i.e., by adding them to the list of free items),
216  so that this memory could be used for subsequently
217  added items. Thus, at each moment the memory allocated
218  is determined by the maximum number of items allocated
219  simultaneously at earlier moments. All memory is
220  deallocated only when the destructor is called.
221 */
222 
223 #ifndef __BLOCK_H__
224 #define __BLOCK_H__
225 
226 #include <stdlib.h>
227 
228 /***********************************************************************/
229 /***********************************************************************/
230 /***********************************************************************/
231 
232 template <class Type> class Block
233 {
234 public:
235  /* Constructor. Arguments are the block size and
236  (optionally) the pointer to the function which
237  will be called if allocation failed; the message
238  passed to this function is "Not enough memory!" */
239  Block(int size, void (*err_function)(char *) = NULL) { first = last = NULL; block_size = size; error_function = err_function; }
240 
241  /* Destructor. Deallocates all items added so far */
242  ~Block() { while (first) { block *next = first -> next; delete[] ((char*)first); first = next; } }
243 
244  /* Allocates 'num' consecutive items; returns pointer
245  to the first item. 'num' cannot be greater than the
246  block size since items must fit in one block */
247  Type *New(int num = 1)
248  {
249  Type *t;
250 
251  if (!last || last->current + num > last->last)
252  {
253  if (last && last->next) last = last -> next;
254  else
255  {
256  block *next = (block *) new char [sizeof(block) + (block_size-1)*sizeof(Type)];
257  if (!next) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
258  if (last) last -> next = next;
259  else first = next;
260  last = next;
261  last -> current = & ( last -> data[0] );
262  last -> last = last -> current + block_size;
263  last -> next = NULL;
264  }
265  }
266 
267  t = last -> current;
268  last -> current += num;
269  return t;
270  }
271 
272  /* Returns the first item (or NULL, if no items were added) */
273  Type *ScanFirst()
274  {
275  for (scan_current_block=first; scan_current_block; scan_current_block = scan_current_block->next)
276  {
277  scan_current_data = & ( scan_current_block -> data[0] );
278  if (scan_current_data < scan_current_block -> current) return scan_current_data ++;
279  }
280  return NULL;
281  }
282 
283  /* Returns the next item (or NULL, if all items have been read)
284  Can be called only if previous ScanFirst() or ScanNext()
285  call returned not NULL. */
286  Type *ScanNext()
287  {
288  while (scan_current_data >= scan_current_block -> current)
289  {
290  scan_current_block = scan_current_block -> next;
291  if (!scan_current_block) return NULL;
292  scan_current_data = & ( scan_current_block -> data[0] );
293  }
294  return scan_current_data ++;
295  }
296 
297  /* Marks all elements as empty */
298  void Reset()
299  {
300  block *b;
301  if (!first) return;
302  for (b=first; ; b=b->next)
303  {
304  b -> current = & ( b -> data[0] );
305  if (b == last) break;
306  }
307  last = first;
308  }
309 
310 /***********************************************************************/
311 
312 private:
313 
314  typedef struct block_st
315  {
316  Type *current, *last;
317  struct block_st *next;
318  Type data[1];
319  } block;
320 
321  int block_size;
322  block *first;
323  block *last;
324 
325  block *scan_current_block;
326  Type *scan_current_data;
327 
328  void (*error_function)(char *);
329 };
330 
331 /***********************************************************************/
332 /***********************************************************************/
333 /***********************************************************************/
334 
335 template <class Type> class DBlock
336 {
337 public:
338  /* Constructor. Arguments are the block size and
339  (optionally) the pointer to the function which
340  will be called if allocation failed; the message
341  passed to this function is "Not enough memory!" */
342  DBlock(int size, void (*err_function)(char *) = NULL) { first = NULL; first_free = NULL; block_size = size; error_function = err_function; }
343 
344  /* Destructor. Deallocates all items added so far */
345  ~DBlock() { while (first) { block *next = first -> next; delete[] ((char*)first); first = next; } }
346 
347  /* Allocates one item */
348  Type *New()
349  {
350  block_item *item;
351 
352  if (!first_free)
353  {
354  block *next = first;
355  first = (block *) new char [sizeof(block) + (block_size-1)*sizeof(block_item)];
356  if (!first) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
357  first_free = & (first -> data[0] );
358  for (item=first_free; item<first_free+block_size-1; item++)
359  item -> next_free = item + 1;
360  item -> next_free = NULL;
361  first -> next = next;
362  }
363 
364  item = first_free;
365  first_free = item -> next_free;
366  return (Type *) item;
367  }
368 
369  /* Deletes an item allocated previously */
370  void Delete(Type *t)
371  {
372  ((block_item *) t) -> next_free = first_free;
373  first_free = (block_item *) t;
374  }
375 
376 /***********************************************************************/
377 
378 private:
379 
380  typedef union block_item_st
381  {
382  Type t;
383  block_item_st *next_free;
384  } block_item;
385 
386  typedef struct block_st
387  {
388  struct block_st *next;
389  block_item data[1];
390  } block;
391 
392  int block_size;
393  block *first;
394  block_item *first_free;
395 
396  void (*error_function)(char *);
397 };
398 
399 
400 #endif
401 
402 
403 
404 /* graph.h */
405 /*
406  This software library implements the maxflow algorithm
407  described in
408 
409  An Experimental Comparison of Min-Cut/Max-Flow Algorithms
410  for Energy Minimization in Vision.
411  Yuri Boykov and Vladimir Kolmogorov.
412  In IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI),
413  September 2004
414 
415  This algorithm was developed by Yuri Boykov and Vladimir Kolmogorov
416  at Siemens Corporate Research. To make it available for public use,
417  it was later reimplemented by Vladimir Kolmogorov based on open publications.
418 
419  If you use this software for research purposes, you should cite
420  the aforementioned paper in any resulting publication.
421 
422  ----------------------------------------------------------------
423 
424  For description, license, example usage, discussion of graph representation and memory usage see README.TXT.
425 */
426 
427 #ifndef __GRAPH_H__
428 #define __GRAPH_H__
429 
430 //#include "block.h"
431 #include <stdio.h>
432 /*
433  Nodes, arcs and pointers to nodes are
434  added in blocks for memory and time efficiency.
435  Below are numbers of items in blocks
436 */
437 #define NODE_BLOCK_SIZE 512
438 #define ARC_BLOCK_SIZE 1024
439 #define NODEPTR_BLOCK_SIZE 128
440 
441 class Graph
442 {
443 public:
444  typedef enum
445  {
446  SOURCE = 0,
447  SINK = 1
448  } termtype; /* terminals */
449 
450  /* Type of edge weights.
451  Can be changed to char, int, float, double, ... */
452  typedef double captype;
453  /* Type of total flow */
454  typedef double flowtype;
455 
456  typedef void * node_id;
457 
458  /* interface functions */
459 
460  /* Constructor. Optional argument is the pointer to the
461  function which will be called if an error occurs;
462  an error message is passed to this function. If this
463  argument is omitted, exit(1) will be called. */
464  Graph(void (*err_function)(char *) = NULL);
465 
466  /* Destructor */
467  ~Graph();
468 
469  /* Adds a node to the graph */
470  node_id add_node();
471 
472  /* Adds a bidirectional edge between 'from' and 'to'
473  with the weights 'cap' and 'rev_cap' */
474  void add_edge(node_id from, node_id to, captype cap, captype rev_cap);
475 
476  /* Sets the weights of the edges 'SOURCE->i' and 'i->SINK'
477  Can be called at most once for each node before any call to 'add_tweights'.
478  Weights can be negative */
479  void set_tweights(node_id i, captype cap_source, captype cap_sink);
480 
481  /* Adds new edges 'SOURCE->i' and 'i->SINK' with corresponding weights
482  Can be called multiple times for each node.
483  Weights can be negative */
484  void add_tweights(node_id i, captype cap_source, captype cap_sink);
485 
486  /* After the maxflow is computed, this function returns to which
487  segment the node 'i' belongs (Graph::SOURCE or Graph::SINK) */
488  termtype what_segment(node_id i);
489 
490  /* Computes the maxflow. Can be called only once. */
491  flowtype maxflow();
492 
493 /***********************************************************************/
494 /***********************************************************************/
495 /***********************************************************************/
496 
497 private:
498  /* internal variables and functions */
499 
500  struct arc_forward_st;
501  struct arc_reverse_st;
502 
503 #define IS_ODD(a) ((int)(a) & 1)
504 #define MAKE_ODD(a) ((arc_forward *) ((int)(a) | 1))
505 #define MAKE_EVEN(a) ((arc_forward *) ((int)(a) & (~1)))
506 #define MAKE_ODD_REV(a) ((arc_reverse *) ((int)(a) | 1))
507 #define MAKE_EVEN_REV(a) ((arc_reverse *) ((int)(a) & (~1)))
508 
509  /* node structure */
510  typedef struct node_st
511  {
512  /*
513  Usually i->first_out is the first outgoing
514  arc, and (i+1)->first_out-1 is the last outgoing arc.
515  However, it is not always possible, since
516  arcs are allocated in blocks, so arcs corresponding
517  to two consecutive nodes may be in different blocks.
518 
519  If outgoing arcs for i are last in the arc block,
520  then a different mechanism is used. i->first_out
521  is odd in this case; the first outgoing arc
522  is (a+1), and the last outgoing arc is
523  ((arc_forward *)(a->shift))-1, where
524  a = (arc_forward *) (((char *)(i->first_out)) + 1);
525 
526  Similar mechanism is used for incoming arcs.
527  */
528  arc_forward_st *first_out; /* first outcoming arc */
529  arc_reverse_st *first_in; /* first incoming arc */
530 
531  arc_forward_st *parent; /* describes node's parent
532  if IS_ODD(parent) then MAKE_EVEN(parent) points to 'arc_reverse',
533  otherwise parent points to 'arc_forward' */
534 
535  node_st *next; /* pointer to the next active node
536  (or to itself if it is the last node in the list) */
537 
538  int TS; /* timestamp showing when DIST was computed */
539  int DIST; /* distance to the terminal */
540  short is_sink; /* flag showing whether the node is in the source or in the sink tree */
541 
542  captype tr_cap; /* if tr_cap > 0 then tr_cap is residual capacity of the arc SOURCE->node
543  otherwise -tr_cap is residual capacity of the arc node->SINK */
544  } node;
545 
546  /* arc structures */
547 #define NEIGHBOR_NODE(i, shift) ((node *) ((char *)(i) + (shift)))
548 #define NEIGHBOR_NODE_REV(i, shift) ((node *) ((char *)(i) - (shift)))
549  typedef struct arc_forward_st
550  {
551  int shift; /* node_to = NEIGHBOR_NODE(node_from, shift) */
552  captype r_cap; /* residual capacity */
553  captype r_rev_cap; /* residual capacity of the reverse arc*/
554  } arc_forward;
555 
556  typedef struct arc_reverse_st
557  {
558  arc_forward *sister; /* reverse arc */
559  } arc_reverse;
560 
561  /* 'pointer to node' structure */
562  typedef struct nodeptr_st
563  {
564  node_st *ptr;
565  nodeptr_st *next;
566  } nodeptr;
567 
568  typedef struct node_block_st
569  {
570  node *current;
571  struct node_block_st *next;
572  node nodes[NODE_BLOCK_SIZE];
573  } node_block;
574 
575 #define last_node LAST_NODE.LAST_NODE
576 
577  typedef struct arc_for_block_st
578  {
579  char *start; /* the actual start address of this block.
580  May be different from 'this' since 'this'
581  must be at an even address. */
582  arc_forward *current;
583  struct arc_for_block_st *next;
584  arc_forward arcs_for[ARC_BLOCK_SIZE]; /* all arcs must be at even addresses */
585  union
586  {
587  arc_forward dummy;
588  node *LAST_NODE; /* used in graph consruction */
589  } LAST_NODE;
590  } arc_for_block;
591 
592  typedef struct arc_rev_block_st
593  {
594  char *start; /* the actual start address of this block.
595  May be different from 'this' since 'this'
596  must be at an even address. */
597  arc_reverse *current;
598  struct arc_rev_block_st *next;
599  arc_reverse arcs_rev[ARC_BLOCK_SIZE]; /* all arcs must be at even addresses */
600  union
601  {
602  arc_reverse dummy;
603  node *LAST_NODE; /* used in graph consruction */
604  } LAST_NODE;
605  } arc_rev_block;
606 
607  node_block *node_block_first;
608  arc_for_block *arc_for_block_first;
609  arc_rev_block *arc_rev_block_first;
610  DBlock<nodeptr> *nodeptr_block;
611 
612  void (*error_function)(char *); /* this function is called if a error occurs,
613  with a corresponding error message
614  (or exit(1) is called if it's NULL) */
615 
616  flowtype flow; /* total flow */
617 
618 /***********************************************************************/
619 
620  node *queue_first[2], *queue_last[2]; /* list of active nodes */
621  nodeptr *orphan_first, *orphan_last; /* list of pointers to orphans */
622  int TIME; /* monotonically increasing global counter */
623 
624 /***********************************************************************/
625 
626  /* functions for processing active list */
627  void set_active(node *i);
628  node *next_active();
629 
630  void prepare_graph();
631  void maxflow_init();
632  void augment(node *s_start, node *t_start, captype *cap_middle, captype *rev_cap_middle);
633  void process_source_orphan(node *i);
634  void process_sink_orphan(node *i);
635 };
636 /* graph.cpp */
637 
638 
639 //#include <stdio.h>
640 //#include "graph.h"
641 
642 inline Graph::Graph(void (*err_function)(char *))
643 {
644  error_function = err_function;
645  node_block_first = NULL;
646  arc_for_block_first = NULL;
647  arc_rev_block_first = NULL;
648  flow = 0;
649 }
650 
651 inline Graph::~Graph()
652 {
653  while (node_block_first)
654  {
655  node_block *next = node_block_first -> next;
656  delete node_block_first;
657  node_block_first = next;
658  }
659 
660  while (arc_for_block_first)
661  {
662  arc_for_block *next = arc_for_block_first -> next;
663  delete arc_for_block_first -> start;
664  arc_for_block_first = next;
665  }
666 
667  while (arc_rev_block_first)
668  {
669  arc_rev_block *next = arc_rev_block_first -> next;
670  delete arc_rev_block_first -> start;
671  arc_rev_block_first = next;
672  }
673 }
674 
675 inline Graph::node_id Graph::add_node()
676 {
677  node *i;
678 
679  if (!node_block_first || node_block_first->current+1 > &node_block_first->nodes[NODE_BLOCK_SIZE-1])
680  {
681  node_block *next = node_block_first;
682  node_block_first = (node_block *) new node_block;
683  if (!node_block_first) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
684  node_block_first -> current = & ( node_block_first -> nodes[0] );
685  node_block_first -> next = next;
686  }
687 
688  i = node_block_first -> current ++;
689  i -> first_out = (arc_forward *) 0;
690  i -> first_in = (arc_reverse *) 0;
691 
692  i -> tr_cap = 0;
693 
694  return (node_id) i;
695 }
696 
697 inline void Graph::add_edge(node_id from, node_id to, captype cap, captype rev_cap)
698 {
699  arc_forward *a_for;
700  arc_reverse *a_rev;
701 
702  if (!arc_for_block_first || arc_for_block_first->current+1 > &arc_for_block_first->arcs_for[ARC_BLOCK_SIZE])
703  {
704  arc_for_block *next = arc_for_block_first;
705  char *ptr = new char[sizeof(arc_for_block)+1];
706  if (!ptr) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
707  if ((int)ptr & 1) arc_for_block_first = (arc_for_block *) (ptr + 1);
708  else arc_for_block_first = (arc_for_block *) ptr;
709  arc_for_block_first -> start = ptr;
710  arc_for_block_first -> current = & ( arc_for_block_first -> arcs_for[0] );
711  arc_for_block_first -> next = next;
712  }
713 
714  if (!arc_rev_block_first || arc_rev_block_first->current+1 > &arc_rev_block_first->arcs_rev[ARC_BLOCK_SIZE])
715  {
716  arc_rev_block *next = arc_rev_block_first;
717  char *ptr = new char[sizeof(arc_rev_block)+1];
718  if (!ptr) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
719  if ((int)ptr & 1) arc_rev_block_first = (arc_rev_block *) (ptr + 1);
720  else arc_rev_block_first = (arc_rev_block *) ptr;
721  arc_rev_block_first -> start = ptr;
722  arc_rev_block_first -> current = & ( arc_rev_block_first -> arcs_rev[0] );
723  arc_rev_block_first -> next = next;
724  }
725 
726  a_for = arc_for_block_first -> current ++;
727  a_rev = arc_rev_block_first -> current ++;
728 
729  a_rev -> sister = (arc_forward *) from;
730  a_for -> shift = (int) to;
731  a_for -> r_cap = cap;
732  a_for -> r_rev_cap = rev_cap;
733 
734  ((node *)from) -> first_out =
735  (arc_forward *) ((int)(((node *)from) -> first_out) + 1);
736  ((node *)to) -> first_in =
737  (arc_reverse *) ((int)(((node *)to) -> first_in) + 1);
738 }
739 
740 inline void Graph::set_tweights(node_id i, captype cap_source, captype cap_sink)
741 {
742  flow += (cap_source < cap_sink) ? cap_source : cap_sink;
743  ((node*)i) -> tr_cap = cap_source - cap_sink;
744 }
745 
746 inline void Graph::add_tweights(node_id i, captype cap_source, captype cap_sink)
747 {
748  register captype delta = ((node*)i) -> tr_cap;
749  if (delta > 0) cap_source += delta;
750  else cap_sink -= delta;
751  flow += (cap_source < cap_sink) ? cap_source : cap_sink;
752  ((node*)i) -> tr_cap = cap_source - cap_sink;
753 }
754 
755 /*
756  Converts arcs added by 'add_edge()' calls
757  to a forward star graph representation.
758 
759  Linear time algorithm.
760  No or little additional memory is allocated
761  during this process
762  (it may be necessary to allocate additional
763  arc blocks, since arcs corresponding to the
764  same node must be contiguous, i.e. be in one
765  arc block.)
766 */
767 inline void Graph::prepare_graph()
768 {
769  node *i;
770  arc_for_block *ab_for, *ab_for_first;
771  arc_rev_block *ab_rev, *ab_rev_first, *ab_rev_scan;
772  arc_forward *a_for;
773  arc_reverse *a_rev, *a_rev_scan, a_rev_tmp;
774  node_block *nb;
775  bool for_flag = false, rev_flag = false;
776  int k;
777 
778  if (!arc_rev_block_first)
779  {
780  node_id from = add_node(), to = add_node();
781  add_edge(from, to, 1, 0);
782  }
783 
784  /* FIRST STAGE */
785  a_rev_tmp.sister = NULL;
786  for (a_rev=arc_rev_block_first->current; a_rev<&arc_rev_block_first->arcs_rev[ARC_BLOCK_SIZE]; a_rev++)
787  {
788  a_rev -> sister = NULL;
789  }
790 
791  ab_for = ab_for_first = arc_for_block_first;
792  ab_rev = ab_rev_first = ab_rev_scan = arc_rev_block_first;
793  a_for = &ab_for->arcs_for[0];
794  a_rev = a_rev_scan = &ab_rev->arcs_rev[0];
795 
796  for (nb=node_block_first; nb; nb=nb->next)
797  {
798  for (i=&nb->nodes[0]; i<nb->current; i++)
799  {
800  /* outgoing arcs */
801  k = (int) i -> first_out;
802  if (a_for + k > &ab_for->arcs_for[ARC_BLOCK_SIZE])
803  {
804  if (k > ARC_BLOCK_SIZE) { if (error_function) (*error_function)("# of arcs per node exceeds block size!"); exit(1); }
805  if (for_flag) ab_for = NULL;
806  else { ab_for = ab_for -> next; ab_rev_scan = ab_rev_scan -> next; }
807  if (ab_for == NULL)
808  {
809  arc_for_block *next = arc_for_block_first;
810  char *ptr = new char[sizeof(arc_for_block)+1];
811  if (!ptr) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
812  if ((int)ptr & 1) arc_for_block_first = (arc_for_block *) (ptr + 1);
813  else arc_for_block_first = (arc_for_block *) ptr;
814  arc_for_block_first -> start = ptr;
815  arc_for_block_first -> current = & ( arc_for_block_first -> arcs_for[0] );
816  arc_for_block_first -> next = next;
817  ab_for = arc_for_block_first;
818  for_flag = true;
819  }
820  else a_rev_scan = &ab_rev_scan->arcs_rev[0];
821  a_for = &ab_for->arcs_for[0];
822  }
823  if (ab_rev_scan)
824  {
825  a_rev_scan += k;
826  i -> parent = (arc_forward *) a_rev_scan;
827  }
828  else i -> parent = (arc_forward *) &a_rev_tmp;
829  a_for += k;
830  i -> first_out = a_for;
831  ab_for -> last_node = i;
832 
833  /* incoming arcs */
834  k = (int) i -> first_in;
835  if (a_rev + k > &ab_rev->arcs_rev[ARC_BLOCK_SIZE])
836  {
837  if (k > ARC_BLOCK_SIZE) { if (error_function) (*error_function)("# of arcs per node exceeds block size!"); exit(1); }
838  if (rev_flag) ab_rev = NULL;
839  else ab_rev = ab_rev -> next;
840  if (ab_rev == NULL)
841  {
842  arc_rev_block *next = arc_rev_block_first;
843  char *ptr = new char[sizeof(arc_rev_block)+1];
844  if (!ptr) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
845  if ((int)ptr & 1) arc_rev_block_first = (arc_rev_block *) (ptr + 1);
846  else arc_rev_block_first = (arc_rev_block *) ptr;
847  arc_rev_block_first -> start = ptr;
848  arc_rev_block_first -> current = & ( arc_rev_block_first -> arcs_rev[0] );
849  arc_rev_block_first -> next = next;
850  ab_rev = arc_rev_block_first;
851  rev_flag = true;
852  }
853  a_rev = &ab_rev->arcs_rev[0];
854  }
855  a_rev += k;
856  i -> first_in = a_rev;
857  ab_rev -> last_node = i;
858  }
859  /* i is the last node in block */
860  i -> first_out = a_for;
861  i -> first_in = a_rev;
862  }
863 
864  /* SECOND STAGE */
865  for (ab_for=arc_for_block_first; ab_for; ab_for=ab_for->next)
866  {
867  ab_for -> current = ab_for -> last_node -> first_out;
868  }
869 
870  for ( ab_for=ab_for_first, ab_rev=ab_rev_first;
871  ab_for;
872  ab_for=ab_for->next, ab_rev=ab_rev->next )
873  for ( a_for=&ab_for->arcs_for[0], a_rev=&ab_rev->arcs_rev[0];
874  a_for<&ab_for->arcs_for[ARC_BLOCK_SIZE];
875  a_for++, a_rev++ )
876  {
877  arc_forward *af;
878  arc_reverse *ar;
879  node *from;
880  int shift = 0, shift_new;
881  captype r_cap, r_rev_cap, r_cap_new, r_rev_cap_new;
882 
883  if (!(from=(node *)(a_rev->sister))) continue;
884  af = a_for;
885  ar = a_rev;
886 
887  do
888  {
889  ar -> sister = NULL;
890 
891  shift_new = ((char *)(af->shift)) - (char *)from;
892  r_cap_new = af -> r_cap;
893  r_rev_cap_new = af -> r_rev_cap;
894  if (shift)
895  {
896  af -> shift = shift;
897  af -> r_cap = r_cap;
898  af -> r_rev_cap = r_rev_cap;
899  }
900  shift = shift_new;
901  r_cap = r_cap_new;
902  r_rev_cap = r_rev_cap_new;
903 
904  af = -- from -> first_out;
905  if ((arc_reverse *)(from->parent) != &a_rev_tmp)
906  {
907  from -> parent = (arc_forward *)(((arc_reverse *)(from -> parent)) - 1);
908  ar = (arc_reverse *)(from -> parent);
909  }
910  } while (from=(node *)(ar->sister));
911 
912  af -> shift = shift;
913  af -> r_cap = r_cap;
914  af -> r_rev_cap = r_rev_cap;
915  }
916 
917  for (ab_for=arc_for_block_first; ab_for; ab_for=ab_for->next)
918  {
919  i = ab_for -> last_node;
920  a_for = i -> first_out;
921  ab_for -> current -> shift = a_for -> shift;
922  ab_for -> current -> r_cap = a_for -> r_cap;
923  ab_for -> current -> r_rev_cap = a_for -> r_rev_cap;
924  a_for -> shift = (int) (ab_for -> current + 1);
925  i -> first_out = (arc_forward *) (((char *)a_for) - 1);
926  }
927 
928  /* THIRD STAGE */
929  for (ab_rev=arc_rev_block_first; ab_rev; ab_rev=ab_rev->next)
930  {
931  ab_rev -> current = ab_rev -> last_node -> first_in;
932  }
933 
934  for (nb=node_block_first; nb; nb=nb->next)
935  for (i=&nb->nodes[0]; i<nb->current; i++)
936  {
937  arc_forward *a_for_first, *a_for_last;
938 
939  a_for_first = i -> first_out;
940  if (IS_ODD(a_for_first))
941  {
942  a_for_first = (arc_forward *) (((char *)a_for_first) + 1);
943  a_for_last = (arc_forward *) ((a_for_first ++) -> shift);
944  }
945  else a_for_last = (i + 1) -> first_out;
946 
947  for (a_for=a_for_first; a_for<a_for_last; a_for++)
948  {
949  node *to = NEIGHBOR_NODE(i, a_for -> shift);
950  a_rev = -- to -> first_in;
951  a_rev -> sister = a_for;
952  }
953  }
954 
955  for (ab_rev=arc_rev_block_first; ab_rev; ab_rev=ab_rev->next)
956  {
957  i = ab_rev -> last_node;
958  a_rev = i -> first_in;
959  ab_rev -> current -> sister = a_rev -> sister;
960  a_rev -> sister = (arc_forward *) (ab_rev -> current + 1);
961  i -> first_in = (arc_reverse *) (((char *)a_rev) - 1);
962  }
963 }
964 
965 /* maxflow.cpp */
966 
967 //#include <stdio.h>
968 //#include "graph.h"
969 
970 /*
971  special constants for node->parent
972 */
973 #define TERMINAL ( (arc_forward *) 1 ) /* to terminal */
974 #define ORPHAN ( (arc_forward *) 2 ) /* orphan */
975 
976 #define INFINITE_D 1000000000 /* infinite distance to the terminal */
977 
978 /***********************************************************************/
979 
980 /*
981  Functions for processing active list.
982  i->next points to the next node in the list
983  (or to i, if i is the last node in the list).
984  If i->next is NULL iff i is not in the list.
985 
986  There are two queues. Active nodes are added
987  to the end of the second queue and read from
988  the front of the first queue. If the first queue
989  is empty, it is replaced by the second queue
990  (and the second queue becomes empty).
991 */
992 
993 inline void Graph::set_active(node *i)
994 {
995  if (!i->next)
996  {
997  /* it's not in the list yet */
998  if (queue_last[1]) queue_last[1] -> next = i;
999  else queue_first[1] = i;
1000  queue_last[1] = i;
1001  i -> next = i;
1002  }
1003 }
1004 
1005 /*
1006  Returns the next active node.
1007  If it is connected to the sink, it stays in the list,
1008  otherwise it is removed from the list
1009 */
1010 inline Graph::node * Graph::next_active()
1011 {
1012  node *i;
1013 
1014  while ( 1 )
1015  {
1016  if (!(i=queue_first[0]))
1017  {
1018  queue_first[0] = i = queue_first[1];
1019  queue_last[0] = queue_last[1];
1020  queue_first[1] = NULL;
1021  queue_last[1] = NULL;
1022  if (!i) return NULL;
1023  }
1024 
1025  /* remove it from the active list */
1026  if (i->next == i) queue_first[0] = queue_last[0] = NULL;
1027  else queue_first[0] = i -> next;
1028  i -> next = NULL;
1029 
1030  /* a node in the list is active iff it has a parent */
1031  if (i->parent) return i;
1032  }
1033 }
1034 
1035 /***********************************************************************/
1036 
1037 inline void Graph::maxflow_init()
1038 {
1039  node *i;
1040  node_block *nb;
1041 
1042  queue_first[0] = queue_last[0] = NULL;
1043  queue_first[1] = queue_last[1] = NULL;
1044  orphan_first = NULL;
1045 
1046  for (nb=node_block_first; nb; nb=nb->next)
1047  for (i=&nb->nodes[0]; i<nb->current; i++)
1048  {
1049  i -> next = NULL;
1050  i -> TS = 0;
1051  if (i->tr_cap > 0)
1052  {
1053  /* i is connected to the source */
1054  i -> is_sink = 0;
1055  i -> parent = TERMINAL;
1056  set_active(i);
1057  i -> TS = 0;
1058  i -> DIST = 1;
1059  }
1060  else if (i->tr_cap < 0)
1061  {
1062  /* i is connected to the sink */
1063  i -> is_sink = 1;
1064  i -> parent = TERMINAL;
1065  set_active(i);
1066  i -> TS = 0;
1067  i -> DIST = 1;
1068  }
1069  else
1070  {
1071  i -> parent = NULL;
1072  }
1073  }
1074  TIME = 0;
1075 }
1076 
1077 /***********************************************************************/
1078 
1079 inline void Graph::augment(node *s_start, node *t_start, captype *cap_middle, captype *rev_cap_middle)
1080 {
1081  node *i;
1082  arc_forward *a;
1083  captype bottleneck;
1084  nodeptr *np;
1085 
1086 
1087  /* 1. Finding bottleneck capacity */
1088  /* 1a - the source tree */
1089  bottleneck = *cap_middle;
1090  for (i=s_start; ; )
1091  {
1092  a = i -> parent;
1093  if (a == TERMINAL) break;
1094  if (IS_ODD(a))
1095  {
1096  a = MAKE_EVEN(a);
1097  if (bottleneck > a->r_cap) bottleneck = a -> r_cap;
1098  i = NEIGHBOR_NODE_REV(i, a -> shift);
1099  }
1100  else
1101  {
1102  if (bottleneck > a->r_rev_cap) bottleneck = a -> r_rev_cap;
1103  i = NEIGHBOR_NODE(i, a -> shift);
1104  }
1105  }
1106  if (bottleneck > i->tr_cap) bottleneck = i -> tr_cap;
1107  /* 1b - the sink tree */
1108  for (i=t_start; ; )
1109  {
1110  a = i -> parent;
1111  if (a == TERMINAL) break;
1112  if (IS_ODD(a))
1113  {
1114  a = MAKE_EVEN(a);
1115  if (bottleneck > a->r_rev_cap) bottleneck = a -> r_rev_cap;
1116  i = NEIGHBOR_NODE_REV(i, a -> shift);
1117  }
1118  else
1119  {
1120  if (bottleneck > a->r_cap) bottleneck = a -> r_cap;
1121  i = NEIGHBOR_NODE(i, a -> shift);
1122  }
1123  }
1124  if (bottleneck > - i->tr_cap) bottleneck = - i -> tr_cap;
1125 
1126 
1127  /* 2. Augmenting */
1128  /* 2a - the source tree */
1129  *rev_cap_middle += bottleneck;
1130  *cap_middle -= bottleneck;
1131  for (i=s_start; ; )
1132  {
1133  a = i -> parent;
1134  if (a == TERMINAL) break;
1135  if (IS_ODD(a))
1136  {
1137  a = MAKE_EVEN(a);
1138  a -> r_rev_cap += bottleneck;
1139  a -> r_cap -= bottleneck;
1140  if (!a->r_cap)
1141  {
1142  /* add i to the adoption list */
1143  i -> parent = ORPHAN;
1144  np = nodeptr_block -> New();
1145  np -> ptr = i;
1146  np -> next = orphan_first;
1147  orphan_first = np;
1148  }
1149  i = NEIGHBOR_NODE_REV(i, a -> shift);
1150  }
1151  else
1152  {
1153  a -> r_cap += bottleneck;
1154  a -> r_rev_cap -= bottleneck;
1155  if (!a->r_rev_cap)
1156  {
1157  /* add i to the adoption list */
1158  i -> parent = ORPHAN;
1159  np = nodeptr_block -> New();
1160  np -> ptr = i;
1161  np -> next = orphan_first;
1162  orphan_first = np;
1163  }
1164  i = NEIGHBOR_NODE(i, a -> shift);
1165  }
1166  }
1167  i -> tr_cap -= bottleneck;
1168  if (!i->tr_cap)
1169  {
1170  /* add i to the adoption list */
1171  i -> parent = ORPHAN;
1172  np = nodeptr_block -> New();
1173  np -> ptr = i;
1174  np -> next = orphan_first;
1175  orphan_first = np;
1176  }
1177  /* 2b - the sink tree */
1178  for (i=t_start; ; )
1179  {
1180  a = i -> parent;
1181  if (a == TERMINAL) break;
1182  if (IS_ODD(a))
1183  {
1184  a = MAKE_EVEN(a);
1185  a -> r_cap += bottleneck;
1186  a -> r_rev_cap -= bottleneck;
1187  if (!a->r_rev_cap)
1188  {
1189  /* add i to the adoption list */
1190  i -> parent = ORPHAN;
1191  np = nodeptr_block -> New();
1192  np -> ptr = i;
1193  np -> next = orphan_first;
1194  orphan_first = np;
1195  }
1196  i = NEIGHBOR_NODE_REV(i, a -> shift);
1197  }
1198  else
1199  {
1200  a -> r_rev_cap += bottleneck;
1201  a -> r_cap -= bottleneck;
1202  if (!a->r_cap)
1203  {
1204  /* add i to the adoption list */
1205  i -> parent = ORPHAN;
1206  np = nodeptr_block -> New();
1207  np -> ptr = i;
1208  np -> next = orphan_first;
1209  orphan_first = np;
1210  }
1211  i = NEIGHBOR_NODE(i, a -> shift);
1212  }
1213  }
1214  i -> tr_cap += bottleneck;
1215  if (!i->tr_cap)
1216  {
1217  /* add i to the adoption list */
1218  i -> parent = ORPHAN;
1219  np = nodeptr_block -> New();
1220  np -> ptr = i;
1221  np -> next = orphan_first;
1222  orphan_first = np;
1223  }
1224 
1225 
1226  flow += bottleneck;
1227 }
1228 
1229 /***********************************************************************/
1230 
1231 inline void Graph::process_source_orphan(node *i)
1232 {
1233  node *j;
1234  arc_forward *a0_for, *a0_for_first, *a0_for_last;
1235  arc_reverse *a0_rev, *a0_rev_first, *a0_rev_last;
1236  arc_forward *a0_min = NULL, *a;
1237  nodeptr *np;
1238  int d, d_min = INFINITE_D;
1239 
1240  /* trying to find a new parent */
1241  a0_for_first = i -> first_out;
1242  if (IS_ODD(a0_for_first))
1243  {
1244  a0_for_first = (arc_forward *) (((char *)a0_for_first) + 1);
1245  a0_for_last = (arc_forward *) ((a0_for_first ++) -> shift);
1246  }
1247  else a0_for_last = (i + 1) -> first_out;
1248  a0_rev_first = i -> first_in;
1249  if (IS_ODD(a0_rev_first))
1250  {
1251  a0_rev_first = (arc_reverse *) (((char *)a0_rev_first) + 1);
1252  a0_rev_last = (arc_reverse *) ((a0_rev_first ++) -> sister);
1253  }
1254  else a0_rev_last = (i + 1) -> first_in;
1255 
1256 
1257  for (a0_for=a0_for_first; a0_for<a0_for_last; a0_for++)
1258  if (a0_for->r_rev_cap)
1259  {
1260  j = NEIGHBOR_NODE(i, a0_for -> shift);
1261  if (!j->is_sink && (a=j->parent))
1262  {
1263  /* checking the origin of j */
1264  d = 0;
1265  while ( 1 )
1266  {
1267  if (j->TS == TIME)
1268  {
1269  d += j -> DIST;
1270  break;
1271  }
1272  a = j -> parent;
1273  d ++;
1274  if (a==TERMINAL)
1275  {
1276  j -> TS = TIME;
1277  j -> DIST = 1;
1278  break;
1279  }
1280  if (a==ORPHAN) { d = INFINITE_D; break; }
1281  if (IS_ODD(a))
1282  j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
1283  else
1284  j = NEIGHBOR_NODE(j, a -> shift);
1285  }
1286  if (d<INFINITE_D) /* j originates from the source - done */
1287  {
1288  if (d<d_min)
1289  {
1290  a0_min = a0_for;
1291  d_min = d;
1292  }
1293  /* set marks along the path */
1294  for (j=NEIGHBOR_NODE(i, a0_for->shift); j->TS!=TIME; )
1295  {
1296  j -> TS = TIME;
1297  j -> DIST = d --;
1298  a = j->parent;
1299  if (IS_ODD(a))
1300  j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
1301  else
1302  j = NEIGHBOR_NODE(j, a -> shift);
1303  }
1304  }
1305  }
1306  }
1307  for (a0_rev=a0_rev_first; a0_rev<a0_rev_last; a0_rev++)
1308  {
1309  a0_for = a0_rev -> sister;
1310  if (a0_for->r_cap)
1311  {
1312  j = NEIGHBOR_NODE_REV(i, a0_for -> shift);
1313  if (!j->is_sink && (a=j->parent))
1314  {
1315  /* checking the origin of j */
1316  d = 0;
1317  while ( 1 )
1318  {
1319  if (j->TS == TIME)
1320  {
1321  d += j -> DIST;
1322  break;
1323  }
1324  a = j -> parent;
1325  d ++;
1326  if (a==TERMINAL)
1327  {
1328  j -> TS = TIME;
1329  j -> DIST = 1;
1330  break;
1331  }
1332  if (a==ORPHAN) { d = INFINITE_D; break; }
1333  if (IS_ODD(a))
1334  j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
1335  else
1336  j = NEIGHBOR_NODE(j, a -> shift);
1337  }
1338  if (d<INFINITE_D) /* j originates from the source - done */
1339  {
1340  if (d<d_min)
1341  {
1342  a0_min = MAKE_ODD(a0_for);
1343  d_min = d;
1344  }
1345  /* set marks along the path */
1346  for (j=NEIGHBOR_NODE_REV(i,a0_for->shift); j->TS!=TIME; )
1347  {
1348  j -> TS = TIME;
1349  j -> DIST = d --;
1350  a = j->parent;
1351  if (IS_ODD(a))
1352  j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
1353  else
1354  j = NEIGHBOR_NODE(j, a -> shift);
1355  }
1356  }
1357  }
1358  }
1359  }
1360 
1361  if (i->parent = a0_min)
1362  {
1363  i -> TS = TIME;
1364  i -> DIST = d_min + 1;
1365  }
1366  else
1367  {
1368  /* no parent is found */
1369  i -> TS = 0;
1370 
1371  /* process neighbors */
1372  for (a0_for=a0_for_first; a0_for<a0_for_last; a0_for++)
1373  {
1374  j = NEIGHBOR_NODE(i, a0_for -> shift);
1375  if (!j->is_sink && (a=j->parent))
1376  {
1377  if (a0_for->r_rev_cap) set_active(j);
1378  if (a!=TERMINAL && a!=ORPHAN && IS_ODD(a) && NEIGHBOR_NODE_REV(j, MAKE_EVEN(a)->shift)==i)
1379  {
1380  /* add j to the adoption list */
1381  j -> parent = ORPHAN;
1382  np = nodeptr_block -> New();
1383  np -> ptr = j;
1384  if (orphan_last) orphan_last -> next = np;
1385  else orphan_first = np;
1386  orphan_last = np;
1387  np -> next = NULL;
1388  }
1389  }
1390  }
1391  for (a0_rev=a0_rev_first; a0_rev<a0_rev_last; a0_rev++)
1392  {
1393  a0_for = a0_rev -> sister;
1394  j = NEIGHBOR_NODE_REV(i, a0_for -> shift);
1395  if (!j->is_sink && (a=j->parent))
1396  {
1397  if (a0_for->r_cap) set_active(j);
1398  if (a!=TERMINAL && a!=ORPHAN && !IS_ODD(a) && NEIGHBOR_NODE(j, a->shift)==i)
1399  {
1400  /* add j to the adoption list */
1401  j -> parent = ORPHAN;
1402  np = nodeptr_block -> New();
1403  np -> ptr = j;
1404  if (orphan_last) orphan_last -> next = np;
1405  else orphan_first = np;
1406  orphan_last = np;
1407  np -> next = NULL;
1408  }
1409  }
1410  }
1411  }
1412 }
1413 
1414 inline void Graph::process_sink_orphan(node *i)
1415 {
1416  node *j;
1417  arc_forward *a0_for, *a0_for_first, *a0_for_last;
1418  arc_reverse *a0_rev, *a0_rev_first, *a0_rev_last;
1419  arc_forward *a0_min = NULL, *a;
1420  nodeptr *np;
1421  int d, d_min = INFINITE_D;
1422 
1423  /* trying to find a new parent */
1424  a0_for_first = i -> first_out;
1425  if (IS_ODD(a0_for_first))
1426  {
1427  a0_for_first = (arc_forward *) (((char *)a0_for_first) + 1);
1428  a0_for_last = (arc_forward *) ((a0_for_first ++) -> shift);
1429  }
1430  else a0_for_last = (i + 1) -> first_out;
1431  a0_rev_first = i -> first_in;
1432  if (IS_ODD(a0_rev_first))
1433  {
1434  a0_rev_first = (arc_reverse *) (((char *)a0_rev_first) + 1);
1435  a0_rev_last = (arc_reverse *) ((a0_rev_first ++) -> sister);
1436  }
1437  else a0_rev_last = (i + 1) -> first_in;
1438 
1439 
1440  for (a0_for=a0_for_first; a0_for<a0_for_last; a0_for++)
1441  if (a0_for->r_cap)
1442  {
1443  j = NEIGHBOR_NODE(i, a0_for -> shift);
1444  if (j->is_sink && (a=j->parent))
1445  {
1446  /* checking the origin of j */
1447  d = 0;
1448  while ( 1 )
1449  {
1450  if (j->TS == TIME)
1451  {
1452  d += j -> DIST;
1453  break;
1454  }
1455  a = j -> parent;
1456  d ++;
1457  if (a==TERMINAL)
1458  {
1459  j -> TS = TIME;
1460  j -> DIST = 1;
1461  break;
1462  }
1463  if (a==ORPHAN) { d = INFINITE_D; break; }
1464  if (IS_ODD(a))
1465  j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
1466  else
1467  j = NEIGHBOR_NODE(j, a -> shift);
1468  }
1469  if (d<INFINITE_D) /* j originates from the sink - done */
1470  {
1471  if (d<d_min)
1472  {
1473  a0_min = a0_for;
1474  d_min = d;
1475  }
1476  /* set marks along the path */
1477  for (j=NEIGHBOR_NODE(i, a0_for->shift); j->TS!=TIME; )
1478  {
1479  j -> TS = TIME;
1480  j -> DIST = d --;
1481  a = j->parent;
1482  if (IS_ODD(a))
1483  j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
1484  else
1485  j = NEIGHBOR_NODE(j, a -> shift);
1486  }
1487  }
1488  }
1489  }
1490  for (a0_rev=a0_rev_first; a0_rev<a0_rev_last; a0_rev++)
1491  {
1492  a0_for = a0_rev -> sister;
1493  if (a0_for->r_rev_cap)
1494  {
1495  j = NEIGHBOR_NODE_REV(i, a0_for -> shift);
1496  if (j->is_sink && (a=j->parent))
1497  {
1498  /* checking the origin of j */
1499  d = 0;
1500  while ( 1 )
1501  {
1502  if (j->TS == TIME)
1503  {
1504  d += j -> DIST;
1505  break;
1506  }
1507  a = j -> parent;
1508  d ++;
1509  if (a==TERMINAL)
1510  {
1511  j -> TS = TIME;
1512  j -> DIST = 1;
1513  break;
1514  }
1515  if (a==ORPHAN) { d = INFINITE_D; break; }
1516  if (IS_ODD(a))
1517  j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
1518  else
1519  j = NEIGHBOR_NODE(j, a -> shift);
1520  }
1521  if (d<INFINITE_D) /* j originates from the sink - done */
1522  {
1523  if (d<d_min)
1524  {
1525  a0_min = MAKE_ODD(a0_for);
1526  d_min = d;
1527  }
1528  /* set marks along the path */
1529  for (j=NEIGHBOR_NODE_REV(i,a0_for->shift); j->TS!=TIME; )
1530  {
1531  j -> TS = TIME;
1532  j -> DIST = d --;
1533  a = j->parent;
1534  if (IS_ODD(a))
1535  j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
1536  else
1537  j = NEIGHBOR_NODE(j, a -> shift);
1538  }
1539  }
1540  }
1541  }
1542  }
1543 
1544  if (i->parent = a0_min)
1545  {
1546  i -> TS = TIME;
1547  i -> DIST = d_min + 1;
1548  }
1549  else
1550  {
1551  /* no parent is found */
1552  i -> TS = 0;
1553 
1554  /* process neighbors */
1555  for (a0_for=a0_for_first; a0_for<a0_for_last; a0_for++)
1556  {
1557  j = NEIGHBOR_NODE(i, a0_for -> shift);
1558  if (j->is_sink && (a=j->parent))
1559  {
1560  if (a0_for->r_cap) set_active(j);
1561  if (a!=TERMINAL && a!=ORPHAN && IS_ODD(a) && NEIGHBOR_NODE_REV(j, MAKE_EVEN(a)->shift)==i)
1562  {
1563  /* add j to the adoption list */
1564  j -> parent = ORPHAN;
1565  np = nodeptr_block -> New();
1566  np -> ptr = j;
1567  if (orphan_last) orphan_last -> next = np;
1568  else orphan_first = np;
1569  orphan_last = np;
1570  np -> next = NULL;
1571  }
1572  }
1573  }
1574  for (a0_rev=a0_rev_first; a0_rev<a0_rev_last; a0_rev++)
1575  {
1576  a0_for = a0_rev -> sister;
1577  j = NEIGHBOR_NODE_REV(i, a0_for -> shift);
1578  if (j->is_sink && (a=j->parent))
1579  {
1580  if (a0_for->r_rev_cap) set_active(j);
1581  if (a!=TERMINAL && a!=ORPHAN && !IS_ODD(a) && NEIGHBOR_NODE(j, a->shift)==i)
1582  {
1583  /* add j to the adoption list */
1584  j -> parent = ORPHAN;
1585  np = nodeptr_block -> New();
1586  np -> ptr = j;
1587  if (orphan_last) orphan_last -> next = np;
1588  else orphan_first = np;
1589  orphan_last = np;
1590  np -> next = NULL;
1591  }
1592  }
1593  }
1594  }
1595 }
1596 
1597 /***********************************************************************/
1598 
1599 inline Graph::flowtype Graph::maxflow()
1600 {
1601  node *i, *j, *current_node = NULL, *s_start, *t_start;
1602  captype *cap_middle, *rev_cap_middle;
1603  arc_forward *a_for, *a_for_first, *a_for_last;
1604  arc_reverse *a_rev, *a_rev_first, *a_rev_last;
1605  nodeptr *np, *np_next;
1606 
1607  prepare_graph();
1608  maxflow_init();
1609  nodeptr_block = new DBlock<nodeptr>(NODEPTR_BLOCK_SIZE, error_function);
1610 
1611  while ( 1 )
1612  {
1613  if (i=current_node)
1614  {
1615  i -> next = NULL; /* remove active flag */
1616  if (!i->parent) i = NULL;
1617  }
1618  if (!i)
1619  {
1620  if (!(i = next_active())) break;
1621  }
1622 
1623  /* growth */
1624  s_start = NULL;
1625 
1626  a_for_first = i -> first_out;
1627  if (IS_ODD(a_for_first))
1628  {
1629  a_for_first = (arc_forward *) (((char *)a_for_first) + 1);
1630  a_for_last = (arc_forward *) ((a_for_first ++) -> shift);
1631  }
1632  else a_for_last = (i + 1) -> first_out;
1633  a_rev_first = i -> first_in;
1634  if (IS_ODD(a_rev_first))
1635  {
1636  a_rev_first = (arc_reverse *) (((char *)a_rev_first) + 1);
1637  a_rev_last = (arc_reverse *) ((a_rev_first ++) -> sister);
1638  }
1639  else a_rev_last = (i + 1) -> first_in;
1640 
1641  if (!i->is_sink)
1642  {
1643  /* grow source tree */
1644  for (a_for=a_for_first; a_for<a_for_last; a_for++)
1645  if (a_for->r_cap)
1646  {
1647  j = NEIGHBOR_NODE(i, a_for -> shift);
1648  if (!j->parent)
1649  {
1650  j -> is_sink = 0;
1651  j -> parent = MAKE_ODD(a_for);
1652  j -> TS = i -> TS;
1653  j -> DIST = i -> DIST + 1;
1654  set_active(j);
1655  }
1656  else if (j->is_sink)
1657  {
1658  s_start = i;
1659  t_start = j;
1660  cap_middle = & ( a_for -> r_cap );
1661  rev_cap_middle = & ( a_for -> r_rev_cap );
1662  break;
1663  }
1664  else if (j->TS <= i->TS &&
1665  j->DIST > i->DIST)
1666  {
1667  /* heuristic - trying to make the distance from j to the source shorter */
1668  j -> parent = MAKE_ODD(a_for);
1669  j -> TS = i -> TS;
1670  j -> DIST = i -> DIST + 1;
1671  }
1672  }
1673  if (!s_start)
1674  for (a_rev=a_rev_first; a_rev<a_rev_last; a_rev++)
1675  {
1676  a_for = a_rev -> sister;
1677  if (a_for->r_rev_cap)
1678  {
1679  j = NEIGHBOR_NODE_REV(i, a_for -> shift);
1680  if (!j->parent)
1681  {
1682  j -> is_sink = 0;
1683  j -> parent = a_for;
1684  j -> TS = i -> TS;
1685  j -> DIST = i -> DIST + 1;
1686  set_active(j);
1687  }
1688  else if (j->is_sink)
1689  {
1690  s_start = i;
1691  t_start = j;
1692  cap_middle = & ( a_for -> r_rev_cap );
1693  rev_cap_middle = & ( a_for -> r_cap );
1694  break;
1695  }
1696  else if (j->TS <= i->TS &&
1697  j->DIST > i->DIST)
1698  {
1699  /* heuristic - trying to make the distance from j to the source shorter */
1700  j -> parent = a_for;
1701  j -> TS = i -> TS;
1702  j -> DIST = i -> DIST + 1;
1703  }
1704  }
1705  }
1706  }
1707  else
1708  {
1709  /* grow sink tree */
1710  for (a_for=a_for_first; a_for<a_for_last; a_for++)
1711  if (a_for->r_rev_cap)
1712  {
1713  j = NEIGHBOR_NODE(i, a_for -> shift);
1714  if (!j->parent)
1715  {
1716  j -> is_sink = 1;
1717  j -> parent = MAKE_ODD(a_for);
1718  j -> TS = i -> TS;
1719  j -> DIST = i -> DIST + 1;
1720  set_active(j);
1721  }
1722  else if (!j->is_sink)
1723  {
1724  s_start = j;
1725  t_start = i;
1726  cap_middle = & ( a_for -> r_rev_cap );
1727  rev_cap_middle = & ( a_for -> r_cap );
1728  break;
1729  }
1730  else if (j->TS <= i->TS &&
1731  j->DIST > i->DIST)
1732  {
1733  /* heuristic - trying to make the distance from j to the sink shorter */
1734  j -> parent = MAKE_ODD(a_for);
1735  j -> TS = i -> TS;
1736  j -> DIST = i -> DIST + 1;
1737  }
1738  }
1739  for (a_rev=a_rev_first; a_rev<a_rev_last; a_rev++)
1740  {
1741  a_for = a_rev -> sister;
1742  if (a_for->r_cap)
1743  {
1744  j = NEIGHBOR_NODE_REV(i, a_for -> shift);
1745  if (!j->parent)
1746  {
1747  j -> is_sink = 1;
1748  j -> parent = a_for;
1749  j -> TS = i -> TS;
1750  j -> DIST = i -> DIST + 1;
1751  set_active(j);
1752  }
1753  else if (!j->is_sink)
1754  {
1755  s_start = j;
1756  t_start = i;
1757  cap_middle = & ( a_for -> r_cap );
1758  rev_cap_middle = & ( a_for -> r_rev_cap );
1759  break;
1760  }
1761  else if (j->TS <= i->TS &&
1762  j->DIST > i->DIST)
1763  {
1764  /* heuristic - trying to make the distance from j to the sink shorter */
1765  j -> parent = a_for;
1766  j -> TS = i -> TS;
1767  j -> DIST = i -> DIST + 1;
1768  }
1769  }
1770  }
1771  }
1772 
1773  TIME ++;
1774 
1775  if (s_start)
1776  {
1777  i -> next = i; /* set active flag */
1778  current_node = i;
1779 
1780  /* augmentation */
1781  augment(s_start, t_start, cap_middle, rev_cap_middle);
1782  /* augmentation end */
1783 
1784  /* adoption */
1785  while (np=orphan_first)
1786  {
1787  np_next = np -> next;
1788  np -> next = NULL;
1789 
1790  while (np=orphan_first)
1791  {
1792  orphan_first = np -> next;
1793  i = np -> ptr;
1794  nodeptr_block -> Delete(np);
1795  if (!orphan_first) orphan_last = NULL;
1796  if (i->is_sink) process_sink_orphan(i);
1797  else process_source_orphan(i);
1798  }
1799 
1800  orphan_first = np_next;
1801  }
1802  /* adoption end */
1803  }
1804  else current_node = NULL;
1805  }
1806 
1807  delete nodeptr_block;
1808 
1809  return flow;
1810 }
1811 
1812 /***********************************************************************/
1813 
1814 inline Graph::termtype Graph::what_segment(node_id i)
1815 {
1816  if (((node*)i)->parent && !((node*)i)->is_sink) return SOURCE;
1817  return SINK;
1818 }
1819 
1820 #endif
Generated on Sun Jan 13 2013 04:06:21 for Surface Mesh Segmentation by   doxygen 1.8.1.2