diff --git a/Demo.m b/Demo.m
index 13b9f588d7778d3c73b5c5272800ae12439b5fce..63640af522c4d4dc9b458836ae034494d0454c5c 100644
--- a/Demo.m
+++ b/Demo.m
@@ -1,10 +1,20 @@
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
diff --git a/Export2inp.m b/Export2inp.m
index 440cdcd6f6b2e72a602a8bfb06a0f9f29ff949d1..a24a2763a3e9ff77481e2179faab2ed7b1c86936 100644
--- a/Export2inp.m
+++ b/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
diff --git a/Functions/functions_GraphCut/GraphCutMex.cpp b/Functions/functions_GraphCut/GraphCutMex.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0cd28e75fd2a0f9a132a276576c3706d2e36c25c
--- /dev/null
+++ b/Functions/functions_GraphCut/GraphCutMex.cpp
@@ -0,0 +1,100 @@
+#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
diff --git a/Functions/functions_GraphCut/GraphCutMex.mexa64 b/Functions/functions_GraphCut/GraphCutMex.mexa64
new file mode 100644
index 0000000000000000000000000000000000000000..a9df40a4455e24c1d68dffa4196f96e4cd9f291e
Binary files /dev/null and b/Functions/functions_GraphCut/GraphCutMex.mexa64 differ
diff --git a/Functions/functions_GraphCut/GraphCutMex.mexmaci64 b/Functions/functions_GraphCut/GraphCutMex.mexmaci64
new file mode 100644
index 0000000000000000000000000000000000000000..75cdb2c698cc1256beeb3c1826286d56f9dceda8
Binary files /dev/null and b/Functions/functions_GraphCut/GraphCutMex.mexmaci64 differ
diff --git a/Functions/functions_GraphCut/GraphCutMex.mexw64 b/Functions/functions_GraphCut/GraphCutMex.mexw64
new file mode 100644
index 0000000000000000000000000000000000000000..91d726b3416ff2f451c7efe36a90c2c845a2e55f
Binary files /dev/null and b/Functions/functions_GraphCut/GraphCutMex.mexw64 differ
diff --git a/Functions/functions_GraphCut/GraphCutMex.mexw64.pdb b/Functions/functions_GraphCut/GraphCutMex.mexw64.pdb
new file mode 100644
index 0000000000000000000000000000000000000000..6fd6cb3465dbb977dc57a697b6e9a9237d3bb484
Binary files /dev/null and b/Functions/functions_GraphCut/GraphCutMex.mexw64.pdb differ
diff --git a/Functions/functions_GraphCut/block.h b/Functions/functions_GraphCut/block.h
new file mode 100644
index 0000000000000000000000000000000000000000..93675f2f6ae8b8eb4dea5692e4f31302626eacca
--- /dev/null
+++ b/Functions/functions_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/Functions/functions_GraphCut/compile_mex_functions.m b/Functions/functions_GraphCut/compile_mex_functions.m
new file mode 100644
index 0000000000000000000000000000000000000000..798ba7880444c4e1d99fa9167c337f0cbbd2a122
--- /dev/null
+++ b/Functions/functions_GraphCut/compile_mex_functions.m
@@ -0,0 +1,7 @@
+% 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
diff --git a/Functions/functions_GraphCut/graph.cpp b/Functions/functions_GraphCut/graph.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..19b0be67b958de6ca151c2ef3224040c6fc7d8d2
--- /dev/null
+++ b/Functions/functions_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/Functions/functions_GraphCut/graph.h b/Functions/functions_GraphCut/graph.h
new file mode 100644
index 0000000000000000000000000000000000000000..87efe7ef57a36269bed7725233e994650d8949dd
--- /dev/null
+++ b/Functions/functions_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/Functions/functions_GraphCut/instances.inc b/Functions/functions_GraphCut/instances.inc
new file mode 100644
index 0000000000000000000000000000000000000000..9c9ee37c27d6e3f0434bd277f309c17b2caa8c8e
--- /dev/null
+++ b/Functions/functions_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/Functions/functions_GraphCut/maxflow.cpp b/Functions/functions_GraphCut/maxflow.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..38b5dcc53794ea37776755674593a138bb294f30
--- /dev/null
+++ b/Functions/functions_GraphCut/maxflow.cpp
@@ -0,0 +1,702 @@
+/* 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
diff --git a/Functions/functions_auxiliary/click_points.m b/Functions/functions_auxiliary/click_points.m
new file mode 100644
index 0000000000000000000000000000000000000000..8e10f2013c1fd623dadc6c8d1ab820c640b345af
--- /dev/null
+++ b/Functions/functions_auxiliary/click_points.m
@@ -0,0 +1,13 @@
+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
diff --git a/Functions/functions_auxiliary/compute_normals.m b/Functions/functions_auxiliary/compute_normals.m
new file mode 100644
index 0000000000000000000000000000000000000000..741500e0946b4e2d24730e2ae25af669e1885e95
--- /dev/null
+++ b/Functions/functions_auxiliary/compute_normals.m
@@ -0,0 +1,23 @@
+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
diff --git a/Functions/functions_auxiliary/distribute_points.m b/Functions/functions_auxiliary/distribute_points.m
new file mode 100644
index 0000000000000000000000000000000000000000..df0df4631c4ed1600279e50321b3715c1932f198
--- /dev/null
+++ b/Functions/functions_auxiliary/distribute_points.m
@@ -0,0 +1,42 @@
+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
diff --git a/Functions/functions_auxiliary/ginputWhite.m b/Functions/functions_auxiliary/ginputWhite.m
new file mode 100644
index 0000000000000000000000000000000000000000..c9c1b220a25f861a2aac1791e128a814eb429343
--- /dev/null
+++ b/Functions/functions_auxiliary/ginputWhite.m
@@ -0,0 +1,319 @@
+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
diff --git a/Functions/functions_auxiliary/registerImages.m b/Functions/functions_auxiliary/registerImages.m
new file mode 100644
index 0000000000000000000000000000000000000000..8d225e4361c4a390ac81efd9510b6e7b9f5c59d0
--- /dev/null
+++ b/Functions/functions_auxiliary/registerImages.m
@@ -0,0 +1,54 @@
+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
+
diff --git a/Functions/functions_auxiliary/regularization_matrix_open.m b/Functions/functions_auxiliary/regularization_matrix_open.m
new file mode 100644
index 0000000000000000000000000000000000000000..85816be0ec223068c0f5c15f29d277a517414fc2
--- /dev/null
+++ b/Functions/functions_auxiliary/regularization_matrix_open.m
@@ -0,0 +1,14 @@
+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
diff --git a/Functions/functions_auxiliary/slice_points_interpolation.m b/Functions/functions_auxiliary/slice_points_interpolation.m
new file mode 100644
index 0000000000000000000000000000000000000000..e841338ab4ef382b471331f8b6b0e69951c80a7f
--- /dev/null
+++ b/Functions/functions_auxiliary/slice_points_interpolation.m
@@ -0,0 +1,48 @@
+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
+
+
diff --git a/Functions/functions_cost/smoothsurface_cost.m b/Functions/functions_cost/smoothsurface_cost.m
new file mode 100644
index 0000000000000000000000000000000000000000..eb93a6736624c64e9af9d86b8d0d53564461e4ff
--- /dev/null
+++ b/Functions/functions_cost/smoothsurface_cost.m
@@ -0,0 +1,28 @@
+function cost = smoothsurface_cost(V,mu,std,s1,s2)
+% abda@dtu.dk
+
+if nargin == 1
+ mu = 23000;
+ std = 4000;
+ s1 = 5;
+ s2 = 1.5;
+elseif nargin == 3
+ s1 = 5;
+ s2 = 1.5;
+end
+
+V_exp = 1/(std*sqrt(2*pi))*exp(-(V-mu).^2/(2*std^2));
+V_exp = V_exp/(max(V_exp(:)));
+g = get_gauss(s1);
+[~,dg] = get_gauss(s2);
+cost = imfilter(imfilter(V_exp, dg','replicate'),g,'replicate');
+
+cost = (cost-min(cost(:)))/(max(cost(:))-min(cost(:)));
+
+function [g,dg] = get_gauss(s)
+
+x = -ceil(3*s):ceil(3*s);
+g = exp(-x.^2/(2*s*s));
+g = g/sum(g);
+dg = -x/(s^2).*g;
+%ddg = -g/(s^2) - x/(s^2).*dg;
\ No newline at end of file
diff --git a/Functions/functions_curvature/surfature.m b/Functions/functions_curvature/surfature.m
new file mode 100644
index 0000000000000000000000000000000000000000..93ef859710d991a642b48153fe43208f0b8fce3b
--- /dev/null
+++ b/Functions/functions_curvature/surfature.m
@@ -0,0 +1,70 @@
+function [K,H,Pmax,Pmin] = surfature(X,Y,Z)
+% SURFATURE - COMPUTE GAUSSIAN AND MEAN CURVATURES OF A SURFACE
+% [K,H] = SURFATURE(X,Y,Z), WHERE X,Y,Z ARE 2D ARRAYS OF POINTS ON THE
+% SURFACE. K AND H ARE THE GAUSSIAN AND MEAN CURVATURES, RESPECTIVELY.
+% SURFATURE RETURNS 2 ADDITIONAL ARGUEMENTS,
+% [K,H,Pmax,Pmin] = SURFATURE(...), WHERE Pmax AND Pmin ARE THE MINIMUM
+% AND MAXIMUM CURVATURES AT EACH POINT, RESPECTIVELY.
+%
+% modified by bepi@dtu.dk 2020.03.17
+%
+% Example
+% [X,Y,Z] = peaks;
+% [K,H,P1,P2] = surfature(X,Y,Z);
+% surf(X,Y,Z,H,'facecolor','interp');
+% set(gca,'clim',[-1,1])
+
+% First Derivatives
+[Xu,Xv] = gradient(X);
+[Yu,Yv] = gradient(Y);
+[Zu,Zv] = gradient(Z);
+
+% Second Derivatives
+[Xuu,~] = gradient(Xu);
+[Yuu,~] = gradient(Yu);
+[Zuu,~] = gradient(Zu);
+
+[Xuv,Xvv] = gradient(Xv);
+[Yuv,Yvv] = gradient(Yv);
+[Zuv,Zvv] = gradient(Zv);
+
+% Reshape 2D Arrays into Vectors
+Xu = Xu(:); Yu = Yu(:); Zu = Zu(:);
+Xv = Xv(:); Yv = Yv(:); Zv = Zv(:);
+Xuu = Xuu(:); Yuu = Yuu(:); Zuu = Zuu(:);
+Xuv = Xuv(:); Yuv = Yuv(:); Zuv = Zuv(:);
+Xvv = Xvv(:); Yvv = Yvv(:); Zvv = Zvv(:);
+
+Xu = [Xu Yu Zu];
+Xv = [Xv Yv Zv];
+Xuu = [Xuu Yuu Zuu];
+Xuv = [Xuv Yuv Zuv];
+Xvv = [Xvv Yvv Zvv];
+
+% First fundamental Coeffecients of the surface (E,F,G)
+E = dot(Xu,Xu,2);
+F = dot(Xu,Xv,2);
+G = dot(Xv,Xv,2);
+
+m = cross(Xu,Xv,2);
+p = sqrt(dot(m,m,2));
+n = m./[p p p];
+
+% Second fundamental Coeffecients of the surface (L,M,N)
+L = dot(Xuu,n,2);
+M = dot(Xuv,n,2);
+N = dot(Xvv,n,2);
+
+[s,t] = size(Z);
+
+% Gaussian Curvature
+K = (L.*N - M.^2)./(E.*G - F.^2);
+K = reshape(K,s,t);
+
+% Mean Curvature
+H = (E.*N + G.*L - 2.*F.*M)./(2*(E.*G - F.^2));
+H = reshape(H,s,t);
+
+% Principal Curvatures
+Pmax = H + sqrt(H.^2 - K);
+Pmin = H - sqrt(H.^2 - K);
diff --git a/Functions/functions_gridcut/grid_cut.m b/Functions/functions_gridcut/grid_cut.m
new file mode 100644
index 0000000000000000000000000000000000000000..9635300bddbcf1a3721fe34dd0d295d7c71b664d
--- /dev/null
+++ b/Functions/functions_gridcut/grid_cut.m
@@ -0,0 +1,252 @@
+function s = grid_cut(s_cost, r_cost, delta_xy, wrap_xy, delta_lu)
+%GRID_CUT detects terrain-like surfaces in a volumetric data.
+%
+% Detected surfaces are terrain-like, i.e. height z is a function of the
+% position (x,y). Only on-surface costs, only in-region costs or both type
+% of costs may be used. When using in-region costs for multiple surfaces,
+% surfaces are non-intersecting and ordered, so that the first surface has
+% the smallest z values, and the last surface has the largest z values.
+%
+% S = GRID_CUT(S_COST,R_COST,DELTA_XY,WRAP_XY,DELTA_LU)
+%
+% S_COST is a size (X,Y,Z,K) surface cost function. Cost (x,y,z,k) is an
+% inverse to likelihood of a voxel (x,y,z) containing surface k.
+% K is the number of surfaces and Z is up. When using only in-region
+% costs, S_COST should either contain only zeros or be empty.
+% R_COST is a size (X,Y,Z,K+1) region cost function. Cost (x,y,z,k) is an
+% inverse to likelihood of a voxel (x,y,z) being in a region between
+% the surface k and k+1, where the first and the last region are
+% boundedend with the volume boundary. When using only on-surface
+% costs, R_COST should either contain only zeros or be empty.
+% DELTA_XY is a size (K,2) array of stiffness parameters for x and y. If a
+% size (1,2) array is given, the same stiffness parameteres are used
+% for all surfaces. If size (K,1) array is given, the same stiffness
+% is used for x and y.
+% WRAP_XY is a length 2 array of boolean wrap options for x and y. If not
+% given defaults to false.
+% DELTA_LU is a size (K-1,2) array of lower and upper surface overlaps,
+% [dl du]. The constraint is: dl <= surf_k - surf_k+1 <= du. When
+% using in-region costs, the surfaces are non-intersecting and
+% ordered which requires dl > 0 and du >= dl. If a size (1,2) array
+% is given, the same overlap parameteres are used for all surfaces.
+% If not given, defaults to "no constraint" (dl=-Z, du=Z) when using
+% only on-surface costs and to "no overlap" (dl=1, du=Z) otherwise.
+% S is a size (X,Y,K) matrix of z coordinates for K segmented surfaces.
+%
+% Based on:
+% "Optimal surface segmentation in volumetric images-a graph-theoretic
+% approach." Li, Wu, Chen and Sonka. PAMI 2006.
+% "Incorporation of Regional Information in Optimal 3-D Graph Search with
+% Application for Intraretinal Layer Segmentation of Optical Coherence
+% Tomography Images". Haeker, Wu, Abramoff, Kardon and Sonka. IPMI 2007.
+%
+% Author: Vedrana Andersen Dahl, vand@dtu.dk, 2013
+%
+% NOTE: needs Kolmogorovs implementation of the algorithm from
+% "An Experimental Comparison of Min-Cut/Max-Flow Algorithms for
+% Energy Minimization in Vision."
+% Yuri Boykov and Vladimir Kolmogorov. PAMI 2004
+
+% assigning an unique index to each vertex in the first 3D graph (first surface)
+
+if nargin<5
+ delta_lu = [];
+end
+if nargin<4
+ wrap_xy = [];
+end
+[s_cost, r_cost, delta_xy, wrap_xy, delta_lu] = assign_defaults...
+ (s_cost, r_cost, delta_xy, wrap_xy, delta_lu);
+[X,Y,Z,K] = size(s_cost);
+dimension = [X Y Z];
+indices = reshape(1:prod(dimension),dimension);
+base = indices(:,:,1);
+next_surface = prod(dimension); % adding to indices is a shift to next surface
+layers = (0:(K-1))*next_surface;
+
+wrap = [wrap_xy 0]; % we do not wrap in z direction
+
+% EDGES WHICH ARE THE SAME FOR ALL SURFACES:
+% intracolumn arcs pointing down, Equation (2)
+Ea1 = displacement_to_edges(indices, [0 0 -1], wrap); % first surface
+Ea = repmat(Ea1,[K 1]) + kron(layers(:),ones(size(Ea1))); % all surfaces
+Ea = [Ea,ones(size(Ea,1),1)*[inf,0]]; % assigning [inf 0] weight
+% base edges, part of intercolumn arcs, Equation (3)
+erxb1 = displacement_to_edges(base, [1 0 0], wrap);
+eryb1 = displacement_to_edges(base, [0 1 0], wrap);
+Erb1 = [erxb1; eryb1]; % first surface
+Erb = repmat(Erb1,[K 1]) + kron(layers(:),ones(size(Erb1))); % all surfaces
+% assigning [inf inf] weight, base edges are in both directions
+Erb = [Erb,ones(size(Erb,1),1)*[inf,inf]];
+
+% EDGES WHICH DEPEND ON SURFACE STIFFNESS:
+% slanted edges, part of intercolumn arcs, Equation (3)
+Erpm = []; % preallocation length could be computed from X,Y,Z,delta and wrap
+for k = 1:K
+ erxp = displacement_to_edges(indices, [1 0 -delta_xy(k,1)], wrap);
+ erxm = displacement_to_edges(indices, [-1 0 -delta_xy(k,1)], wrap);
+ eryp = displacement_to_edges(indices, [0 1 -delta_xy(k,2)], wrap);
+ erym = displacement_to_edges(indices, [0 -1 -delta_xy(k,2)], wrap);
+ erpm = [erxp; erxm; eryp; erym] + layers(k);
+ Erpm = [Erpm;erpm];
+end
+% assigning [inf 0] weight, slanted edges are in one direction
+Erpm = [Erpm,ones(size(Erpm,1),1)*[inf,0]];
+Er = [Erpm; Erb]; % all inter edges, for all surfaces
+E = [Ea;Er]; % all intracolumn and intercolumn edges, for all surfaces
+
+% EDGES WHICH DEPEND ON DISTANCE BETWEEN SURFACES:
+% intersurface arc, Equation (4)
+if K>1 % only if we have more than 1 surface
+ Es = []; % preallocation length could be computed
+ for k = 1:K-1
+ esl = displacement_to_edges(indices, [0 0 delta_lu(k,1)], wrap);
+ esl = esl + ones(size(esl,1),1)*layers([k,k+1]);
+ esu = displacement_to_edges(indices, [0 0 -delta_lu(k,2)], wrap);
+ esu = esu + ones(size(esu,1),1)*layers([k+1,k]);
+ es = [esl;esu];
+ Es = [Es;es];
+ end
+ Es = [Es,ones(size(Es,1),1)*[inf,0]]; % assigning [inf 0] weight
+
+ % intersurface base edges
+ esb = [layers(1:end-1)',layers(2:end)']+1; %first vertex in all surfaces
+ Esb = [esb,ones(K-1,1)*[inf,inf]]; % assigning [inf inf]
+ E = [E;Es;Esb]; % all intracolumn, intercolumn and intersurface edges, for all surfaces
+end
+
+% EDGES WHICH DEPEND ON SURFACE LIKELIHOOD:
+% up to here we do not use s_cost -- consider efficient algorithm if solving
+% multiple problems of the same size
+% vertex s_cost, Equation (1), done simulatniously for all surfaces
+w_on = -1*ones([dimension,K]); % to prevent empty solution, see second half of section 4.1
+w_on(:,:,2:end,:) = double(s_cost(:,:,2:end,:))-double(s_cost(:,:,1:end-1,:));
+% In case of layered surfaces vertices which can't be realized should be removed,
+% i.e. topmost vertices of lower surface and lowest vertices of higher surface.
+% Instead of removing vertices, I assign inf weight to topmost vertices.
+% And this only when delta_lu are both positive (both negative can be
+% avoided by proper ordering of surfaces).
+if K>1 % only if we have more than 1 surface
+ for k=1:K-1 % assigning inf on-surface weights instead of removing vertices
+ if delta_lu(k,1)>0 && delta_lu(k,1)<= delta_lu(k,2) % bit clumsy
+ w_on(:,:,end-delta_lu(k,1)+1:end,k) = inf;
+ end
+ end
+end
+
+% In-region cost, converting to double to avoid problems when using images
+w_in = double(r_cost(:,:,:,1:end-1))-double(r_cost(:,:,:,2:end));
+% There are issues concerning in-region cost e.g. a better way of preventing
+% empty solutions or topmost solutions.
+w_in(:,:,1,:) = -inf; % preventing empty solution
+
+w = w_on + w_in;
+
+Vp_ix = find(w(:)>=0); % positive vertices, to be connected to sink
+Vm_ix = find(w(:)<0); % negative vertices, to be connected to source
+Es = [Vm_ix, -w(Vm_ix), zeros(length(Vm_ix),1)]; % source edges
+Et = [Vp_ix, zeros(length(Vp_ix),1), w(Vp_ix)]; % sink edges
+Est = [Es;Et]; % all terminal edges
+
+% FINDING GRAPH CUT USING MAGIC
+Scut = GraphCutMex(prod(dimension)*K,Est,E);
+
+% retreiving surfaces as the upper envelope of Scut
+S = zeros(X,Y,Z,K);
+S(Scut) = 1;
+s = zeros(X,Y,K);
+for ki = 1:K
+ for x = 1:X
+ for y = 1:Y
+ s(x,y,ki) = find(S(x,y,:,ki),1,'last');
+ end
+ end
+end
+end
+
+% NOTE: when some of the intersurface edges point upwards, i.e. when
+% delta_l is positive, the lowest delta_l rows of a graph for higher
+% surface can not appear in any feasible solution and can be removed.
+% This has NOT been implemented here.
+
+function edges = displacement_to_edges(indices, disp, wrap)
+% indices -- 3D volume of indices
+% disp -- length 3 vector of displacements, given as [x y,z]
+% wrapping -- length 3 vector of boolean options, given as [wrapx wrapy wrapz]
+% edges -- two columns of indices, each line is an edge.
+
+[x_from,x_into] = displace_1D_indices(size(indices,1),disp(1),wrap(1));
+[y_from,y_into] = displace_1D_indices(size(indices,2),disp(2),wrap(2));
+[z_from,z_into] = displace_1D_indices(size(indices,3),disp(3),wrap(3));
+
+indices_from = indices(x_from,y_from,z_from);
+indices_into = indices(x_into,y_into,z_into);
+
+edges = [indices_from(:),indices_into(:)];
+end
+
+function [from,into] = displace_1D_indices(dim,disp,wrap)
+% dim -- length of the indices vector
+% disp -- length of the displacement
+% wrap -- boolean, indicating trimming or wrapping
+
+indices = 1:dim;
+if wrap
+ from = indices([1-min(disp,0):end,1:-min(disp,0)]);
+ into = indices([1+max(disp,0):end,1:max(disp,0)]);
+else
+ from = indices(1-min(disp,0):end-max(disp,0));
+ into = indices(1+max(disp,0):end+min(disp,0));
+end
+end
+
+function [s_cost, r_cost, delta_xy, wrap_xy, delta_lu] = assign_defaults...
+ (s_cost, r_cost, delta_xy, wrap_xy, delta_lu)
+% either cost_s or cost_r need to be given (not empty and not all zeros)
+if isempty(s_cost)
+ s_cost = zeros(size(r_cost)-[0,0,0,1]);
+end
+if isempty(r_cost)
+ r_cost = zeros([size(s_cost,1),size(s_cost,2),size(s_cost,3),...
+ size(s_cost,4)]+[0,0,0,1]); % to allow s_cost to be 3D
+ regional = false;
+elseif all(r_cost(:)==0)
+ regional = false;
+else
+ regional = true;
+end
+
+[X,Y,Z,K] = size(s_cost);
+[Xr,Yr,Zr,Krplus1] = size(r_cost);
+if ~X==Xr || ~Y==Yr || ~Z==Zr || ~(K+1)==Krplus1
+ error('Error using grid_cut. Dimensions of s_cost and r_cost must agree.')
+end
+
+% smoothness constraint may be given once for all surfaces
+if size(delta_xy,1)==1
+ delta_xy = ones(K,1)*delta_xy;
+end
+% smoothness constraint may be given once for both directions
+if size(delta_xy,2)==1
+ delta_xy = delta_xy*ones(1,2);
+end
+
+if isempty(wrap_xy) || (numel(wrap_xy)==1 && wrap_xy==0)
+ wrap_xy = [0 0];
+end
+
+% overlap constraint defaults to no overlap in region case
+% and no constraint in surface case
+if isempty(delta_lu)
+ if regional % no overlap
+ delta_lu = [1 Z];
+ else % no constraint
+ delta_lu = [-Z Z];
+ end
+end
+
+% overlap constraint may be given once for all surface pairs
+if all(size(delta_lu)==[1,2])
+ delta_lu = ones(K-1,1)*delta_lu;
+end
+end
\ No newline at end of file
diff --git a/Functions/functions_mesh/do_mesh_example.m b/Functions/functions_mesh/do_mesh_example.m
new file mode 100644
index 0000000000000000000000000000000000000000..6166d1ab9f0dd78280818b522fbfda17c15a0719
--- /dev/null
+++ b/Functions/functions_mesh/do_mesh_example.m
@@ -0,0 +1,87 @@
+clear all
+
+Lx = 2;
+Ly = 1.3;
+Lz = 0.5;
+
+nx = 2;%10
+ny = 3;%7
+nz = 4;%5
+
+cnt = 1;
+for i=1:nx
+ for j=1:ny
+ for k=1:nz
+ nodes(cnt,:) = [cnt i*Lx/nx j*Ly/ny k*Lz/nz];
+ cnt = cnt + 1;
+ end
+ end
+end
+
+number_of_nodes = cnt-1;
+
+clf,plot3(nodes(:,2),nodes(:,3),nodes(:,4),'.')
+
+elements = [];
+cnt = 1;
+for i=1:nx-1
+ for j=1:ny-1
+ for k=1:nz-1
+ x0 = i*Lx/nx;
+ y0 = j*Ly/ny;
+ z0 = k*Lz/nz;
+ x1 = (i+1)*Lx/nx;
+ y1 = (j+1)*Ly/ny;
+ z1 = (k+1)*Lz/nz;
+
+ dist_to_nodes = ((nodes(:,2)-x0).^2+(nodes(:,3)-y0).^2+(nodes(:,4)-z0).^2);
+ [minval,minidx] = min(dist_to_nodes);
+ corner(1) = minidx;
+ dist_to_nodes = ((nodes(:,2)-x1).^2+(nodes(:,3)-y0).^2+(nodes(:,4)-z0).^2);
+ [minval,minidx] = min(dist_to_nodes);
+ corner(2) = minidx;
+ dist_to_nodes = ((nodes(:,2)-x1).^2+(nodes(:,3)-y1).^2+(nodes(:,4)-z0).^2);
+ [minval,minidx] = min(dist_to_nodes);
+ corner(3) = minidx;
+ dist_to_nodes = ((nodes(:,2)-x0).^2+(nodes(:,3)-y1).^2+(nodes(:,4)-z0).^2);
+ [minval,minidx] = min(dist_to_nodes);
+ corner(4) = minidx;
+
+ dist_to_nodes = ((nodes(:,2)-x0).^2+(nodes(:,3)-y0).^2+(nodes(:,4)-z1).^2);
+ [minval,minidx] = min(dist_to_nodes);
+ corner(5) = minidx;
+
+ dist_to_nodes = ((nodes(:,2)-x1).^2+(nodes(:,3)-y0).^2+(nodes(:,4)-z1).^2);
+ [minval,minidx] = min(dist_to_nodes);
+ corner(6) = minidx;
+
+ dist_to_nodes = ((nodes(:,2)-x1).^2+(nodes(:,3)-y1).^2+(nodes(:,4)-z1).^2);
+ [minval,minidx] = min(dist_to_nodes);
+ corner(7) = minidx;
+ dist_to_nodes = ((nodes(:,2)-x0).^2+(nodes(:,3)-y1).^2+(nodes(:,4)-z1).^2);
+ [minval,minidx] = min(dist_to_nodes);
+ corner(8) = minidx;
+
+ elements(cnt,:)=[cnt corner([1 2 3 4 5 6 7 8 ])];
+ cnt=cnt+1;
+
+ end
+ end
+end
+
+figure,patch('Faces',elements(:,2:end),'Vertices',nodes(:,2:4))
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% fid = fopen('output_mesh_example.inp','w');
+% fprintf(fid,'*Part, name=my_mesh \n');
+% fprintf(fid,'*Node \n');
+% for i=1:size(nodes,1)
+% fprintf(fid,'%d, %f, %f, %f \n',nodes(i,:));
+% end
+% fprintf(fid,'*Element, type=C3D8R \n');
+% for i = 1:size(elements,1)
+% fprintf(fid,'%d, %d, %d, %d, %d, %d, %d, %d, %d \n',elements(i,:));
+% end
+%
+% fprintf(fid,'*End Part \n');
\ No newline at end of file
diff --git a/Functions/functions_mesh/elements_mesh.m b/Functions/functions_mesh/elements_mesh.m
new file mode 100644
index 0000000000000000000000000000000000000000..f62f26a20d0c53f42c160fd2bf29251b8c16900e
--- /dev/null
+++ b/Functions/functions_mesh/elements_mesh.m
@@ -0,0 +1,19 @@
+function ele = elements_mesh(mesh)
+
+sz0 = size(mesh);
+ny = sz0(1);
+nz = sz0(2);
+cnt = 1;
+ele = zeros((nz-1)*(ny-1),9);
+%h = waitbar(0,'1','Name',inputname(mesh));
+for j = 1:ny-1
+ for k = 1:nz-1
+ i = (j-1)*nz+k;
+ ele(cnt,:)=[cnt i ny*nz+i ny*nz+nz+i nz+i i+1 ny*nz+i+1 ny*nz+nz+i+1 nz+i+1];
+ cnt=cnt+1;
+ end
+ %waitbar(j/(ny-1),h,sprintf('j: %d',j))
+end
+%close(h)
+
+end
\ No newline at end of file
diff --git a/Functions/functions_mesh/write_mesh.m b/Functions/functions_mesh/write_mesh.m
new file mode 100644
index 0000000000000000000000000000000000000000..9495f455914a59a7cdee8f15784c05d714fa7c90
--- /dev/null
+++ b/Functions/functions_mesh/write_mesh.m
@@ -0,0 +1,18 @@
+function write_mesh(nodes,elements,filename)
+
+fid = fopen([filename,'.inp'],'w');
+[~,name,~] = fileparts(filename);
+fprintf(fid,'*Part, name=%s \n',name);
+fprintf(fid,'*Node \n');
+for i=1:size(nodes,1)
+ fprintf(fid,'%d, %f, %f, %f \n',nodes(i,:));
+end
+
+fprintf(fid,'*Element, type=C3D8R \n');
+for i = 1:size(elements,1)
+ fprintf(fid,'%d, %d, %d, %d, %d, %d, %d, %d, %d \n',elements(i,:));
+end
+
+fprintf(fid,'*End Part \n');
+fclose(fid);
+end
\ No newline at end of file
diff --git a/Functions/functions_visualization/show_matvol.m b/Functions/functions_visualization/show_matvol.m
new file mode 100644
index 0000000000000000000000000000000000000000..24992ce1a504149b4dd875305d1dae26aa0a9de4
--- /dev/null
+++ b/Functions/functions_visualization/show_matvol.m
@@ -0,0 +1,48 @@
+function show_matvol(V,clim)
+% SHOW_MATVOL Shows volumetric data
+% show_matvol(V,clim)
+% Inputs: volume, grayscale color limits (optional)
+% Copyright 2018 Vedrana Andersen Dahl, vand@dtu.dk
+
+figure('Units','Normalized','Position',[0.1 0.3 0.5 0.6],...
+ 'KeyPressFcn',@key_press);
+dim = size(V);
+if nargin<2 || isempty(clim)
+ clim = [min(double(V(:))),max(double(V(:)))];
+end
+z = round(0.5*(dim(3)+1));
+%update_drawing
+i = imagesc(V(:,:,z),clim);
+title(['slice ',num2str(z),'/',num2str(dim(3))]), axis image ij
+colormap gray
+drawnow
+
+%%%%%%%%%% CALLBACK FUNCTIONS %%%%%%%%%%
+ function key_press(~,object)
+ % keyboard commands
+ switch object.Key
+ case 'uparrow'
+ z = min(z+1,dim(3));
+ case 'downarrow'
+ z = max(z-1,1);
+ case 'rightarrow'
+ z = min(z+10,dim(3));
+ case 'leftarrow'
+ z = max(z-10,1);
+ case 'pagedown'
+ z = min(z+50,dim(3));
+ case 'pageup'
+ z = max(z-50,1);
+ end
+ update_drawing
+ end
+
+%%%%%%%%%% HELPING FUNCTIONS %%%%%%%%%%
+ function update_drawing
+ set(i,'CData',V(:,:,z))
+ title(['slice ',num2str(z),'/',num2str(dim(3))])
+ drawnow
+ end
+
+end
+
diff --git a/Functions/functions_visualization/suptitle_modified.m b/Functions/functions_visualization/suptitle_modified.m
new file mode 100644
index 0000000000000000000000000000000000000000..3fafe385b8903535b811ec06fa7cd647510c3af7
--- /dev/null
+++ b/Functions/functions_visualization/suptitle_modified.m
@@ -0,0 +1,109 @@
+function hout=suptitle_modified(str,interpreter)
+%SUPTITLE puts a title above all subplots.
+%
+% SUPTITLE('text') adds text to the top of the figure
+% above all subplots (a "super title"). Use this function
+% after all subplot commands.
+%
+% SUPTITLE is a helper function for yeastdemo.
+
+% Copyright 2003-2014 The MathWorks, Inc.
+
+
+% Warning: If the figure or axis units are non-default, this
+% function will temporarily change the units.
+
+% Parameters used to position the supertitle.
+
+% Amount of the figure window devoted to subplots
+plotregion = .92;
+
+% Y position of title in normalized coordinates
+titleypos = .95;
+
+% Fontsize for supertitle
+fs = get(gcf,'defaultaxesfontsize')+4;
+
+% Fudge factor to adjust y spacing between subplots
+fudge=1;
+
+haold = gca;
+figunits = get(gcf,'units');
+
+% Get the (approximate) difference between full height (plot + title
+% + xlabel) and bounding rectangle.
+
+if ~strcmp(figunits,'pixels')
+ set(gcf,'units','pixels');
+ pos = get(gcf,'position');
+ set(gcf,'units',figunits);
+else
+ pos = get(gcf,'position');
+end
+ff = (fs-4)*1.27*5/pos(4)*fudge;
+
+% The 5 here reflects about 3 characters of height below
+% an axis and 2 above. 1.27 is pixels per point.
+
+% Determine the bounding rectangle for all the plots
+
+h = findobj(gcf,'Type','axes');
+
+oldUnits = get(h, {'Units'});
+if ~all(strcmp(oldUnits, 'normalized'))
+ % This code is based on normalized units, so we need to temporarily
+ % change the axes to normalized units.
+ set(h, 'Units', 'normalized');
+ cleanup = onCleanup(@()resetUnits(h, oldUnits));
+end
+
+max_y=0;
+min_y=1;
+oldtitle = [];
+numAxes = length(h);
+thePositions = zeros(numAxes,4);
+for i=1:numAxes
+ pos=get(h(i),'pos');
+ thePositions(i,:) = pos;
+ if ~strcmp(get(h(i),'Tag'),'suptitle')
+ if pos(2) < min_y
+ min_y=pos(2)-ff/5*3;
+ end
+ if pos(4)+pos(2) > max_y
+ max_y=pos(4)+pos(2)+ff/5*2;
+ end
+ else
+ oldtitle = h(i);
+ end
+end
+
+if max_y > plotregion
+ scale = (plotregion-min_y)/(max_y-min_y);
+ for i=1:numAxes
+ pos = thePositions(i,:);
+ pos(2) = (pos(2)-min_y)*scale+min_y;
+ pos(4) = pos(4)*scale-(1-scale)*ff/5*3;
+ set(h(i),'position',pos);
+ end
+end
+
+np = get(gcf,'nextplot');
+set(gcf,'nextplot','add');
+if ~isempty(oldtitle)
+ delete(oldtitle);
+end
+axes('pos',[0 1 1 1],'visible','off','Tag','suptitle');
+ht=text(.5,titleypos-1,str,'Interpreter',interpreter);set(ht,'horizontalalignment','center','fontsize',fs);
+set(gcf,'nextplot',np);
+axes(haold);
+if nargout
+ hout=ht;
+end
+end
+
+function resetUnits(h, oldUnits)
+ % Reset units on axes object. Note that one of these objects could have
+ % been an old supertitle that has since been deleted.
+ valid = isgraphics(h);
+ set(h(valid), {'Units'}, oldUnits(valid));
+end
\ No newline at end of file
diff --git a/LayerDetection.m b/LayerDetection.m
index 6db4c491bf13de32d5ec8f0acf9eb26b352e8b9d..81252ca57cce20e8d22b08c552a60c8358323c36 100644
--- a/LayerDetection.m
+++ b/LayerDetection.m
@@ -1,7 +1,4 @@
%% detecting layers
-flag_save = 1;
-flag_fig = 0; % if you want to check the results
-
% create result folder
save_folder = 'Results\layer\';
if ~exist(save_folder,'dir')
diff --git a/PointSelection.m b/PointSelection.m
index 0d60119e889a36cf68d0a753aef487c5b97146e5..7eb0e6fc980b731f18b7e1a4e7846dbc4ea80864 100644
--- a/PointSelection.m
+++ b/PointSelection.m
@@ -1,189 +1,185 @@
%% move from pipeline to here to select the points
-flag_run = 1;
save_folder = 'Results\point\';
if ~exist(save_folder,'dir')
mkdir(save_folder)
end
-if flag_run
- %% contrast enhancement
- ce = [.7 1.5];
-
- points = [];
- points.contrast_enhansed_range = ce;%contrast enhansed range
- %% initail points layer
-
- %f = msgbox({'Suggestions:';'0 & 45 dg: first & last slice';'90 & 180 dg: first, mid & last or first, 1/3, 2/3 & last'},'Choose Slices','warn');
- fprintf('Suggestions:\n0 & 45 dg: first & last slice\n90 & 180 dg: first, middle & last slice or\n first, 1/3, 2/3 & last slice\n')
- show_matvol(Vol)
- % move to desire slice and run this line to get the points for that slice
- fprintf('Controls\n--------\n up-arrow : Next slice in volume \n down-arrow : Previous slice in volume \n right-arrow : +10 slices \n left-arrow : -10 slices \n PgUp : -50 slices \n PgDown : +50 slices \n')
- fprintf('Choose the slices by scrolling in the volume and press Enter \n DO NOT CLOSE THE FIGURE!')
- input('');
-
- answer = inputdlg('Enter space-separated slice numbers:','Slices',[1 50],{'100 2550'});
- points.slice_no = str2num(answer{1});
-
- num_layers = length(points.slice_no);
-
- %% select points, measure the length and distribute points
- for l = 1:2
- flag_layer = layer_name{l};
- fprintf('%s layer:\n',flag_layer)
- switch length(points.slice_no)
- case 2
- fprintf('Go to slice %d and select points on the %s layer\n',points.slice_no(1),flag_layer)
- input('');
- [x1,y1] = click_points;
- fprintf('Go to slice %d and select %d points on the %s layer\n',points.slice_no(2),length(x1),flag_layer)
- input('');
- [x2,y2] = click_points;
-
- pts1 = round([x1',y1']);
- pts2 = round([x2',y2']);
- points.two_slice = cat(3,pts1,pts2);
-
- if l == 1 %compute the length of the layers. this gives us number of px per layers and we keep it same for inner and outer layer
- dx = [];
- dy = [];
- len0 = [];
- for s = 1:2
- dx = points.two_slice(1:end-1,1,s)-points.two_slice(2:end,1,s);
- dy = points.two_slice(1:end-1,2,s)-points.two_slice(2:end,2,s);
- len0(1,s) = sum(sqrt(dx.^2+dy.^2));
- end
- m_l = round(mean(len0));
- end
-
- p_new1 = distribute_points(pts1,'number',m_l,1);
- p_new2 = distribute_points(pts2,'number',m_l,1);
-
- points.all_slice = slice_points_interpolation(p_new2,p_new1,...
- points.slice_no(2),points.slice_no(1),0);
+%% initail points layer
+points = [];
+points.contrast_enhansed_range = ce;%contrast enhansed range
+
+fprintf('Suggestions:\n0 & 45 dg: first & last slice\n90 & 180 dg: first, middle & last slice or\n first, 1/3, 2/3 & last slice\n')
+show_matvol(Vol)
+
+% move to desire slice and run this line to get the points for that slice
+fprintf('Controls\n--------\n up-arrow : Next slice in volume \n down-arrow : Previous slice in volume \n right-arrow : +10 slices \n left-arrow : -10 slices \n PgUp : -50 slices \n PgDown : +50 slices \n')
+fprintf('Choose the slices by scrolling in the volume and press Enter \n DO NOT CLOSE THE FIGURE!')
+input('');
+
+answer = inputdlg('Enter space-separated slice numbers:','Slices',[1 50],{'100 2550'});
+points.slice_no = str2num(answer{1});
+
+num_layers = length(points.slice_no);
+
+%% select points, measure the length and distribute points
+for l = 1:2
+ flag_layer = layer_name{l};
+ fprintf('%s layer:\n',flag_layer)
+ switch length(points.slice_no)
+ case 2
+ fprintf('Go to slice %d and select points on the %s layer\n',points.slice_no(1),flag_layer)
+ input('');
+ [x1,y1] = click_points;
+ fprintf('Go to slice %d and select %d points on the %s layer\n',points.slice_no(2),length(x1),flag_layer)
+ input('');
+ [x2,y2] = click_points;
+
+ pts1 = round([x1',y1']);
+ pts2 = round([x2',y2']);
+ points.two_slice = cat(3,pts1,pts2);
+
+ if l == 1 %compute the length of the layers. this gives us number of px per layers and we keep it same for inner and outer layer
dx = [];
dy = [];
- len1 = [];
- for s = 1:size(points.all_slice,1)
- dx = points.all_slice(s,1,1:end-1)-points.all_slice(s,1,2:end);
- dy = points.all_slice(s,2,1:end-1)-points.all_slice(s,2,2:end);
- len1(1,s) = sum(sqrt(dx.^2+dy.^2));
- end
- fprintf('Approximated length ~ %d or %d \n',m_l,round(mean(len1)));
-
- case 3
-
- fprintf('Go to slice %d and select points on the %s layer\n',points.slice_no(1),flag_layer)
- input('');
- [x1,y1] = click_points;
- fprintf('Go to slice %d and select %d points on the %s layer\n',points.slice_no(2),length(x1),flag_layer)
- input('');
- [x2,y2] = click_points;
- fprintf('Go to slice %d and select %d points on the %s layer\n',points.slice_no(3),length(x1),flag_layer)
- input('');
- [x3,y3] = click_points;
-
- pts1 = round([x1',y1']);
- pts2 = round([x2',y2']);
- pts3 = round([x3',y3']);
-
- points.three_slice = cat(3,pts1,pts2,pts3);
-
- if l == 1 %compute the length of the layers. this gives us number of px per layers and we keep it same for inner and outer layer
- dx = [];
- dy = [];
- len0 = [];
- for s = 1:3
- dx = points.three_slice(1:end-1,1,s)-points.three_slice(2:end,1,s);
- dy = points.three_slice(1:end-1,2,s)-points.three_slice(2:end,2,s);
- len0(1,s) = sum(sqrt(dx.^2+dy.^2));
- end
- m_l = round(mean(len0));
+ len0 = [];
+ for s = 1:2
+ dx = points.two_slice(1:end-1,1,s)-points.two_slice(2:end,1,s);
+ dy = points.two_slice(1:end-1,2,s)-points.two_slice(2:end,2,s);
+ len0(1,s) = sum(sqrt(dx.^2+dy.^2));
end
-
- p_new1 = distribute_points(pts1,'number',m_l,1);
- p_new2 = distribute_points(pts2,'number',m_l,1);
- p_new3 = distribute_points(pts3,'number',m_l,1);
-
- slice_23 = slice_points_interpolation(p_new3,p_new2,...
- points.slice_no(3),points.slice_no(2),0);
- slice_12 = slice_points_interpolation(p_new2,p_new1,...
- points.slice_no(2)-1,points.slice_no(1),0);
- points.all_slice = [slice_12;slice_23];
-
+ m_l = round(mean(len0));
+ end
+
+ p_new1 = distribute_points(pts1,'number',m_l,1);
+ p_new2 = distribute_points(pts2,'number',m_l,1);
+
+ points.all_slice = slice_points_interpolation(p_new2,p_new1,...
+ points.slice_no(2),points.slice_no(1),0);
+ dx = [];
+ dy = [];
+ len1 = [];
+ for s = 1:size(points.all_slice,1)
+ dx = points.all_slice(s,1,1:end-1)-points.all_slice(s,1,2:end);
+ dy = points.all_slice(s,2,1:end-1)-points.all_slice(s,2,2:end);
+ len1(1,s) = sum(sqrt(dx.^2+dy.^2));
+ end
+ fprintf('Approximated length ~ %d or %d \n',m_l,round(mean(len1)));
+
+ case 3
+
+ fprintf('Go to slice %d and select points on the %s layer\n',points.slice_no(1),flag_layer)
+ input('');
+ [x1,y1] = click_points;
+ fprintf('Go to slice %d and select %d points on the %s layer\n',points.slice_no(2),length(x1),flag_layer)
+ input('');
+ [x2,y2] = click_points;
+ fprintf('Go to slice %d and select %d points on the %s layer\n',points.slice_no(3),length(x1),flag_layer)
+ input('');
+ [x3,y3] = click_points;
+
+ pts1 = round([x1',y1']);
+ pts2 = round([x2',y2']);
+ pts3 = round([x3',y3']);
+
+ points.three_slice = cat(3,pts1,pts2,pts3);
+
+ if l == 1 %compute the length of the layers. this gives us number of px per layers and we keep it same for inner and outer layer
dx = [];
dy = [];
- len1 = [];
- for s = 1:size(points.all_slice,1)
- dx = points.all_slice(s,1,1:end-1)-points.all_slice(s,1,2:end);
- dy = points.all_slice(s,2,1:end-1)-points.all_slice(s,2,2:end);
- len1(1,s) = sum(sqrt(dx.^2+dy.^2));
+ len0 = [];
+ for s = 1:3
+ dx = points.three_slice(1:end-1,1,s)-points.three_slice(2:end,1,s);
+ dy = points.three_slice(1:end-1,2,s)-points.three_slice(2:end,2,s);
+ len0(1,s) = sum(sqrt(dx.^2+dy.^2));
end
- fprintf('Approximated length ~ %d or %d \n',m_l,round(mean(len1)));
-
- case 4
-
- fprintf('Go to slice %d and select points on the %s layer\n',points.slice_no(1),flag_layer)
- input('');
- [x1,y1] = click_points;
- fprintf('Go to slice %d and select %d points on the %s layer\n',points.slice_no(2),length(x1),flag_layer)
- input('');
- [x2,y2] = click_points;
- fprintf('Go to slice %d and select %d points on the %s layer\n',points.slice_no(3),length(x1),flag_layer)
- input('');
- [x3,y3] = click_points;
- fprintf('Go to slice %d and select %d points on the %s layer\n',points.slice_no(4),length(x1),flag_layer)
- input('');
- [x4,y4] = click_points;
-
- pts1 = round([x1',y1']);
- pts2 = round([x2',y2']);
- pts3 = round([x3',y3']);
- pts4 = round([x4',y4']);
-
- points.four_slice = cat(3,pts1,pts2,pts3,pts4);
-
- if l == 1 %compute the length of the layers. this gives us number of px per layers and we keep it same for inner and outer layer
- dx = [];
- dy = [];
- len0 = [];
- for s = 1:4
- dx = points.four_slice(1:end-1,1,s)-points.four_slice(2:end,1,s);
- dy = points.four_slice(1:end-1,2,s)-points.four_slice(2:end,2,s);
- len0(1,s) = sum(sqrt(dx.^2+dy.^2));
- end
- m_l = round(mean(len0));
- end
-
- p_new1 = distribute_points(pts1,'number',m_l,1);
- p_new2 = distribute_points(pts2,'number',m_l,1);
- p_new3 = distribute_points(pts3,'number',m_l,1);
- p_new4 = distribute_points(pts4,'number',m_l,1);
-
- slice_34 = slice_points_interpolation(p_new4,p_new3,...
- points.slice_no(4),points.slice_no(3),0);
- slice_23 = slice_points_interpolation(p_new3,p_new2,...
- points.slice_no(3),points.slice_no(2),0);
- slice_12 = slice_points_interpolation(p_new2,p_new1,...
- points.slice_no(2)-1,points.slice_no(1),0);
- points.all_slice = [slice_12;slice_23;slice_34];
-
+ m_l = round(mean(len0));
+ end
+
+ p_new1 = distribute_points(pts1,'number',m_l,1);
+ p_new2 = distribute_points(pts2,'number',m_l,1);
+ p_new3 = distribute_points(pts3,'number',m_l,1);
+
+ slice_23 = slice_points_interpolation(p_new3,p_new2,...
+ points.slice_no(3),points.slice_no(2),0);
+ slice_12 = slice_points_interpolation(p_new2,p_new1,...
+ points.slice_no(2)-1,points.slice_no(1),0);
+ points.all_slice = [slice_12;slice_23];
+
+ dx = [];
+ dy = [];
+ len1 = [];
+ for s = 1:size(points.all_slice,1)
+ dx = points.all_slice(s,1,1:end-1)-points.all_slice(s,1,2:end);
+ dy = points.all_slice(s,2,1:end-1)-points.all_slice(s,2,2:end);
+ len1(1,s) = sum(sqrt(dx.^2+dy.^2));
+ end
+ fprintf('Approximated length ~ %d or %d \n',m_l,round(mean(len1)));
+
+ case 4
+
+ fprintf('Go to slice %d and select points on the %s layer\n',points.slice_no(1),flag_layer)
+ input('');
+ [x1,y1] = click_points;
+ fprintf('Go to slice %d and select %d points on the %s layer\n',points.slice_no(2),length(x1),flag_layer)
+ input('');
+ [x2,y2] = click_points;
+ fprintf('Go to slice %d and select %d points on the %s layer\n',points.slice_no(3),length(x1),flag_layer)
+ input('');
+ [x3,y3] = click_points;
+ fprintf('Go to slice %d and select %d points on the %s layer\n',points.slice_no(4),length(x1),flag_layer)
+ input('');
+ [x4,y4] = click_points;
+
+ pts1 = round([x1',y1']);
+ pts2 = round([x2',y2']);
+ pts3 = round([x3',y3']);
+ pts4 = round([x4',y4']);
+
+ points.four_slice = cat(3,pts1,pts2,pts3,pts4);
+
+ if l == 1 %compute the length of the layers. this gives us number of px per layers and we keep it same for inner and outer layer
dx = [];
dy = [];
- len1 = [];
- for s = 1:size(points.all_slice,1)
- dx = points.all_slice(s,1,1:end-1)-points.all_slice(s,1,2:end);
- dy = points.all_slice(s,2,1:end-1)-points.all_slice(s,2,2:end);
- len1(1,s) = sum(sqrt(dx.^2+dy.^2));
+ len0 = [];
+ for s = 1:4
+ dx = points.four_slice(1:end-1,1,s)-points.four_slice(2:end,1,s);
+ dy = points.four_slice(1:end-1,2,s)-points.four_slice(2:end,2,s);
+ len0(1,s) = sum(sqrt(dx.^2+dy.^2));
end
- fprintf('Approximated length ~ %d or %d \n',m_l,round(mean(len1)));
- end
+ m_l = round(mean(len0));
+ end
+
+ p_new1 = distribute_points(pts1,'number',m_l,1);
+ p_new2 = distribute_points(pts2,'number',m_l,1);
+ p_new3 = distribute_points(pts3,'number',m_l,1);
+ p_new4 = distribute_points(pts4,'number',m_l,1);
+
+ slice_34 = slice_points_interpolation(p_new4,p_new3,...
+ points.slice_no(4),points.slice_no(3),0);
+ slice_23 = slice_points_interpolation(p_new3,p_new2,...
+ points.slice_no(3),points.slice_no(2),0);
+ slice_12 = slice_points_interpolation(p_new2,p_new1,...
+ points.slice_no(2)-1,points.slice_no(1),0);
+ points.all_slice = [slice_12;slice_23;slice_34];
+
+ dx = [];
+ dy = [];
+ len1 = [];
+ for s = 1:size(points.all_slice,1)
+ dx = points.all_slice(s,1,1:end-1)-points.all_slice(s,1,2:end);
+ dy = points.all_slice(s,2,1:end-1)-points.all_slice(s,2,2:end);
+ len1(1,s) = sum(sqrt(dx.^2+dy.^2));
+ end
+ fprintf('Approximated length ~ %d or %d \n',m_l,round(mean(len1)));
+ end
- %% save points
+ %% save points
+ if flag_save
save([save_folder,flag_layer,'_points_',sample,'_',dg],'points')
end
- close all
end
+close all
%% test the points by visualization
% for s = 1:points.slice_no(2)-points.slice_no(1)
diff --git a/PrepMeasure.m b/PrepMeasure.m
index 1370eef2b7f28d6ebfb441d97621e4083ecfd4ff..dbd7bad1e48ad2c54f8ae030713f2378fc93f6c3 100644
--- a/PrepMeasure.m
+++ b/PrepMeasure.m
@@ -1,7 +1,5 @@
%% measuring the distance betwee two sruface and length of each surface
-addpath(genpath('Functions'))
-
-flag_save = 1;
+addpath(genpath('.\Functions'))
save_folder = 'Results\measure\';
if ~exist(save_folder,'dir')
diff --git a/PrepRegistration.m b/PrepRegistration.m
index 211d329215aea0355f9430e7f95050a058461d17..0ae1eeb3f3572c52eaeafd6c466bfc27b7f8dd8a 100644
--- a/PrepRegistration.m
+++ b/PrepRegistration.m
@@ -1,7 +1,5 @@
%% summing over 10 pixels along normal direction of the smoothed surface
-addpath(genpath('Functions'))
-
-flag_save = 1;
+addpath(genpath('.\Functions'))
save_folder = 'Results\registration\';
if ~exist(save_folder,'dir')
diff --git a/Registration.m b/Registration.m
index 403980ce921c0215822481c06de776cbb98e0cd5..7460ebd479aac485f2c1a5a6fb3e47762cd0c603 100644
--- a/Registration.m
+++ b/Registration.m
@@ -36,7 +36,9 @@ end
Reg_dg45_outer = registerImages(dg45_outer,dg180_outer,30);
Reg_dg90_outer = registerImages(dg90_outer,dg180_outer,30);
- save([save_folder,'reg_tfrom_outer_',sample],'*_outer')
+ if flag_save
+ save([save_folder,'reg_tfrom_outer_',sample],'*_outer')
+ end
%end
@@ -86,8 +88,10 @@ end
subplot(223)
imshowpair(inner_90dg,inner.dg180_inner),title('90 to 180')
suptitle_modified(['inner using inner T ',sample],'none')
- saveas(gcf,[save_folder,'inner_using_inner_T_',sample,'.fig'])
- saveas(gcf,[save_folder,'inner_using_inner_T_',sample,'.png'])
+ if flag_save
+ saveas(gcf,[save_folder,'inner_using_inner_T_',sample,'.fig'])
+ saveas(gcf,[save_folder,'inner_using_inner_T_',sample,'.png'])
+ end
% check all outer agains 180 using inner transform T
MOVING = outer.dg0_outer;%imrotate(layers(:,:,1),180);
@@ -119,8 +123,10 @@ end
subplot(223)
imshowpair(outer_90dg,outer.dg180_outer),title('90 to 180')
suptitle_modified(['outer using inner T ',sample],'none')
- saveas(gcf,[save_folder,'outer_using_inner_T_',sample,'.fig'])
- saveas(gcf,[save_folder,'outer_using_inner_T_',sample,'.png'])
+ if flag_save
+ saveas(gcf,[save_folder,'outer_using_inner_T_',sample,'.fig'])
+ saveas(gcf,[save_folder,'outer_using_inner_T_',sample,'.png'])
+ end
%end
close all
@@ -238,34 +244,36 @@ end
subplot(236)
imshowpair(outer_90dg(:,:,s),outer.dg180_outer),title('90 to 180')
suptitle_modified(sample,'none')
-
- saveas(gcf,sprintf('%sReg_s%d_%s.png',save_folder,s,sample))
+ if flag_save
+ saveas(gcf,sprintf('%sReg_s%d_%s.png',save_folder,s,sample))
+ end
+ end
+ if flag_save
+ disp('Save results ...')
+ layers = inner_0dg;
+ Ref = inner_0dg_Ref;
+ save([save_folder,'Reg_inner_layers_',sample,'_0'],'layers','Ref')
+ layers = inner_45dg;
+ Ref = inner_45dg_Ref;
+ save([save_folder,'Reg_inner_layers_',sample,'_45'],'layers','Ref')
+ layers = inner_90dg;
+ Ref = inner_90dg_Ref;
+ save([save_folder,'Reg_inner_layers_',sample,'_90'],'layers','Ref')
+ layers = inner_180dg;
+ save([save_folder,'Reg_inner_layers_',sample,'_180'],'layers')
+
+ layers = outer_0dg;
+ Ref = outer_0dg_Ref;
+ save([save_folder,'Reg_outer_layers_',sample,'_0'],'layers','Ref')
+ layers = outer_45dg;
+ Ref = outer_45dg_Ref;
+ save([save_folder,'Reg_outer_layers_',sample,'_45'],'layers','Ref')
+ layers = outer_90dg;
+ Ref = outer_90dg_Ref;
+ save([save_folder,'Reg_outer_layers_',sample,'_90'],'layers','Ref')
+ layers = outer_180dg;
+ save([save_folder,'Reg_outer_layers_',sample,'_180'],'layers')
end
- disp('Save results ...')
- layers = inner_0dg;
- Ref = inner_0dg_Ref;
- save([save_folder,'Reg_inner_layers_',sample,'_0'],'layers','Ref')
- layers = inner_45dg;
- Ref = inner_45dg_Ref;
- save([save_folder,'Reg_inner_layers_',sample,'_45'],'layers','Ref')
- layers = inner_90dg;
- Ref = inner_90dg_Ref;
- save([save_folder,'Reg_inner_layers_',sample,'_90'],'layers','Ref')
- layers = inner_180dg;
- save([save_folder,'Reg_inner_layers_',sample,'_180'],'layers')
-
- layers = outer_0dg;
- Ref = outer_0dg_Ref;
- save([save_folder,'Reg_outer_layers_',sample,'_0'],'layers','Ref')
- layers = outer_45dg;
- Ref = outer_45dg_Ref;
- save([save_folder,'Reg_outer_layers_',sample,'_45'],'layers','Ref')
- layers = outer_90dg;
- Ref = outer_90dg_Ref;
- save([save_folder,'Reg_outer_layers_',sample,'_90'],'layers','Ref')
- layers = outer_180dg;
- save([save_folder,'Reg_outer_layers_',sample,'_180'],'layers')
-
%end
close all
@@ -363,25 +371,26 @@ end
clear tmp Reg
end
- disp('Save results ...')
- points = inner_0dg;
- save([save_folder,'Reg_inner_points_',sample,'_0'],'points')
- points = inner_45dg;
- save([save_folder,'Reg_inner_points_',sample,'_45'],'points')
- points = inner_90dg;
- save([save_folder,'Reg_inner_points_',sample,'_90'],'points')
- points = inner_180dg;
- save([save_folder,'Reg_inner_points_',sample,'_180'],'points')
-
- points = outer_0dg;
- save([save_folder,'Reg_outer_points_',sample,'_0'],'points')
- points = outer_45dg;
- save([save_folder,'Reg_outer_points_',sample,'_45'],'points')
- points = outer_90dg;
- save([save_folder,'Reg_outer_points_',sample,'_90'],'points')
- points = outer_180dg;
- save([save_folder,'Reg_outer_points_',sample,'_180'],'points')
-
+ if flag_save
+ disp('Save results ...')
+ points = inner_0dg;
+ save([save_folder,'Reg_inner_points_',sample,'_0'],'points')
+ points = inner_45dg;
+ save([save_folder,'Reg_inner_points_',sample,'_45'],'points')
+ points = inner_90dg;
+ save([save_folder,'Reg_inner_points_',sample,'_90'],'points')
+ points = inner_180dg;
+ save([save_folder,'Reg_inner_points_',sample,'_180'],'points')
+
+ points = outer_0dg;
+ save([save_folder,'Reg_outer_points_',sample,'_0'],'points')
+ points = outer_45dg;
+ save([save_folder,'Reg_outer_points_',sample,'_45'],'points')
+ points = outer_90dg;
+ save([save_folder,'Reg_outer_points_',sample,'_90'],'points')
+ points = outer_180dg;
+ save([save_folder,'Reg_outer_points_',sample,'_180'],'points')
+ end
%end
close all
diff --git a/Transformation.m b/Transformation.m
index 351e9a23474ff6e4c97d60924c74220ec7931364..0779d1c5f018665a510435f2ad04311c1895f2e0 100644
--- a/Transformation.m
+++ b/Transformation.m
@@ -1,4 +1,4 @@
-addpath(genpath('..\Functions'))
+addpath(genpath('.\Functions'))
samples = {'out_out_CD','out_out_MD','in_in_CD','in_in_MD'};
layer_name = {'outer','inner'};
@@ -17,9 +17,6 @@ if ~exist(save_folder,'dir')
mkdir(save_folder)
end
-flag_save = 1;
-flag_fig = 1;
-
% if you want to run the measure separately you can use the for loop
%for i = 1:4
sample = samples{i};
@@ -119,9 +116,6 @@ if ~exist(save_folder,'dir')
mkdir(save_folder)
end
-flag_save = 1;
-flag_fig = 1;
-
% if you want to run the measure separately you can use the for loop
%for i = 1:4
sample = samples{i};
@@ -193,9 +187,6 @@ flag_fig = 1;
%% move all the bend to the concave part for inner
-flag_save = 1;
-flag_fig = 1;
-
load_folder_layers = 'Results\layer\transfer\0dg_rest\';
save_folder = 'Results\layer\transfer\align_0dg_rest\';