add the functions and modified the codes

added the relevant functions and modified the codes accordingly

Demo.m

 close all
 clear
 
-addpath(genpath('../Functions'))
+addpath(genpath('./Functions'))
 
 %% Defaults
-Data_path = '..\..\..\..\..\..\LINX FP08.001\Activity 7 - bending - visualizations and tiff stacks\mat_file\';
+% if you know the path otherwise choose the path using the 
+answer = questdlg('Write data path or browse to the data folder?', 'Data Path','Write','Browse','Browse');
+switch answer
+    case 'Write'
+        Data_path = '..\..\..\HCP Anywhere\LINX FP08.001\Activity 7 - bending - visualizations and tiff stacks\mat_file\';
+    case 'Browse'
+        Data_path = uigetdir;
+        Data_path = fullfile(Data_path,'\');
+        fprintf('Data Path: %s \n', Data_path)
+end
+
 samples = {'out_out_CD','out_out_MD','in_in_CD','in_in_MD'};
 degree = {'0','45','90','180'};
 layer_name = {'outer','inner'};
@@ -12,12 +22,17 @@ layer_name = {'outer','inner'};
 
 %% extract poitns for the pipeline
 clc
-flag_points = 0; % 1 if no point has been extracted and first time use of the code
+flag_points = 1; % 1 if no point has been extracted and first time use of the code
+flag_save = 0; % 0: to not save the results; 1: to save the results
+
+answer = inputdlg('Contrast Enhansed value:','Contrast',[1 50],{'.7 1.5'}); % this is default based on the data scanned at DTU
+ce = str2num(answer{1});
+
 
 if flag_points
     for i = 1
         sample = samples{i};
-        for j = 1:4
+        for j = 1
             dg = degree{j};
 
             %% load data
@@ -35,7 +50,7 @@ clearvars -except Data_path samples degree layer_name
 clc
 answer = questdlg('Do you want to enhance the contrast? ','Contrast:','No','Yes','No');
 
-for i = 1
+for i = 1:4
     sample = samples{i};
     for j = 1:4
         %% load data
@@ -67,40 +82,48 @@ for i = 1
                         range = 20:-1:-20;
                 end
             end
+            flag_save = 1;
+            flag_fig = 0; % if you want to check the results
             LayerDetection
         end
         
          %% measure
         disp('Measure ...')
+        flag_save = 1;
         PrepMeasure
         
          %% registration
         disp('Registration ...')
+        flag_save = 1;
         PrepRegistration
         
     end
 end
 
 %% visualize the measure
-for i = 1%:4
+for i = 1:4
+    flag_save = 1;
     Measure
 end
 close all
 
 %% apply the registration
-for i = 1%:4
+for i = 1:4
+    flag_save = 1;
     Registration
 end
 close all
 
 %% apply the alignment and transformation
-for i = 1%:4
+for i = 1:4
+    flag_save = 1;
+    flag_fig = 1;
     Transformation
 end
 close all
 
 %% convert the meshes
-for i = 1%:4
+for i = 1:4
     Export2inp
 end
 close all

Export2inp.m

@@ -3,8 +3,6 @@
 samples = {'out_out_CD','out_out_MD','in_in_CD','in_in_MD'};
 degree = {'0','45','90','180'};
 
-flag_save = 1;
-
 load_folder_layers = 'Results\layer\transfer\align_0dg_rest\';
 
 save_folder =  'Results\export\';
@@ -29,7 +27,7 @@ end
           reshape(permute(aligned_inner_dg0(:,:,3,:),[2,1,3,4]),[sz0(1)*sz0(2)*sz0(4),1])];
     no = [(1:size(no,1))',no];
     
-    disp('Write mesh ...')
+    disp('Write mesh inner 0 ...')
     parfeval(@write_mesh,0,no,ele,[save_folder,'inner_',sample,'_0']);
     
     no = [reshape(permute(aligned_outer_dg0(:,:,1,:),[2,1,3,4]),[sz0(1)*sz0(2)*sz0(4),1]), ...
@@ -37,7 +35,7 @@ end
           reshape(permute(aligned_outer_dg0(:,:,3,:),[2,1,3,4]),[sz0(1)*sz0(2)*sz0(4),1])];
     no = [(1:size(no,1))',no];
 
-    disp('Write mesh ...')
+    disp('Write mesh outer 0 ...')
     parfeval(@write_mesh,0,no,ele,[save_folder,'outer_',sample,'_0']);
     
     % 45
@@ -46,7 +44,7 @@ end
           reshape(permute(aligned_inner_dg45(:,:,3,:),[2,1,3,4]),[sz0(1)*sz0(2)*sz0(4),1])];
     no = [(1:size(no,1))',no];
 
-    disp('Write mesh ...')
+    disp('Write mesh inner 45 ...')
     parfeval(@write_mesh,0,no,ele,[save_folder,'inner_',sample,'_45']);
     
     no = [reshape(permute(aligned_outer_dg45(:,:,1,:),[2,1,3,4]),[sz0(1)*sz0(2)*sz0(4),1]), ...
@@ -54,7 +52,7 @@ end
           reshape(permute(aligned_outer_dg45(:,:,3,:),[2,1,3,4]),[sz0(1)*sz0(2)*sz0(4),1])];
     no = [(1:size(no,1))',no];
 
-    disp('Write mesh ...')
+    disp('Write mesh outer 45 ...')
     parfeval(@write_mesh,0,no,ele,[save_folder,'outer_',sample,'_45']);
     
     % 90
@@ -63,7 +61,7 @@ end
           reshape(permute(aligned_inner_dg90(:,:,3,:),[2,1,3,4]),[sz0(1)*sz0(2)*sz0(4),1])];
     no = [(1:size(no,1))',no];
 
-    disp('Write mesh ...')
+    disp('Write mesh inner 90 ...')
     parfeval(@write_mesh,0,no,ele,[save_folder,'inner_',sample,'_90']);
     
     no = [reshape(permute(aligned_outer_dg90(:,:,1,:),[2,1,3,4]),[sz0(1)*sz0(2)*sz0(4),1]), ...
@@ -71,7 +69,7 @@ end
           reshape(permute(aligned_outer_dg90(:,:,3,:),[2,1,3,4]),[sz0(1)*sz0(2)*sz0(4),1])];
     no = [(1:size(no,1))',no];
 
-    disp('Write mesh ...')
+    disp('Write mesh outer 90 ...')
     parfeval(@write_mesh,0,no,ele,[save_folder,'outer_',sample,'_90']);
     
     % 180
@@ -80,7 +78,7 @@ end
           reshape(permute(aligned_inner_dg180(:,:,3,:),[2,1,3,4]),[sz0(1)*sz0(2)*sz0(4),1])];
     no = [(1:size(no,1))',no];
 
-    disp('Write mesh ...')
+    disp('Write mesh inner 180 ...')
     parfeval(@write_mesh,0,no,ele,[save_folder,'inner_',sample,'_180']);
     
     no = [reshape(permute(aligned_outer_dg180(:,:,1,:),[2,1,3,4]),[sz0(1)*sz0(2)*sz0(4),1]), ...
@@ -88,7 +86,7 @@ end
           reshape(permute(aligned_outer_dg180(:,:,3,:),[2,1,3,4]),[sz0(1)*sz0(2)*sz0(4),1])];
     no = [(1:size(no,1))',no];
 
-    disp('Write mesh ...')
+    disp('Write mesh outer 180 s...')
     parfeval(@write_mesh,0,no,ele,[save_folder,'outer_',sample,'_180']);
     
 %end
\ No newline at end of file

Functions/functions_GraphCut/GraphCutMex.cpp

+#include <math.h>
+#include "mex.h"
+#include "graph.h"
+#include <vector>
+
+void
+mexFunction(int nout, mxArray *out[], 
+            int nin, const mxArray *in[])
+{
+    if (nin != 3) 
+        mexErrMsgTxt("Three arguments are required (nNodes,TerminalWeights,EdgeWeights)") ;
+    if (nout > 2) 
+        mexErrMsgTxt("Too many output arguments.");
+  
+
+    int nNodes = (int) *mxGetPr(in[0]);
+    
+
+    const int* twDim = mxGetDimensions(in[1]) ;
+	// mexPrintf("Dim is:  %d\n", twDim[0]);
+	// mexPrintf("Dim is:  %d\n", twDim[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]);
+		//mexPrintf("The loop for TerminalWeight is:  %d\n", cTw);
+    }
+
+	// mexPrintf("ewRows is:  %d\n", ewRows);
+    // mexEvalString("drawnow");
+
+    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]);
+		//mexPrintf("The loop for Edge Weight is:  %d\n", cEw);
+    }
+	
+    double flow=G.maxflow();
+    
+    std::vector<int> SourceCut;
+
+	// mexPrintf("nNodes is:  %d\n", nNodes);
+	// mexEvalString("drawnow");
+
+    for(int cNode=0;cNode<nNodes;cNode++)
+    {
+        if(G.what_segment(cNode)== GraphType::SOURCE)
+            SourceCut.push_back(cNode+1);
+		//mexPrintf("The loop for graph is:  %d\n", cNode);
+    }
+
+
+    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;
+		//mexPrintf("The loop for pOut:  %d\n", pOut);
+    }
+
+
+    if(nout==2)
+    {
+        out[1] = mxCreateDoubleMatrix(1, 1, mxREAL) ;
+        *mxGetPr(out[1])=flow;
+    }
+ 
+}
\ No newline at end of file

Functions/functions_GraphCut/block.h

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

Functions/functions_GraphCut/compile_mex_functions.m

+% compiling mex functions
+% using libut makes it possible to use (ctrl + C) to stop C++ from inside
+% MATLAB
+
+mex -g -O LIBPATH="C:/Program Files/MATLAB/R2019b/extern/lib/win64/microsoft" -llibut -compatibleArrayDims GraphCutMex.cpp MaxFlow.cpp graph.cpp 
+
+mex -g -O -llibut -compatibleArrayDims GraphCutMex.cpp MaxFlow.cpp graph.cpp 

Functions/functions_GraphCut/graph.cpp

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

Functions/functions_GraphCut/graph.h

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

Functions/functions_GraphCut/instances.inc

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

Functions/functions_GraphCut/maxflow.cpp

+/* maxflow.cpp */
+
+
+#include <stdio.h>
+#include "graph.h"
+#include "instances.inc"
+#include "mex.h"
+
+
+/*
+	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 */
+#ifdef __cplusplus /*function for stopping from MATLAB*/
+extern "C" bool utIsInterruptPending();
+#else
+extern bool utIsInterruptPending();
+#endif
+
+/***********************************************************************/
+
+/*
+	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 )
+	{	
+		if (utIsInterruptPending()) {        /* check for a Ctrl-C event */
+			mexPrintf("Ctrl-C Detected. END\n\n");
+			mexEvalString("drawnow");
+			return -1;
+		}
+		// mexEvalString("drawnow");
+		// 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

Functions/functions_auxiliary/click_points.m

+function [xp,yp] = click_points
+
+hold on
+[x,y,button] = ginputWhite(1);
+xp = [];
+yp = [];
+while button==1
+    xp = [xp,x];
+    yp = [yp,y];
+    plot(x,y,'co')
+    [x,y,button] = ginputWhite(1);    
+end
+end
\ No newline at end of file

Functions/functions_auxiliary/compute_normals.m

+function N = compute_normals(S,open)
+% normals for a closed (or open!) curve
+
+if nargin<2
+    open = false;
+end
+
+X = S([end,1:end,1],:); % extended S
+
+dX = X(1:end-1,:)-X(2:end,:); % dX
+Ne = [dX(:,2),-dX(:,1)]; % edge normals orthogonal to dX
+de = sum(Ne.^2,2).^0.5; % for normalization
+nonzero = de~=0;
+Ne(nonzero,:) = Ne(nonzero,:)./de(nonzero,[1 1]); % edge normals
+N = 0.5*(Ne(1:end-1,:)+Ne(2:end,:)); % vertices normals
+d = sum(N.^2,2).^0.5; % for normalization
+nonzero = d~=0;
+N(nonzero,:) = N(nonzero,:)./d(nonzero,[1 1]);
+
+if open
+    N(1,:) = Ne(2,:);
+    N(end,:) = Ne(end-1,:);
+end
\ No newline at end of file

Functions/functions_auxiliary/distribute_points.m

+function p_new = distribute_points(p,type,value,open)
+% DISTRIBUTE_POINTS   Distributes snakes points equidistantly
+% (NOTE: this version can also handle open snake)
+%
+%   DISTRIBUTE_POINTS(P) keeps the number of points
+%   DISTRIBUTE_POINTS(P,'number',N) returns N points
+%   DISTRIBUTE_POINTS(P,'ael',d) returns average edge length d
+%   Author: vand@dtu.dk
+
+if nargin<4
+    open = false;
+end
+
+if ~open
+    p = p([1:end,1],:); % closing the curve
+end
+N = size(p,1); % number of points (+ 1 for closed curve)
+dist = sqrt(sum(diff(p).^2,2)); % edge segment lengths
+t = [0;cumsum(dist)]; % current positions
+% if we want the fixed edge length then the new N could be computed 
+% from the total length of the curve which is t(end)
+
+if nargin<2 
+    N_new = N;
+else
+    switch lower(type)
+        case 'number'
+            N_new = value+~open;
+        case 'ael'
+           N_new = round(t(end)/value+~open);
+        otherwise
+           error('Unknown distribution type.')   
+    end    
+end
+
+tq = linspace(0,t(end),N_new)'; % equidistant positions
+%p_new(:,1) = interp1(t,p(:,1),tq); % distributed x
+%p_new(:,2) = interp1(t,p(:,2),tq); % distributed y
+p_new = interp1(t,p,tq); 
+if ~open
+    p_new = p_new(1:end-1,:); % opening the curve again
+end

Functions/functions_auxiliary/ginputWhite.m

+function [out1,out2,out3] = ginputWhite(arg1)
+%GINPUT Graphical input from mouse.
+%   [X,Y] = GINPUT(N) gets N points from the current axes and returns
+%   the X- and Y-coordinates in length N vectors X and Y.  The cursor
+%   can be positioned using a mouse.  Data points are entered by pressing
+%   a mouse button or any key on the keyboard except carriage return,
+%   which terminates the input before N points are entered. If the current
+%   axes is a geographic axes, the coordinates returned are latitude and
+%   longitude instead of X and Y.
+%
+%   [X,Y] = GINPUT gathers an unlimited number of points until the
+%   return key is pressed.
+%
+%   [X,Y,BUTTON] = GINPUT(N) returns a third result, BUTTON, that
+%   contains a vector of integers specifying which mouse button was
+%   used (1,2,3 from left) or ASCII numbers if a key on the keyboard
+%   was used.
+%
+%   Examples:
+%       [x,y] = ginput;
+%
+%       [x,y] = ginput(5);
+%
+%       [x, y, button] = ginput(1);
+%
+%   See also GTEXT, WAITFORBUTTONPRESS.
+
+%   Copyright 1984-2018 The MathWorks, Inc.
+
+out1 = []; out2 = []; out3 = []; y = [];
+
+if ~matlab.ui.internal.isFigureShowEnabled
+    error(message('MATLAB:hg:NoDisplayNoFigureSupport', 'ginput'))
+end
+    
+    % Check Inputs
+    if nargin == 0
+        how_many = -1;
+        b = [];
+    else
+        how_many = arg1;
+        b = [];
+        if  ~isPositiveScalarIntegerNumber(how_many) 
+            error(message('MATLAB:ginput:NeedPositiveInt'))
+        end
+        if how_many == 0
+            % If input argument is equal to zero points,
+            % give a warning and return empty for the outputs.            
+            warning (message('MATLAB:ginput:InputArgumentZero'));
+        end
+    end
+    
+    % Get figure
+    fig = gcf;
+    drawnow;
+    figure(gcf);
+    
+    % Make sure the figure has an axes
+    gca(fig);    
+    
+    % Setup the figure to disable interactive modes and activate pointers. 
+    initialState = setupFcn(fig);
+    
+    % onCleanup object to restore everything to original state in event of
+    % completion, closing of figure errors or ctrl+c. 
+    c = onCleanup(@() restoreFcn(initialState));
+    
+    drawnow
+    char = 0;
+    
+    while how_many ~= 0
+        waserr = 0;
+        try
+            keydown = wfbp;
+        catch %#ok<CTCH>
+            waserr = 1;
+        end
+        if(waserr == 1)
+            if(ishghandle(fig))
+                cleanup(c);
+                error(message('MATLAB:ginput:Interrupted'));
+            else
+                cleanup(c);
+                error(message('MATLAB:ginput:FigureDeletionPause'));
+            end
+        end
+        % g467403 - ginput failed to discern clicks/keypresses on the figure it was
+        % registered to operate on and any other open figures whose handle
+        % visibility were set to off
+        figchildren = allchild(0);
+        if ~isempty(figchildren)
+            ptr_fig = figchildren(1);
+        else
+            error(message('MATLAB:ginput:FigureUnavailable'));
+        end
+        %         old code -> ptr_fig = get(0,'CurrentFigure'); Fails when the
+        %         clicked figure has handlevisibility set to callback
+        if(ptr_fig == fig)
+            if keydown
+                char = get(fig, 'CurrentCharacter');
+                button = abs(get(fig, 'CurrentCharacter'));
+            else
+                button = get(fig, 'SelectionType');
+                if strcmp(button,'open')
+                    button = 1;
+                elseif strcmp(button,'normal')
+                    button = 1;
+                elseif strcmp(button,'extend')
+                    button = 2;
+                elseif strcmp(button,'alt')
+                    button = 3;
+                else
+                    error(message('MATLAB:ginput:InvalidSelection'))
+                end
+            end
+            
+            if(char == 13) % & how_many ~= 0)
+                % if the return key was pressed, char will == 13,
+                % and that's our signal to break out of here whether
+                % or not we have collected all the requested data
+                % points.
+                % If this was an early breakout, don't include
+                % the <Return> key info in the return arrays.
+                % We will no longer count it if it's the last input.
+                break;
+            end
+            
+            axes_handle = gca;            
+            if ~(isa(axes_handle,'matlab.graphics.axis.Axes') ...
+                    || isa(axes_handle,'matlab.graphics.axis.GeographicAxes'))
+                % If gca is not an axes, warn but keep listening for clicks. 
+                % (There may still be other subplots with valid axes)
+                warning(message('MATLAB:Chart:UnsupportedConvenienceFunction', 'ginput', axes_handle.Type));
+                continue            
+            end
+            
+            drawnow;
+            pt = get(axes_handle, 'CurrentPoint');            
+            how_many = how_many - 1;
+            
+
+            
+            out1 = [out1;pt(1,1)]; %#ok<AGROW>
+            y = [y;pt(1,2)]; %#ok<AGROW>
+            b = [b;button]; %#ok<AGROW>
+        end
+    end
+    
+    % Cleanup and Restore 
+    cleanup(c);
+    
+    if nargout > 1
+        out2 = y;
+        if nargout > 2
+            out3 = b;
+        end
+    else
+        out1 = [out1 y];
+    end
+    
+end
+
+function valid = isPositiveScalarIntegerNumber(how_many)
+valid = ~isa(how_many, 'matlab.graphics.Graphics') && ... % not a graphics handle
+        ~ischar(how_many) && ...            % is numeric
+        isscalar(how_many) && ...           % is scalar
+        (fix(how_many) == how_many) && ...  % is integer in value
+        how_many >= 0;                      % is positive
+end
+
+function key = wfbp
+%WFBP   Replacement for WAITFORBUTTONPRESS that has no side effects.
+
+fig = gcf;
+current_char = []; %#ok<NASGU>
+
+% Now wait for that buttonpress, and check for error conditions
+waserr = 0;
+try
+    h=findall(fig,'Type','uimenu','Accelerator','C');   % Disabling ^C for edit menu so the only ^C is for
+    set(h,'Accelerator','');                            % interrupting the function.
+    keydown = waitforbuttonpress;
+    current_char = double(get(fig,'CurrentCharacter')); % Capturing the character.
+    if~isempty(current_char) && (keydown == 1)          % If the character was generated by the
+        if(current_char == 3)                           % current keypress AND is ^C, set 'waserr'to 1
+            waserr = 1;                                 % so that it errors out.
+        end
+    end
+    
+    set(h,'Accelerator','C');                           % Set back the accelerator for edit menu.
+catch %#ok<CTCH>
+    waserr = 1;
+end
+drawnow;
+if(waserr == 1)
+    set(h,'Accelerator','C');                          % Set back the accelerator if it errored out.
+    error(message('MATLAB:ginput:Interrupted'));
+end
+
+if nargout>0, key = keydown; end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+end
+
+function initialState = setupFcn(fig)
+
+% Store Figure Handle. 
+initialState.figureHandle = fig; 
+
+% Suspend figure functions
+initialState.uisuspendState = uisuspend(fig);
+
+% Disable Plottools Buttons
+initialState.toolbar = findobj(allchild(fig),'flat','Type','uitoolbar');
+if ~isempty(initialState.toolbar)
+    initialState.ptButtons = [uigettool(initialState.toolbar,'Plottools.PlottoolsOff'), ...
+        uigettool(initialState.toolbar,'Plottools.PlottoolsOn')];
+    initialState.ptState = get (initialState.ptButtons,'Enable');
+    set (initialState.ptButtons,'Enable','off');
+end
+
+% Disable AxesToolbar
+initialState.axes = findobj(allchild(fig),'-isa','matlab.graphics.axis.AbstractAxes');
+tb = get(initialState.axes, 'Toolbar');
+if ~isempty(tb) && ~iscell(tb)
+    initialState.toolbarVisible{1} = tb.Visible;
+    tb.Visible = 'off';
+else
+    for i=1:numel(tb)
+        if ~isempty(tb{i})
+            initialState.toolbarVisible{i} = tb{i}.Visible;
+            tb{i}.Visible = 'off';
+        end
+    end
+end
+
+%Setup empty pointer
+cdata = NaN(16,16);
+hotspot = [8,8];
+set(gcf,'Pointer','custom','PointerShapeCData',cdata,'PointerShapeHotSpot',hotspot)
+
+% Create uicontrols to simulate fullcrosshair pointer.
+initialState.CrossHair = createCrossHair(fig);
+
+% Adding this to enable automatic updating of currentpoint on the figure 
+% This function is also used to update the display of the fullcrosshair
+% pointer and make them track the currentpoint.
+set(fig,'WindowButtonMotionFcn',@(o,e) dummy()); % Add dummy so that the CurrentPoint is constantly updated
+initialState.MouseListener = addlistener(fig,'WindowMouseMotion', @(o,e) updateCrossHair(o,initialState.CrossHair));
+
+% Get the initial Figure Units
+initialState.fig_units = get(fig,'Units');
+end
+
+function restoreFcn(initialState)
+if ishghandle(initialState.figureHandle)
+    delete(initialState.CrossHair);
+    
+    % Figure Units
+    set(initialState.figureHandle,'Units',initialState.fig_units);
+    
+    set(initialState.figureHandle,'WindowButtonMotionFcn','');
+    delete(initialState.MouseListener);
+    
+    % Plottools Icons
+    if ~isempty(initialState.toolbar) && ~isempty(initialState.ptButtons)
+        set (initialState.ptButtons(1),'Enable',initialState.ptState{1});
+        set (initialState.ptButtons(2),'Enable',initialState.ptState{2});
+    end
+    
+    % Restore axestoolbar
+    for i=1:numel(initialState.axes)
+        if ~isempty(initialState.axes(i).Toolbar)
+            initialState.axes(i).Toolbar.Visible_I = initialState.toolbarVisible{i};
+        end
+    end    
+    
+    % UISUSPEND
+    uirestore(initialState.uisuspendState);    
+end
+end
+
+function updateCrossHair(fig, crossHair)
+% update cross hair for figure.
+gap = 3; % 3 pixel view port between the crosshairs
+cp = hgconvertunits(fig, [fig.CurrentPoint 0 0], fig.Units, 'pixels', fig);
+cp = cp(1:2);
+figPos = hgconvertunits(fig, fig.Position, fig.Units, 'pixels', fig.Parent);
+figWidth = figPos(3);
+figHeight = figPos(4);
+
+% Early return if point is outside the figure
+if cp(1) < gap || cp(2) < gap || cp(1)>figWidth-gap || cp(2)>figHeight-gap
+    return
+end
+
+set(crossHair, 'Visible', 'on');
+thickness = 1; % 1 Pixel thin lines. 
+set(crossHair(1), 'Position', [0 cp(2) cp(1)-gap thickness]);
+set(crossHair(2), 'Position', [cp(1)+gap cp(2) figWidth-cp(1)-gap thickness]);
+set(crossHair(3), 'Position', [cp(1) 0 thickness cp(2)-gap]);
+set(crossHair(4), 'Position', [cp(1) cp(2)+gap thickness figHeight-cp(2)-gap]);
+end
+
+function crossHair = createCrossHair(fig)
+% Create thin uicontrols with black backgrounds to simulate fullcrosshair pointer.
+% 1: horizontal left, 2: horizontal right, 3: vertical bottom, 4: vertical top
+for k = 1:4
+    crossHair(k) = uicontrol(fig, 'Style', 'text', 'Visible', 'off', 'Units', 'pixels', 'BackgroundColor', [1 1 1], 'HandleVisibility', 'off', 'HitTest', 'off'); %#ok<AGROW>
+end
+end
+
+function cleanup(c)
+if isvalid(c)
+    delete(c);
+end
+end
+
+function dummy(~,~) 
+end

Functions/functions_auxiliary/registerImages.m

+function [MOVINGREG] = registerImages(MOVING,FIXED,MatchThreshold)
+%registerImages  Register grayscale images using auto-generated code from Registration Estimator app.
+%  [MOVINGREG] = registerImages(MOVING,FIXED) Register grayscale images
+%  MOVING and FIXED using auto-generated code from the Registration
+%  Estimator app. The values for all registration parameters were set
+%  interactively in the app and result in the registered image stored in the
+%  structure array MOVINGREG.
+
+% Auto-generated by registrationEstimator app on 29-Apr-2020
+%-----------------------------------------------------------
+
+
+% Feature-based techniques require license to Computer Vision Toolbox
+checkLicense()
+
+% Default spatial referencing objects
+fixedRefObj = imref2d(size(FIXED));
+movingRefObj = imref2d(size(MOVING));
+
+% Detect SURF features
+fixedPoints = detectSURFFeatures(FIXED,'MetricThreshold',750.000000,'NumOctaves',3,'NumScaleLevels',5);
+movingPoints = detectSURFFeatures(MOVING,'MetricThreshold',750.000000,'NumOctaves',3,'NumScaleLevels',5);
+
+% Extract features
+[fixedFeatures,fixedValidPoints] = extractFeatures(FIXED,fixedPoints,'Upright',false);
+[movingFeatures,movingValidPoints] = extractFeatures(MOVING,movingPoints,'Upright',false);
+
+% Match features
+indexPairs = matchFeatures(fixedFeatures,movingFeatures,'MatchThreshold',MatchThreshold,'MaxRatio',0.01*MatchThreshold);
+fixedMatchedPoints = fixedValidPoints(indexPairs(:,1));
+movingMatchedPoints = movingValidPoints(indexPairs(:,2));
+MOVINGREG.FixedMatchedFeatures = fixedMatchedPoints;
+MOVINGREG.MovingMatchedFeatures = movingMatchedPoints;
+
+% Apply transformation - Results may not be identical between runs because of the randomized nature of the algorithm
+tform = estimateGeometricTransform(movingMatchedPoints,fixedMatchedPoints,'affine');
+MOVINGREG.Transformation = tform;
+MOVINGREG.RegisteredImage = imwarp(MOVING, movingRefObj, tform, 'OutputView', fixedRefObj, 'SmoothEdges', true);
+
+% Store spatial referencing object
+MOVINGREG.SpatialRefObj = fixedRefObj;
+
+end
+
+function checkLicense()
+
+% Check for license to Computer Vision Toolbox
+CVSTStatus = license('test','Video_and_Image_Blockset');
+if ~CVSTStatus
+    error(message('images:imageRegistration:CVSTRequired'));
+end
+
+end
+

Functions/functions_auxiliary/regularization_matrix_open.m

+function B = regularization_matrix_open(N,alpha,beta)
+% B is an NxN matrix for imposing elasticity and rigidity to snakes
+% alpha is weigth for second derivative (elasticity)
+% beta is weigth for (-)fourth derivative (rigidity)
+
+r = zeros(1,N);
+r(1:3) = alpha*[-2 1 0]/4 + beta*[-6 4 -1]/16;
+r(end-1:end) = alpha*[0 1]/4 + beta*[-1 4]/16;
+A = toeplitz(r);
+A([1,2,end-1,end],:) = 0;
+A(2,1:3) = alpha*[1 -2 1]/4;
+A(end-1,end-2:end) = alpha*[1 -2 1]/4;
+
+B = inv(eye(N)-A);
\ No newline at end of file

Functions/functions_auxiliary/slice_points_interpolation.m

+function pk = slice_points_interpolation(p1,p2,n,m,flag_fig) 
+% SLICE_POINTS_INTERPOLATION takes two sets of points and the range of 
+% interpolition n & m and return pk which is the interpolants points 
+% [n:(n-m+1):m] 
+% 
+% the sampling is done by computing
+%  v = p1-p2
+% pk = p2 + d/norm(v)*v where d is the increment of the distances
+%
+% written by bepi@dtu.dk first version 2019.12.02
+
+%default no figure
+if nargout < 3
+    flag_fig = 0;
+end
+
+d_p = sqrt(sum((p1-p2).^2,2));
+
+pk = nan(n-m+1,3,length(d_p));
+for j = 1:length(d_p)
+    sim = linspace(0,d_p(j),n-m+1);
+    for i = 1:length(sim)
+        if sim(i) == 0
+            pk(i,1:2,j) = p2(j,:);
+        else
+            pk(i,1:2,j) = p2(j,:) + (sim(i)/d_p(j)).*(p1(j,:)-p2(j,:));
+        end
+        pk(i,3,j) = m + i-1;%n,m = (220:800)';
+    end
+end
+
+if flag_fig 
+    figure,
+    plot3([p2(:,1),p2(:,1)]',...
+          [p1(:,2),p1(:,2)]',...
+          [ones(size(p1,1),1)*n,ones(size(p1,1),1)*m]')
+    xlabel('x') ; ylabel('y') ; zlabel('z')
+    hold on
+
+    for s = 1:size(pk,3)
+        tmp = squeeze(pk(:,:,s));
+        plot3(tmp(:,1),tmp(:,2),tmp(:,3),'-xm')
+    end
+end
+
+end
+
+