diff --git a/Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.cpp b/Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.cpp new file mode 100755 index 0000000000000000000000000000000000000000..0806d9e5510e4246601fcabe5fd5a19f950afccd --- /dev/null +++ b/Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.cpp @@ -0,0 +1,87 @@ +#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 diff --git a/Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.mexa64 b/Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.mexa64 new file mode 100755 index 0000000000000000000000000000000000000000..a9df40a4455e24c1d68dffa4196f96e4cd9f291e Binary files /dev/null and b/Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.mexa64 differ diff --git a/Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.mexmaci64 b/Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.mexmaci64 new file mode 100755 index 0000000000000000000000000000000000000000..75cdb2c698cc1256beeb3c1826286d56f9dceda8 Binary files /dev/null and b/Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.mexmaci64 differ diff --git a/Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.mexw64 b/Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.mexw64 new file mode 100755 index 0000000000000000000000000000000000000000..d5b12fb8fa67766a0392f5226a075f06885a5b75 Binary files /dev/null and b/Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.mexw64 differ diff --git a/Chapter07/Chapter07_mcode/GraphCut/MaxFlow.cpp b/Chapter07/Chapter07_mcode/GraphCut/MaxFlow.cpp new file mode 100755 index 0000000000000000000000000000000000000000..b91856eed8105b0b053250396bd529656e996bc3 --- /dev/null +++ b/Chapter07/Chapter07_mcode/GraphCut/MaxFlow.cpp @@ -0,0 +1,690 @@ +/* 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 diff --git a/Chapter07/Chapter07_mcode/GraphCut/block.h b/Chapter07/Chapter07_mcode/GraphCut/block.h new file mode 100755 index 0000000000000000000000000000000000000000..14346efc32d94f7b3ebad384d6c166fb173c21d7 --- /dev/null +++ b/Chapter07/Chapter07_mcode/GraphCut/block.h @@ -0,0 +1,268 @@ +/* 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 + diff --git a/Chapter07/Chapter07_mcode/GraphCut/compile_mex_functions.m b/Chapter07/Chapter07_mcode/GraphCut/compile_mex_functions.m new file mode 100755 index 0000000000000000000000000000000000000000..7bcfafe6b146e3b0fda7d00638e2160479ef6959 --- /dev/null +++ b/Chapter07/Chapter07_mcode/GraphCut/compile_mex_functions.m @@ -0,0 +1,3 @@ +% compiling mex functions + +mex GraphCutMex.cpp MaxFlow.cpp graph.cpp diff --git a/Chapter07/Chapter07_mcode/GraphCut/graph.cpp b/Chapter07/Chapter07_mcode/GraphCut/graph.cpp new file mode 100755 index 0000000000000000000000000000000000000000..8ba3d5a7f1c9a173fc09626abf348db9c43255e6 --- /dev/null +++ b/Chapter07/Chapter07_mcode/GraphCut/graph.cpp @@ -0,0 +1,114 @@ +/* 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))); + } + } +} diff --git a/Chapter07/Chapter07_mcode/GraphCut/graph.h b/Chapter07/Chapter07_mcode/GraphCut/graph.h new file mode 100755 index 0000000000000000000000000000000000000000..b45dd369ffbe244a365b8b81392e63479ffd424b --- /dev/null +++ b/Chapter07/Chapter07_mcode/GraphCut/graph.h @@ -0,0 +1,506 @@ +/* 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 diff --git a/Chapter07/Chapter07_mcode/GraphCut/instances.inc b/Chapter07/Chapter07_mcode/GraphCut/instances.inc new file mode 100755 index 0000000000000000000000000000000000000000..0bcbd57d426ea2143299bef0407ce4a91744fcce --- /dev/null +++ b/Chapter07/Chapter07_mcode/GraphCut/instances.inc @@ -0,0 +1,16 @@ +#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>; + diff --git a/Chapter07/Chapter07_mcode/functions/grid_cut.m b/Chapter07/Chapter07_mcode/functions/grid_cut.m new file mode 100755 index 0000000000000000000000000000000000000000..cabd67a756b0db2d03e19a167f7379ae07eabb4a --- /dev/null +++ b/Chapter07/Chapter07_mcode/functions/grid_cut.m @@ -0,0 +1,252 @@ +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 diff --git a/Chapter07/Chapter07_mcode/in_region_cost_example.m b/Chapter07/Chapter07_mcode/in_region_cost_example.m new file mode 100644 index 0000000000000000000000000000000000000000..d5f725ae3192ae25996d00124c8e617dc6a32748 --- /dev/null +++ b/Chapter07/Chapter07_mcode/in_region_cost_example.m @@ -0,0 +1,25 @@ +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) + diff --git a/Chapter07/Chapter07_mcode/on_surface_cost_example.m b/Chapter07/Chapter07_mcode/on_surface_cost_example.m new file mode 100755 index 0000000000000000000000000000000000000000..fe8e858890a3dc9cf6c36b67c69c235cc6618e8b --- /dev/null +++ b/Chapter07/Chapter07_mcode/on_surface_cost_example.m @@ -0,0 +1,44 @@ +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 diff --git a/Chapter07/Chapter07_mcode/w1_unwrapping_image_exercise.m b/Chapter07/Chapter07_mcode/w1_unwrapping_image_exercise.m new file mode 100644 index 0000000000000000000000000000000000000000..9692cd5f8944ecbfd218465265d553eae0301466 --- /dev/null +++ b/Chapter07/Chapter07_mcode/w1_unwrapping_image_exercise.m @@ -0,0 +1,18 @@ +% 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 diff --git a/Chapter07/in_region_cost_example.py b/Chapter07/in_region_cost_example.py new file mode 100755 index 0000000000000000000000000000000000000000..c1ac7ce163472686335b664349e2d3d57e42d08d --- /dev/null +++ b/Chapter07/in_region_cost_example.py @@ -0,0 +1,38 @@ +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') + + + diff --git a/Chapter07/on_surface_cost_example.py b/Chapter07/on_surface_cost_example.py new file mode 100755 index 0000000000000000000000000000000000000000..eeab4b9b6b4776db4c5ab83cfbe45dae7f18f313 --- /dev/null +++ b/Chapter07/on_surface_cost_example.py @@ -0,0 +1,66 @@ +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') diff --git a/Chapter07/w1_unwrapping_image_exercise.py b/Chapter07/w1_unwrapping_image_exercise.py new file mode 100644 index 0000000000000000000000000000000000000000..c67fffe91ef1310c56cac640045228192d08eda3 --- /dev/null +++ b/Chapter07/w1_unwrapping_image_exercise.py @@ -0,0 +1,35 @@ +#!/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')