diff --git a/Chapter05/Chapter05_mcode/GraphCut/GraphCutMex.cpp b/Chapter05/Chapter05_mcode/GraphCut/GraphCutMex.cpp new file mode 100755 index 0000000000000000000000000000000000000000..0806d9e5510e4246601fcabe5fd5a19f950afccd --- /dev/null +++ b/Chapter05/Chapter05_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/Chapter05/Chapter05_mcode/GraphCut/GraphCutMex.mexa64 b/Chapter05/Chapter05_mcode/GraphCut/GraphCutMex.mexa64 new file mode 100755 index 0000000000000000000000000000000000000000..a9df40a4455e24c1d68dffa4196f96e4cd9f291e Binary files /dev/null and b/Chapter05/Chapter05_mcode/GraphCut/GraphCutMex.mexa64 differ diff --git a/Chapter05/Chapter05_mcode/GraphCut/GraphCutMex.mexmaci64 b/Chapter05/Chapter05_mcode/GraphCut/GraphCutMex.mexmaci64 new file mode 100755 index 0000000000000000000000000000000000000000..75cdb2c698cc1256beeb3c1826286d56f9dceda8 Binary files /dev/null and b/Chapter05/Chapter05_mcode/GraphCut/GraphCutMex.mexmaci64 differ diff --git a/Chapter05/Chapter05_mcode/GraphCut/GraphCutMex.mexw64 b/Chapter05/Chapter05_mcode/GraphCut/GraphCutMex.mexw64 new file mode 100755 index 0000000000000000000000000000000000000000..d5b12fb8fa67766a0392f5226a075f06885a5b75 Binary files /dev/null and b/Chapter05/Chapter05_mcode/GraphCut/GraphCutMex.mexw64 differ diff --git a/Chapter05/Chapter05_mcode/GraphCut/MaxFlow.cpp b/Chapter05/Chapter05_mcode/GraphCut/MaxFlow.cpp new file mode 100755 index 0000000000000000000000000000000000000000..b91856eed8105b0b053250396bd529656e996bc3 --- /dev/null +++ b/Chapter05/Chapter05_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/Chapter05/Chapter05_mcode/GraphCut/block.h b/Chapter05/Chapter05_mcode/GraphCut/block.h new file mode 100755 index 0000000000000000000000000000000000000000..14346efc32d94f7b3ebad384d6c166fb173c21d7 --- /dev/null +++ b/Chapter05/Chapter05_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/Chapter05/Chapter05_mcode/GraphCut/compile_mex_functions.m b/Chapter05/Chapter05_mcode/GraphCut/compile_mex_functions.m new file mode 100755 index 0000000000000000000000000000000000000000..7bcfafe6b146e3b0fda7d00638e2160479ef6959 --- /dev/null +++ b/Chapter05/Chapter05_mcode/GraphCut/compile_mex_functions.m @@ -0,0 +1,3 @@ +% compiling mex functions + +mex GraphCutMex.cpp MaxFlow.cpp graph.cpp diff --git a/Chapter05/Chapter05_mcode/GraphCut/graph.cpp b/Chapter05/Chapter05_mcode/GraphCut/graph.cpp new file mode 100755 index 0000000000000000000000000000000000000000..8ba3d5a7f1c9a173fc09626abf348db9c43255e6 --- /dev/null +++ b/Chapter05/Chapter05_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/Chapter05/Chapter05_mcode/GraphCut/graph.h b/Chapter05/Chapter05_mcode/GraphCut/graph.h new file mode 100755 index 0000000000000000000000000000000000000000..b45dd369ffbe244a365b8b81392e63479ffd424b --- /dev/null +++ b/Chapter05/Chapter05_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/Chapter05/Chapter05_mcode/GraphCut/instances.inc b/Chapter05/Chapter05_mcode/GraphCut/instances.inc new file mode 100755 index 0000000000000000000000000000000000000000..0bcbd57d426ea2143299bef0407ce4a91744fcce --- /dev/null +++ b/Chapter05/Chapter05_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/Chapter05/Chapter05_mcode/circles_modeling_EMPTY.m b/Chapter05/Chapter05_mcode/circles_modeling_EMPTY.m new file mode 100644 index 0000000000000000000000000000000000000000..b1169bf5deda5e47ceb78020987d20bf2ea50e43 --- /dev/null +++ b/Chapter05/Chapter05_mcode/circles_modeling_EMPTY.m @@ -0,0 +1,66 @@ +clear +close all +addpath('functions') + +%% loading data +path = '../../../../Data/week5/'; +I = double(imread([path, 'noisy_circles.png'])); + +segmentations = {}; % cell-array where I'll place diifferent segmentations +GT = double(imread([path, 'noise_free_circles.png'])); + +[mu, ~, S_gt] = unique(GT); +S_gt = reshape(S_gt, size(I)); +segmentations{1} = S_gt; + +%% finding some configurations (segmentations) using conventional methods +S_t = ones(size(I))+(I>100)+(I>160); % thresholded +segmentations{2} = S_t; + +D_s = imfilter(I, fspecial('gaussian', 7 ,1), 'replicate'); +S_g = ones(size(I))+(D_s>100)+(D_s>160); % gaussian fitered and thresholded +segmentations{3} = S_g; + +D_m = medfilt2(I,[5 5],'symmetric'); +S_m = ones(size(I))+(D_m>100)+(D_m>160); % median fitered and thresholded +segmentations{4} = S_m; + +%% visualization +figure +imagesc(I, [0,255]), axis image, colormap gray + +figure +N = length(segmentations); +beta = 100; + +for i = 1:N + + subplot(3,N,i) + imagesc(segmentations{i}), axis image + [V1,V2] = segmentation_energy(segmentations{i}, I, mu, beta); + title({['likelihood: ', num2str(round(V1))],['prior: ', num2str(V2)],... + ['posterior: ', num2str(round(V1+V2))]}) + + subplot(3, N, i+N) + segmentation_histogram(I, segmentations{i}) + xlabel('intensity'), ylabel('count'), %legend('all','1','2','3') + xlim([0,255]), box on + + subplot(3, N, i+2*N) + err = S_gt - segmentations{i}; + imagesc(err, [-2,2]), axis image + title(['error: ',num2str(sum(err(:)>0))]) + colormap(gca,[0,0,0.5; 0,0.25,1; 1,1,1; 1,0,0; 0.5,0,0]) +end + +function [U1,U2] = segmentation_energy(S, I, mu, beta) +% TODO -- add your code here + +% likelihood energy +U1 = 0; + +% prior energy +U2 = 0; + +end + diff --git a/Chapter05/Chapter05_mcode/gender_labeling.m b/Chapter05/Chapter05_mcode/gender_labeling.m new file mode 100755 index 0000000000000000000000000000000000000000..047d991b1ce95a3710ea8c286eb99130978e9c57 --- /dev/null +++ b/Chapter05/Chapter05_mcode/gender_labeling.m @@ -0,0 +1,24 @@ +clear, close all +addpath GraphCut + +d = [179 174 182 162 175 165]; % heights (data) +mu = [181, 165]; % means of two classes +beta = 100; % weight of the prior term +w_s = (d(:)-mu(1)).^2; % source weight +w_t = (d(:)-mu(2)).^2; % sink weights +N = numel(d); % number of graph nodes +indices = (1:N)'; % an index for each person + +% terminal and internal edge matrix +E_terminal = [indices,[w_s,w_t]]; +E_internal = [indices(1:end-1),indices(2:end),beta*ones(N-1,2)]; + +[Scut, flow] = GraphCutMex(N,E_terminal,E_internal); % here it happens +disp(['Maximum flow: ',num2str(flow)]) + +% displaying results +S = 'MMMMMM'; +S(Scut) = 'F'; +for i=1:6 + disp(['Person ',num2str(i),' is estimated as ',S(i)]) +end diff --git a/Chapter05/Chapter05_mcode/multilabel_MRF.m b/Chapter05/Chapter05_mcode/multilabel_MRF.m new file mode 100644 index 0000000000000000000000000000000000000000..ab3d4538db2388464fe2164cc64e8b4c4c69d3b7 --- /dev/null +++ b/Chapter05/Chapter05_mcode/multilabel_MRF.m @@ -0,0 +1,91 @@ +function [S, iter] = multilabel_MRF(U, dim, beta, max_iter) +%MULTILABEL_MRF Iterative alpha expansion. +% [S, iter] = MULTILABEL_MRF(U, dim, beta, max_iter) finds a solution for +% the Potts model, with the 2-clique potentials in {0,beta}, and the +% 1-clique potentials for each vertex and label given in input U. +% +% U is a N-times-K matrix with energy of each pixel/voxel for each label, +% where N is a total number of pixels/voxels and K is a number of labels. +% dim is a vector indicating the dimensions (size) of problem. It must be +% such that prod(dim)=N. +% beta is a scalar weight of the neighorhood edges (2-clique potental for +% different labels). +% max_iter is a maximum number of alpha expansion cycles (a cycle includes +% K alpha expansion steps), defaults to 10. +% S is a length N vector of segmentation labels form {1,...,K}. S is of +% type uint8, so up til 255 labels are suported. +% iter is a number of actually performed cycles. +% +% Using Kolmorogov+Zabih PAMI 2004 formulation (with directed edges) of +% alpha expansion from Boykov et al. PAMI 2001 (with auxiliary vertices), +% but with source being not-alpha and sink as alpha (opposite to Boykov +% formulation). +% +% vand@dtu.dk, 2014, fixed small bug 2021 + +if nargin<4 + max_iter = 10; +end + +K = size(U,2); % number of classes +[~,S] = min(U,[],2); % initial segmentation +S = uint8(S); % assumes no more than 256 labels + +% preparing for n-links, 6-neigborhood of each vertex (4-neigbhorhood for 2D) +indices = reshape(1:prod(dim),dim); +edge_x = indices(1:end-1,:,:); +edge_y = indices(:,1:end-1,:); +edge_z = indices(:,:,1:end-1); % empty for 2D segmentation +edge_n = zeros(numel(edge_x)+numel(edge_y)+numel(edge_z),4); +edge_n(:,[1,2]) = [edge_x(:),edge_x(:)+1; edge_y(:),edge_y(:)+dim(1); ... + edge_z(:),edge_z(:)+dim(1)*dim(2)]; + +% preparing for t-links +edge_t = zeros(prod(dim),3); +edge_t(:,1) = indices(:); + +iter = 0; +changed = true; + +while changed && (iter<max_iter) + + S_it = S; +% for alpha = randperm(K) % expanding classes in random order + for alpha = 1:K + + % Each pair with different labels requires an adjustment in a + % n-link (-beta) and an adjustment in a t-lint (+beta), so that + % only configuration with 0 potential is (alpha,alpha), i.e. + % (relabeled,relabeled). Applies also to pairs with alpha, as alpha + % is always relabeled due to inf, so only other label matters. + edge_n(:,[3,4]) = beta; + diff_pairs = S(edge_n(:,1))~=S(edge_n(:,2)); % neighboors with different labels + edge_n(diff_pairs,4) = 0; % adjusting n-link (minus beta) + + % Finding vertices that need to have t-link adjusted, + % some should be adjusted multiple times, i.e. E(diff_pairs,1) is a + % non-unique list which gets transformed to vector of occurences. + [pos,temp,~] = unique(sort(edge_n(diff_pairs,1))); + from_pair = zeros(prod(dim),1); + % So much adjusting each vertex needs: + from_pair(pos) = [temp(2:end); numel(temp)+1] - temp; + + edge_t(:,3) = U(:,alpha) ; + % First term in edge_t(:,2) is equivalent to U(sub2lin(indices(:),S(:))) + % i.e. U(current label) while the second term is adjustment (+beta). + edge_t(:,2) = U((double(S(:))-1)*prod(dim)+indices(:))+ beta*from_pair; + edge_t(indices(S==alpha),2) = inf; % to force re-labeling with alpha + + % solving + Scut = GraphCutMex(prod(dim),edge_t,edge_n); + S(Scut) = alpha; % source is not-alpha, source set changes to alpha + + end + changed = any(S(:)~=S_it(:)); % change in labeling + iter = iter+1; +end +if changed + warning(['Alpha expansion not converged after ',num2str(max_iter),... + ' iterations.']) +end +S = reshape(S,dim); \ No newline at end of file