added chapter 7

Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.cpp

+#include <math.h>
+#include "mex.h"
+
+#include "graph.h"
+#include <vector>
+
+void
+mexFunction(int nout, mxArray *out[], 
+            int nin, const mxArray *in[])
+{
+
+    if (nin != 3) 
+        mexErrMsgTxt("Three arguments are required (nNodes,TerminalWeights,EdgeWeights)") ;
+    if (nout > 2) 
+        mexErrMsgTxt("Too many output arguments.");
+  
+
+    int nNodes = (int) *mxGetPr(in[0]);
+    
+
+    const int* twDim = mxGetDimensions(in[1]) ;
+    int twRows = twDim[0] ;
+    int twCols = twDim[1] ;
+    double* twPt = mxGetPr(in[1]) ;
+    if(twCols!=3)
+        mexErrMsgTxt("The Terminal Weight matix should have 3 columns, (Node,sink,source).");
+
+    
+    const int* ewDim = mxGetDimensions(in[2]) ;
+    int ewRows = ewDim[0] ;
+    int ewCols = ewDim[1] ;
+    double* ewPt = mxGetPr(in[2]) ;
+    if(ewCols!=4)
+        mexErrMsgTxt("The Terminal Weight matix should have 4 columns, (From,To,Capacity,Rev_Capacity).");
+    
+  
+    typedef Graph<double,double,double> GraphType;
+	GraphType G(static_cast<int>(nNodes), static_cast<int>(ewRows+twRows)); 
+    G.add_node(nNodes);  
+  
+    for(int cTw=0;cTw<twRows;cTw++)
+    {
+        //Test for nodes in range
+        int node=(int)twPt[cTw]-1;
+        if(node<0 || node>=nNodes)
+            mexErrMsgTxt("index out of bounds in TerminalWeight Matrix.");
+        G.add_tweights(node,twPt[cTw+twRows],twPt[cTw+2*twRows]);
+    }
+    
+    for(int cEw=0;cEw<ewRows;cEw++)
+    {
+        //Test for nodes in range
+        int From=(int)ewPt[cEw]-1;
+        int To=(int)ewPt[cEw+ewRows]-1;
+        if(From<0 || From>=nNodes)
+            mexErrMsgTxt("From index out of bounds in Edge Weight Matrix.");
+        if(To<0 || To>=nNodes)
+            mexErrMsgTxt("To index out of bounds in Edge Weight Matrix.");
+
+        G.add_edge(From,To,ewPt[cEw+2*ewRows],ewPt[cEw+3*ewRows]);
+    }
+
+    double flow=G.maxflow();
+    
+    std::vector<int> SourceCut;
+    for(int cNode=0;cNode<nNodes;cNode++)
+    {
+        if(G.what_segment(cNode)== GraphType::SOURCE)
+            SourceCut.push_back(cNode+1);
+    }
+       
+    out[0] = mxCreateDoubleMatrix(SourceCut.size(), 1, mxREAL) ;
+    double* pOut=mxGetPr(out[0]);
+    std::vector<int>::const_iterator Itt(SourceCut.begin());
+    for(;Itt!=SourceCut.end();++Itt)
+    {
+        *pOut=*Itt;
+        ++pOut;
+    }
+    
+    if(nout==2)
+    {
+        out[1] = mxCreateDoubleMatrix(1, 1, mxREAL) ;
+        *mxGetPr(out[1])=flow;
+    }
+    
+}
\ No newline at end of file

Chapter07/Chapter07_mcode/GraphCut/MaxFlow.cpp

+/* maxflow.cpp */
+
+
+#include <stdio.h>
+#include "graph.h"
+#include "instances.inc"
+
+
+/*
+	special constants for node->parent
+*/
+#define TERMINAL ( (arc *) 1 )		/* to terminal */
+#define ORPHAN   ( (arc *) 2 )		/* orphan */
+
+
+#define INFINITE_D ((int)(((unsigned)-1)/2))		/* infinite distance to the terminal */
+
+/***********************************************************************/
+
+/*
+	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).
+*/
+
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	inline void Graph<captype,tcaptype,flowtype>::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
+*/
+template <typename captype, typename tcaptype, typename flowtype> 
+	inline typename Graph<captype,tcaptype,flowtype>::node* Graph<captype,tcaptype,flowtype>::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;
+	}
+}
+
+/***********************************************************************/
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	inline void Graph<captype,tcaptype,flowtype>::set_orphan_front(node *i)
+{
+	nodeptr *np;
+	i -> parent = ORPHAN;
+	np = nodeptr_block -> New();
+	np -> ptr = i;
+	np -> next = orphan_first;
+	orphan_first = np;
+}
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	inline void Graph<captype,tcaptype,flowtype>::set_orphan_rear(node *i)
+{
+	nodeptr *np;
+	i -> parent = ORPHAN;
+	np = nodeptr_block -> New();
+	np -> ptr = i;
+	if (orphan_last) orphan_last -> next = np;
+	else             orphan_first        = np;
+	orphan_last = np;
+	np -> next = NULL;
+}
+
+/***********************************************************************/
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	inline void Graph<captype,tcaptype,flowtype>::add_to_changed_list(node *i)
+{
+	if (changed_list && !i->is_in_changed_list)
+	{
+		node_id* ptr = changed_list->New();
+		*ptr = (node_id)(i - nodes);
+		i->is_in_changed_list = true;
+	}
+}
+
+/***********************************************************************/
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	void Graph<captype,tcaptype,flowtype>::maxflow_init()
+{
+	node *i;
+
+	queue_first[0] = queue_last[0] = NULL;
+	queue_first[1] = queue_last[1] = NULL;
+	orphan_first = NULL;
+
+	TIME = 0;
+
+	for (i=nodes; i<node_last; i++)
+	{
+		i -> next = NULL;
+		i -> is_marked = 0;
+		i -> is_in_changed_list = 0;
+		i -> TS = TIME;
+		if (i->tr_cap > 0)
+		{
+			/* i is connected to the source */
+			i -> is_sink = 0;
+			i -> parent = TERMINAL;
+			set_active(i);
+			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 -> DIST = 1;
+		}
+		else
+		{
+			i -> parent = NULL;
+		}
+	}
+}
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	void Graph<captype,tcaptype,flowtype>::maxflow_reuse_trees_init()
+{
+	node* i;
+	node* j;
+	node* queue = queue_first[1];
+	arc* a;
+	nodeptr* np;
+
+	queue_first[0] = queue_last[0] = NULL;
+	queue_first[1] = queue_last[1] = NULL;
+	orphan_first = orphan_last = NULL;
+
+	TIME ++;
+
+	while ((i=queue))
+	{
+		queue = i->next;
+		if (queue == i) queue = NULL;
+		i->next = NULL;
+		i->is_marked = 0;
+		set_active(i);
+
+		if (i->tr_cap == 0)
+		{
+			if (i->parent) set_orphan_rear(i);
+			continue;
+		}
+
+		if (i->tr_cap > 0)
+		{
+			if (!i->parent || i->is_sink)
+			{
+				i->is_sink = 0;
+				for (a=i->first; a; a=a->next)
+				{
+					j = a->head;
+					if (!j->is_marked)
+					{
+						if (j->parent == a->sister) set_orphan_rear(j);
+						if (j->parent && j->is_sink && a->r_cap > 0) set_active(j);
+					}
+				}
+				add_to_changed_list(i);
+			}
+		}
+		else
+		{
+			if (!i->parent || !i->is_sink)
+			{
+				i->is_sink = 1;
+				for (a=i->first; a; a=a->next)
+				{
+					j = a->head;
+					if (!j->is_marked)
+					{
+						if (j->parent == a->sister) set_orphan_rear(j);
+						if (j->parent && !j->is_sink && a->sister->r_cap > 0) set_active(j);
+					}
+				}
+				add_to_changed_list(i);
+			}
+		}
+		i->parent = TERMINAL;
+		i -> TS = TIME;
+		i -> DIST = 1;
+	}
+
+	//test_consistency();
+
+	/* adoption */
+	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);
+	}
+	/* adoption end */
+
+	//test_consistency();
+}
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	void Graph<captype,tcaptype,flowtype>::augment(arc *middle_arc)
+{
+	node *i;
+	arc *a;
+	tcaptype bottleneck;
+
+
+	/* 1. Finding bottleneck capacity */
+	/* 1a - the source tree */
+	bottleneck = middle_arc -> r_cap;
+	for (i=middle_arc->sister->head; ; i=a->head)
+	{
+		a = i -> parent;
+		if (a == TERMINAL) break;
+		if (bottleneck > a->sister->r_cap) bottleneck = a -> sister -> r_cap;
+	}
+	if (bottleneck > i->tr_cap) bottleneck = i -> tr_cap;
+	/* 1b - the sink tree */
+	for (i=middle_arc->head; ; i=a->head)
+	{
+		a = i -> parent;
+		if (a == TERMINAL) break;
+		if (bottleneck > a->r_cap) bottleneck = a -> r_cap;
+	}
+	if (bottleneck > - i->tr_cap) bottleneck = - i -> tr_cap;
+
+
+	/* 2. Augmenting */
+	/* 2a - the source tree */
+	middle_arc -> sister -> r_cap += bottleneck;
+	middle_arc -> r_cap -= bottleneck;
+	for (i=middle_arc->sister->head; ; i=a->head)
+	{
+		a = i -> parent;
+		if (a == TERMINAL) break;
+		a -> r_cap += bottleneck;
+		a -> sister -> r_cap -= bottleneck;
+		if (!a->sister->r_cap)
+		{
+			set_orphan_front(i); // add i to the beginning of the adoption list
+		}
+	}
+	i -> tr_cap -= bottleneck;
+	if (!i->tr_cap)
+	{
+		set_orphan_front(i); // add i to the beginning of the adoption list
+	}
+	/* 2b - the sink tree */
+	for (i=middle_arc->head; ; i=a->head)
+	{
+		a = i -> parent;
+		if (a == TERMINAL) break;
+		a -> sister -> r_cap += bottleneck;
+		a -> r_cap -= bottleneck;
+		if (!a->r_cap)
+		{
+			set_orphan_front(i); // add i to the beginning of the adoption list
+		}
+	}
+	i -> tr_cap += bottleneck;
+	if (!i->tr_cap)
+	{
+		set_orphan_front(i); // add i to the beginning of the adoption list
+	}
+
+
+	flow += bottleneck;
+}
+
+/***********************************************************************/
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	void Graph<captype,tcaptype,flowtype>::process_source_orphan(node *i)
+{
+	node *j;
+	arc *a0, *a0_min = NULL, *a;
+	int d, d_min = INFINITE_D;
+
+	/* trying to find a new parent */
+	for (a0=i->first; a0; a0=a0->next)
+	if (a0->sister->r_cap)
+	{
+		j = a0 -> head;
+		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; }
+				j = a -> head;
+			}
+			if (d<INFINITE_D) /* j originates from the source - done */
+			{
+				if (d<d_min)
+				{
+					a0_min = a0;
+					d_min = d;
+				}
+				/* set marks along the path */
+				for (j=a0->head; j->TS!=TIME; j=j->parent->head)
+				{
+					j -> TS = TIME;
+					j -> DIST = d --;
+				}
+			}
+		}
+	}
+
+	if (i->parent = a0_min)
+	{
+		i -> TS = TIME;
+		i -> DIST = d_min + 1;
+	}
+	else
+	{
+		/* no parent is found */
+		add_to_changed_list(i);
+
+		/* process neighbors */
+		for (a0=i->first; a0; a0=a0->next)
+		{
+			j = a0 -> head;
+			if (!j->is_sink && (a=j->parent))
+			{
+				if (a0->sister->r_cap) set_active(j);
+				if (a!=TERMINAL && a!=ORPHAN && a->head==i)
+				{
+					set_orphan_rear(j); // add j to the end of the adoption list
+				}
+			}
+		}
+	}
+}
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	void Graph<captype,tcaptype,flowtype>::process_sink_orphan(node *i)
+{
+	node *j;
+	arc *a0, *a0_min = NULL, *a;
+	int d, d_min = INFINITE_D;
+
+	/* trying to find a new parent */
+	for (a0=i->first; a0; a0=a0->next)
+	if (a0->r_cap)
+	{
+		j = a0 -> head;
+		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; }
+				j = a -> head;
+			}
+			if (d<INFINITE_D) /* j originates from the sink - done */
+			{
+				if (d<d_min)
+				{
+					a0_min = a0;
+					d_min = d;
+				}
+				/* set marks along the path */
+				for (j=a0->head; j->TS!=TIME; j=j->parent->head)
+				{
+					j -> TS = TIME;
+					j -> DIST = d --;
+				}
+			}
+		}
+	}
+
+	if (i->parent = a0_min)
+	{
+		i -> TS = TIME;
+		i -> DIST = d_min + 1;
+	}
+	else
+	{
+		/* no parent is found */
+		add_to_changed_list(i);
+
+		/* process neighbors */
+		for (a0=i->first; a0; a0=a0->next)
+		{
+			j = a0 -> head;
+			if (j->is_sink && (a=j->parent))
+			{
+				if (a0->r_cap) set_active(j);
+				if (a!=TERMINAL && a!=ORPHAN && a->head==i)
+				{
+					set_orphan_rear(j); // add j to the end of the adoption list
+				}
+			}
+		}
+	}
+}
+
+/***********************************************************************/
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	flowtype Graph<captype,tcaptype,flowtype>::maxflow(bool reuse_trees, Block<node_id>* _changed_list)
+{
+	node *i, *j, *current_node = NULL;
+	arc *a;
+	nodeptr *np, *np_next;
+
+	if (!nodeptr_block)
+	{
+		nodeptr_block = new DBlock<nodeptr>(NODEPTR_BLOCK_SIZE, error_function);
+	}
+
+	changed_list = _changed_list;
+	if (maxflow_iteration == 0 && reuse_trees) { if (error_function) (*error_function)("reuse_trees cannot be used in the first call to maxflow()!"); exit(1); }
+	if (changed_list && !reuse_trees) { if (error_function) (*error_function)("changed_list cannot be used without reuse_trees!"); exit(1); }
+
+	if (reuse_trees) maxflow_reuse_trees_init();
+	else             maxflow_init();
+
+	// main loop
+	while ( 1 )
+	{
+		// test_consistency(current_node);
+
+		if ((i=current_node))
+		{
+			i -> next = NULL; /* remove active flag */
+			if (!i->parent) i = NULL;
+		}
+		if (!i)
+		{
+			if (!(i = next_active())) break;
+		}
+
+		/* growth */
+		if (!i->is_sink)
+		{
+			/* grow source tree */
+			for (a=i->first; a; a=a->next)
+			if (a->r_cap)
+			{
+				j = a -> head;
+				if (!j->parent)
+				{
+					j -> is_sink = 0;
+					j -> parent = a -> sister;
+					j -> TS = i -> TS;
+					j -> DIST = i -> DIST + 1;
+					set_active(j);
+					add_to_changed_list(j);
+				}
+				else if (j->is_sink) 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 -> sister;
+					j -> TS = i -> TS;
+					j -> DIST = i -> DIST + 1;
+				}
+			}
+		}
+		else
+		{
+			/* grow sink tree */
+			for (a=i->first; a; a=a->next)
+			if (a->sister->r_cap)
+			{
+				j = a -> head;
+				if (!j->parent)
+				{
+					j -> is_sink = 1;
+					j -> parent = a -> sister;
+					j -> TS = i -> TS;
+					j -> DIST = i -> DIST + 1;
+					set_active(j);
+					add_to_changed_list(j);
+				}
+				else if (!j->is_sink) { a = a -> sister; 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 -> sister;
+					j -> TS = i -> TS;
+					j -> DIST = i -> DIST + 1;
+				}
+			}
+		}
+
+		TIME ++;
+
+		if (a)
+		{
+			i -> next = i; /* set active flag */
+			current_node = i;
+
+			/* augmentation */
+			augment(a);
+			/* 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;
+	}
+	// test_consistency();
+
+	if (!reuse_trees || (maxflow_iteration % 64) == 0)
+	{
+		delete nodeptr_block; 
+		nodeptr_block = NULL; 
+	}
+
+	maxflow_iteration ++;
+	return flow;
+}
+
+/***********************************************************************/
+
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	void Graph<captype,tcaptype,flowtype>::test_consistency(node* current_node)
+{
+	node *i;
+	arc *a;
+	int r;
+	int num1 = 0, num2 = 0;
+
+	// test whether all nodes i with i->next!=NULL are indeed in the queue
+	for (i=nodes; i<node_last; i++)
+	{
+		if (i->next || i==current_node) num1 ++;
+	}
+	for (r=0; r<3; r++)
+	{
+		i = (r == 2) ? current_node : queue_first[r];
+		if (i)
+		for ( ; ; i=i->next)
+		{
+			num2 ++;
+			if (i->next == i)
+			{
+				if (r<2) assert(i == queue_last[r]);
+				else     assert(i == current_node);
+				break;
+			}
+		}
+	}
+	assert(num1 == num2);
+
+	for (i=nodes; i<node_last; i++)
+	{
+		// test whether all edges in seach trees are non-saturated
+		if (i->parent == NULL) {}
+		else if (i->parent == ORPHAN) {}
+		else if (i->parent == TERMINAL)
+		{
+			if (!i->is_sink) assert(i->tr_cap > 0);
+			else             assert(i->tr_cap < 0);
+		}
+		else
+		{
+			if (!i->is_sink) assert (i->parent->sister->r_cap > 0);
+			else             assert (i->parent->r_cap > 0);
+		}
+		// test whether passive nodes in search trees have neighbors in
+		// a different tree through non-saturated edges
+		if (i->parent && !i->next)
+		{
+			if (!i->is_sink)
+			{
+				assert(i->tr_cap >= 0);
+				for (a=i->first; a; a=a->next)
+				{
+					if (a->r_cap > 0) assert(a->head->parent && !a->head->is_sink);
+				}
+			}
+			else
+			{
+				assert(i->tr_cap <= 0);
+				for (a=i->first; a; a=a->next)
+				{
+					if (a->sister->r_cap > 0) assert(a->head->parent && a->head->is_sink);
+				}
+			}
+		}
+		// test marking invariants
+		if (i->parent && i->parent!=ORPHAN && i->parent!=TERMINAL)
+		{
+			assert(i->TS <= i->parent->head->TS);
+			if (i->TS == i->parent->head->TS) assert(i->DIST > i->parent->head->DIST);
+		}
+	}
+}
+
+
+	void Inst()
+	{
+		Graph<int,int,int> G(1,1);
+		G.maxflow();
+	}
\ No newline at end of file

Chapter07/Chapter07_mcode/GraphCut/block.h

+/* block.h */
+/*
+	Template classes Block and DBlock
+	Implement adding and deleting items of the same type in blocks.
+
+	If there there are many items then using Block or DBlock
+	is more efficient than using 'new' and 'delete' both in terms
+	of memory and time since
+	(1) On some systems there is some minimum amount of memory
+	    that 'new' can allocate (e.g., 64), so if items are
+	    small that a lot of memory is wasted.
+	(2) 'new' and 'delete' are designed for items of varying size.
+	    If all items has the same size, then an algorithm for
+	    adding and deleting can be made more efficient.
+	(3) All Block and DBlock functions are inline, so there are
+	    no extra function calls.
+
+	Differences between Block and DBlock:
+	(1) DBlock allows both adding and deleting items,
+	    whereas Block allows only adding items.
+	(2) Block has an additional operation of scanning
+	    items added so far (in the order in which they were added).
+	(3) Block allows to allocate several consecutive
+	    items at a time, whereas DBlock can add only a single item.
+
+	Note that no constructors or destructors are called for items.
+
+	Example usage for items of type 'MyType':
+
+	///////////////////////////////////////////////////
+	#include "block.h"
+	#define BLOCK_SIZE 1024
+	typedef struct { int a, b; } MyType;
+	MyType *ptr, *array[10000];
+
+	...
+
+	Block<MyType> *block = new Block<MyType>(BLOCK_SIZE);
+
+	// adding items
+	for (int i=0; i<sizeof(array); i++)
+	{
+		ptr = block -> New();
+		ptr -> a = ptr -> b = rand();
+	}
+
+	// reading items
+	for (ptr=block->ScanFirst(); ptr; ptr=block->ScanNext())
+	{
+		printf("%d %d\n", ptr->a, ptr->b);
+	}
+
+	delete block;
+
+	...
+
+	DBlock<MyType> *dblock = new DBlock<MyType>(BLOCK_SIZE);
+	
+	// adding items
+	for (int i=0; i<sizeof(array); i++)
+	{
+		array[i] = dblock -> New();
+	}
+
+	// deleting items
+	for (int i=0; i<sizeof(array); i+=2)
+	{
+		dblock -> Delete(array[i]);
+	}
+
+	// adding items
+	for (int i=0; i<sizeof(array); i++)
+	{
+		array[i] = dblock -> New();
+	}
+
+	delete dblock;
+
+	///////////////////////////////////////////////////
+
+	Note that DBlock deletes items by marking them as
+	empty (i.e., by adding them to the list of free items),
+	so that this memory could be used for subsequently
+	added items. Thus, at each moment the memory allocated
+	is determined by the maximum number of items allocated
+	simultaneously at earlier moments. All memory is
+	deallocated only when the destructor is called.
+*/
+
+#ifndef __BLOCK_H__
+#define __BLOCK_H__
+
+#include <stdlib.h>
+
+/***********************************************************************/
+/***********************************************************************/
+/***********************************************************************/
+
+template <class Type> class Block
+{
+public:
+	/* Constructor. Arguments are the block size and
+	   (optionally) the pointer to the function which
+	   will be called if allocation failed; the message
+	   passed to this function is "Not enough memory!" */
+	Block(int size, void (*err_function)(char *) = NULL) { first = last = NULL; block_size = size; error_function = err_function; }
+
+	/* Destructor. Deallocates all items added so far */
+	~Block() { while (first) { block *next = first -> next; delete first; first = next; } }
+
+	/* Allocates 'num' consecutive items; returns pointer
+	   to the first item. 'num' cannot be greater than the
+	   block size since items must fit in one block */
+	Type *New(int num = 1)
+	{
+		Type *t;
+
+		if (!last || last->current + num > last->last)
+		{
+			if (last && last->next) last = last -> next;
+			else
+			{
+				block *next = (block *) new char [sizeof(block) + (block_size-1)*sizeof(Type)];
+				if (!next) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
+				if (last) last -> next = next;
+				else first = next;
+				last = next;
+				last -> current = & ( last -> data[0] );
+				last -> last = last -> current + block_size;
+				last -> next = NULL;
+			}
+		}
+
+		t = last -> current;
+		last -> current += num;
+		return t;
+	}
+
+	/* Returns the first item (or NULL, if no items were added) */
+	Type *ScanFirst()
+	{
+		for (scan_current_block=first; scan_current_block; scan_current_block = scan_current_block->next)
+		{
+			scan_current_data = & ( scan_current_block -> data[0] );
+			if (scan_current_data < scan_current_block -> current) return scan_current_data ++;
+		}
+		return NULL;
+	}
+
+	/* Returns the next item (or NULL, if all items have been read)
+	   Can be called only if previous ScanFirst() or ScanNext()
+	   call returned not NULL. */
+	Type *ScanNext()
+	{
+		while (scan_current_data >= scan_current_block -> current)
+		{
+			scan_current_block = scan_current_block -> next;
+			if (!scan_current_block) return NULL;
+			scan_current_data = & ( scan_current_block -> data[0] );
+		}
+		return scan_current_data ++;
+	}
+
+	/* Marks all elements as empty */
+	void Reset()
+	{
+		block *b;
+		if (!first) return;
+		for (b=first; ; b=b->next)
+		{
+			b -> current = & ( b -> data[0] );
+			if (b == last) break;
+		}
+		last = first;
+	}
+
+/***********************************************************************/
+
+private:
+
+	typedef struct block_st
+	{
+		Type					*current, *last;
+		struct block_st			*next;
+		Type					data[1];
+	} block;
+
+	int		block_size;
+	block	*first;
+	block	*last;
+
+	block	*scan_current_block;
+	Type	*scan_current_data;
+
+	void	(*error_function)(char *);
+};
+
+/***********************************************************************/
+/***********************************************************************/
+/***********************************************************************/
+
+template <class Type> class DBlock
+{
+public:
+	/* Constructor. Arguments are the block size and
+	   (optionally) the pointer to the function which
+	   will be called if allocation failed; the message
+	   passed to this function is "Not enough memory!" */
+	DBlock(int size, void (*err_function)(char *) = NULL) { first = NULL; first_free = NULL; block_size = size; error_function = err_function; }
+
+	/* Destructor. Deallocates all items added so far */
+	~DBlock() { while (first) { block *next = first -> next; delete first; first = next; } }
+
+	/* Allocates one item */
+	Type *New()
+	{
+		block_item *item;
+
+		if (!first_free)
+		{
+			block *next = first;
+			first = (block *) new char [sizeof(block) + (block_size-1)*sizeof(block_item)];
+			if (!first) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
+			first_free = & (first -> data[0] );
+			for (item=first_free; item<first_free+block_size-1; item++)
+				item -> next_free = item + 1;
+			item -> next_free = NULL;
+			first -> next = next;
+		}
+
+		item = first_free;
+		first_free = item -> next_free;
+		return (Type *) item;
+	}
+
+	/* Deletes an item allocated previously */
+	void Delete(Type *t)
+	{
+		((block_item *) t) -> next_free = first_free;
+		first_free = (block_item *) t;
+	}
+
+/***********************************************************************/
+
+private:
+
+	typedef union block_item_st
+	{
+		Type			t;
+		block_item_st	*next_free;
+	} block_item;
+
+	typedef struct block_st
+	{
+		struct block_st			*next;
+		block_item				data[1];
+	} block;
+
+	int			block_size;
+	block		*first;
+	block_item	*first_free;
+
+	void	(*error_function)(char *);
+};
+
+
+#endif
+

Chapter07/Chapter07_mcode/GraphCut/compile_mex_functions.m

+% compiling mex functions
+
+mex  GraphCutMex.cpp MaxFlow.cpp graph.cpp 

Chapter07/Chapter07_mcode/GraphCut/graph.cpp

+/* graph.cpp */
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "graph.h"
+#include "instances.inc"
+
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	Graph<captype, tcaptype, flowtype>::Graph(int node_num_max, int edge_num_max, void (*err_function)(char *))
+	: node_num(0),
+	  nodeptr_block(NULL),
+	  error_function(err_function)
+{
+	if (node_num_max < 16) node_num_max = 16;
+	if (edge_num_max < 16) edge_num_max = 16;
+
+	nodes = (node*) malloc(node_num_max*sizeof(node));
+	arcs = (arc*) malloc(2*edge_num_max*sizeof(arc));
+	if (!nodes || !arcs) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
+
+	node_last = nodes;
+	node_max = nodes + node_num_max;
+	arc_last = arcs;
+	arc_max = arcs + 2*edge_num_max;
+
+	maxflow_iteration = 0;
+	flow = 0;
+}
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	Graph<captype,tcaptype,flowtype>::~Graph()
+{
+	if (nodeptr_block) 
+	{ 
+		delete nodeptr_block; 
+		nodeptr_block = NULL; 
+	}
+	free(nodes);
+	free(arcs);
+}
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	void Graph<captype,tcaptype,flowtype>::reset()
+{
+	node_last = nodes;
+	arc_last = arcs;
+	node_num = 0;
+
+	if (nodeptr_block) 
+	{ 
+		delete nodeptr_block; 
+		nodeptr_block = NULL; 
+	}
+
+	maxflow_iteration = 0;
+	flow = 0;
+}
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	void Graph<captype,tcaptype,flowtype>::reallocate_nodes(int num)
+{
+	int node_num_max = (int)(node_max - nodes);
+	node* nodes_old = nodes;
+
+	node_num_max += node_num_max / 2;
+	if (node_num_max < node_num + num) node_num_max = node_num + num;
+	nodes = (node*) realloc(nodes_old, node_num_max*sizeof(node));
+	if (!nodes) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
+
+	node_last = nodes + node_num;
+	node_max = nodes + node_num_max;
+
+	if (nodes != nodes_old)
+	{
+		arc* a;
+		for (a=arcs; a<arc_last; a++)
+		{
+			a->head = (node*) ((char*)a->head + (((char*) nodes) - ((char*) nodes_old)));
+		}
+	}
+}
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	void Graph<captype,tcaptype,flowtype>::reallocate_arcs()
+{
+	int arc_num_max = (int)(arc_max - arcs);
+	int arc_num = (int)(arc_last - arcs);
+	arc* arcs_old = arcs;
+
+	arc_num_max += arc_num_max / 2; if (arc_num_max & 1) arc_num_max ++;
+	arcs = (arc*) realloc(arcs_old, arc_num_max*sizeof(arc));
+	if (!arcs) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
+
+	arc_last = arcs + arc_num;
+	arc_max = arcs + arc_num_max;
+
+	if (arcs != arcs_old)
+	{
+		node* i;
+		arc* a;
+		for (i=nodes; i<node_last; i++)
+		{
+			if (i->first) i->first = (arc*) ((char*)i->first + (((char*) arcs) - ((char*) arcs_old)));
+		}
+		for (a=arcs; a<arc_last; a++)
+		{
+			if (a->next) a->next = (arc*) ((char*)a->next + (((char*) arcs) - ((char*) arcs_old)));
+			a->sister = (arc*) ((char*)a->sister + (((char*) arcs) - ((char*) arcs_old)));
+		}
+	}
+}

Chapter07/Chapter07_mcode/GraphCut/graph.h

+/* graph.h */
+/*
+	This software library implements the maxflow algorithm
+	described in
+
+		"An Experimental Comparison of Min-Cut/Max-Flow Algorithms for Energy Minimization in Vision."
+		Yuri Boykov and Vladimir Kolmogorov.
+		In IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI), 
+		September 2004
+
+	This algorithm was developed by Yuri Boykov and Vladimir Kolmogorov
+	at Siemens Corporate Research. To make it available for public use,
+	it was later reimplemented by Vladimir Kolmogorov based on open publications.
+
+	If you use this software for research purposes, you should cite
+	the aforementioned paper in any resulting publication.
+
+	----------------------------------------------------------------------
+
+	REUSING TREES:
+
+	Starting with version 3.0, there is a also an option of reusing search
+	trees from one maxflow computation to the next, as described in
+
+		"Efficiently Solving Dynamic Markov Random Fields Using Graph Cuts."
+		Pushmeet Kohli and Philip H.S. Torr
+		International Conference on Computer Vision (ICCV), 2005
+
+	If you use this option, you should cite
+	the aforementioned paper in any resulting publication.
+*/
+	
+
+
+/*
+	For description, license, example usage see README.TXT.
+*/
+
+#ifndef __GRAPH_H__
+#define __GRAPH_H__
+
+#include <string.h>
+#include "block.h"
+
+#include <assert.h>
+// NOTE: in UNIX you need to use -DNDEBUG preprocessor option to supress assert's!!!
+
+
+
+// captype: type of edge capacities (excluding t-links)
+// tcaptype: type of t-links (edges between nodes and terminals)
+// flowtype: type of total flow
+//
+// Current instantiations are in instances.inc
+template <typename captype, typename tcaptype, typename flowtype> class Graph
+{
+public:
+	typedef enum
+	{
+		SOURCE	= 0,
+		SINK	= 1
+	} termtype; // terminals 
+	typedef int node_id;
+
+	/////////////////////////////////////////////////////////////////////////
+	//                     BASIC INTERFACE FUNCTIONS                       //
+    //              (should be enough for most applications)               //
+	/////////////////////////////////////////////////////////////////////////
+
+	// Constructor. 
+	// The first argument gives an estimate of the maximum number of nodes that can be added
+	// to the graph, and the second argument is an estimate of the maximum number of edges.
+	// The last (optional) argument is the pointer to the function which will be called 
+	// if an error occurs; an error message is passed to this function. 
+	// If this argument is omitted, exit(1) will be called.
+	//
+	// IMPORTANT: It is possible to add more nodes to the graph than node_num_max 
+	// (and node_num_max can be zero). However, if the count is exceeded, then 
+	// the internal memory is reallocated (increased by 50%) which is expensive. 
+	// Also, temporarily the amount of allocated memory would be more than twice than needed.
+	// Similarly for edges.
+	// If you wish to avoid this overhead, you can download version 2.2, where nodes and edges are stored in blocks.
+	Graph(int node_num_max, int edge_num_max, void (*err_function)(char *) = NULL);
+
+	// Destructor
+	~Graph();
+
+	// Adds node(s) to the graph. By default, one node is added (num=1); then first call returns 0, second call returns 1, and so on. 
+	// If num>1, then several nodes are added, and node_id of the first one is returned.
+	// IMPORTANT: see note about the constructor 
+	node_id add_node(int num = 1);
+
+	// Adds a bidirectional edge between 'i' and 'j' with the weights 'cap' and 'rev_cap'.
+	// IMPORTANT: see note about the constructor 
+	void add_edge(node_id i, node_id j, captype cap, captype rev_cap);
+
+	// Adds new edges 'SOURCE->i' and 'i->SINK' with corresponding weights.
+	// Can be called multiple times for each node.
+	// Weights can be negative.
+	// NOTE: the number of such edges is not counted in edge_num_max.
+	//       No internal memory is allocated by this call.
+	void add_tweights(node_id i, tcaptype cap_source, tcaptype cap_sink);
+
+
+	// Computes the maxflow. Can be called several times.
+	// FOR DESCRIPTION OF reuse_trees, SEE mark_node().
+	// FOR DESCRIPTION OF changed_list, SEE remove_from_changed_list().
+	flowtype maxflow(bool reuse_trees = false, Block<node_id>* changed_list = NULL);
+
+	// After the maxflow is computed, this function returns to which
+	// segment the node 'i' belongs (Graph<captype,tcaptype,flowtype>::SOURCE or Graph<captype,tcaptype,flowtype>::SINK).
+	//
+	// Occasionally there may be several minimum cuts. If a node can be assigned
+	// to both the source and the sink, then default_segm is returned.
+	termtype what_segment(node_id i, termtype default_segm = SOURCE);
+
+
+
+	//////////////////////////////////////////////
+	//       ADVANCED INTERFACE FUNCTIONS       //
+	//      (provide access to the graph)       //
+	//////////////////////////////////////////////
+
+private:
+	struct node;
+	struct arc;
+
+public:
+
+	////////////////////////////
+	// 1. Reallocating graph. //
+	////////////////////////////
+
+	// Removes all nodes and edges. 
+	// After that functions add_node() and add_edge() must be called again. 
+	//
+	// Advantage compared to deleting Graph and allocating it again:
+	// no calls to delete/new (which could be quite slow).
+	//
+	// If the graph structure stays the same, then an alternative
+	// is to go through all nodes/edges and set new residual capacities
+	// (see functions below).
+	void reset();
+
+	////////////////////////////////////////////////////////////////////////////////
+	// 2. Functions for getting pointers to arcs and for reading graph structure. //
+	//    NOTE: adding new arcs may invalidate these pointers (if reallocation    //
+	//    happens). So it's best not to add arcs while reading graph structure.   //
+	////////////////////////////////////////////////////////////////////////////////
+
+	// The following two functions return arcs in the same order that they
+	// were added to the graph. NOTE: for each call add_edge(i,j,cap,cap_rev)
+	// the first arc returned will be i->j, and the second j->i.
+	// If there are no more arcs, then the function can still be called, but
+	// the returned arc_id is undetermined.
+	typedef arc* arc_id;
+	arc_id get_first_arc();
+	arc_id get_next_arc(arc_id a);
+
+	// other functions for reading graph structure
+	int get_node_num() { return node_num; }
+	int get_arc_num() { return (int)(arc_last - arcs); }
+	void get_arc_ends(arc_id a, node_id& i, node_id& j); // returns i,j to that a = i->j
+
+	///////////////////////////////////////////////////
+	// 3. Functions for reading residual capacities. //
+	///////////////////////////////////////////////////
+
+	// returns residual capacity of SOURCE->i minus residual capacity of i->SINK
+	tcaptype get_trcap(node_id i); 
+	// returns residual capacity of arc a
+	captype get_rcap(arc* a);
+
+	/////////////////////////////////////////////////////////////////
+	// 4. Functions for setting residual capacities.               //
+	//    NOTE: If these functions are used, the value of the flow //
+	//    returned by maxflow() will not be valid!                 //
+	/////////////////////////////////////////////////////////////////
+
+	void set_trcap(node_id i, tcaptype trcap); 
+	void set_rcap(arc* a, captype rcap);
+
+	////////////////////////////////////////////////////////////////////
+	// 5. Functions related to reusing trees & list of changed nodes. //
+	////////////////////////////////////////////////////////////////////
+
+	// If flag reuse_trees is true while calling maxflow(), then search trees
+	// are reused from previous maxflow computation. 
+	// In this case before calling maxflow() the user must
+	// specify which parts of the graph have changed by calling mark_node():
+	//   add_tweights(i),set_trcap(i)    => call mark_node(i)
+	//   add_edge(i,j),set_rcap(a)       => call mark_node(i); mark_node(j)
+	//
+	// This option makes sense only if a small part of the graph is changed.
+	// The initialization procedure goes only through marked nodes then.
+	// 
+	// mark_node(i) can either be called before or after graph modification.
+	// Can be called more than once per node, but calls after the first one
+	// do not have any effect.
+	// 
+	// NOTE: 
+	//   - This option cannot be used in the first call to maxflow().
+	//   - It is not necessary to call mark_node() if the change is ``not essential'',
+	//     i.e. sign(trcap) is preserved for a node and zero/nonzero status is preserved for an arc.
+	//   - To check that you marked all necessary nodes, you can call maxflow(false) after calling maxflow(true).
+	//     If everything is correct, the two calls must return the same value of flow. (Useful for debugging).
+	void mark_node(node_id i);
+
+	// If changed_list is not NULL while calling maxflow(), then the algorithm
+	// keeps a list of nodes which could potentially have changed their segmentation label.
+	// Nodes which are not in the list are guaranteed to keep their old segmentation label (SOURCE or SINK).
+	// Example usage:
+	//
+	//		typedef Graph<int,int,int> G;
+	//		G* g = new Graph(nodeNum, edgeNum);
+	//		Block<G::node_id>* changed_list = new Block<G::node_id>(128);
+	//
+	//		... // add nodes and edges
+	//
+	//		g->maxflow(); // first call should be without arguments
+	//		for (int iter=0; iter<10; iter++)
+	//		{
+	//			... // change graph, call mark_node() accordingly
+	//
+	//			g->maxflow(true, changed_list);
+	//			G::node_id* ptr;
+	//			for (ptr=changed_list->ScanFirst(); ptr; ptr=changed_list->ScanNext())
+	//			{
+	//				G::node_id i = *ptr; assert(i>=0 && i<nodeNum);
+	//				g->remove_from_changed_list(i);
+	//				// do something with node i...
+	//				if (g->what_segment(i) == G::SOURCE) { ... }
+	//			}
+	//			changed_list->Reset();
+	//		}
+	//		delete changed_list;
+	//		
+	// NOTE:
+	//  - If changed_list option is used, then reuse_trees must be used as well.
+	//  - In the example above, the user may omit calls g->remove_from_changed_list(i) and changed_list->Reset() in a given iteration.
+	//    Then during the next call to maxflow(true, &changed_list) new nodes will be added to changed_list.
+	//  - If the next call to maxflow() does not use option reuse_trees, then calling remove_from_changed_list()
+	//    is not necessary. ("changed_list->Reset()" or "delete changed_list" should still be called, though).
+	void remove_from_changed_list(node_id i) 
+	{ 
+		assert(i>=0 && i<node_num && nodes[i].is_in_changed_list); 
+		nodes[i].is_in_changed_list = 0;
+	}
+
+
+
+
+
+
+/////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////
+	
+private:
+	// internal variables and functions
+
+	struct node
+	{
+		arc			*first;		// first outcoming arc
+
+		arc			*parent;	// node's parent
+		node		*next;		// pointer to the next active node
+								//   (or to itself if it is the last node in the list)
+		int			TS;			// timestamp showing when DIST was computed
+		int			DIST;		// distance to the terminal
+		int			is_sink : 1;	// flag showing whether the node is in the source or in the sink tree (if parent!=NULL)
+		int			is_marked : 1;	// set by mark_node()
+		int			is_in_changed_list : 1; // set by maxflow if 
+
+		tcaptype	tr_cap;		// if tr_cap > 0 then tr_cap is residual capacity of the arc SOURCE->node
+								// otherwise         -tr_cap is residual capacity of the arc node->SINK 
+
+	};
+
+	struct arc
+	{
+		node		*head;		// node the arc points to
+		arc			*next;		// next arc with the same originating node
+		arc			*sister;	// reverse arc
+
+		captype		r_cap;		// residual capacity
+	};
+
+	struct nodeptr
+	{
+		node    	*ptr;
+		nodeptr		*next;
+	};
+	static const int NODEPTR_BLOCK_SIZE = 128;
+
+	node				*nodes, *node_last, *node_max; // node_last = nodes+node_num, node_max = nodes+node_num_max;
+	arc					*arcs, *arc_last, *arc_max; // arc_last = arcs+2*edge_num, arc_max = arcs+2*edge_num_max;
+
+	int					node_num;
+
+	DBlock<nodeptr>		*nodeptr_block;
+
+	void	(*error_function)(char *);	// this function is called if a error occurs,
+										// with a corresponding error message
+										// (or exit(1) is called if it's NULL)
+
+	flowtype			flow;		// total flow
+
+	// reusing trees & list of changed pixels
+	int					maxflow_iteration; // counter
+	Block<node_id>		*changed_list;
+
+	/////////////////////////////////////////////////////////////////////////
+
+	node				*queue_first[2], *queue_last[2];	// list of active nodes
+	nodeptr				*orphan_first, *orphan_last;		// list of pointers to orphans
+	int					TIME;								// monotonically increasing global counter
+
+	/////////////////////////////////////////////////////////////////////////
+
+	void reallocate_nodes(int num); // num is the number of new nodes
+	void reallocate_arcs();
+
+	// functions for processing active list
+	void set_active(node *i);
+	node *next_active();
+
+	// functions for processing orphans list
+	void set_orphan_front(node* i); // add to the beginning of the list
+	void set_orphan_rear(node* i);  // add to the end of the list
+
+	void add_to_changed_list(node* i);
+
+	void maxflow_init();             // called if reuse_trees == false
+	void maxflow_reuse_trees_init(); // called if reuse_trees == true
+	void augment(arc *middle_arc);
+	void process_source_orphan(node *i);
+	void process_sink_orphan(node *i);
+
+	void test_consistency(node* current_node=NULL); // debug function
+};
+
+
+
+
+
+
+
+
+
+
+
+///////////////////////////////////////
+// Implementation - inline functions //
+///////////////////////////////////////
+
+
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	inline typename Graph<captype,tcaptype,flowtype>::node_id Graph<captype,tcaptype,flowtype>::add_node(int num)
+{
+	assert(num > 0);
+
+	if (node_last + num > node_max) reallocate_nodes(num);
+
+	if (num == 1)
+	{
+		node_last -> first = NULL;
+		node_last -> tr_cap = 0;
+		node_last -> is_marked = 0;
+		node_last -> is_in_changed_list = 0;
+
+		node_last ++;
+		return node_num ++;
+	}
+	else
+	{
+		memset(node_last, 0, num*sizeof(node));
+
+		node_id i = node_num;
+		node_num += num;
+		node_last += num;
+		return i;
+	}
+}
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	inline void Graph<captype,tcaptype,flowtype>::add_tweights(node_id i, tcaptype cap_source, tcaptype cap_sink)
+{
+	assert(i >= 0 && i < node_num);
+
+	tcaptype delta = nodes[i].tr_cap;
+	if (delta > 0) cap_source += delta;
+	else           cap_sink   -= delta;
+	flow += (cap_source < cap_sink) ? cap_source : cap_sink;
+	nodes[i].tr_cap = cap_source - cap_sink;
+}
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	inline void Graph<captype,tcaptype,flowtype>::add_edge(node_id _i, node_id _j, captype cap, captype rev_cap)
+{
+	assert(_i >= 0 && _i < node_num);
+	assert(_j >= 0 && _j < node_num);
+	assert(_i != _j);
+	assert(cap >= 0);
+	assert(rev_cap >= 0);
+
+	if (arc_last == arc_max) reallocate_arcs();
+
+	arc *a = arc_last ++;
+	arc *a_rev = arc_last ++;
+
+	node* i = nodes + _i;
+	node* j = nodes + _j;
+
+	a -> sister = a_rev;
+	a_rev -> sister = a;
+	a -> next = i -> first;
+	i -> first = a;
+	a_rev -> next = j -> first;
+	j -> first = a_rev;
+	a -> head = j;
+	a_rev -> head = i;
+	a -> r_cap = cap;
+	a_rev -> r_cap = rev_cap;
+}
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	inline typename Graph<captype,tcaptype,flowtype>::arc* Graph<captype,tcaptype,flowtype>::get_first_arc()
+{
+	return arcs;
+}
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	inline typename Graph<captype,tcaptype,flowtype>::arc* Graph<captype,tcaptype,flowtype>::get_next_arc(arc* a) 
+{
+	return a + 1; 
+}
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	inline void Graph<captype,tcaptype,flowtype>::get_arc_ends(arc* a, node_id& i, node_id& j)
+{
+	assert(a >= arcs && a < arc_last);
+	i = (node_id) (a->sister->head - nodes);
+	j = (node_id) (a->head - nodes);
+}
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	inline tcaptype Graph<captype,tcaptype,flowtype>::get_trcap(node_id i)
+{
+	assert(i>=0 && i<node_num);
+	return nodes[i].tr_cap;
+}
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	inline captype Graph<captype,tcaptype,flowtype>::get_rcap(arc* a)
+{
+	assert(a >= arcs && a < arc_last);
+	return a->r_cap;
+}
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	inline void Graph<captype,tcaptype,flowtype>::set_trcap(node_id i, tcaptype trcap)
+{
+	assert(i>=0 && i<node_num); 
+	nodes[i].tr_cap = trcap;
+}
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	inline void Graph<captype,tcaptype,flowtype>::set_rcap(arc* a, captype rcap)
+{
+	assert(a >= arcs && a < arc_last);
+	a->r_cap = rcap;
+}
+
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	inline typename Graph<captype,tcaptype,flowtype>::termtype Graph<captype,tcaptype,flowtype>::what_segment(node_id i, termtype default_segm)
+{
+	if (nodes[i].parent)
+	{
+		return (nodes[i].is_sink) ? SINK : SOURCE;
+	}
+	else
+	{
+		return default_segm;
+	}
+}
+
+template <typename captype, typename tcaptype, typename flowtype> 
+	inline void Graph<captype,tcaptype,flowtype>::mark_node(node_id _i)
+{
+	node* i = nodes + _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;
+	}
+	i->is_marked = 1;
+}
+
+
+#endif

Chapter07/Chapter07_mcode/GraphCut/instances.inc

+#include "graph.h"
+
+#ifdef _MSC_VER
+#pragma warning(disable: 4661)
+#endif
+
+// Instantiations: <captype, tcaptype, flowtype>
+// IMPORTANT: 
+//    flowtype should be 'larger' than tcaptype 
+//    tcaptype should be 'larger' than captype
+
+template class Graph<int,int,int>;
+template class Graph<short,int,int>;
+template class Graph<float,float,float>;
+template class Graph<double,double,double>;
+

Chapter07/Chapter07_mcode/functions/grid_cut.m

+function s = grid_cut(s_cost, r_cost, delta_xy, wrap_xy, delta_lu)
+%GRID_CUT detects terrain-like surfaces in a volumetric data.
+%
+% Detected surfaces are terrain-like, i.e. height z is a function of the
+% position (x,y). Only on-surface costs, only in-region costs or both type
+% of costs may be used. When  using in-region costs for multiple surfaces,
+% surfaces are non-intersecting and ordered, so that the first surface has
+% the smallest z values, and the last surface has the largest z values.
+%
+% S = GRID_CUT(S_COST,R_COST,DELTA_XY,WRAP_XY,DELTA_LU)
+%
+% S_COST is a size (X,Y,Z,K) surface cost function. Cost (x,y,z,k) is an
+%       inverse to likelihood of a voxel (x,y,z) containing surface k.
+%       K is the number of surfaces and Z is up. When using only in-region
+%       costs, S_COST should either contain only zeros or be empty.
+% R_COST is a size (X,Y,Z,K+1) region cost function. Cost (x,y,z,k) is an
+%       inverse to likelihood of a voxel (x,y,z) being in a region between
+%       the surface k and k+1, where the first and the last region are
+%       boundedend with the volume boundary. When using only on-surface
+%       costs, R_COST should either contain only zeros or be empty.
+% DELTA_XY is a size (K,2) array of stiffness parameters for x and y. If a
+%       size (1,2) array is given, the same stiffness parameteres are used
+%       for all surfaces. If size (K,1) array is given, the same stiffness
+%       is used for x and y.
+% WRAP_XY is a length 2 array of boolean wrap options for x and y. If not
+%       given defaults to false.
+% DELTA_LU is a size (K-1,2) array of lower and upper surface overlaps, 
+%       [dl du]. The constraint is: dl <= surf_k - surf_k+1 <= du. When 
+%       using in-region costs, the surfaces are non-intersecting and 
+%       ordered which requires dl > 0 and du >= dl. If a size (1,2) array 
+%       is given, the  same overlap parameteres are used for all surfaces. 
+%       If not given, defaults to "no constraint" (dl=-Z, du=Z) when using 
+%       only on-surface costs and to "no overlap" (dl=1, du=Z) otherwise.
+% S is a size (X,Y,K) matrix of z coordinates for K segmented surfaces.
+%
+% Based on:
+%  "Optimal surface segmentation in volumetric images-a graph-theoretic
+%   approach."  Li, Wu, Chen and Sonka. PAMI 2006.
+%  "Incorporation of Regional Information in Optimal 3-D Graph Search with
+%   Application for Intraretinal Layer Segmentation of Optical Coherence
+%   Tomography Images". Haeker, Wu, Abramoff, Kardon and Sonka. IPMI 2007.
+%
+% Author: Vedrana Andersen Dahl, vand@dtu.dk, 2013
+%
+% NOTE: needs Kolmogorovs implementation of the algorithm from
+%		"An Experimental Comparison of Min-Cut/Max-Flow Algorithms for
+%       Energy Minimization in Vision."
+%		Yuri Boykov and Vladimir Kolmogorov. PAMI 2004
+
+% assigning an unique index to each vertex in the first 3D graph (first surface)
+
+if nargin<5
+    delta_lu = [];
+end
+if nargin<4
+    wrap_xy = [];
+end
+[s_cost, r_cost, delta_xy, wrap_xy, delta_lu] = assign_defaults...
+    (s_cost, r_cost, delta_xy, wrap_xy, delta_lu);
+[X,Y,Z,K] = size(s_cost);
+dimension = [X Y Z];
+indices = reshape(1:prod(dimension),dimension);
+base = indices(:,:,1);
+next_surface = prod(dimension); % adding to indices is a shift to next surface
+layers = (0:(K-1))*next_surface;
+
+wrap = [wrap_xy 0]; % we do not wrap in z direction
+
+% EDGES WHICH ARE THE SAME FOR ALL SURFACES:
+% intracolumn arcs pointing down, Equation (2)
+Ea1 = displacement_to_edges(indices, [0 0 -1], wrap); % first surface
+Ea = repmat(Ea1,[K 1]) + kron(layers(:),ones(size(Ea1))); % all surfaces
+Ea = [Ea,ones(size(Ea,1),1)*[inf,0]]; % assigning [inf 0] weight
+% base edges, part of intercolumn arcs, Equation (3)
+erxb1 = displacement_to_edges(base, [1 0 0], wrap);
+eryb1 = displacement_to_edges(base, [0 1 0], wrap);
+Erb1 = [erxb1; eryb1]; % first surface
+Erb = repmat(Erb1,[K 1]) + kron(layers(:),ones(size(Erb1))); % all surfaces
+% assigning [inf inf] weight, base edges are in both directions
+Erb = [Erb,ones(size(Erb,1),1)*[inf,inf]];
+
+% EDGES WHICH DEPEND ON SURFACE STIFFNESS:
+% slanted edges, part of intercolumn arcs, Equation (3)
+Erpm = []; % preallocation length could be computed from X,Y,Z,delta and wrap
+for k = 1:K
+    erxp = displacement_to_edges(indices, [1 0 -delta_xy(k,1)], wrap);
+    erxm = displacement_to_edges(indices, [-1 0 -delta_xy(k,1)], wrap);
+    eryp = displacement_to_edges(indices, [0 1 -delta_xy(k,2)], wrap);
+    erym = displacement_to_edges(indices, [0 -1 -delta_xy(k,2)], wrap);
+    erpm  = [erxp; erxm; eryp; erym] + layers(k);
+    Erpm = [Erpm;erpm];
+end
+% assigning [inf 0] weight, slanted edges are in one direction
+Erpm = [Erpm,ones(size(Erpm,1),1)*[inf,0]];
+Er = [Erpm; Erb]; % all inter edges, for all surfaces
+E = [Ea;Er]; % all intracolumn and intercolumn edges, for all surfaces
+
+% EDGES WHICH DEPEND ON DISTANCE BETWEEN SURFACES:
+% intersurface arc, Equation (4)
+if K>1 % only if we have more than 1 surface
+    Es = []; % preallocation length could be computed
+    for k = 1:K-1
+        esl = displacement_to_edges(indices, [0 0 delta_lu(k,1)], wrap);
+        esl = esl + ones(size(esl,1),1)*layers([k,k+1]);
+        esu = displacement_to_edges(indices, [0 0 -delta_lu(k,2)], wrap);
+        esu = esu + ones(size(esu,1),1)*layers([k+1,k]);
+        es = [esl;esu];
+        Es = [Es;es];
+    end
+    Es = [Es,ones(size(Es,1),1)*[inf,0]]; % assigning [inf 0] weight
+    
+    % intersurface base edges
+    esb = [layers(1:end-1)',layers(2:end)']+1; %first vertex in all surfaces
+    Esb = [esb,ones(K-1,1)*[inf,inf]]; % assigning [inf inf]
+    E = [E;Es;Esb]; % all intracolumn, intercolumn and intersurface edges, for all surfaces
+end
+
+% EDGES WHICH DEPEND ON SURFACE LIKELIHOOD:
+% up to here we do not use s_cost -- consider efficient algorithm if solving
+% multiple problems of the same size
+% vertex s_cost, Equation (1), done simulatniously for all surfaces
+w_on = -1*ones([dimension,K]); % to prevent empty solution, see second half of section 4.1
+w_on(:,:,2:end,:) = double(s_cost(:,:,2:end,:))-double(s_cost(:,:,1:end-1,:));
+% In case of layered surfaces vertices which can't be realized should be removed, 
+% i.e. topmost vertices of lower surface and lowest vertices of higher surface.
+% Instead of removing vertices, I assign inf weight to topmost vertices. 
+% And this only when delta_lu are both positive (both negative can be 
+% avoided by proper ordering of surfaces).
+if K>1 % only if we have more than 1 surface
+    for k=1:K-1 % assigning inf on-surface weights instead of removing vertices
+        if delta_lu(k,1)>0 && delta_lu(k,1)<= delta_lu(k,2) % bit clumsy
+            w_on(:,:,end-delta_lu(k,1)+1:end,k) = inf;                        
+        end    
+    end
+end
+
+% In-region cost, converting to double to avoid problems when using images
+w_in = double(r_cost(:,:,:,1:end-1))-double(r_cost(:,:,:,2:end));
+% There are issues concerning in-region cost e.g. a better way of preventing
+% empty solutions or topmost solutions. 
+w_in(:,:,1,:) = -inf; % preventing empty solution
+
+w = w_on + w_in;
+
+Vp_ix = find(w(:)>=0); % positive vertices, to be connected to sink
+Vm_ix = find(w(:)<0); % negative vertices, to be connected to source
+Es = [Vm_ix, -w(Vm_ix), zeros(length(Vm_ix),1)]; % source edges
+Et = [Vp_ix, zeros(length(Vp_ix),1), w(Vp_ix)]; % sink edges
+Est = [Es;Et]; % all terminal edges
+
+% FINDING GRAPH CUT USING MAGIC
+[Scut, ~] = GraphCutMex(prod(dimension)*K,Est,E);
+
+% retreiving surfaces as the upper envelope of Scut
+S = zeros(X,Y,Z,K);
+S(Scut) = 1;
+s = zeros(X,Y,K);
+for ki = 1:K
+    for x = 1:X
+        for y = 1:Y
+            s(x,y,ki) = find(S(x,y,:,ki),1,'last');
+        end
+    end
+end
+end
+
+% NOTE: when some of the intersurface edges point upwards, i.e. when
+% delta_l is positive, the lowest delta_l rows of a graph for higher
+% surface can not appear in any feasible solution and can be removed.
+% This has NOT been implemented here.
+
+function edges = displacement_to_edges(indices, disp, wrap)
+% indices -- 3D volume of indices
+% disp -- length 3 vector of displacements, given as [x y,z]
+% wrapping -- length 3 vector of boolean options, given as [wrapx wrapy wrapz]
+% edges -- two columns of indices, each line is an edge.
+
+[x_from,x_into] = displace_1D_indices(size(indices,1),disp(1),wrap(1));
+[y_from,y_into] = displace_1D_indices(size(indices,2),disp(2),wrap(2));
+[z_from,z_into] = displace_1D_indices(size(indices,3),disp(3),wrap(3));
+
+indices_from = indices(x_from,y_from,z_from);
+indices_into = indices(x_into,y_into,z_into);
+
+edges = [indices_from(:),indices_into(:)];
+end
+
+function [from,into] = displace_1D_indices(dim,disp,wrap)
+% dim -- length of the indices vector
+% disp -- length of the displacement
+% wrap -- boolean, indicating trimming or wrapping
+
+indices = 1:dim;
+if wrap
+    from = indices([1-min(disp,0):end,1:-min(disp,0)]);
+    into = indices([1+max(disp,0):end,1:max(disp,0)]);
+else
+    from = indices(1-min(disp,0):end-max(disp,0));
+    into = indices(1+max(disp,0):end+min(disp,0));
+end
+end
+
+function [s_cost, r_cost, delta_xy, wrap_xy, delta_lu] = assign_defaults...
+    (s_cost, r_cost, delta_xy, wrap_xy, delta_lu)
+% either cost_s or cost_r need to be given (not empty and not all zeros)
+if isempty(s_cost)
+    s_cost = zeros(size(r_cost)-[0,0,0,1]);
+end
+if isempty(r_cost)
+    r_cost = zeros([size(s_cost,1),size(s_cost,2),size(s_cost,3),...
+        size(s_cost,4)]+[0,0,0,1]); % to allow s_cost to be 3D
+    regional = false;
+elseif all(r_cost(:)==0)
+    regional = false;
+else
+    regional = true;
+end
+
+[X,Y,Z,K] = size(s_cost);
+[Xr,Yr,Zr,Krplus1] = size(r_cost);
+if ~X==Xr || ~Y==Yr || ~Z==Zr || ~(K+1)==Krplus1
+    error('Error using grid_cut. Dimensions of s_cost and r_cost must agree.')
+end
+
+% smoothness constraint may be given once for all surfaces
+if size(delta_xy,1)==1
+    delta_xy = ones(K,1)*delta_xy;
+end
+% smoothness constraint may be given once for both directions
+if size(delta_xy,2)==1
+    delta_xy = delta_xy*ones(1,2);
+end
+
+if isempty(wrap_xy) || (numel(wrap_xy)==1 && wrap_xy==0)
+    wrap_xy = [0 0];
+end
+
+% overlap constraint defaults to no overlap in region case
+% and no constraint in surface case
+if isempty(delta_lu)
+    if regional % no overlap
+        delta_lu = [1 Z];
+    else % no constraint
+        delta_lu = [-Z Z];
+    end
+end
+
+% overlap constraint may be given once for all surface pairs
+if all(size(delta_lu)==[1,2])
+    delta_lu = ones(K-1,1)*delta_lu;
+end
+end
\ No newline at end of file

Chapter07/Chapter07_mcode/in_region_cost_example.m

+clear, close all
+addpath GraphCut functions
+I = imread('../../../../Data/week7/peaks_image.png');
+figure, 
+subplot(1,2,1)
+imagesc(I), colormap gray, axis image
+
+% COST FUNCTION
+% The region in the middle is bright compared to two darker regions. 
+P = permute(I, [2,3,1]); % making sure that up is the third dimension
+region_cost = cat(4, P, (255-P), P);
+
+% GEOMETRIC CONSTRAINS
+delta_xy = 1; % smoothness very constrained, try also 3 to see less smoothness
+wrap_xy = 0; % 0 for terren-like surfaces
+delta_ul = [1 size(I,1)]; % can be very close, but may not overlap
+
+% CUT
+s = grid_cut([],region_cost,delta_xy,wrap_xy,delta_ul);
+
+% VISUALIZATION
+subplot(1,2,2)
+imagesc(I), axis image ij, colormap gray, hold on
+plot(1:size(I,2),permute(s,[1,3,2]),'r', 'LineWidth',2)
+

Chapter07/Chapter07_mcode/on_surface_cost_example.m

+clear, close all
+addpath GraphCut functions
+
+lw = 2;
+
+I = imread('../../../../Data/week7/layers_A.png');
+surface_cost = permute(I,[2,3,1]);
+dim = size(I);
+
+%% input
+figure
+subplot(1,4,1)
+imagesc(I), axis image, colormap gray
+title('input image'), drawnow
+
+%% one line
+delta = 3;
+wrap =  false;
+s = grid_cut(surface_cost, [], delta, wrap);
+subplot(1,4,2)
+imagesc(I), axis image, colormap gray, hold on
+plot(1:dim(2), s, 'r', 'LineWidth', lw)
+title(['delta = ',num2str(delta)]), drawnow
+
+%% a smoother line
+delta = 1;
+wrap =  false;
+s = grid_cut(surface_cost, [], delta, wrap);
+subplot(1,4,3)
+imagesc(I), axis image, colormap gray, hold on
+plot(1:dim(2), s, 'r', 'LineWidth', lw)
+title(['delta = ',num2str(delta)]), drawnow
+
+%% two lines
+costs = cat(4,surface_cost,surface_cost);
+delta = 3;
+wrap = false;
+overlap = [15, size(I,1)];
+
+s = grid_cut(costs, [], delta, wrap, overlap);
+subplot(1,4,4)
+imagesc(I), axis image, colormap gray, hold on
+plot(1:dim(2),permute(s,[1,3,2]),'r', 'LineWidth',lw)
+title('two dark lines'), drawnow
\ No newline at end of file

Chapter07/Chapter07_mcode/w1_unwrapping_image_exercise.m

+% optional exercise 1.1.5 
+I = imread('../../../../Data/week1/dental/slice100.png');
+
+a = 180; % number of angles for unfolding
+angles = (0 : (a-1)) *2*pi/a; % angular coordinate
+
+center = (1+size(I))/2;
+r = 1:min(size(I)/2); % radial coordinate for unwrapping
+
+X = center(1) + r'*cos(angles);
+Y = center(2) + r'*sin(angles);
+
+F = griddedInterpolant(double(I));
+U = F(X,Y);
+
+figure
+subplot(121), imagesc(I), axis image, colormap gray
+subplot(122), imagesc(U), axis image, colormap gray
\ No newline at end of file

Chapter07/in_region_cost_example.py

+import numpy as np 
+import matplotlib.pyplot as plt
+import skimage.io 
+import slgbuilder
+
+
+I = skimage.io.imread('../../../../Data/week7/peaks_image.png').astype(np.int32)
+
+fig, ax = plt.subplots(1,2)
+ax[0].imshow(I, cmap='gray')
+
+
+layers = [slgbuilder.GraphObject(0*I), slgbuilder.GraphObject(0*I)] # no on-surface cost
+helper = slgbuilder.MaxflowBuilder()
+helper.add_objects(layers)
+
+# Addin regional costs, 
+# the region in the middle is bright compared to two darker regions.
+helper.add_layered_region_cost(layers[0], I, 255-I)
+helper.add_layered_region_cost(layers[1], 255-I, I)
+
+# Adding geometric constrains
+helper.add_layered_boundary_cost()
+helper.add_layered_smoothness(delta=1, wrap=False)  
+helper.add_layered_containment(layers[0], layers[1], min_margin=1)
+
+# Cut
+helper.solve()
+segmentations = [helper.what_segments(l).astype(np.int32) for l in layers]
+segmentation_lines = [s.shape[0] - np.argmax(s[::-1,:], axis=0) - 1 for s in segmentations]
+
+# Visualization
+ax[1].imshow(I, cmap='gray')
+for line in segmentation_lines:
+    ax[1].plot(line, 'r')
+
+
+

Chapter07/on_surface_cost_example.py

+import numpy as np 
+import matplotlib.pyplot as plt
+import skimage.io 
+import slgbuilder
+
+#%% input
+I = skimage.io.imread('../../../../Data/week7/layers_A.png').astype(np.int32)
+
+fig, ax = plt.subplots(1,4)
+ax[0].imshow(I, cmap='gray')
+ax[0].set_title('input image')
+
+#%% one line
+delta = 3
+
+layer = slgbuilder.GraphObject(I)
+helper = slgbuilder.MaxflowBuilder()
+helper.add_object(layer)
+helper.add_layered_boundary_cost()
+helper.add_layered_smoothness(delta=delta, wrap=False)
+
+helper.solve()
+segmentation = helper.what_segments(layer)
+segmentation_line = segmentation.shape[0] - np.argmax(segmentation[::-1,:], axis=0) - 1
+
+ax[1].imshow(I, cmap='gray')
+ax[1].plot(segmentation_line, 'r')
+ax[1].set_title(f'delta = {delta}')
+
+
+#%% a smoother line
+delta = 1
+
+layer = slgbuilder.GraphObject(I)
+helper = slgbuilder.MaxflowBuilder()
+helper.add_object(layer)
+helper.add_layered_boundary_cost()
+helper.add_layered_smoothness(delta=delta, wrap=False)
+
+helper.solve()
+segmentation = helper.what_segments(layer)
+segmentation_line = segmentation.shape[0] - np.argmax(segmentation[::-1,:], axis=0) - 1
+
+ax[2].imshow(I, cmap='gray')
+ax[2].plot(segmentation_line, 'r')
+ax[2].set_title(f'delta = {delta}')
+
+
+#%% tow lines
+layers = [slgbuilder.GraphObject(I), slgbuilder.GraphObject(I)]
+delta = 3
+
+helper = slgbuilder.MaxflowBuilder()
+helper.add_objects(layers)
+helper.add_layered_boundary_cost()
+helper.add_layered_smoothness(delta=delta, wrap=False)  
+helper.add_layered_containment(layers[0], layers[1], min_margin=15)
+
+helper.solve()
+segmentations = [helper.what_segments(l).astype(np.int32) for l in layers]
+segmentation_lines = [s.shape[0] - np.argmax(s[::-1,:], axis=0) - 1 for s in segmentations]
+
+ax[3].imshow(I, cmap='gray')
+for line in segmentation_lines:
+    ax[3].plot(line, 'r')
+ax[3].set_title('two dark lines')

Chapter07/w1_unwrapping_image_exercise.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Feb 27 09:14:57 2020
+
+@author: vand
+"""
+
+# optional exercise 1.1.5 
+
+import skimage.io
+import numpy as np
+import scipy.interpolate
+import matplotlib.pyplot as plt
+
+I = skimage.io.imread('../../../../Data/week1/dental/slice100.png')
+
+
+a = 180 # number of angles for unfolding
+angles = np.arange(a)*2*np.pi/a # angular coordinate
+
+center = (np.array(I.shape)-1)/2
+r = int(min(I.shape)/2)
+radii = np.arange(r) + 1 #radial coordinate for unwrapping
+
+X = center[0] + np.outer(radii,np.cos(angles))
+Y = center[1] + np.outer(radii,np.sin(angles))
+
+F = scipy.interpolate.interp2d(np.arange(I.shape[0]), np.arange(I.shape[1]), I)
+U = np.array([F(p[0],p[1]) for p in np.c_[Y.ravel(),X.ravel()]])
+U = U.reshape((r,a)).astype(np.uint8)
+
+fig, ax = plt.subplots(1,2)
+ax[0].imshow(I, cmap='gray')
+ax[1].imshow(U, cmap='gray')