1
// Copyright (c) 2001 Yuri Boykov
2
// All rights reserved.
4
// This file is part of CGAL (www.cgal.org).
5
// You can redistribute it and/or modify it under the terms of the GNU
6
// General Public License as published by the Free Software Foundation,
7
// either version 3 of the License, or (at your option) any later version.
9
// Licensees holding a valid commercial license may use this file in
10
// accordance with the commercial license agreement provided with the software.
12
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
13
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
15
###################################################################
17
# MAXFLOW - software for computing mincut/maxflow in a graph #
19
# http://www.cs.ucl.ac.uk/staff/V.Kolmogorov/software.html #
21
# Yuri Boykov (yuri@csd.uwo.ca) #
22
# Vladimir Kolmogorov (v.kolmogorov@cs.ucl.ac.uk) #
25
###################################################################
29
This software library implements the maxflow algorithm
32
An Experimental Comparison of Min-Cut/Max-Flow Algorithms
33
for Energy Minimization in Vision.
34
Yuri Boykov and Vladimir Kolmogorov.
35
In IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI),
38
This algorithm was developed by Yuri Boykov and Vladimir Kolmogorov
39
at Siemens Corporate Research. To make it available for public use,
40
it was later reimplemented by Vladimir Kolmogorov based on open publications.
42
If you use this software for research purposes, you should cite
43
the aforementioned paper in any resulting publication.
45
Tested under windows, Visual C++ 6.0 compiler and unix (SunOS 5.8
46
and RedHat Linux 7.0, GNU c++ compiler).
48
##################################################################
52
Copyright UCL Business PLC
54
This program is available under dual licence:
56
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.
57
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.
58
GNU General Public License can be found at http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
60
2) Proprietary Licence from UCL Business PLC.
61
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
63
##################################################################
65
3. Graph representation.
67
There are two versions of the algorithm using different
68
graph representations (adjacency list and forward star).
69
The former one uses more than twice as much memory as the
70
latter one but is 10-20% faster.
72
Memory allocation (assuming that all capacities are 'short' - 2 bytes):
75
------------------------------------------
76
Adjacency list | *24 bytes | *14 bytes
77
Forward star | *28 bytes | 6 bytes
79
(* means that often it should be rounded up to be a multiple of 4
80
- some compilers (e.g. Visual C++) seem to round up elements
81
of arrays unless the are structures containing only char[].)
83
Note that arcs are always added in pairs - in forward and reverse directions.
84
Arcs between nodes and terminals (the source and the sink) are
85
not stored as arcs, but rather as a part of nodes.
87
The assumption for the forward star representation is that
88
the maximum number of arcs per node (except the source
89
and the sink) is much less than ARC_BLOCK_SIZE (1024 by default).
91
Both versions have the same interface.
93
##################################################################
97
This section shows how to use the library to compute
98
a minimum cut on the following graph:
112
///////////////////////////////////////////////////
119
Graph::node_id nodes[2];
120
Graph *g = new Graph();
122
nodes[0] = g -> add_node();
123
nodes[1] = g -> add_node();
124
g -> set_tweights(nodes[0], 1, 5);
125
g -> set_tweights(nodes[1], 2, 6);
126
g -> add_edge(nodes[0], nodes[1], 3, 4);
128
Graph::flowtype flow = g -> maxflow();
130
printf("Flow = %d\n", flow);
131
printf("Minimum cut:\n");
132
if (g->what_segment(nodes[0]) == Graph::SOURCE)
133
printf("node0 is in the SOURCE set\n");
135
printf("node0 is in the SINK set\n");
136
if (g->what_segment(nodes[1]) == Graph::SOURCE)
137
printf("node1 is in the SOURCE set\n");
139
printf("node1 is in the SINK set\n");
144
///////////////////////////////////////////////////
149
Template classes Block and DBlock
150
Implement adding and deleting items of the same type in blocks.
152
If there there are many items then using Block or DBlock
153
is more efficient than using 'new' and 'delete' both in terms
154
of memory and time since
155
(1) On some systems there is some minimum amount of memory
156
that 'new' can allocate (e.g., 64), so if items are
157
small that a lot of memory is wasted.
158
(2) 'new' and 'delete' are designed for items of varying size.
159
If all items has the same size, then an algorithm for
160
adding and deleting can be made more efficient.
161
(3) All Block and DBlock functions are inline, so there are
162
no extra function calls.
164
Differences between Block and DBlock:
165
(1) DBlock allows both adding and deleting items,
166
whereas Block allows only adding items.
167
(2) Block has an additional operation of scanning
168
items added so far (in the order in which they were added).
169
(3) Block allows to allocate several consecutive
170
items at a time, whereas DBlock can add only a single item.
172
Note that no constructors or destructors are called for items.
174
Example usage for items of type 'MyType':
176
///////////////////////////////////////////////////
178
#define BLOCK_SIZE 1024
179
typedef struct { int a, b; } MyType;
180
MyType *ptr, *array[10000];
184
Block<MyType> *block = new Block<MyType>(BLOCK_SIZE);
187
for (int i=0; i<sizeof(array); i++)
189
ptr = block -> New();
190
ptr -> a = ptr -> b = rand();
194
for (ptr=block->ScanFirst(); ptr; ptr=block->ScanNext())
196
printf("%d %d\n", ptr->a, ptr->b);
203
DBlock<MyType> *dblock = new DBlock<MyType>(BLOCK_SIZE);
206
for (int i=0; i<sizeof(array); i++)
208
array[i] = dblock -> New();
212
for (int i=0; i<sizeof(array); i+=2)
214
dblock -> Delete(array[i]);
218
for (int i=0; i<sizeof(array); i++)
220
array[i] = dblock -> New();
225
///////////////////////////////////////////////////
227
Note that DBlock deletes items by marking them as
228
empty (i.e., by adding them to the list of free items),
229
so that this memory could be used for subsequently
230
added items. Thus, at each moment the memory allocated
231
is determined by the maximum number of items allocated
232
simultaneously at earlier moments. All memory is
233
deallocated only when the destructor is called.
241
/***********************************************************************/
242
/***********************************************************************/
243
/***********************************************************************/
245
template <class Type> class Block
248
/* Constructor. Arguments are the block size and
249
(optionally) the pointer to the function which
250
will be called if allocation failed; the message
251
passed to this function is "Not enough memory!" */
252
Block(int size, void (*err_function)(const char *) = NULL) {
255
error_function = err_function;
258
/* Destructor. Deallocates all items added so far */
261
block *next = first -> next;
262
delete[] ((char*)first);
267
/* Allocates 'num' consecutive items; returns pointer
268
to the first item. 'num' cannot be greater than the
269
block size since items must fit in one block */
270
Type *New(int num = 1) {
273
if (!last || last->current + num > last->last) {
274
if (last && last->next) last = last -> next;
276
block *next = (block *) new char [sizeof(block) + (block_size-1)*sizeof(Type)];
278
if (error_function) (*error_function)("Not enough memory!");
281
if (last) last -> next = next;
284
last -> current = & ( last -> data[0] );
285
last -> last = last -> current + block_size;
291
last -> current += num;
295
/* Returns the first item (or NULL, if no items were added) */
297
for (scan_current_block=first; scan_current_block;
298
scan_current_block = scan_current_block->next) {
299
scan_current_data = & ( scan_current_block -> data[0] );
300
if (scan_current_data < scan_current_block -> current) return scan_current_data
306
/* Returns the next item (or NULL, if all items have been read)
307
Can be called only if previous ScanFirst() or ScanNext()
308
call returned not NULL. */
310
while (scan_current_data >= scan_current_block -> current) {
311
scan_current_block = scan_current_block -> next;
312
if (!scan_current_block) return NULL;
313
scan_current_data = & ( scan_current_block -> data[0] );
315
return scan_current_data ++;
318
/* Marks all elements as empty */
322
for (b=first; ; b=b->next) {
323
b -> current = & ( b -> data[0] );
324
if (b == last) break;
329
/***********************************************************************/
333
typedef struct block_st {
334
Type *current, *last;
335
struct block_st *next;
343
block *scan_current_block;
344
Type *scan_current_data;
346
void (*error_function)(const char *);
349
/***********************************************************************/
350
/***********************************************************************/
351
/***********************************************************************/
353
template <class Type> class DBlock
356
/* Constructor. Arguments are the block size and
357
(optionally) the pointer to the function which
358
will be called if allocation failed; the message
359
passed to this function is "Not enough memory!" */
360
DBlock(int size, void (*err_function)(const char *) = NULL) {
364
error_function = err_function;
367
/* Destructor. Deallocates all items added so far */
370
block *next = first -> next;
371
delete[] ((char*)first);
376
/* Allocates one item */
382
first = (block *) new char [sizeof(block) + (block_size-1)*sizeof(block_item)];
384
if (error_function) (*error_function)("Not enough memory!");
387
first_free = & (first -> data[0] );
388
for (item=first_free; item<first_free+block_size-1; item++)
389
item -> next_free = item + 1;
390
item -> next_free = NULL;
391
first -> next = next;
395
first_free = item -> next_free;
396
return (Type *) item;
399
/* Deletes an item allocated previously */
400
void Delete(Type *t) {
401
((block_item *) t) -> next_free = first_free;
402
first_free = (block_item *) t;
405
/***********************************************************************/
409
typedef union block_item_st {
411
block_item_st *next_free;
414
typedef struct block_st {
415
struct block_st *next;
421
block_item *first_free;
423
void (*error_function)(const char *);
433
This software library implements the maxflow algorithm
436
An Experimental Comparison of Min-Cut/Max-Flow Algorithms
437
for Energy Minimization in Vision.
438
Yuri Boykov and Vladimir Kolmogorov.
439
In IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI),
442
This algorithm was developed by Yuri Boykov and Vladimir Kolmogorov
443
at Siemens Corporate Research. To make it available for public use,
444
it was later reimplemented by Vladimir Kolmogorov based on open publications.
446
If you use this software for research purposes, you should cite
447
the aforementioned paper in any resulting publication.
449
----------------------------------------------------------------
451
For description, license, example usage, discussion of graph representation and memory usage see README.TXT.
460
Nodes, arcs and pointers to nodes are
461
added in blocks for memory and time efficiency.
462
Below are numbers of items in blocks
464
#define NODE_BLOCK_SIZE 512
465
#define ARC_BLOCK_SIZE 1024
466
#define NODEPTR_BLOCK_SIZE 128
468
template <std::size_t size>
471
template<> struct Int_to_ptr<sizeof(int)> {
474
#if INT_MAX != LONG_MAX
475
template<> struct Int_to_ptr<sizeof(long)> {
479
template<> struct Int_to_ptr<sizeof(long long)> {
480
typedef long long type;
491
} termtype; /* terminals */
493
/* Type of edge weights.
494
Can be changed to char, int, float, double, ... */
495
typedef double captype;
496
/* Type of total flow */
497
typedef double flowtype;
499
typedef void * node_id;
501
/* interface functions */
503
/* Constructor. Optional argument is the pointer to the
504
function which will be called if an error occurs;
505
an error message is passed to this function. If this
506
argument is omitted, exit(1) will be called. */
507
Graph(void (*err_function)(const char *) = NULL);
512
/* Adds a node to the graph */
515
/* Adds a bidirectional edge between 'from' and 'to'
516
with the weights 'cap' and 'rev_cap' */
517
void add_edge(node_id from, node_id to, captype cap, captype rev_cap);
519
/* Sets the weights of the edges 'SOURCE->i' and 'i->SINK'
520
Can be called at most once for each node before any call to 'add_tweights'.
521
Weights can be negative */
522
void set_tweights(node_id i, captype cap_source, captype cap_sink);
524
/* Adds new edges 'SOURCE->i' and 'i->SINK' with corresponding weights
525
Can be called multiple times for each node.
526
Weights can be negative */
527
void add_tweights(node_id i, captype cap_source, captype cap_sink);
529
/* After the maxflow is computed, this function returns to which
530
segment the node 'i' belongs (Graph::SOURCE or Graph::SINK) */
531
termtype what_segment(node_id i);
533
/* Computes the maxflow. Can be called only once. */
536
/***********************************************************************/
537
/***********************************************************************/
538
/***********************************************************************/
541
/* internal variables and functions */
543
struct arc_forward_st;
544
struct arc_reverse_st;
546
typedef Int_to_ptr< sizeof(void*) >::type INTEGER;
547
#define IS_ODD(a) ((INTEGER)(a) & 1)
548
#define MAKE_ODD(a) ((arc_forward *) ((INTEGER)(a) | 1))
549
#define MAKE_EVEN(a) ((arc_forward *) ((INTEGER)(a) & (~1)))
550
#define MAKE_ODD_REV(a) ((arc_reverse *) ((INTEGER)(a) | 1))
551
#define MAKE_EVEN_REV(a) ((arc_reverse *) ((INTEGER)(a) & (~1)))
552
#define POINTER_TO_INTEGER(ptr) ((INTEGER) ptr)
557
typedef struct node_st {
559
Usually i->first_out is the first outgoing
560
arc, and (i+1)->first_out-1 is the last outgoing arc.
561
However, it is not always possible, since
562
arcs are allocated in blocks, so arcs corresponding
563
to two consecutive nodes may be in different blocks.
565
If outgoing arcs for i are last in the arc block,
566
then a different mechanism is used. i->first_out
567
is odd in this case; the first outgoing arc
568
is (a+1), and the last outgoing arc is
569
((arc_forward *)(a->shift))-1, where
570
a = (arc_forward *) (((char *)(i->first_out)) + 1);
572
Similar mechanism is used for incoming arcs.
574
arc_forward_st *first_out; /* first outcoming arc */
575
arc_reverse_st *first_in; /* first incoming arc */
577
arc_forward_st *parent; /* describes node's parent
578
if IS_ODD(parent) then MAKE_EVEN(parent) points to 'arc_reverse',
579
otherwise parent points to 'arc_forward' */
581
node_st *next; /* pointer to the next active node
582
(or to itself if it is the last node in the list) */
584
int TS; /* timestamp showing when DIST was computed */
585
int DIST; /* distance to the terminal */
586
short is_sink; /* flag showing whether the node is in the source or in the sink tree */
588
captype tr_cap; /* if tr_cap > 0 then tr_cap is residual capacity of the arc SOURCE->node
589
otherwise -tr_cap is residual capacity of the arc node->SINK */
593
#define NEIGHBOR_NODE(i, shift) ((node *) ((char *)(i) + (shift)))
594
#define NEIGHBOR_NODE_REV(i, shift) ((node *) ((char *)(i) - (shift)))
595
typedef struct arc_forward_st {
596
INTEGER shift; /* node_to = NEIGHBOR_NODE(node_from, shift) */
597
captype r_cap; /* residual capacity */
598
captype r_rev_cap; /* residual capacity of the reverse arc*/
601
typedef struct arc_reverse_st {
602
arc_forward *sister; /* reverse arc */
605
/* 'pointer to node' structure */
606
typedef struct nodeptr_st {
611
typedef struct node_block_st {
613
struct node_block_st *next;
614
node nodes[NODE_BLOCK_SIZE];
617
#define last_node LAST_NODE.LAST_NODE
619
typedef struct arc_for_block_st {
620
char *start; /* the actual start address of this block.
621
May be different from 'this' since 'this'
622
must be at an even address. */
623
arc_forward *current;
624
struct arc_for_block_st *next;
626
arcs_for[ARC_BLOCK_SIZE]; /* all arcs must be at even addresses */
629
node *LAST_NODE; /* used in graph consruction */
633
typedef struct arc_rev_block_st {
634
char *start; /* the actual start address of this block.
635
May be different from 'this' since 'this'
636
must be at an even address. */
637
arc_reverse *current;
638
struct arc_rev_block_st *next;
640
arcs_rev[ARC_BLOCK_SIZE]; /* all arcs must be at even addresses */
643
node *LAST_NODE; /* used in graph consruction */
647
node_block *node_block_first;
648
arc_for_block *arc_for_block_first;
649
arc_rev_block *arc_rev_block_first;
650
DBlock<nodeptr> *nodeptr_block;
652
void (*error_function)(const char
653
*); /* this function is called if a error occurs,
654
with a corresponding error message
655
(or exit(1) is called if it's NULL) */
657
flowtype flow; /* total flow */
659
/***********************************************************************/
661
node *queue_first[2], *queue_last[2]; /* list of active nodes */
662
nodeptr *orphan_first, *orphan_last; /* list of pointers to orphans */
663
int TIME; /* monotonically increasing global counter */
665
/***********************************************************************/
667
/* functions for processing active list */
668
void set_active(node *i);
671
void prepare_graph();
673
void augment(node *s_start, node *t_start, captype *cap_middle,
674
captype *rev_cap_middle);
675
void process_source_orphan(node *i);
676
void process_sink_orphan(node *i);
684
inline Graph::Graph(void (*err_function)(const char *))
686
error_function = err_function;
687
node_block_first = NULL;
688
arc_for_block_first = NULL;
689
arc_rev_block_first = NULL;
693
inline Graph::~Graph()
695
while (node_block_first) {
696
node_block *next = node_block_first -> next;
697
delete node_block_first;
698
node_block_first = next;
701
while (arc_for_block_first) {
702
arc_for_block *next = arc_for_block_first -> next;
703
delete[] arc_for_block_first -> start;
704
arc_for_block_first = next;
707
while (arc_rev_block_first) {
708
arc_rev_block *next = arc_rev_block_first -> next;
709
delete[] arc_rev_block_first -> start;
710
arc_rev_block_first = next;
714
inline Graph::node_id Graph::add_node()
718
if (!node_block_first
719
|| node_block_first->current+1 > &node_block_first->nodes[NODE_BLOCK_SIZE-1]) {
720
node_block *next = node_block_first;
721
node_block_first = (node_block *) new node_block;
722
if (!node_block_first) {
723
if (error_function) (*error_function)("Not enough memory!");
726
node_block_first -> current = & ( node_block_first -> nodes[0] );
727
node_block_first -> next = next;
730
i = node_block_first -> current ++;
731
i -> first_out = (arc_forward *) 0;
732
i -> first_in = (arc_reverse *) 0;
739
inline void Graph::add_edge(node_id from, node_id to, captype cap,
745
if (!arc_for_block_first
746
|| arc_for_block_first->current+1 >
747
&arc_for_block_first->arcs_for[ARC_BLOCK_SIZE]) {
748
arc_for_block *next = arc_for_block_first;
749
char *ptr = new char[sizeof(arc_for_block)+1];
751
if (error_function) (*error_function)("Not enough memory!");
754
if (IS_ODD(ptr)) arc_for_block_first = (arc_for_block *) (ptr + 1);
755
else arc_for_block_first = (arc_for_block *) ptr;
756
arc_for_block_first -> start = ptr;
757
arc_for_block_first -> current = & ( arc_for_block_first -> arcs_for[0] );
758
arc_for_block_first -> next = next;
761
if (!arc_rev_block_first
762
|| arc_rev_block_first->current+1 >
763
&arc_rev_block_first->arcs_rev[ARC_BLOCK_SIZE]) {
764
arc_rev_block *next = arc_rev_block_first;
765
char *ptr = new char[sizeof(arc_rev_block)+1];
767
if (error_function) (*error_function)("Not enough memory!");
770
if (IS_ODD(ptr)) arc_rev_block_first = (arc_rev_block *) (ptr + 1);
771
else arc_rev_block_first = (arc_rev_block *) ptr;
772
arc_rev_block_first -> start = ptr;
773
arc_rev_block_first -> current = & ( arc_rev_block_first -> arcs_rev[0] );
774
arc_rev_block_first -> next = next;
777
a_for = arc_for_block_first -> current ++;
778
a_rev = arc_rev_block_first -> current ++;
780
a_rev -> sister = (arc_forward *) from;
781
a_for -> shift = POINTER_TO_INTEGER(to);
782
a_for -> r_cap = cap;
783
a_for -> r_rev_cap = rev_cap;
785
((node *)from) -> first_out =
786
(arc_forward *) (POINTER_TO_INTEGER(((node *)from) -> first_out) + 1);
787
((node *)to) -> first_in =
788
(arc_reverse *) (POINTER_TO_INTEGER(((node *)to) -> first_in) + 1);
791
inline void Graph::set_tweights(node_id i, captype cap_source, captype cap_sink)
793
flow += (cap_source < cap_sink) ? cap_source : cap_sink;
794
((node*)i) -> tr_cap = cap_source - cap_sink;
797
inline void Graph::add_tweights(node_id i, captype cap_source, captype cap_sink)
799
register captype delta = ((node*)i) -> tr_cap;
800
if (delta > 0) cap_source += delta;
801
else cap_sink -= delta;
802
flow += (cap_source < cap_sink) ? cap_source : cap_sink;
803
((node*)i) -> tr_cap = cap_source - cap_sink;
807
Converts arcs added by 'add_edge()' calls
808
to a forward star graph representation.
810
Linear time algorithm.
811
No or little additional memory is allocated
813
(it may be necessary to allocate additional
814
arc blocks, since arcs corresponding to the
815
same node must be contiguous, i.e. be in one
818
inline void Graph::prepare_graph()
821
arc_for_block *ab_for, *ab_for_first;
822
arc_rev_block *ab_rev, *ab_rev_first, *ab_rev_scan;
824
arc_reverse *a_rev, *a_rev_scan, *a_rev_tmp=new arc_reverse;
826
bool for_flag = false, rev_flag = false;
829
if (!arc_rev_block_first) {
830
node_id from = add_node(), to = add_node();
831
add_edge(from, to, 1, 0);
835
a_rev_tmp->sister = NULL;
836
for (a_rev=arc_rev_block_first->current;
837
a_rev<&arc_rev_block_first->arcs_rev[ARC_BLOCK_SIZE]; a_rev++) {
838
a_rev -> sister = NULL;
841
ab_for = ab_for_first = arc_for_block_first;
842
ab_rev = ab_rev_first = ab_rev_scan = arc_rev_block_first;
843
a_for = &ab_for->arcs_for[0];
844
a_rev = a_rev_scan = &ab_rev->arcs_rev[0];
846
for (nb=node_block_first; nb; nb=nb->next) {
847
for (i=&nb->nodes[0]; i<nb->current; i++) {
849
k = POINTER_TO_INTEGER(i -> first_out);
850
if (a_for + k > &ab_for->arcs_for[ARC_BLOCK_SIZE]) {
851
if (k > ARC_BLOCK_SIZE) {
852
if (error_function) (*error_function)("# of arcs per node exceeds block size!");
855
if (for_flag) ab_for = NULL;
857
ab_for = ab_for -> next;
858
ab_rev_scan = ab_rev_scan -> next;
860
if (ab_for == NULL) {
861
arc_for_block *next = arc_for_block_first;
862
char *ptr = new char[sizeof(arc_for_block)+1];
864
if (error_function) (*error_function)("Not enough memory!");
867
if (IS_ODD(ptr)) arc_for_block_first = (arc_for_block *) (ptr + 1);
868
else arc_for_block_first = (arc_for_block *) ptr;
869
arc_for_block_first -> start = ptr;
870
arc_for_block_first -> current = & ( arc_for_block_first -> arcs_for[0] );
871
arc_for_block_first -> next = next;
872
ab_for = arc_for_block_first;
874
} else a_rev_scan = &ab_rev_scan->arcs_rev[0];
875
a_for = &ab_for->arcs_for[0];
879
i -> parent = (arc_forward *) a_rev_scan;
880
} else i -> parent = (arc_forward *) a_rev_tmp;
882
i -> first_out = a_for;
883
ab_for -> last_node = i;
886
k = POINTER_TO_INTEGER(i -> first_in);
887
if (a_rev + k > &ab_rev->arcs_rev[ARC_BLOCK_SIZE]) {
888
if (k > ARC_BLOCK_SIZE) {
889
if (error_function) (*error_function)("# of arcs per node exceeds block size!");
892
if (rev_flag) ab_rev = NULL;
893
else ab_rev = ab_rev -> next;
894
if (ab_rev == NULL) {
895
arc_rev_block *next = arc_rev_block_first;
896
char *ptr = new char[sizeof(arc_rev_block)+1];
898
if (error_function) (*error_function)("Not enough memory!");
901
if (IS_ODD(ptr)) arc_rev_block_first = (arc_rev_block *) (ptr + 1);
902
else arc_rev_block_first = (arc_rev_block *) ptr;
903
arc_rev_block_first -> start = ptr;
904
arc_rev_block_first -> current = & ( arc_rev_block_first -> arcs_rev[0] );
905
arc_rev_block_first -> next = next;
906
ab_rev = arc_rev_block_first;
909
a_rev = &ab_rev->arcs_rev[0];
912
i -> first_in = a_rev;
913
ab_rev -> last_node = i;
915
/* i is the last node in block */
916
i -> first_out = a_for;
917
i -> first_in = a_rev;
921
for (ab_for=arc_for_block_first; ab_for; ab_for=ab_for->next) {
922
ab_for -> current = ab_for -> last_node -> first_out;
925
for ( ab_for=ab_for_first, ab_rev=ab_rev_first;
927
ab_for=ab_for->next, ab_rev=ab_rev->next )
928
for ( a_for=&ab_for->arcs_for[0], a_rev=&ab_rev->arcs_rev[0];
929
a_for<&ab_for->arcs_for[ARC_BLOCK_SIZE];
934
INTEGER shift = 0, shift_new;
935
captype r_cap=0, r_rev_cap=0, r_cap_new, r_rev_cap_new;
937
if (!(from=(node *)(a_rev->sister))) continue;
944
shift_new = ((char *)(af->shift)) - (char *)from;
945
r_cap_new = af -> r_cap;
946
r_rev_cap_new = af -> r_rev_cap;
950
af -> r_rev_cap = r_rev_cap;
954
r_rev_cap = r_rev_cap_new;
956
af = -- from -> first_out;
957
if ((arc_reverse *)(from->parent) != a_rev_tmp) {
958
from -> parent = (arc_forward *)(((arc_reverse *)(from -> parent)) - 1);
959
ar = (arc_reverse *)(from -> parent);
961
} while ( (from=(node *)(ar->sister)) );
965
af -> r_rev_cap = r_rev_cap;
968
for (ab_for=arc_for_block_first; ab_for; ab_for=ab_for->next) {
969
i = ab_for -> last_node;
970
a_for = i -> first_out;
971
ab_for -> current -> shift = a_for -> shift;
972
ab_for -> current -> r_cap = a_for -> r_cap;
973
ab_for -> current -> r_rev_cap = a_for -> r_rev_cap;
974
a_for -> shift = POINTER_TO_INTEGER(ab_for -> current + 1);
975
i -> first_out = (arc_forward *) (((char *)a_for) - 1);
979
for (ab_rev=arc_rev_block_first; ab_rev; ab_rev=ab_rev->next) {
980
ab_rev -> current = ab_rev -> last_node -> first_in;
983
for (nb=node_block_first; nb; nb=nb->next)
984
for (i=&nb->nodes[0]; i<nb->current; i++) {
985
arc_forward *a_for_first, *a_for_last;
987
a_for_first = i -> first_out;
988
if (IS_ODD(a_for_first)) {
989
a_for_first = (arc_forward *) (((char *)a_for_first) + 1);
990
a_for_last = (arc_forward *) ((a_for_first ++) -> shift);
991
} else a_for_last = (i + 1) -> first_out;
993
for (a_for=a_for_first; a_for<a_for_last; a_for++) {
994
node *to = NEIGHBOR_NODE(i, a_for -> shift);
995
a_rev = -- to -> first_in;
996
a_rev -> sister = a_for;
1000
for (ab_rev=arc_rev_block_first; ab_rev; ab_rev=ab_rev->next) {
1001
i = ab_rev -> last_node;
1002
a_rev = i -> first_in;
1003
ab_rev -> current -> sister = a_rev -> sister;
1004
a_rev -> sister = (arc_forward *) (ab_rev -> current + 1);
1005
i -> first_in = (arc_reverse *) (((char *)a_rev) - 1);
1012
//#include <stdio.h>
1013
//#include "graph.h"
1016
special constants for node->parent
1018
#define TERMINAL ( (arc_forward *) 1 ) /* to terminal */
1019
#define ORPHAN ( (arc_forward *) 2 ) /* orphan */
1021
#define INFINITE_D 1000000000 /* infinite distance to the terminal */
1023
/***********************************************************************/
1026
Functions for processing active list.
1027
i->next points to the next node in the list
1028
(or to i, if i is the last node in the list).
1029
If i->next is NULL iff i is not in the list.
1031
There are two queues. Active nodes are added
1032
to the end of the second queue and read from
1033
the front of the first queue. If the first queue
1034
is empty, it is replaced by the second queue
1035
(and the second queue becomes empty).
1038
inline void Graph::set_active(node *i)
1041
/* it's not in the list yet */
1042
if (queue_last[1]) queue_last[1] -> next = i;
1043
else queue_first[1] = i;
1050
Returns the next active node.
1051
If it is connected to the sink, it stays in the list,
1052
otherwise it is removed from the list
1054
inline Graph::node * Graph::next_active()
1059
if (!(i=queue_first[0])) {
1060
queue_first[0] = i = queue_first[1];
1061
queue_last[0] = queue_last[1];
1062
queue_first[1] = NULL;
1063
queue_last[1] = NULL;
1064
if (!i) return NULL;
1067
/* remove it from the active list */
1068
if (i->next == i) queue_first[0] = queue_last[0] = NULL;
1069
else queue_first[0] = i -> next;
1072
/* a node in the list is active iff it has a parent */
1073
if (i->parent) return i;
1077
/***********************************************************************/
1079
inline void Graph::maxflow_init()
1084
queue_first[0] = queue_last[0] = NULL;
1085
queue_first[1] = queue_last[1] = NULL;
1086
orphan_first = NULL;
1088
for (nb=node_block_first; nb; nb=nb->next)
1089
for (i=&nb->nodes[0]; i<nb->current; i++) {
1092
if (i->tr_cap > 0) {
1093
/* i is connected to the source */
1095
i -> parent = TERMINAL;
1099
} else if (i->tr_cap < 0) {
1100
/* i is connected to the sink */
1102
i -> parent = TERMINAL;
1113
/***********************************************************************/
1115
inline void Graph::augment(node *s_start, node *t_start, captype *cap_middle,
1116
captype *rev_cap_middle)
1124
/* 1. Finding bottleneck capacity */
1125
/* 1a - the source tree */
1126
bottleneck = *cap_middle;
1127
for (i=s_start; ; ) {
1129
if (a == TERMINAL) break;
1132
if (bottleneck > a->r_cap) bottleneck = a -> r_cap;
1133
i = NEIGHBOR_NODE_REV(i, a -> shift);
1135
if (bottleneck > a->r_rev_cap) bottleneck = a -> r_rev_cap;
1136
i = NEIGHBOR_NODE(i, a -> shift);
1139
if (bottleneck > i->tr_cap) bottleneck = i -> tr_cap;
1140
/* 1b - the sink tree */
1141
for (i=t_start; ; ) {
1143
if (a == TERMINAL) break;
1146
if (bottleneck > a->r_rev_cap) bottleneck = a -> r_rev_cap;
1147
i = NEIGHBOR_NODE_REV(i, a -> shift);
1149
if (bottleneck > a->r_cap) bottleneck = a -> r_cap;
1150
i = NEIGHBOR_NODE(i, a -> shift);
1153
if (bottleneck > - i->tr_cap) bottleneck = - i -> tr_cap;
1157
/* 2a - the source tree */
1158
*rev_cap_middle += bottleneck;
1159
*cap_middle -= bottleneck;
1160
for (i=s_start; ; ) {
1162
if (a == TERMINAL) break;
1165
a -> r_rev_cap += bottleneck;
1166
a -> r_cap -= bottleneck;
1168
/* add i to the adoption list */
1169
i -> parent = ORPHAN;
1170
np = nodeptr_block -> New();
1172
np -> next = orphan_first;
1175
i = NEIGHBOR_NODE_REV(i, a -> shift);
1177
a -> r_cap += bottleneck;
1178
a -> r_rev_cap -= bottleneck;
1179
if (!a->r_rev_cap) {
1180
/* add i to the adoption list */
1181
i -> parent = ORPHAN;
1182
np = nodeptr_block -> New();
1184
np -> next = orphan_first;
1187
i = NEIGHBOR_NODE(i, a -> shift);
1190
i -> tr_cap -= bottleneck;
1192
/* add i to the adoption list */
1193
i -> parent = ORPHAN;
1194
np = nodeptr_block -> New();
1196
np -> next = orphan_first;
1199
/* 2b - the sink tree */
1200
for (i=t_start; ; ) {
1202
if (a == TERMINAL) break;
1205
a -> r_cap += bottleneck;
1206
a -> r_rev_cap -= bottleneck;
1207
if (!a->r_rev_cap) {
1208
/* add i to the adoption list */
1209
i -> parent = ORPHAN;
1210
np = nodeptr_block -> New();
1212
np -> next = orphan_first;
1215
i = NEIGHBOR_NODE_REV(i, a -> shift);
1217
a -> r_rev_cap += bottleneck;
1218
a -> r_cap -= bottleneck;
1220
/* add i to the adoption list */
1221
i -> parent = ORPHAN;
1222
np = nodeptr_block -> New();
1224
np -> next = orphan_first;
1227
i = NEIGHBOR_NODE(i, a -> shift);
1230
i -> tr_cap += bottleneck;
1232
/* add i to the adoption list */
1233
i -> parent = ORPHAN;
1234
np = nodeptr_block -> New();
1236
np -> next = orphan_first;
1244
/***********************************************************************/
1246
inline void Graph::process_source_orphan(node *i)
1249
arc_forward *a0_for, *a0_for_first, *a0_for_last;
1250
arc_reverse *a0_rev, *a0_rev_first, *a0_rev_last;
1251
arc_forward *a0_min = NULL, *a;
1253
int d, d_min = INFINITE_D;
1255
/* trying to find a new parent */
1256
a0_for_first = i -> first_out;
1257
if (IS_ODD(a0_for_first)) {
1258
a0_for_first = (arc_forward *) (((char *)a0_for_first) + 1);
1259
a0_for_last = (arc_forward *) ((a0_for_first ++) -> shift);
1260
} else a0_for_last = (i + 1) -> first_out;
1261
a0_rev_first = i -> first_in;
1262
if (IS_ODD(a0_rev_first)) {
1263
a0_rev_first = (arc_reverse *) (((char *)a0_rev_first) + 1);
1264
a0_rev_last = (arc_reverse *) ((a0_rev_first ++) -> sister);
1265
} else a0_rev_last = (i + 1) -> first_in;
1268
for (a0_for=a0_for_first; a0_for<a0_for_last; a0_for++)
1269
if (a0_for->r_rev_cap) {
1270
j = NEIGHBOR_NODE(i, a0_for -> shift);
1271
if (!j->is_sink && (a=j->parent)) {
1272
/* checking the origin of j */
1275
if (j->TS == TIME) {
1291
j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
1293
j = NEIGHBOR_NODE(j, a -> shift);
1295
if (d<INFINITE_D) { /* j originates from the source - done */
1300
/* set marks along the path */
1301
for (j=NEIGHBOR_NODE(i, a0_for->shift); j->TS!=TIME; ) {
1306
j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
1308
j = NEIGHBOR_NODE(j, a -> shift);
1313
for (a0_rev=a0_rev_first; a0_rev<a0_rev_last; a0_rev++) {
1314
a0_for = a0_rev -> sister;
1315
if (a0_for->r_cap) {
1316
j = NEIGHBOR_NODE_REV(i, a0_for -> shift);
1317
if (!j->is_sink && (a=j->parent)) {
1318
/* checking the origin of j */
1321
if (j->TS == TIME) {
1337
j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
1339
j = NEIGHBOR_NODE(j, a -> shift);
1341
if (d<INFINITE_D) { /* j originates from the source - done */
1343
a0_min = MAKE_ODD(a0_for);
1346
/* set marks along the path */
1347
for (j=NEIGHBOR_NODE_REV(i,a0_for->shift); j->TS!=TIME; ) {
1352
j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
1354
j = NEIGHBOR_NODE(j, a -> shift);
1361
if ( (i->parent = a0_min) ) {
1363
i -> DIST = d_min + 1;
1365
/* no parent is found */
1368
/* process neighbors */
1369
for (a0_for=a0_for_first; a0_for<a0_for_last; a0_for++) {
1370
j = NEIGHBOR_NODE(i, a0_for -> shift);
1371
if (!j->is_sink && (a=j->parent)) {
1372
if (a0_for->r_rev_cap) set_active(j);
1373
if (a!=TERMINAL && a!=ORPHAN && IS_ODD(a)
1374
&& NEIGHBOR_NODE_REV(j, MAKE_EVEN(a)->shift)==i) {
1375
/* add j to the adoption list */
1376
j -> parent = ORPHAN;
1377
np = nodeptr_block -> New();
1379
if (orphan_last) orphan_last -> next = np;
1380
else orphan_first = np;
1386
for (a0_rev=a0_rev_first; a0_rev<a0_rev_last; a0_rev++) {
1387
a0_for = a0_rev -> sister;
1388
j = NEIGHBOR_NODE_REV(i, a0_for -> shift);
1389
if (!j->is_sink && (a=j->parent)) {
1390
if (a0_for->r_cap) set_active(j);
1391
if (a!=TERMINAL && a!=ORPHAN && !IS_ODD(a) && NEIGHBOR_NODE(j, a->shift)==i) {
1392
/* add j to the adoption list */
1393
j -> parent = ORPHAN;
1394
np = nodeptr_block -> New();
1396
if (orphan_last) orphan_last -> next = np;
1397
else orphan_first = np;
1406
inline void Graph::process_sink_orphan(node *i)
1409
arc_forward *a0_for, *a0_for_first, *a0_for_last;
1410
arc_reverse *a0_rev, *a0_rev_first, *a0_rev_last;
1411
arc_forward *a0_min = NULL, *a;
1413
int d, d_min = INFINITE_D;
1415
/* trying to find a new parent */
1416
a0_for_first = i -> first_out;
1417
if (IS_ODD(a0_for_first)) {
1418
a0_for_first = (arc_forward *) (((char *)a0_for_first) + 1);
1419
a0_for_last = (arc_forward *) ((a0_for_first ++) -> shift);
1420
} else a0_for_last = (i + 1) -> first_out;
1421
a0_rev_first = i -> first_in;
1422
if (IS_ODD(a0_rev_first)) {
1423
a0_rev_first = (arc_reverse *) (((char *)a0_rev_first) + 1);
1424
a0_rev_last = (arc_reverse *) ((a0_rev_first ++) -> sister);
1425
} else a0_rev_last = (i + 1) -> first_in;
1428
for (a0_for=a0_for_first; a0_for<a0_for_last; a0_for++)
1429
if (a0_for->r_cap) {
1430
j = NEIGHBOR_NODE(i, a0_for -> shift);
1431
if (j->is_sink && (a=j->parent)) {
1432
/* checking the origin of j */
1435
if (j->TS == TIME) {
1451
j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
1453
j = NEIGHBOR_NODE(j, a -> shift);
1455
if (d<INFINITE_D) { /* j originates from the sink - done */
1460
/* set marks along the path */
1461
for (j=NEIGHBOR_NODE(i, a0_for->shift); j->TS!=TIME; ) {
1466
j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
1468
j = NEIGHBOR_NODE(j, a -> shift);
1473
for (a0_rev=a0_rev_first; a0_rev<a0_rev_last; a0_rev++) {
1474
a0_for = a0_rev -> sister;
1475
if (a0_for->r_rev_cap) {
1476
j = NEIGHBOR_NODE_REV(i, a0_for -> shift);
1477
if (j->is_sink && (a=j->parent)) {
1478
/* checking the origin of j */
1481
if (j->TS == TIME) {
1497
j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
1499
j = NEIGHBOR_NODE(j, a -> shift);
1501
if (d<INFINITE_D) { /* j originates from the sink - done */
1503
a0_min = MAKE_ODD(a0_for);
1506
/* set marks along the path */
1507
for (j=NEIGHBOR_NODE_REV(i,a0_for->shift); j->TS!=TIME; ) {
1512
j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
1514
j = NEIGHBOR_NODE(j, a -> shift);
1521
if ( (i->parent = a0_min) ) {
1523
i -> DIST = d_min + 1;
1525
/* no parent is found */
1528
/* process neighbors */
1529
for (a0_for=a0_for_first; a0_for<a0_for_last; a0_for++) {
1530
j = NEIGHBOR_NODE(i, a0_for -> shift);
1531
if (j->is_sink && (a=j->parent)) {
1532
if (a0_for->r_cap) set_active(j);
1533
if (a!=TERMINAL && a!=ORPHAN && IS_ODD(a)
1534
&& NEIGHBOR_NODE_REV(j, MAKE_EVEN(a)->shift)==i) {
1535
/* add j to the adoption list */
1536
j -> parent = ORPHAN;
1537
np = nodeptr_block -> New();
1539
if (orphan_last) orphan_last -> next = np;
1540
else orphan_first = np;
1546
for (a0_rev=a0_rev_first; a0_rev<a0_rev_last; a0_rev++) {
1547
a0_for = a0_rev -> sister;
1548
j = NEIGHBOR_NODE_REV(i, a0_for -> shift);
1549
if (j->is_sink && (a=j->parent)) {
1550
if (a0_for->r_rev_cap) set_active(j);
1551
if (a!=TERMINAL && a!=ORPHAN && !IS_ODD(a) && NEIGHBOR_NODE(j, a->shift)==i) {
1552
/* add j to the adoption list */
1553
j -> parent = ORPHAN;
1554
np = nodeptr_block -> New();
1556
if (orphan_last) orphan_last -> next = np;
1557
else orphan_first = np;
1566
/***********************************************************************/
1568
inline Graph::flowtype Graph::maxflow()
1570
node *i, *j, *current_node = NULL, *s_start, *t_start=NULL;
1571
captype *cap_middle=NULL, *rev_cap_middle=NULL;
1572
arc_forward *a_for, *a_for_first, *a_for_last;
1573
arc_reverse *a_rev, *a_rev_first, *a_rev_last;
1574
nodeptr *np, *np_next;
1578
nodeptr_block = new DBlock<nodeptr>(NODEPTR_BLOCK_SIZE, error_function);
1581
if ( (i=current_node) ) {
1582
i -> next = NULL; /* remove active flag */
1583
if (!i->parent) i = NULL;
1586
if (!(i = next_active())) break;
1592
a_for_first = i -> first_out;
1593
if (IS_ODD(a_for_first)) {
1594
a_for_first = (arc_forward *) (((char *)a_for_first) + 1);
1595
a_for_last = (arc_forward *) ((a_for_first ++) -> shift);
1596
} else a_for_last = (i + 1) -> first_out;
1597
a_rev_first = i -> first_in;
1598
if (IS_ODD(a_rev_first)) {
1599
a_rev_first = (arc_reverse *) (((char *)a_rev_first) + 1);
1600
a_rev_last = (arc_reverse *) ((a_rev_first ++) -> sister);
1601
} else a_rev_last = (i + 1) -> first_in;
1604
/* grow source tree */
1605
for (a_for=a_for_first; a_for<a_for_last; a_for++)
1607
j = NEIGHBOR_NODE(i, a_for -> shift);
1610
j -> parent = MAKE_ODD(a_for);
1612
j -> DIST = i -> DIST + 1;
1614
} else if (j->is_sink) {
1617
cap_middle = & ( a_for -> r_cap );
1618
rev_cap_middle = & ( a_for -> r_rev_cap );
1620
} else if (j->TS <= i->TS &&
1621
j->DIST > i->DIST) {
1622
/* heuristic - trying to make the distance from j to the source shorter */
1623
j -> parent = MAKE_ODD(a_for);
1625
j -> DIST = i -> DIST + 1;
1629
for (a_rev=a_rev_first; a_rev<a_rev_last; a_rev++) {
1630
a_for = a_rev -> sister;
1631
if (a_for->r_rev_cap) {
1632
j = NEIGHBOR_NODE_REV(i, a_for -> shift);
1635
j -> parent = a_for;
1637
j -> DIST = i -> DIST + 1;
1639
} else if (j->is_sink) {
1642
cap_middle = & ( a_for -> r_rev_cap );
1643
rev_cap_middle = & ( a_for -> r_cap );
1645
} else if (j->TS <= i->TS &&
1646
j->DIST > i->DIST) {
1647
/* heuristic - trying to make the distance from j to the source shorter */
1648
j -> parent = a_for;
1650
j -> DIST = i -> DIST + 1;
1655
/* grow sink tree */
1656
for (a_for=a_for_first; a_for<a_for_last; a_for++)
1657
if (a_for->r_rev_cap) {
1658
j = NEIGHBOR_NODE(i, a_for -> shift);
1661
j -> parent = MAKE_ODD(a_for);
1663
j -> DIST = i -> DIST + 1;
1665
} else if (!j->is_sink) {
1668
cap_middle = & ( a_for -> r_rev_cap );
1669
rev_cap_middle = & ( a_for -> r_cap );
1671
} else if (j->TS <= i->TS &&
1672
j->DIST > i->DIST) {
1673
/* heuristic - trying to make the distance from j to the sink shorter */
1674
j -> parent = MAKE_ODD(a_for);
1676
j -> DIST = i -> DIST + 1;
1679
for (a_rev=a_rev_first; a_rev<a_rev_last; a_rev++) {
1680
a_for = a_rev -> sister;
1682
j = NEIGHBOR_NODE_REV(i, a_for -> shift);
1685
j -> parent = a_for;
1687
j -> DIST = i -> DIST + 1;
1689
} else if (!j->is_sink) {
1692
cap_middle = & ( a_for -> r_cap );
1693
rev_cap_middle = & ( a_for -> r_rev_cap );
1695
} else if (j->TS <= i->TS &&
1696
j->DIST > i->DIST) {
1697
/* heuristic - trying to make the distance from j to the sink shorter */
1698
j -> parent = a_for;
1700
j -> DIST = i -> DIST + 1;
1709
i -> next = i; /* set active flag */
1713
augment(s_start, t_start, cap_middle, rev_cap_middle);
1714
/* augmentation end */
1717
while ( (np=orphan_first) ) {
1718
np_next = np -> next;
1721
while ( (np=orphan_first) ) {
1722
orphan_first = np -> next;
1724
nodeptr_block -> Delete(np);
1725
if (!orphan_first) orphan_last = NULL;
1726
if (i->is_sink) process_sink_orphan(i);
1727
else process_source_orphan(i);
1730
orphan_first = np_next;
1733
} else current_node = NULL;
1736
delete nodeptr_block;
1741
/***********************************************************************/
1743
inline Graph::termtype Graph::what_segment(node_id i)
1745
if (((node*)i)->parent && !((node*)i)->is_sink) return SOURCE;
added
removed
include/CGAL/internal/auxiliary/graph.h
