www.pudn.com > GraphCut.zip > Graph.cpp
// Graph.cpp: implementation of the Graph class.
//
//////////////////////////////////////////////////////////////////////
/* graph.cpp */
/*
Copyright 2001 Vladimir Kolmogorov (vnk@cs.cornell.edu), Yuri Boykov (yuri@csd.uwo.ca).
This program is free software; you can redistribute it and/or modify
it 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.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
#include
#include "stdafx.h"
#include "GraphCut.h"
#include "Graph.h"
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
/*
special constants for node->parent
*/
#define TERMINAL ( (arc_forward *) 1 ) /* to terminal */
#define ORPHAN ( (arc_forward *) 2 ) /* orphan */
#define INFINITE_D 1000000000 /* infinite distance to the terminal */
/***********************************************************************/
Graph::Graph(void (*err_function)(char *))
{
error_function = err_function;
node_block_first = NULL;
arc_for_block_first = NULL;
arc_rev_block_first = NULL;
flow = 0;
}
Graph::~Graph()
{
while (node_block_first)
{
node_block *next = node_block_first -> next;
delete node_block_first;
node_block_first = next;
}
while (arc_for_block_first)
{
arc_for_block *next = arc_for_block_first -> next;
delete arc_for_block_first -> start;
arc_for_block_first = next;
}
while (arc_rev_block_first)
{
arc_rev_block *next = arc_rev_block_first -> next;
delete arc_rev_block_first -> start;
arc_rev_block_first = next;
}
}
Graph::node_id Graph::add_node()
{
node *i;
if (!node_block_first || node_block_first->current+1 > &node_block_first->nodes[NODE_BLOCK_SIZE-1])
{
node_block *next = node_block_first;
node_block_first = (node_block *) new node_block;
if (!node_block_first) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
node_block_first -> current = & ( node_block_first -> nodes[0] );
node_block_first -> next = next;
}
i = node_block_first -> current ++;
i -> first_out = (arc_forward *) 0;
i -> first_in = (arc_reverse *) 0;
i -> tr_cap = 0;
return (node_id) i;
}
void Graph::add_edge(node_id from, node_id to, captype cap, captype rev_cap)
{
arc_forward *a_for;
arc_reverse *a_rev;
if (!arc_for_block_first || arc_for_block_first->current+1 > &arc_for_block_first->arcs_for[ARC_BLOCK_SIZE])
{
arc_for_block *next = arc_for_block_first;
char *ptr = new char[sizeof(arc_for_block)+1];
if (!ptr) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
if ((int)ptr & 1) arc_for_block_first = (arc_for_block *) (ptr + 1);
else arc_for_block_first = (arc_for_block *) ptr;
arc_for_block_first -> start = ptr;
arc_for_block_first -> current = & ( arc_for_block_first -> arcs_for[0] );
arc_for_block_first -> next = next;
}
if (!arc_rev_block_first || arc_rev_block_first->current+1 > &arc_rev_block_first->arcs_rev[ARC_BLOCK_SIZE])
{
arc_rev_block *next = arc_rev_block_first;
char *ptr = new char[sizeof(arc_rev_block)+1];
if (!ptr) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
if ((int)ptr & 1) arc_rev_block_first = (arc_rev_block *) (ptr + 1);
else arc_rev_block_first = (arc_rev_block *) ptr;
arc_rev_block_first -> start = ptr;
arc_rev_block_first -> current = & ( arc_rev_block_first -> arcs_rev[0] );
arc_rev_block_first -> next = next;
}
a_for = arc_for_block_first -> current ++;
a_rev = arc_rev_block_first -> current ++;
a_rev -> sister = (arc_forward *) from;
a_for -> shift = (int) to;
a_for -> r_cap = cap;
a_for -> r_rev_cap = rev_cap;
((node *)from) -> first_out =
(arc_forward *) ((int)(((node *)from) -> first_out) + 1);
((node *)to) -> first_in =
(arc_reverse *) ((int)(((node *)to) -> first_in) + 1);
}
void Graph::set_tweights(node_id i, captype cap_source, captype cap_sink)
{
flow += (cap_source < cap_sink) ? cap_source : cap_sink;
((node*)i) -> tr_cap = cap_source - cap_sink;
}
void Graph::add_tweights(node_id i, captype cap_source, captype cap_sink)
{
register captype delta = ((node*)i) -> tr_cap;
if (delta > 0) cap_source += delta;
else cap_sink -= delta;
flow += (cap_source < cap_sink) ? cap_source : cap_sink;
((node*)i) -> tr_cap = cap_source - cap_sink;
}
/*
Converts arcs added by 'add_edge()' calls
to a forward star graph representation.
Linear time algorithm.
No or little additional memory is allocated
during this process
(it may be necessary to allocate additional
arc blocks, since arcs corresponding to the
same node must be contiguous, i.e. be in one
arc block.)
*/
void Graph::prepare_graph()
{
node *i;
arc_for_block *ab_for, *ab_for_first;
arc_rev_block *ab_rev, *ab_rev_first, *ab_rev_scan;
arc_forward *a_for;
arc_reverse *a_rev, *a_rev_scan, a_rev_tmp;
node_block *nb;
bool for_flag = false, rev_flag = false;
int k;
if (!arc_rev_block_first)
{
node_id from = add_node(), to = add_node();
add_edge(from, to, 1, 0);
}
/* FIRST STAGE */
a_rev_tmp.sister = NULL;
for (a_rev=arc_rev_block_first->current; a_rev<&arc_rev_block_first->arcs_rev[ARC_BLOCK_SIZE]; a_rev++)
{
a_rev -> sister = NULL;
}
ab_for = ab_for_first = arc_for_block_first;
ab_rev = ab_rev_first = ab_rev_scan = arc_rev_block_first;
a_for = &ab_for->arcs_for[0];
a_rev = a_rev_scan = &ab_rev->arcs_rev[0];
for (nb=node_block_first; nb; nb=nb->next)
{
for (i=&nb->nodes[0]; icurrent; i++)
{
/* outgoing arcs */
k = (int) i -> first_out;
if (a_for + k > &ab_for->arcs_for[ARC_BLOCK_SIZE])
{
if (k > ARC_BLOCK_SIZE) { if (error_function) (*error_function)("# of arcs per node exceeds block size!"); exit(1); }
if (for_flag) ab_for = NULL;
else { ab_for = ab_for -> next; ab_rev_scan = ab_rev_scan -> next; }
if (ab_for == NULL)
{
arc_for_block *next = arc_for_block_first;
char *ptr = new char[sizeof(arc_for_block)+1];
if (!ptr) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
if ((int)ptr & 1) arc_for_block_first = (arc_for_block *) (ptr + 1);
else arc_for_block_first = (arc_for_block *) ptr;
arc_for_block_first -> start = ptr;
arc_for_block_first -> current = & ( arc_for_block_first -> arcs_for[0] );
arc_for_block_first -> next = next;
ab_for = arc_for_block_first;
for_flag = true;
}
else a_rev_scan = &ab_rev_scan->arcs_rev[0];
a_for = &ab_for->arcs_for[0];
}
if (ab_rev_scan)
{
a_rev_scan += k;
i -> parent = (arc_forward *) a_rev_scan;
}
else i -> parent = (arc_forward *) &a_rev_tmp;
a_for += k;
i -> first_out = a_for;
ab_for -> last_node = i;
/* incoming arcs */
k = (int) i -> first_in;
if (a_rev + k > &ab_rev->arcs_rev[ARC_BLOCK_SIZE])
{
if (k > ARC_BLOCK_SIZE) { if (error_function) (*error_function)("# of arcs per node exceeds block size!"); exit(1); }
if (rev_flag) ab_rev = NULL;
else ab_rev = ab_rev -> next;
if (ab_rev == NULL)
{
arc_rev_block *next = arc_rev_block_first;
char *ptr = new char[sizeof(arc_rev_block)+1];
if (!ptr) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
if ((int)ptr & 1) arc_rev_block_first = (arc_rev_block *) (ptr + 1);
else arc_rev_block_first = (arc_rev_block *) ptr;
arc_rev_block_first -> start = ptr;
arc_rev_block_first -> current = & ( arc_rev_block_first -> arcs_rev[0] );
arc_rev_block_first -> next = next;
ab_rev = arc_rev_block_first;
rev_flag = true;
}
a_rev = &ab_rev->arcs_rev[0];
}
a_rev += k;
i -> first_in = a_rev;
ab_rev -> last_node = i;
}
/* i is the last node in block */
i -> first_out = a_for;
i -> first_in = a_rev;
}
/* SECOND STAGE */
for (ab_for=arc_for_block_first; ab_for; ab_for=ab_for->next)
{
ab_for -> current = ab_for -> last_node -> first_out;
}
for ( ab_for=ab_for_first, ab_rev=ab_rev_first;
ab_for;
ab_for=ab_for->next, ab_rev=ab_rev->next )
for ( a_for=&ab_for->arcs_for[0], a_rev=&ab_rev->arcs_rev[0];
a_for<&ab_for->arcs_for[ARC_BLOCK_SIZE];
a_for++, a_rev++ )
{
arc_forward *af;
arc_reverse *ar;
node *from;
int shift = 0, shift_new;
captype r_cap, r_rev_cap, r_cap_new, r_rev_cap_new;
if (!(from=(node *)(a_rev->sister))) continue;
af = a_for;
ar = a_rev;
do
{
ar -> sister = NULL;
shift_new = ((char *)(af->shift)) - (char *)from;
r_cap_new = af -> r_cap;
r_rev_cap_new = af -> r_rev_cap;
if (shift)
{
af -> shift = shift;
af -> r_cap = r_cap;
af -> r_rev_cap = r_rev_cap;
}
shift = shift_new;
r_cap = r_cap_new;
r_rev_cap = r_rev_cap_new;
af = -- from -> first_out;
if ((arc_reverse *)(from->parent) != &a_rev_tmp)
{
from -> parent = (arc_forward *)(((arc_reverse *)(from -> parent)) - 1);
ar = (arc_reverse *)(from -> parent);
}
} while (from=(node *)(ar->sister));
af -> shift = shift;
af -> r_cap = r_cap;
af -> r_rev_cap = r_rev_cap;
}
for (ab_for=arc_for_block_first; ab_for; ab_for=ab_for->next)
{
i = ab_for -> last_node;
a_for = i -> first_out;
ab_for -> current -> shift = a_for -> shift;
ab_for -> current -> r_cap = a_for -> r_cap;
ab_for -> current -> r_rev_cap = a_for -> r_rev_cap;
a_for -> shift = (int) (ab_for -> current + 1);
i -> first_out = (arc_forward *) (((char *)a_for) - 1);
}
/* THIRD STAGE */
for (ab_rev=arc_rev_block_first; ab_rev; ab_rev=ab_rev->next)
{
ab_rev -> current = ab_rev -> last_node -> first_in;
}
for (nb=node_block_first; nb; nb=nb->next)
for (i=&nb->nodes[0]; icurrent; i++)
{
arc_forward *a_for_first, *a_for_last;
a_for_first = i -> first_out;
if (IS_ODD(a_for_first))
{
a_for_first = (arc_forward *) (((char *)a_for_first) + 1);
a_for_last = (arc_forward *) ((a_for_first ++) -> shift);
}
else a_for_last = (i + 1) -> first_out;
for (a_for=a_for_first; a_for shift);
a_rev = -- to -> first_in;
a_rev -> sister = a_for;
}
}
for (ab_rev=arc_rev_block_first; ab_rev; ab_rev=ab_rev->next)
{
i = ab_rev -> last_node;
a_rev = i -> first_in;
ab_rev -> current -> sister = a_rev -> sister;
a_rev -> sister = (arc_forward *) (ab_rev -> current + 1);
i -> first_in = (arc_reverse *) (((char *)a_rev) - 1);
}
}
/* maxflow.cpp */
/*
Copyright 2001 Vladimir Kolmogorov (vnk@cs.cornell.edu), Yuri Boykov (yuri@csd.uwo.ca).
This program is free software; you can redistribute it and/or modify
it 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.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
/*
Functions for processing active list.
i->next points to the next node in the list
(or to i, if i is the last node in the list).
If i->next is NULL iff i is not in the list.
There are two queues. Active nodes are added
to the end of the second queue and read from
the front of the first queue. If the first queue
is empty, it is replaced by the second queue
(and the second queue becomes empty).
*/
inline void Graph::set_active(node *i)
{
if (!i->next)
{
/* it's not in the list yet */
if (queue_last[1]) queue_last[1] -> next = i;
else queue_first[1] = i;
queue_last[1] = i;
i -> next = i;
}
}
/*
Returns the next active node.
If it is connected to the sink, it stays in the list,
otherwise it is removed from the list
*/
inline Graph::node * Graph::next_active()
{
node *i;
while ( 1 )
{
if (!(i=queue_first[0]))
{
queue_first[0] = i = queue_first[1];
queue_last[0] = queue_last[1];
queue_first[1] = NULL;
queue_last[1] = NULL;
if (!i) return NULL;
}
/* remove it from the active list */
if (i->next == i) queue_first[0] = queue_last[0] = NULL;
else queue_first[0] = i -> next;
i -> next = NULL;
/* a node in the list is active iff it has a parent */
if (i->parent) return i;
}
}
/***********************************************************************/
void Graph::maxflow_init()
{
node *i;
node_block *nb;
queue_first[0] = queue_last[0] = NULL;
queue_first[1] = queue_last[1] = NULL;
orphan_first = NULL;
for (nb=node_block_first; nb; nb=nb->next)
for (i=&nb->nodes[0]; icurrent; i++)
{
i -> next = NULL;
i -> TS = 0;
if (i->tr_cap > 0)
{
/* i is connected to the source */
i -> is_sink = 0;
i -> parent = TERMINAL;
set_active(i);
i -> TS = 0;
i -> DIST = 1;
}
else if (i->tr_cap < 0)
{
/* i is connected to the sink */
i -> is_sink = 1;
i -> parent = TERMINAL;
set_active(i);
i -> TS = 0;
i -> DIST = 1;
}
else
{
i -> parent = NULL;
}
}
TIME = 0;
}
/***********************************************************************/
void Graph::augment(node *s_start, node *t_start, captype *cap_middle, captype *rev_cap_middle)
{
node *i;
arc_forward *a;
captype bottleneck;
nodeptr *np;
/* 1. Finding bottleneck capacity */
/* 1a - the source tree */
bottleneck = *cap_middle;
for (i=s_start; ; )
{
a = i -> parent;
if (a == TERMINAL) break;
if (IS_ODD(a))
{
a = MAKE_EVEN(a);
if (bottleneck > a->r_cap) bottleneck = a -> r_cap;
i = NEIGHBOR_NODE_REV(i, a -> shift);
}
else
{
if (bottleneck > a->r_rev_cap) bottleneck = a -> r_rev_cap;
i = NEIGHBOR_NODE(i, a -> shift);
}
}
if (bottleneck > i->tr_cap) bottleneck = i -> tr_cap;
/* 1b - the sink tree */
for (i=t_start; ; )
{
a = i -> parent;
if (a == TERMINAL) break;
if (IS_ODD(a))
{
a = MAKE_EVEN(a);
if (bottleneck > a->r_rev_cap) bottleneck = a -> r_rev_cap;
i = NEIGHBOR_NODE_REV(i, a -> shift);
}
else
{
if (bottleneck > a->r_cap) bottleneck = a -> r_cap;
i = NEIGHBOR_NODE(i, a -> shift);
}
}
if (bottleneck > - i->tr_cap) bottleneck = - i -> tr_cap;
/* 2. Augmenting */
/* 2a - the source tree */
*rev_cap_middle += bottleneck;
*cap_middle -= bottleneck;
for (i=s_start; ; )
{
a = i -> parent;
if (a == TERMINAL) break;
if (IS_ODD(a))
{
a = MAKE_EVEN(a);
a -> r_rev_cap += bottleneck;
a -> r_cap -= bottleneck;
if (!a->r_cap)
{
/* add i to the adoption list */
i -> parent = ORPHAN;
np = nodeptr_block -> New();
np -> ptr = i;
np -> next = orphan_first;
orphan_first = np;
}
i = NEIGHBOR_NODE_REV(i, a -> shift);
}
else
{
a -> r_cap += bottleneck;
a -> r_rev_cap -= bottleneck;
if (!a->r_rev_cap)
{
/* add i to the adoption list */
i -> parent = ORPHAN;
np = nodeptr_block -> New();
np -> ptr = i;
np -> next = orphan_first;
orphan_first = np;
}
i = NEIGHBOR_NODE(i, a -> shift);
}
}
i -> tr_cap -= bottleneck;
if (!i->tr_cap)
{
/* add i to the adoption list */
i -> parent = ORPHAN;
np = nodeptr_block -> New();
np -> ptr = i;
np -> next = orphan_first;
orphan_first = np;
}
/* 2b - the sink tree */
for (i=t_start; ; )
{
a = i -> parent;
if (a == TERMINAL) break;
if (IS_ODD(a))
{
a = MAKE_EVEN(a);
a -> r_cap += bottleneck;
a -> r_rev_cap -= bottleneck;
if (!a->r_rev_cap)
{
/* add i to the adoption list */
i -> parent = ORPHAN;
np = nodeptr_block -> New();
np -> ptr = i;
np -> next = orphan_first;
orphan_first = np;
}
i = NEIGHBOR_NODE_REV(i, a -> shift);
}
else
{
a -> r_rev_cap += bottleneck;
a -> r_cap -= bottleneck;
if (!a->r_cap)
{
/* add i to the adoption list */
i -> parent = ORPHAN;
np = nodeptr_block -> New();
np -> ptr = i;
np -> next = orphan_first;
orphan_first = np;
}
i = NEIGHBOR_NODE(i, a -> shift);
}
}
i -> tr_cap += bottleneck;
if (!i->tr_cap)
{
/* add i to the adoption list */
i -> parent = ORPHAN;
np = nodeptr_block -> New();
np -> ptr = i;
np -> next = orphan_first;
orphan_first = np;
}
flow += bottleneck;
}
/***********************************************************************/
void Graph::process_source_orphan(node *i)
{
node *j;
arc_forward *a0_for, *a0_for_first, *a0_for_last;
arc_reverse *a0_rev, *a0_rev_first, *a0_rev_last;
arc_forward *a0_min = NULL, *a;
nodeptr *np;
int d, d_min = INFINITE_D;
/* trying to find a new parent */
a0_for_first = i -> first_out;
if (IS_ODD(a0_for_first))
{
a0_for_first = (arc_forward *) (((char *)a0_for_first) + 1);
a0_for_last = (arc_forward *) ((a0_for_first ++) -> shift);
}
else a0_for_last = (i + 1) -> first_out;
a0_rev_first = i -> first_in;
if (IS_ODD(a0_rev_first))
{
a0_rev_first = (arc_reverse *) (((char *)a0_rev_first) + 1);
a0_rev_last = (arc_reverse *) ((a0_rev_first ++) -> sister);
}
else a0_rev_last = (i + 1) -> first_in;
for (a0_for=a0_for_first; a0_forr_rev_cap)
{
j = NEIGHBOR_NODE(i, a0_for -> shift);
if (!j->is_sink && (a=j->parent))
{
/* checking the origin of j */
d = 0;
while ( 1 )
{
if (j->TS == TIME)
{
d += j -> DIST;
break;
}
a = j -> parent;
d ++;
if (a==TERMINAL)
{
j -> TS = TIME;
j -> DIST = 1;
break;
}
if (a==ORPHAN) { d = INFINITE_D; break; }
if (IS_ODD(a))
j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
else
j = NEIGHBOR_NODE(j, a -> shift);
}
if (dshift); j->TS!=TIME; )
{
j -> TS = TIME;
j -> DIST = d --;
a = j->parent;
if (IS_ODD(a))
j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
else
j = NEIGHBOR_NODE(j, a -> shift);
}
}
}
}
for (a0_rev=a0_rev_first; a0_rev sister;
if (a0_for->r_cap)
{
j = NEIGHBOR_NODE_REV(i, a0_for -> shift);
if (!j->is_sink && (a=j->parent))
{
/* checking the origin of j */
d = 0;
while ( 1 )
{
if (j->TS == TIME)
{
d += j -> DIST;
break;
}
a = j -> parent;
d ++;
if (a==TERMINAL)
{
j -> TS = TIME;
j -> DIST = 1;
break;
}
if (a==ORPHAN) { d = INFINITE_D; break; }
if (IS_ODD(a))
j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
else
j = NEIGHBOR_NODE(j, a -> shift);
}
if (dshift); j->TS!=TIME; )
{
j -> TS = TIME;
j -> DIST = d --;
a = j->parent;
if (IS_ODD(a))
j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
else
j = NEIGHBOR_NODE(j, a -> shift);
}
}
}
}
}
if (i->parent = a0_min)
{
i -> TS = TIME;
i -> DIST = d_min + 1;
}
else
{
/* no parent is found */
i -> TS = 0;
/* process neighbors */
for (a0_for=a0_for_first; a0_for shift);
if (!j->is_sink && (a=j->parent))
{
if (a0_for->r_rev_cap) set_active(j);
if (a!=TERMINAL && a!=ORPHAN && IS_ODD(a) && NEIGHBOR_NODE_REV(j, MAKE_EVEN(a)->shift)==i)
{
/* add j to the adoption list */
j -> parent = ORPHAN;
np = nodeptr_block -> New();
np -> ptr = j;
if (orphan_last) orphan_last -> next = np;
else orphan_first = np;
orphan_last = np;
np -> next = NULL;
}
}
}
for (a0_rev=a0_rev_first; a0_rev sister;
j = NEIGHBOR_NODE_REV(i, a0_for -> shift);
if (!j->is_sink && (a=j->parent))
{
if (a0_for->r_cap) set_active(j);
if (a!=TERMINAL && a!=ORPHAN && !IS_ODD(a) && NEIGHBOR_NODE(j, a->shift)==i)
{
/* add j to the adoption list */
j -> parent = ORPHAN;
np = nodeptr_block -> New();
np -> ptr = j;
if (orphan_last) orphan_last -> next = np;
else orphan_first = np;
orphan_last = np;
np -> next = NULL;
}
}
}
}
}
void Graph::process_sink_orphan(node *i)
{
node *j;
arc_forward *a0_for, *a0_for_first, *a0_for_last;
arc_reverse *a0_rev, *a0_rev_first, *a0_rev_last;
arc_forward *a0_min = NULL, *a;
nodeptr *np;
int d, d_min = INFINITE_D;
/* trying to find a new parent */
a0_for_first = i -> first_out;
if (IS_ODD(a0_for_first))
{
a0_for_first = (arc_forward *) (((char *)a0_for_first) + 1);
a0_for_last = (arc_forward *) ((a0_for_first ++) -> shift);
}
else a0_for_last = (i + 1) -> first_out;
a0_rev_first = i -> first_in;
if (IS_ODD(a0_rev_first))
{
a0_rev_first = (arc_reverse *) (((char *)a0_rev_first) + 1);
a0_rev_last = (arc_reverse *) ((a0_rev_first ++) -> sister);
}
else a0_rev_last = (i + 1) -> first_in;
for (a0_for=a0_for_first; a0_forr_cap)
{
j = NEIGHBOR_NODE(i, a0_for -> shift);
if (j->is_sink && (a=j->parent))
{
/* checking the origin of j */
d = 0;
while ( 1 )
{
if (j->TS == TIME)
{
d += j -> DIST;
break;
}
a = j -> parent;
d ++;
if (a==TERMINAL)
{
j -> TS = TIME;
j -> DIST = 1;
break;
}
if (a==ORPHAN) { d = INFINITE_D; break; }
if (IS_ODD(a))
j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
else
j = NEIGHBOR_NODE(j, a -> shift);
}
if (dshift); j->TS!=TIME; )
{
j -> TS = TIME;
j -> DIST = d --;
a = j->parent;
if (IS_ODD(a))
j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
else
j = NEIGHBOR_NODE(j, a -> shift);
}
}
}
}
for (a0_rev=a0_rev_first; a0_rev sister;
if (a0_for->r_rev_cap)
{
j = NEIGHBOR_NODE_REV(i, a0_for -> shift);
if (j->is_sink && (a=j->parent))
{
/* checking the origin of j */
d = 0;
while ( 1 )
{
if (j->TS == TIME)
{
d += j -> DIST;
break;
}
a = j -> parent;
d ++;
if (a==TERMINAL)
{
j -> TS = TIME;
j -> DIST = 1;
break;
}
if (a==ORPHAN) { d = INFINITE_D; break; }
if (IS_ODD(a))
j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
else
j = NEIGHBOR_NODE(j, a -> shift);
}
if (dshift); j->TS!=TIME; )
{
j -> TS = TIME;
j -> DIST = d --;
a = j->parent;
if (IS_ODD(a))
j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
else
j = NEIGHBOR_NODE(j, a -> shift);
}
}
}
}
}
if (i->parent = a0_min)
{
i -> TS = TIME;
i -> DIST = d_min + 1;
}
else
{
/* no parent is found */
i -> TS = 0;
/* process neighbors */
for (a0_for=a0_for_first; a0_for shift);
if (j->is_sink && (a=j->parent))
{
if (a0_for->r_cap) set_active(j);
if (a!=TERMINAL && a!=ORPHAN && IS_ODD(a) && NEIGHBOR_NODE_REV(j, MAKE_EVEN(a)->shift)==i)
{
/* add j to the adoption list */
j -> parent = ORPHAN;
np = nodeptr_block -> New();
np -> ptr = j;
if (orphan_last) orphan_last -> next = np;
else orphan_first = np;
orphan_last = np;
np -> next = NULL;
}
}
}
for (a0_rev=a0_rev_first; a0_rev sister;
j = NEIGHBOR_NODE_REV(i, a0_for -> shift);
if (j->is_sink && (a=j->parent))
{
if (a0_for->r_rev_cap) set_active(j);
if (a!=TERMINAL && a!=ORPHAN && !IS_ODD(a) && NEIGHBOR_NODE(j, a->shift)==i)
{
/* add j to the adoption list */
j -> parent = ORPHAN;
np = nodeptr_block -> New();
np -> ptr = j;
if (orphan_last) orphan_last -> next = np;
else orphan_first = np;
orphan_last = np;
np -> next = NULL;
}
}
}
}
}
/***********************************************************************/
Graph::flowtype Graph::maxflow()
{
node *i, *j, *current_node = NULL, *s_start, *t_start;
captype *cap_middle, *rev_cap_middle;
arc_forward *a_for, *a_for_first, *a_for_last;
arc_reverse *a_rev, *a_rev_first, *a_rev_last;
nodeptr *np, *np_next;
prepare_graph();
maxflow_init();
nodeptr_block = new DBlock(NODEPTR_BLOCK_SIZE, error_function);
while ( 1 )
{
if (i=current_node)
{
i -> next = NULL; /* remove active flag */
if (!i->parent) i = NULL;
}
if (!i)
{
if (!(i = next_active())) break;
}
/* growth */
s_start = NULL;
a_for_first = i -> first_out;
if (IS_ODD(a_for_first))
{
a_for_first = (arc_forward *) (((char *)a_for_first) + 1);
a_for_last = (arc_forward *) ((a_for_first ++) -> shift);
}
else a_for_last = (i + 1) -> first_out;
a_rev_first = i -> first_in;
if (IS_ODD(a_rev_first))
{
a_rev_first = (arc_reverse *) (((char *)a_rev_first) + 1);
a_rev_last = (arc_reverse *) ((a_rev_first ++) -> sister);
}
else a_rev_last = (i + 1) -> first_in;
if (!i->is_sink)
{
/* grow source tree */
for (a_for=a_for_first; a_forr_cap)
{
j = NEIGHBOR_NODE(i, a_for -> shift);
if (!j->parent)
{
j -> is_sink = 0;
j -> parent = MAKE_ODD(a_for);
j -> TS = i -> TS;
j -> DIST = i -> DIST + 1;
set_active(j);
}
else if (j->is_sink)
{
s_start = i;
t_start = j;
cap_middle = & ( a_for -> r_cap );
rev_cap_middle = & ( a_for -> r_rev_cap );
break;
}
else if (j->TS <= i->TS &&
j->DIST > i->DIST)
{
/* heuristic - trying to make the distance from j to the source shorter */
j -> parent = MAKE_ODD(a_for);
j -> TS = i -> TS;
j -> DIST = i -> DIST + 1;
}
}
if (!s_start)
for (a_rev=a_rev_first; a_rev sister;
if (a_for->r_rev_cap)
{
j = NEIGHBOR_NODE_REV(i, a_for -> shift);
if (!j->parent)
{
j -> is_sink = 0;
j -> parent = a_for;
j -> TS = i -> TS;
j -> DIST = i -> DIST + 1;
set_active(j);
}
else if (j->is_sink)
{
s_start = i;
t_start = j;
cap_middle = & ( a_for -> r_rev_cap );
rev_cap_middle = & ( a_for -> r_cap );
break;
}
else if (j->TS <= i->TS &&
j->DIST > i->DIST)
{
/* heuristic - trying to make the distance from j to the source shorter */
j -> parent = a_for;
j -> TS = i -> TS;
j -> DIST = i -> DIST + 1;
}
}
}
}
else
{
/* grow sink tree */
for (a_for=a_for_first; a_forr_rev_cap)
{
j = NEIGHBOR_NODE(i, a_for -> shift);
if (!j->parent)
{
j -> is_sink = 1;
j -> parent = MAKE_ODD(a_for);
j -> TS = i -> TS;
j -> DIST = i -> DIST + 1;
set_active(j);
}
else if (!j->is_sink)
{
s_start = j;
t_start = i;
cap_middle = & ( a_for -> r_rev_cap );
rev_cap_middle = & ( a_for -> r_cap );
break;
}
else if (j->TS <= i->TS &&
j->DIST > i->DIST)
{
/* heuristic - trying to make the distance from j to the sink shorter */
j -> parent = MAKE_ODD(a_for);
j -> TS = i -> TS;
j -> DIST = i -> DIST + 1;
}
}
for (a_rev=a_rev_first; a_rev sister;
if (a_for->r_cap)
{
j = NEIGHBOR_NODE_REV(i, a_for -> shift);
if (!j->parent)
{
j -> is_sink = 1;
j -> parent = a_for;
j -> TS = i -> TS;
j -> DIST = i -> DIST + 1;
set_active(j);
}
else if (!j->is_sink)
{
s_start = j;
t_start = i;
cap_middle = & ( a_for -> r_cap );
rev_cap_middle = & ( a_for -> r_rev_cap );
break;
}
else if (j->TS <= i->TS &&
j->DIST > i->DIST)
{
/* heuristic - trying to make the distance from j to the sink shorter */
j -> parent = a_for;
j -> TS = i -> TS;
j -> DIST = i -> DIST + 1;
}
}
}
}
TIME ++;
if (s_start)
{
i -> next = i; /* set active flag */
current_node = i;
/* augmentation */
augment(s_start, t_start, cap_middle, rev_cap_middle);
/* augmentation end */
/* adoption */
while (np=orphan_first)
{
np_next = np -> next;
np -> next = NULL;
while (np=orphan_first)
{
orphan_first = np -> next;
i = np -> ptr;
nodeptr_block -> Delete(np);
if (!orphan_first) orphan_last = NULL;
if (i->is_sink) process_sink_orphan(i);
else process_source_orphan(i);
}
orphan_first = np_next;
}
/* adoption end */
}
else current_node = NULL;
}
delete nodeptr_block;
return flow;
}
/***********************************************************************/
Graph::termtype Graph::what_segment(node_id i)
{
if (((node*)i)->parent && !((node*)i)->is_sink) return SOURCE;
return SINK;
}