Add current code for surfseg lib

CMakeLists.txt

+cmake_minimum_required (VERSION 3.5)
+
+project ("surfseg")
+
+set (GEL_LIB_ROOT_DIR ${PROJECT_SOURCE_DIR}/../GEL)
+
+set (CMAKE_CXX_STANDARD 14)
+
+if (UNIX AND NOT APPLE)
+	set (LINUX true)
+endif ()
+
+if (WIN32)
+	# For Windows we manually specify the location of the header- and lib directories
+    set (EXT_LIB_ROOT_DIR ${PROJECT_SOURCE_DIR}/../../../lib)
+	if (CMAKE_SIZEOF_VOID_P EQUAL 8)
+		# 64 bits
+		link_directories (${EXT_LIB_ROOT_DIR}/freeglut/lib/x64 ${EXT_LIB_ROOT_DIR}/glew-2.1.0/lib/Release/x64 ${GEL_LIB_ROOT_DIR}/GEL_WIN/x64/${CMAKE_BUILD_TYPE})
+	elseif (CMAKE_SIZEOF_VOID_P EQUAL 4)
+		# 32 bits
+		link_directories (${EXT_LIB_ROOT_DIR}/freeglut/lib ${EXT_LIB_ROOT_DIR}/glew-2.1.0/lib/Release/Win32 ${GEL_LIB_ROOT_DIR}/GEL_WIN/${CMAKE_BUILD_TYPE})
+	endif ()
+    include_directories (${EXT_LIB_ROOT_DIR}/freeglut/include ${EXT_LIB_ROOT_DIR}/glew-2.1.0/include)
+endif ()
+include_directories (${PROJECT_SOURCE_DIR}/surfseg/include ${PROJECT_SOURCE_DIR}/surfseg/src ${GEL_LIB_ROOT_DIR}/src)
+
+# Include sub-projects.
+add_subdirectory ("surfseg")

CMakeSettings.json

+{
+  "configurations": [
+    {
+      "name": "x64-Debug",
+      "generator": "Ninja",
+      "configurationType": "Debug",
+      "inheritEnvironments": [
+        "msvc_x64"
+      ],
+      "buildRoot": "${projectDir}\\surfseg\\build\\x64-Debug",
+      "cmakeCommandArgs": "",
+      "buildCommandArgs": "-v",
+      "ctestCommandArgs": ""
+    },
+    {
+      "name": "x64-Release",
+      "generator": "Ninja",
+      "configurationType": "Release",
+      "inheritEnvironments": [
+        "msvc_x64"
+      ],
+      "buildRoot": "${projectDir}\\surfseg\\build\\x64-Release",
+      "cmakeCommandArgs": "",
+      "buildCommandArgs": "-v",
+      "ctestCommandArgs": ""
+    }
+  ]
+}

surfseg/CMakeLists.txt

+add_subdirectory ("src")
+add_subdirectory ("matlab")
\ No newline at end of file

surfseg/include/block.h

+/* block.h */
+/* Vladimir Kolmogorov vnk@ist.ac.at */
+/* Last modified 08/05/2012 */
+/*
+	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)(const 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[] ((char*)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 ++;
+	}
+
+	struct iterator; // for overlapping scans
+	Type *ScanFirst(iterator& i)
+	{
+		for (i.scan_current_block=first; i.scan_current_block; i.scan_current_block = i.scan_current_block->next)
+		{
+			i.scan_current_data = & ( i.scan_current_block -> data[0] );
+			if (i.scan_current_data < i.scan_current_block -> current) return i.scan_current_data ++;
+		}
+		return NULL;
+	}
+	Type *ScanNext(iterator& i)
+	{
+		while (i.scan_current_data >= i.scan_current_block -> current)
+		{
+			i.scan_current_block = i.scan_current_block -> next;
+			if (!i.scan_current_block) return NULL;
+			i.scan_current_data = & ( i.scan_current_block -> data[0] );
+		}
+		return i.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;
+public:
+	struct iterator
+	{
+		block	*scan_current_block;
+		Type	*scan_current_data;
+	};
+private:
+	block	*scan_current_block;
+	Type	*scan_current_data;
+
+	void	(*error_function)(const 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)(const 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[] ((char*)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)(const char *);
+};
+
+
+#endif

surfseg/include/dsc.h

+#ifndef DSC_H__
+#define DSC_H__
+
+#include <vector>
+
+#include <GEL/CGLA/Vec3f.h>
+
+#include "simplex_mesh.h"
+
+using namespace CGLA;
+
+/** Deformable Simplicial Complex */
+class DSC : public SimplexMesh {
+public:
+    explicit DSC() = default;
+
+    float tetQuality(TetKey tk) const;
+    std::vector<float> allTetQuality() const;
+
+    void smartLaplacianSmooth(bool interfaceFlagsSet = true);
+    void smartLaplacianSmooth(std::vector<float>& tetQualities,
+        bool interfaceFlagsSet = true);
+};
+
+#endif // DSC_H__

surfseg/include/graph.h

+/* graph.h */
+/*
+    Copyright Vladimir Kolmogorov (vnk@ist.ac.at), Yuri Boykov (yuri@csd.uwo.ca) 
+	Modifications Copyright 2018 Niels Jeppesen (niejep@dtu.dk).
+
+    This file is part of MAXFLOW.
+
+    MAXFLOW is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    MAXFLOW is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with MAXFLOW.  If not, see <http://www.gnu.org/licenses/>.
+
+========================
+
+	version 3.04
+
+	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!!!
+
+
+enum termtype
+{
+	SOURCE	= 0,
+	SINK	= 1
+}; // terminals 
+typedef int node_id;
+
+// 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:
+
+	/////////////////////////////////////////////////////////////////////////
+	//                     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, long long edge_num_max, void (*err_function)(const 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) const;
+
+
+
+	//////////////////////////////////////////////
+	//       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; }
+	long long get_arc_num() { return (long long)(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 = false;
+	}
+
+
+
+
+
+
+/////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////
+	
+private:
+	// internal variables and functions
+
+# pragma pack (1)
+	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)
+		long long	TS;			// timestamp showing when DIST was computed
+		int			DIST;		// distance to the terminal
+
+		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 
+
+		bool		is_sink : true;	// flag showing whether the node is in the source or in the sink tree (if parent!=NULL)
+		bool		is_marked : true;	// set by mark_node()
+		bool		is_in_changed_list : true; // set by maxflow if 
+	};
+
+	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
+	};
+# pragma pack ()
+
+	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)(const 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
+	long long			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 node_id Graph<captype,tcaptype,flowtype>::add_node(int num)
+{
+	assert(num > 0);
+
+	if (node_last + num > node_max) reallocate_nodes(num);
+
+	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 termtype Graph<captype,tcaptype,flowtype>::what_segment(node_id i, termtype default_segm) const
+{
+	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 = true;
+}
+
+
+#endif
\ No newline at end of file

surfseg/include/manifold_mesh.h

+#ifndef MANIFOLD_MESH_H__
+#define MANIFOLD_MESH_H__
+
+#include <vector>
+#include <string>
+
+#include <GEL/CGLA/Vec3f.h>
+
+using namespace CGLA;
+
+/** Half-edge data structure for storing a manifold mesh */
+class ManifoldMesh {
+public:
+	typedef int VertKey;
+	typedef int FaceKey;
+	typedef int EdgeKey;
+
+	static const VertKey INVALID_VERT = -1;
+	static const FaceKey INVALID_FACE = -1;
+	static const EdgeKey INVALID_EDGE = -1;
+
+	struct Vertex {
+		VertKey self;
+		EdgeKey edge; // Edge pointing out from this vertex
+
+		Vertex() :
+			self(INVALID_VERT),
+			edge(INVALID_EDGE) {}
+	};
+
+	struct Face {
+		FaceKey self;
+		EdgeKey edge;
+
+		Face() :
+			self(INVALID_FACE),
+			edge(INVALID_EDGE) {}
+	};
+
+	struct HalfEdge {
+		EdgeKey self;
+		EdgeKey prev;
+		EdgeKey next;
+		EdgeKey twin;
+
+		VertKey vert; // Vertex this edge points out of
+		FaceKey face;
+
+		HalfEdge() :
+			self(INVALID_EDGE),
+			prev(INVALID_EDGE),
+			next(INVALID_EDGE),
+			twin(INVALID_EDGE),
+			vert(INVALID_VERT),
+			face(INVALID_FACE) {}
+	};
+
+	std::vector<Vertex> vertices;
+	std::vector<Face> faces;
+	std::vector<HalfEdge> edges;
+
+	std::vector<Vec3f> vertPositions;
+	std::vector<Vec3f> vertNormals;
+	std::vector<Vec3f> faceNormals;
+
+	explicit ManifoldMesh() = default;
+	ManifoldMesh(const ManifoldMesh& other);
+	ManifoldMesh(ManifoldMesh&& other);
+
+	ManifoldMesh& operator=(const ManifoldMesh& other);
+	ManifoldMesh& operator=(ManifoldMesh&& other);
+
+	void clear();
+	void build(std::vector<Vec3f> vertexPositions,
+		const std::vector<std::vector<int>>& faceIdxLists,
+		size_t expectedEdgeCount = 0);
+
+	void laplaceSmooth(float lambda, int niter = 1);
+	void taubinSmooth(float lambda, float mu, int niter = 1);
+
+	void splitEdge(EdgeKey ek, float a = 0.5f);
+	void flipEdge(EdgeKey ek);
+	void removeEdge(EdgeKey ek);
+
+	void removeInvalid();
+
+	float computeBoundingSphere() const noexcept;
+	float computeBoundingSphere(Vec3f& center) const noexcept;
+
+	std::vector<float> gaussCurvatures() const;
+
+	Vec3f faceNormal(FaceKey fk) const;
+	Vec3f faceNormal(const Face& f) const;
+	void computeFaceNormals(bool normalize = true);
+
+	Vec3f vertexNormal(VertKey vk, bool faceNormalsComputed = false) const;
+	Vec3f vertexNormal(const Vertex& v, bool faceNormalsComputed = false) const;
+	void computeVertexNormals(bool faceNormalsComputed = false, bool normalize = true);
+
+	void saveToObj(const std::string& fname) const;
+	void loadFromObj(const std::string& fname);
+
+#ifndef NO_OPENGL
+	void render(bool vertexNormals) const;
+#endif
+
+	EdgeKey next(EdgeKey ek) const;
+	EdgeKey prev(EdgeKey ek) const;
+	EdgeKey twin(EdgeKey ek) const;
+
+	Vec3f& vpos(const Vertex& v);
+	const Vec3f& vpos(const Vertex& v) const;
+	Vec3f& vpos(VertKey vk);
+	const Vec3f& vpos(VertKey vk) const;
+
+	Vec3f& vnormal(const Vertex& v);
+	const Vec3f& vnormal(const Vertex& v) const;
+	Vec3f& vnormal(VertKey vk);
+	const Vec3f& vnormal(VertKey vk) const;
+
+	Vec3f& fnormal(const Face& f);
+	const Vec3f& fnormal(const Face& f) const;
+	Vec3f& fnormal(FaceKey fk);
+	const Vec3f& fnormal(FaceKey fk) const;
+
+	void assertConsistent() const;
+};
+
+bool operator==(const ManifoldMesh& lhs, const ManifoldMesh& rhs);
+bool operator!=(const ManifoldMesh& lhs, const ManifoldMesh& rhs);
+
+bool operator==(const ManifoldMesh::Vertex& lhs, const ManifoldMesh::Vertex& rhs);
+bool operator!=(const ManifoldMesh::Vertex& lhs, const ManifoldMesh::Vertex& rhs);
+
+bool operator==(const ManifoldMesh::Face& lhs, const ManifoldMesh::Face& rhs);
+bool operator!=(const ManifoldMesh::Face& lhs, const ManifoldMesh::Face& rhs);
+
+bool operator==(const ManifoldMesh::HalfEdge& lhs, const ManifoldMesh::HalfEdge& rhs);
+bool operator!=(const ManifoldMesh::HalfEdge& lhs, const ManifoldMesh::HalfEdge& rhs);
+
+#endif // MANIFOLD_MESH_H__

surfseg/include/simplex_mesh.h

+#ifndef SIMPLEX_MESH_H__
+#define SIMPLEX_MESH_H__
+
+#include <vector>
+#include <unordered_map>
+#include <utility>
+#include <tuple>
+#include <string>
+
+#include <GEL/CGLA/Vec3f.h>
+
+using namespace CGLA;
+
+/** Incidence simplicial data structure for storing 3D simplicial complexes */
+class SimplexMesh {
+public:
+	typedef int VertKey;
+	typedef int EdgeKey;
+	typedef int TriKey;
+	typedef int TetKey;
+
+	static const VertKey INVALID_VERTEX = -1;
+	static const EdgeKey INVALID_EDGE = -1;
+	static const TriKey INVALID_TRIANGLE = -1;
+	static const TetKey INVALID_TETRAHEDRON = -1;
+
+	struct Vertex {
+		bool interf;
+		bool flag;
+		Vec3f pos;
+		VertKey self;
+		EdgeKey edge;
+
+		Vertex() :
+			interf(false),
+			flag(false),
+			pos(0.0f),
+			self(INVALID_VERTEX),
+			edge(INVALID_EDGE) {}
+	};
+
+	struct Edge {
+		bool interf;
+		bool flag;
+		EdgeKey self;
+		VertKey v0;
+		VertKey v1;
+		TriKey triangle;
+
+		Edge() :
+			interf(false),
+			flag(false),
+			self(INVALID_EDGE),
+			v0(INVALID_VERTEX),
+			v1(INVALID_VERTEX),
+			triangle(INVALID_TRIANGLE) {}
+	};
+
+	struct Triangle {
+		bool interf;
+		bool flag;
+		TriKey self;
+		std::array<EdgeKey, 3> edges;
+		TetKey t0;
+		TetKey t1;
+
+		Triangle() :
+			interf(false),
+			flag(false),
+			self(INVALID_TRIANGLE),
+			edges({ INVALID_EDGE, INVALID_EDGE, INVALID_EDGE }),
+			t0(INVALID_TETRAHEDRON),
+			t1(INVALID_TETRAHEDRON) {}
+	};
+
+	struct Tetrahedron {
+		int value;
+		bool flag;
+		TetKey self;
+		std::array<TriKey, 4> tris;
+
+		Tetrahedron() :
+			value(0),
+			flag(false),
+			self(INVALID_TETRAHEDRON),
+			tris({ INVALID_TRIANGLE, INVALID_TRIANGLE, INVALID_TRIANGLE, INVALID_TRIANGLE }) {}
+	};
+
+	std::vector<Vertex> vertices;
+	std::vector<Edge> edges;
+	std::vector<Triangle> triangles;
+	std::vector<Tetrahedron> tets;
+
+	explicit SimplexMesh() = default;
+	explicit SimplexMesh(const SimplexMesh& other);
+	explicit SimplexMesh(SimplexMesh&& other);
+
+	SimplexMesh& operator=(const SimplexMesh& other);
+	SimplexMesh& operator=(SimplexMesh&& other);
+
+	void clear();
+	void build(const std::vector<Vec3f>& vertexPositions,
+		const std::vector<std::array<int, 4>>& tetIdxs,
+		const std::vector<int> *tetValues = nullptr);
+
+	float computeBoundingSphere() const noexcept;
+	float computeBoundingSphere(Vec3f& center) const noexcept;
+
+	bool isTriInterface(TriKey tk, int insideVal = 1) const;
+	bool isTriInterface(const Triangle& t, int insideVal = 1) const;
+	void setInterfaceFlags();
+
+	void clearFlags();
+
+	void saveToDsc(const std::string& fname) const;
+	void loadFromDsc(const std::string& fname);
+
+	void renderInterface(int insideVal = 1) const;
+	void renderWireframe() const;
+
+	std::array<VertKey, 3> getTriVerts(TriKey tk) const;
+	std::array<VertKey, 4> getTetVerts(TetKey tk) const;
+
+	std::array<EdgeKey, 6> getTetEdges(TetKey tk) const;
+
+	std::tuple<std::vector<EdgeKey>, std::vector<TriKey>, std::vector<TetKey>>
+		getVertCoboundary(VertKey vk);
+	std::tuple<std::vector<VertKey>, std::vector<EdgeKey>, std::vector<TriKey>>
+		getVertLink(VertKey vk);
+
+	std::array<TetKey, 4> getAdjTets(TetKey tk) const;
+
+private:
+	EdgeKey addEdge_(int i1, int i2,
+		std::unordered_map<std::pair<int, int>, EdgeKey>& edgeMap);
+
+	TriKey addTriangle_(int i1, int i2, int i3,
+		std::unordered_map<std::array<int, 3>, TriKey>& triMap,
+		const std::unordered_map<std::pair<int, int>, EdgeKey>& edgeMap);
+
+	TetKey addTetrahedron_(int i1, int i2, int i3, int i4,
+		std::unordered_map<std::array<int, 4>, TetKey>& tetMap,
+		const std::unordered_map<std::array<int, 3>, TriKey>& triMap);
+};
+
+bool operator==(const SimplexMesh& lhs, const SimplexMesh& rhs);
+bool operator!=(const SimplexMesh& lhs, const SimplexMesh& rhs);
+
+bool operator==(const SimplexMesh::Vertex& lhs, const SimplexMesh::Vertex& rhs);
+bool operator!=(const SimplexMesh::Vertex& lhs, const SimplexMesh::Vertex& rhs);
+
+bool operator==(const SimplexMesh::Edge& lhs, const SimplexMesh::Edge& rhs);
+bool operator!=(const SimplexMesh::Edge& lhs, const SimplexMesh::Edge& rhs);
+
+bool operator==(const SimplexMesh::Triangle& lhs, const SimplexMesh::Triangle& rhs);
+bool operator!=(const SimplexMesh::Triangle& lhs, const SimplexMesh::Triangle& rhs);
+
+bool operator==(const SimplexMesh::Tetrahedron& lhs, const SimplexMesh::Tetrahedron& rhs);
+bool operator!=(const SimplexMesh::Tetrahedron& lhs, const SimplexMesh::Tetrahedron& rhs);
+
+#endif // SIMPLEX_MESH_H__

surfseg/include/subdivided_icosahedron.h

+#ifndef SUBDIVIDED_ICOSAHEDRON_H__
+#define SUBDIVIDED_ICOSAHEDRON_H__
+
+#include <GEL/CGLA/Vec3f.h>
+#include <GEL/CGLA/Mat3x3f.h>
+
+#include "manifold_mesh.h"
+
+using namespace CGLA;
+
+class SubdividedIcosahedron : public ManifoldMesh {
+    Vec3f center_;
+    float r_;
+    int divisionLevel_;
+
+public:
+    explicit SubdividedIcosahedron() :
+        ManifoldMesh(),
+        center_(0.0f),
+        r_(0.0f),
+        divisionLevel_(1) {};
+    explicit SubdividedIcosahedron(float r);
+    explicit SubdividedIcosahedron(const Vec3f& center, float r, int divLvl = 1);
+    explicit SubdividedIcosahedron(const Mat3x3f& R);
+    explicit SubdividedIcosahedron(const Vec3f& center, const Mat3x3f& R, int divLvl = 1);
+
+    void buildIcosahedron();
+    void buildIcosahedron(const Vec3f& center, float r);
+
+    void rescale(float newR) noexcept;
+    void move(const Vec3f& newCenter) noexcept;
+    void moveAndRescale(const Vec3f& newCenter, float newR) noexcept;
+
+    float r() const noexcept;
+    const Vec3f& center() const noexcept;
+
+    void subdivide(int numDivide = 1);
+
+private:
+    void singleSubdivide_();
+};
+
+#endif // SUBDIVIDED_ICOSAHEDRON_H__

surfseg/include/surface_segment.h

+#ifndef SURFACE_SEGMENT_H__
+#define SURFACE_SEGMENT_H__
+
+#include <vector>
+
+#include <GEL/CGLA/Vec3f.h>
+
+#include "manifold_mesh.h"
+#include "subdivided_icosahedron.h"
+#include "tet_mesh_4d.h"
+#include "volume.h"
+#include "volume4d.h"
+#include "graph.h"
+
+using namespace CGLA;
+
+typedef Graph<float, float, float> FloatGraph;
+
+enum CostType : int {
+	ON_SURFACE = 0,
+	IN_SURFACE = 1
+};
+
+ManifoldMesh surfaceCut(const Volume<float>& cost, const ManifoldMesh& init,
+    int numSamples, float sampleStep, int maxDiff, CostType costType);
+ManifoldMesh surfaceCut(const Volume<float>& cost, const ManifoldMesh&& init,
+    int numSamples, float sampleStep, int maxDiff, CostType costType);
+
+TetMesh4d surfaceCut4d(const Volume4d<float>& cost, TetMesh4d mesh,
+	int numSamples, float sampleStep, int maxDiff, CostType costType, bool bend = false);
+
+ManifoldMesh surfaceCutPlaneSep(const Volume<float>& cost, const ManifoldMesh& init,
+	int numSamples, float sampleStep, int maxDiff, CostType costType, std::vector<Vec3f> planeNormals,
+	std::vector<float> planeDists);
+ManifoldMesh surfaceCutPlaneSep(const Volume<float>& cost, const ManifoldMesh&& init,
+	int numSamples, float sampleStep, int maxDiff, CostType costType, std::vector<Vec3f> planeNormals,
+	std::vector<float> planeDists);
+
+std::vector<ManifoldMesh> surfaceCutPlaneSepQPBO(const Volume<float>& cost, std::vector<ManifoldMesh> meshes,
+	int numSamples, float sampleStep, int maxDiff, CostType costType, const std::vector<Vec3f>& centers,
+	const std::vector<std::vector<size_t>>& connections);
+
+std::vector<ManifoldMesh> surfaceCutPlaneSepDual(const Volume<float>& cost, std::vector<ManifoldMesh> meshes,
+	int numSamples, float sampleStep, int maxDiff, CostType costType, const std::vector<Vec3f>& centers,
+	const std::vector<std::vector<size_t>>& connections);
+std::vector<ManifoldMesh> surfaceCutPlaneSepDualParallel(const Volume<float>& cost, 
+	std::vector<ManifoldMesh> meshes, int numSamples, float sampleStep, int maxDiff, CostType costType,
+	const std::vector<Vec3f>& centers, const std::vector<std::vector<size_t>>& connections, int numThreads);
+
+std::vector<ManifoldMesh> kSurfaceCutNoOverlap(const std::vector<Volume<float>>& costs,
+	const ManifoldMesh& init, int numSamples, float sampleStep, int maxDiff, int minSep, int maxSep,
+	CostType costType);
+
+#endif // SURFACE_SEGMENT_H__

surfseg/include/tet_mesh_4d.h

+#ifndef TET_MESH_4D_H__
+#define TET_MESH_4D_H__
+
+#include <vector>
+#include <unordered_map>
+#include <utility>
+#include <tuple>
+#include <string>
+
+#include <GEL/CGLA/Vec4f.h>
+
+using namespace CGLA;
+
+/** Incidence simplicial data structure for storing 3D simplicial complexes in 4D space.
+  * NOTE: This code was adapted from the SimplexMesh class.
+  */
+class TetMesh4d {
+public:
+	typedef int VertKey;
+	typedef int EdgeKey;
+	typedef int TriKey;
+	typedef int TetKey;
+
+	static const VertKey INVALID_VERTEX = -1;
+	static const EdgeKey INVALID_EDGE = -1;
+	static const TriKey INVALID_TRIANGLE = -1;
+	static const TetKey INVALID_TETRAHEDRON = -1;
+
+	struct Vertex {
+		bool flag;
+		Vec4f pos;
+		Vec4f normal;
+		VertKey self;
+		EdgeKey edge;
+
+		Vertex() :
+			flag(false),
+			pos(0.0f),
+			normal(0.0f),
+			self(INVALID_VERTEX),
+			edge(INVALID_EDGE) {}
+	};
+
+	struct Edge {
+		bool flag;
+		EdgeKey self;
+		VertKey v0;
+		VertKey v1;
+		TriKey triangle;
+
+		Edge() :
+			flag(false),
+			self(INVALID_EDGE),
+			v0(INVALID_VERTEX),
+			v1(INVALID_VERTEX),
+			triangle(INVALID_TRIANGLE) {}
+	};
+
+	struct Triangle {
+		bool flag;
+		TriKey self;
+		std::array<EdgeKey, 3> edges;
+		TetKey t0;
+		TetKey t1;
+
+		Triangle() :
+			flag(false),
+			self(INVALID_TRIANGLE),
+			edges({ INVALID_EDGE, INVALID_EDGE, INVALID_EDGE }),
+			t0(INVALID_TETRAHEDRON),
+			t1(INVALID_TETRAHEDRON) {}
+	};
+
+	struct Tetrahedron {
+		bool flag;
+		TetKey self;
+		std::array<TriKey, 4> tris;
+		std::array<VertKey, 4> verts;
+		Vec4f normal;
+
+		Tetrahedron() :
+			flag(false),
+			self(INVALID_TETRAHEDRON),
+			tris({ INVALID_TRIANGLE, INVALID_TRIANGLE, INVALID_TRIANGLE, INVALID_TRIANGLE }),
+			verts({ INVALID_VERTEX, INVALID_VERTEX, INVALID_VERTEX, INVALID_VERTEX }),
+			normal(0.0f) {}
+	};
+
+	// TODO: Store flags in seperate arrays so we can have proper const qualifiers
+	// TODO: Store positions and normals in seperate vectors like ManifoldMesh
+	// TODO: Profile speed-up of above suggestion
+	std::vector<Vertex> vertices;
+	std::vector<Edge> edges;
+	std::vector<Triangle> triangles;
+	std::vector<Tetrahedron> tets;
+
+	explicit TetMesh4d() = default;
+	explicit TetMesh4d(const std::vector<Vec4f>& vertexPositions,
+		const std::vector<std::array<int, 4>>& tetIdxs);
+	TetMesh4d(const TetMesh4d& other);
+	TetMesh4d(TetMesh4d&& other);
+
+	TetMesh4d& operator=(const TetMesh4d& other);
+	TetMesh4d& operator=(TetMesh4d&& other);
+
+	void clear();
+	void build(const std::vector<Vec4f>& vertexPositions,
+		const std::vector<std::array<int, 4>>& tetIdxs);
+
+	void laplaceSmooth(float lambda, int niter = 1);
+
+	Vec4f computeCenter() const noexcept;
+	float computeBoundingSphere() const noexcept;
+	float computeBoundingSphere(Vec4f& center) const noexcept;
+
+	void clearFlags();
+
+	void saveToDsc(const std::string& fname) const;
+	void loadFromDsc(const std::string& fname);
+
+	Vec4f tetNormal(TetKey tk) const;
+	Vec4f tetNormal(const Tetrahedron& t) const;
+	void computeTetNormals(bool normalize = true);
+
+	Vec4f vertexNormal(VertKey vk, bool tetNormalsComputed = false);
+	Vec4f vertexNormal(const Vertex& v, bool tetNormalsComputed = false);
+	void computeVertexNormals(bool normalize = true, bool tetNormalsComputed = false);
+
+	std::array<VertKey, 3> getTriVerts(TriKey tk) const;
+
+	std::array<EdgeKey, 6> getTetEdges(TetKey tk) const;
+
+	std::tuple<std::vector<EdgeKey>, std::vector<TriKey>, std::vector<TetKey>>
+		getVertCoboundary(VertKey vk);
+	std::tuple<std::vector<VertKey>, std::vector<EdgeKey>, std::vector<TriKey>>
+		getVertLink(VertKey vk);
+
+	std::array<TetKey, 4> getAdjTets(TetKey tk) const;
+
+private:
+	EdgeKey addEdge_(int i1, int i2,
+		std::unordered_map<std::pair<int, int>, EdgeKey>& edgeMap);
+
+	TriKey addTriangle_(int i1, int i2, int i3,
+		std::unordered_map<std::array<int, 3>, TriKey>& triMap,
+		const std::unordered_map<std::pair<int, int>, EdgeKey>& edgeMap);
+
+	TetKey addTetrahedron_(int i1, int i2, int i3, int i4,
+		std::unordered_map<std::array<int, 4>, TetKey>& tetMap,
+		const std::unordered_map<std::array<int, 3>, TriKey>& triMap);
+};
+
+bool operator==(const TetMesh4d& lhs, const TetMesh4d& rhs);
+bool operator!=(const TetMesh4d& lhs, const TetMesh4d& rhs);
+
+bool operator==(const TetMesh4d::Vertex& lhs, const TetMesh4d::Vertex& rhs);
+bool operator!=(const TetMesh4d::Vertex& lhs, const TetMesh4d::Vertex& rhs);
+
+bool operator==(const TetMesh4d::Edge& lhs, const TetMesh4d::Edge& rhs);
+bool operator!=(const TetMesh4d::Edge& lhs, const TetMesh4d::Edge& rhs);
+
+bool operator==(const TetMesh4d::Triangle& lhs, const TetMesh4d::Triangle& rhs);
+bool operator!=(const TetMesh4d::Triangle& lhs, const TetMesh4d::Triangle& rhs);
+
+bool operator==(const TetMesh4d::Tetrahedron& lhs, const TetMesh4d::Tetrahedron& rhs);
+bool operator!=(const TetMesh4d::Tetrahedron& lhs, const TetMesh4d::Tetrahedron& rhs);
+
+#endif // TET_MESH_4D_H__

surfseg/include/util.h

+#ifndef UTIL_H__
+#define UTIL_H__
+
+#include <utility>
+#include <functional>
+#include <array>
+#include <limits>
+#include <type_traits>
+
+namespace std {
+    // We need to make our hash specialization since std::pair does not have one
+    template <>
+    struct hash<std::pair<int, int>> {
+        size_t operator()(const std::pair<int, int>& key) const;
+    };
+
+    template <size_t N>
+    struct hash<std::array<int, N>> {
+        size_t operator()(const std::array<int, N>& key) const
+        {
+            // Code from https://codereview.stackexchange.com/q/171999
+            static std::hash<int> hasher;
+            size_t result = 0;
+            for (int k : key) {
+                result = result * 31 + hasher(k);
+            }
+            return result;
+        }
+    };
+} // std
+
+void graphErrFunc(const char *msg);
+
+template <class Ty>
+constexpr std::enable_if_t<std::is_floating_point<Ty>::value, Ty> infOrMaxF() {
+    return std::numeric_limits<Ty>::infinity();
+}
+
+template <class Ty>
+constexpr std::enable_if_t<!std::is_floating_point<Ty>::value, Ty> infOrMaxF() {
+    return std::numeric_limits<Ty>::max();
+}
+
+template <class Ty>
+constexpr Ty infOrMax = infOrMaxF<Ty>();
+
+#endif // UTIL_H__

surfseg/include/volume.h

+#ifndef VOLUME_H__
+#define VOLUME_H__
+
+#include <cstdint>
+#include <memory>
+
+#include <GEL/CGLA/Vec3i.h>
+#include <GEL/CGLA/Vec3f.h>
+
+using namespace CGLA;
+
+/** Simple structure to hold volumetric data of type Ty */
+template <class Ty>
+struct Volume {
+    typedef Ty DataType;
+
+    std::shared_ptr<Ty> data;
+    size_t nx;
+    size_t ny;
+    size_t nz;
+
+    explicit Volume() :
+        data(),
+        nx(0),
+        ny(0),
+        nz(0) {}
+    explicit Volume(size_t nx, size_t ny, size_t nz) noexcept;
+    Volume(const Volume&) noexcept;
+    Volume(Volume&& other) noexcept;
+
+    Volume& operator=(const Volume&) noexcept;
+    Volume& operator=(Volume&& other) noexcept;
+
+    void alloc();
+    void alloc(size_t nelem);
+    void alloc(size_t nx, size_t ny, size_t nz);
+
+    void clear();
+
+	void copy(const Volume<Ty>& src);
+	void copy(Volume<Ty>&& src) noexcept;
+
+    size_t numElem() const noexcept;
+
+    Ty interp(float x, float y, float z) const;
+    Ty interp(const Vec3f& p) const;
+
+    size_t idx(size_t x, size_t y, size_t z) const noexcept;
+    Vec3i pos(size_t idx) const noexcept;
+	bool contains(const Vec3f& p) const noexcept;
+
+    Ty& at(size_t i);
+    const Ty& at(size_t i) const;
+    Ty& at(size_t x, size_t y, size_t z);
+    const Ty& at(size_t x, size_t y, size_t z) const;
+    Ty& at(const Vec3i& p);
+    const Ty& at(const Vec3i& p) const;
+
+    Ty& operator[](size_t i);
+    const Ty& operator[](size_t i) const;
+    Ty& operator[](const Vec3i& p);
+    const Ty& operator[](const Vec3i& p) const;
+
+    void print() const;
+    void saveToBin(const std::string& fname) const;
+    void loadFromBin(const std::string& fname);
+};
+
+template <class Ty>
+constexpr uint32_t getTypeCode();
+
+#include <volume.inl>
+
+#endif // VOLUME_H__

surfseg/include/volume4d.h

+#ifndef VOLUME_4D_H__
+#define VOLUME_4D_H__
+
+#include <memory>
+
+#include <GEL/CGLA/Vec4i.h>
+#include <GEL/CGLA/Vec4f.h>
+#include <GEL/CGLA/Vec3i.h>
+
+using namespace CGLA;
+
+/** Simple structure to hold time series volumetric data of type Ty.
+ * Adapted from the Volume class.
+ */
+template <class Ty>
+struct Volume4d {
+	typedef Ty DataType;
+
+	std::shared_ptr<Ty> data;
+	size_t nx;
+	size_t ny;
+	size_t nz;
+	size_t nt;
+
+	explicit Volume4d() :
+		data(),
+		nx(0),
+		ny(0),
+		nz(0) {}
+	explicit Volume4d(size_t nx, size_t ny, size_t nz, size_t nt) noexcept;
+	Volume4d(const Volume4d&) noexcept;
+	Volume4d(Volume4d&& other) noexcept;
+
+	Volume4d& operator=(const Volume4d&) noexcept;
+	Volume4d& operator=(Volume4d&& other) noexcept;
+
+	void alloc();
+	void alloc(size_t nelem);
+	void alloc(size_t nx, size_t ny, size_t nz, size_t nt);
+
+	void clear();
+
+	void copy(const Volume4d<Ty>& src);
+	void copy(Volume4d<Ty>&& src) noexcept;
+
+	size_t numElem() const noexcept;
+
+	Ty interp3d(float x, float y, float z, size_t t) const;
+	Ty interp3d(const Vec3f& p, size_t t) const;
+
+	Ty interp(float x, float y, float z, float t) const;
+	Ty interp(const Vec4f& p) const;
+	
+	size_t idx(size_t x, size_t y, size_t z, size_t t) const noexcept;
+	Vec4i pos(size_t idx) const noexcept;
+	bool contains(const Vec4f& p) const noexcept;
+
+	Ty& at(size_t i);
+	const Ty& at(size_t i) const;
+	Ty& at(size_t x, size_t y, size_t z, size_t t);
+	const Ty& at(size_t x, size_t y, size_t z, size_t t) const;
+	Ty& at(const Vec4i& p);
+	const Ty& at(const Vec4i& p) const;
+
+	Ty& operator[](size_t i);
+	const Ty& operator[](size_t i) const;
+	Ty& operator[](const Vec4i& p);
+	const Ty& operator[](const Vec4i& p) const;
+};
+
+#include <volume4d.inl>
+
+#endif // VOLUME_H__

surfseg/matlab/CMakeLists.txt

+find_package(Matlab)
+
+set (MESH_LIB_ROOT_DIR ${PROJECT_SOURCE_DIR}/../MESH)
+
+if (Matlab_FOUND)
+	include_directories (${Matlab_INCLUDE_DIRS})
+	include_directories (${PROJECT_SOURCE_DIR}/Demo/include ${PROJECT_SOURCE_DIR}/Demo/src)
+	include_directories (${GEL_LIB_ROOT_DIR}/src ${MESH_LIB_ROOT_DIR} ${MESH_LIB_ROOT_DIR}/lib3d/include)
+	if (WIN32)
+		link_directories (${MESH_LIB_ROOT_DIR}/build/x64-${CMAKE_BUILD_TYPE})
+	else ()
+		link_directories (${MESH_LIB_ROOT_DIR}/build)
+	endif()
+
+	set (CMAKE_POSITION_INDEPENDENT_CODE ON)
+
+	# 3D Surface cuts
+	matlab_add_mex (NAME mex_surfcut SRC mex_surfcut.cpp matlab_util.cpp LINK_TO surfseg)
+	matlab_add_mex (NAME mex_ksurfcut SRC mex_ksurfcut.cpp matlab_util.cpp LINK_TO surfseg)
+	matlab_add_mex (NAME mex_surfcut_planesep SRC mex_surfcut_planesep.cpp matlab_util.cpp LINK_TO surfseg)
+	matlab_add_mex (NAME mex_surfcut_planesep_dual SRC mex_surfcut_planesep_dual.cpp matlab_util.cpp LINK_TO surfseg)
+	matlab_add_mex (NAME mex_surfcut_planesep_qpbo SRC mex_surfcut_planesep_qpbo.cpp matlab_util.cpp LINK_TO surfseg)
+	
+	# 4D surface cuts
+	matlab_add_mex (NAME mex_surfcut_4d SRC mex_surfcut_4d.cpp matlab_util.cpp LINK_TO surfseg)
+
+	# Mesh utils
+	matlab_add_mex (NAME mex_gauss_curvature SRC mex_gauss_curvature.cpp matlab_util.cpp LINK_TO surfseg)
+	matlab_add_mex (NAME mex_star_intersect SRC mex_star_intersect.cpp matlab_util.cpp LINK_TO surfseg)
+	matlab_add_mex (NAME mex_hausdorff SRC mex_hausdorff.cpp matlab_util.cpp LINK_TO surfseg mesh)
+	matlab_add_mex (NAME mex_subdiv_icosahedron SRC mex_subdiv_icosahedron.cpp matlab_util.cpp LINK_TO surfseg mesh)
+	
+	if (WIN32)
+		# For some reason we have to specify this manually on Windows
+		foreach (EXE mex_surfcut mex_ksurfcut mex_surfcut_planesep mex_surfcut_planesep_dual mex_surfcut_planesep_qpbo mex_surfcut_4d mex_gauss_curvature mex_star_intersect mex_hausdorff mex_subdiv_icosahedron)
+			target_link_libraries (${EXE} ${Matlab_ROOT_DIR}/extern/lib/win64/microsoft/libmat.lib ${Matlab_ROOT_DIR}/extern/lib/win64/microsoft/libmx.lib ${Matlab_ROOT_DIR}/extern/lib/win64/microsoft/libut.lib)
+			set_target_properties (${EXE} PROPERTIES LINK_FLAGS "/export:mexFunction /debug")
+		endforeach ()
+	endif ()
+endif ()
\ No newline at end of file

surfseg/matlab/compile_mex.m

+mex mex_surfcut.cpp matlab_util.cpp -I../include -I../src -I../../../GEL/src CXXFLAGS="-std=c++14 -fPIC" -L../../build/Demo/src -lsurfseg
+mex mex_ksurfcut.cpp matlab_util.cpp -I../include -I../src -I../../../GEL/src CXXFLAGS="-std=c++14 -fPIC" -L../../build/Demo/src -lsurfseg
+mex mex_surfcut_planesep.cpp matlab_util.cpp -I../include -I../src -I../../../GEL/src CXXFLAGS="-std=c++14 -fPIC" -L../../build/Demo/src -lsurfseg
+mex mex_surfcut_planesep_dual.cpp matlab_util.cpp -I../include -I../src -I../../../GEL/src CXXFLAGS="-std=c++14 -fPIC" -L../../build/Demo/src -lsurfseg
+mex mex_surfcut_planesep_qpbo.cpp matlab_util.cpp -I../include -I../src -I../../../GEL/src CXXFLAGS="-std=c++14 -fPIC" -L../../build/Demo/src -lsurfseg
+
+mex mex_surfcut_4d.cpp matlab_util.cpp -I../include -I../src -I../../../GEL/src CXXFLAGS="-std=c++14 -fPIC" -L../../build/Demo/src -lsurfseg
+
+mex mex_gauss_curvature.cpp matlab_util.cpp -I../include -I../src -I../../../GEL/src CXXFLAGS="-std=c++14 -fPIC" -L../../build/Demo/src -lsurfseg
+mex mex_star_intersect.cpp matlab_util.cpp -I../include -I../src -I../../../GEL/src CXXFLAGS="-std=c++14 -fPIC" -L../../build/Demo/src -lsurfseg
+mex mex_subdiv_icosahedron.cpp matlab_util.cpp -I../include -I../src -I../../../GEL/src CXXFLAGS="-std=c++14 -fPIC" -L../../build/Demo/src -lsurfseg
+mex mex_hausdorff.cpp matlab_util.cpp -I../include -I../src -I../../../GEL/src -I../../../MESH -I../../../MESH/lib3d/include CXXFLAGS="-std=c++14 -fPIC" -LDFLAGS="-fPIC" -L../../build/Demo/src -L../../../MESH/build -lsurfseg -lmesh

surfseg/matlab/matlab_util.cpp

+#include "matlab_util.h"
+#include <algorithm>
+
+bool isSize(const mxArray *a, std::initializer_list<int> size)
+{
+	size_t dimA = mxGetNumberOfDimensions(a);
+	int i = 0;
+	for (int s : size) {
+		if (s != -1 && (i >= dimA || s != mxGetDimensions(a)[i])) {
+			return false;
+		}
+		++i;
+	}
+	return true;
+}
+
+int getMaxCompThreads()
+{
+	mxArray *matlabCallOut[1] = { 0 };
+	mxArray *matlabCallIn[1] = { 0 };
+	mexCallMATLAB(1, matlabCallOut, 0, matlabCallIn, "maxNumCompThreads");
+	double *Nthreadsd = mxGetPr(matlabCallOut[0]);
+	return (int)Nthreadsd[0];
+}

surfseg/matlab/matlab_util.h

+#ifndef MATLAB_UTIL_H__
+#define MATLAB_UTIL_H__
+
+#include <string>
+#include <memory>
+#include <typeinfo>
+#include <cstdint>
+#include <initializer_list>
+#include <type_traits>
+
+#include "mex.h"
+#include "matrix.h"
+
+#include <GEL/CGLA/ArithVec.h>
+
+#include "volume.h"
+#include "volume4d.h"
+
+#define ensureOrError(expr, ...) do { \
+	if (!(expr)) { \
+		mexErrMsgIdAndTxt("surfseg:internal", __VA_ARGS__); \
+	} \
+} while(false)
+
+#ifdef __cplusplus
+extern "C" bool utIsInterruptPending();
+#else
+extern bool utIsInterruptPending();
+#endif
+
+template <class Ty>
+constexpr mxClassID getTypeClassId()
+{
+	// We cannot use switch so have to revert to this
+	if (std::is_same<Ty, mxLogical>::value) {
+		return mxLOGICAL_CLASS;
+	} else if (std::is_same<Ty, mxInt8>::value) {
+		return mxINT8_CLASS;
+	} else if (std::is_same<Ty, mxUint8>::value) {
+		return mxUINT8_CLASS;
+	} else if (std::is_same<Ty, mxInt16>::value) {
+		return mxINT16_CLASS;
+	} else if (std::is_same<Ty, mxUint16>::value) {
+		return mxUINT16_CLASS;
+	} else if (std::is_same<Ty, mxInt32>::value) {
+		return mxINT32_CLASS;
+	} else if (std::is_same<Ty, mxUint32>::value) {
+		return mxUINT32_CLASS;
+	} else if (std::is_same<Ty, mxSingle>::value) {
+		return mxSINGLE_CLASS;
+	} else if (std::is_same<Ty, mxDouble>::value) {
+		return mxDOUBLE_CLASS;
+	} else if (std::is_same<Ty, mxInt64>::value) {
+		return mxINT64_CLASS;
+	} else if (std::is_same<Ty, mxUint64>::value) {
+		return mxUINT64_CLASS;
+	} else {
+		return mxUNKNOWN_CLASS;
+	}
+}
+
+constexpr const char *getClassNameFromId(mxClassID id)
+{
+	switch (id) {
+	case mxLOGICAL_CLASS:
+		return "logical";
+	case mxINT8_CLASS:
+		return "int8";
+	case mxUINT8_CLASS:
+		return "uint8";
+	case mxINT16_CLASS:
+		return "int16";
+	case mxUINT16_CLASS:
+		return "uint16";
+	case mxINT32_CLASS:
+		return "int32";
+	case mxUINT32_CLASS:
+		return "uint32";
+	case mxSINGLE_CLASS:
+		return "single";
+	case mxDOUBLE_CLASS:
+		return "double";
+	case mxINT64_CLASS:
+		return "int64";
+	case mxUINT64_CLASS:
+		return "uint64";
+	default:
+		return "unkown";
+	}
+}
+
+template <class Ty>
+inline Ty derefAndCast(const void *ptr, mxClassID id)
+{
+	switch (id) {
+	case mxLOGICAL_CLASS:
+		return static_cast<Ty>(*reinterpret_cast<const mxLogical *>(ptr));
+	case mxINT8_CLASS:
+		return static_cast<Ty>(*reinterpret_cast<const mxInt8 *>(ptr));
+	case mxUINT8_CLASS:
+		return static_cast<Ty>(*reinterpret_cast<const mxUint8 *>(ptr));
+	case mxINT16_CLASS:
+		return static_cast<Ty>(*reinterpret_cast<const mxInt16 *>(ptr));
+	case mxUINT16_CLASS:
+		return static_cast<Ty>(*reinterpret_cast<const mxUint16 *>(ptr));
+	case mxINT32_CLASS:
+		return static_cast<Ty>(*reinterpret_cast<const mxInt32 *>(ptr));
+	case mxUINT32_CLASS:
+		return static_cast<Ty>(*reinterpret_cast<const mxUint32 *>(ptr));
+	case mxSINGLE_CLASS:
+		return static_cast<Ty>(*reinterpret_cast<const mxSingle *>(ptr));
+	case mxDOUBLE_CLASS:
+		return static_cast<Ty>(*reinterpret_cast<const mxDouble *>(ptr));
+	case mxINT64_CLASS:
+		return static_cast<Ty>(*reinterpret_cast<const mxInt64 *>(ptr));
+	case mxUINT64_CLASS:
+		return static_cast<Ty>(*reinterpret_cast<const mxUint64 *>(ptr));
+	default:
+		mexErrMsgIdAndTxt("surfseg:internal", "derefAndCast: Unknown class ID");
+		return Ty(); // NOTE: This will never run as mexErrMsgIdAndTxt aborts execution
+	}
+}
+
+inline bool isVector(const mxArray *a)
+{
+	return mxGetNumberOfDimensions(a) == 2 && (mxGetN(a) == 1 || mxGetM(a) == 1);
+}
+
+inline bool isMatrix(const mxArray *a)
+{
+	return mxGetNumberOfDimensions(a) == 2;
+}
+
+bool isSize(const mxArray *a, std::initializer_list<int> size);
+
+template <class Ty>
+inline void ensureMatchingClass(const mxArray *a, const std::string& errPrefixStr = "Value")
+{
+	const char *errPrefix = errPrefixStr.c_str();
+	ensureOrError(mxGetClassID(a) == getTypeClassId<Ty>(), "%s must be of class %s (was: %s)",
+		errPrefix, getClassNameFromId(getTypeClassId<Ty>()), mxGetClassName(a));
+}
+
+template <class Ty>
+inline Ty getCastScalarChecked(const mxArray *a, const std::string& errPrefixStr = "Value")
+{
+	const char *errPrefix = errPrefixStr.c_str();
+
+	ensureOrError(!mxIsComplex(a), "%s must be real", errPrefix);
+	ensureOrError(mxIsScalar(a), "%s must be a scalar", errPrefix);
+
+	return static_cast<Ty>(mxGetScalar(a));
+}
+
+template <class Ty, unsigned int N>
+inline std::array<Ty, N> getVectorChecked(const mxArray *a, const std::string& errPrefixStr = "Value")
+{
+	mxClassID classId = mxGetClassID(a);
+	const char *errPrefix = errPrefixStr.c_str();
+
+	ensureOrError(!mxIsComplex(a), "%s must be real", errPrefix);
+	ensureMatchingClass<Ty>(a, errPrefixStr);
+	ensureOrError(isVector(a), "%s must be a vector", errPrefix);
+	ensureOrError(mxGetNumberOfElements(a) == N, "%s must have length %d", errPrefix, N);
+
+	std::array<Ty, N> out;
+	const Ty *data = static_cast<Ty *>(mxGetData(a));
+	for (int i = 0; i < N; ++i) {
+		out[i] = data[i];
+	}
+	return out;
+}
+
+template <class VecTy>
+inline VecTy getVectorChecked(const mxArray *a, const std::string& errPrefixStr = "Value")
+{
+	typedef typename VecTy::ScalarType Ty;
+	const size_t N = VecTy::get_dim();
+	mxClassID classId = mxGetClassID(a);
+	const char *errPrefix = errPrefixStr.c_str();
+
+	ensureOrError(!mxIsComplex(a), "%s must be real", errPrefix);
+	ensureMatchingClass<Ty>(a, errPrefixStr);
+	ensureOrError(isVector(a), "%s must be a vector", errPrefix);
+	ensureOrError(mxGetNumberOfElements(a) == N, "%s must have length %d", errPrefix, N);
+
+	VecTy out;
+	const Ty *data = static_cast<Ty *>(mxGetData(a));
+	for (int i = 0; i < N; ++i) {
+		out[i] = data[i];
+	}
+	return out;
+}
+
+template <class VecTy>
+inline VecTy getCastVectorChecked(const mxArray *a, const std::string& errPrefixStr = "Value")
+{
+	typedef typename VecTy::ScalarType Ty;
+	const size_t N = VecTy::get_dim();
+	mxClassID classId = mxGetClassID(a);
+	const char *errPrefix = errPrefixStr.c_str();
+
+	ensureOrError(mxIsNumeric(a), "%s must be numeric", errPrefix);
+	ensureOrError(!mxIsComplex(a), "%s must be real", errPrefix);
+	ensureOrError(isVector(a), "%s must be a vector", errPrefix);
+	ensureOrError(mxGetNumberOfElements(a) == N, "%s must have length %d", errPrefix, N);
+
+	VecTy out;
+	const std::uint8_t *data = static_cast<const std::uint8_t *>(mxGetData(a));
+	for (int i = 0; i < N; ++i) {
+		out[i] = derefAndCast<Ty>(data, classId);
+		data += mxGetElementSize(a);
+	}
+	return out;
+}
+
+template <class Ty>
+inline Volume<Ty> getVolumeChecked(const mxArray *a, const std::string& errPrefixStr = "Value")
+{
+	mxClassID classId = mxGetClassID(a);
+	const char *errPrefix = errPrefixStr.c_str();
+
+	ensureOrError(!mxIsComplex(a), "%s must be real", errPrefix);
+	ensureMatchingClass<Ty>(a, errPrefixStr);
+	ensureOrError(mxGetNumberOfDimensions(a) <= 3, "%s have must at most 3 dimensions", errPrefix);
+
+	const mwSize *dims = mxGetDimensions(a);
+	Volume<Ty> out;
+	out.nx = dims[0];
+	out.ny = dims[1];
+	out.nz = mxGetNumberOfDimensions(a) > 2 ? dims[2] : 1;
+	// Make a shared pointer with no deleter, since the memory comes from MATLAB
+	out.data = std::shared_ptr<Ty>(static_cast<Ty *>(mxGetData(a)), [](auto) {});
+	return out;
+}
+
+template <class Ty>
+inline Volume4d<Ty> getVolume4dChecked(const mxArray *a, const std::string& errPrefixStr = "Value")
+{
+	// TODO: Merge with getVolumeChecked
+	mxClassID classId = mxGetClassID(a);
+	const char *errPrefix = errPrefixStr.c_str();
+
+	ensureOrError(!mxIsComplex(a), "%s must be real", errPrefix);
+	ensureMatchingClass<Ty>(a, errPrefixStr);
+	ensureOrError(mxGetNumberOfDimensions(a) <= 4, "%s have must at most 4 dimensions", errPrefix);
+
+	const mwSize *dims = mxGetDimensions(a);
+	Volume4d<Ty> out;
+	out.nx = dims[0];
+	out.ny = dims[1];
+	out.nz = mxGetNumberOfDimensions(a) > 2 ? dims[2] : 1;
+	out.nt = mxGetNumberOfDimensions(a) > 3 ? dims[3] : 1;
+	// Make a shared pointer with no deleter, since the memory comes from MATLAB
+	out.data = std::shared_ptr<Ty>(static_cast<Ty *>(mxGetData(a)), [](auto) {});
+	return out;
+}
+
+template <class Ty>
+inline std::shared_ptr<Ty> getCastSharedPtrChecked(const mxArray *a, const std::string& errPrefixStr = "Value")
+{
+	mxClassID classId = mxGetClassID(a);
+	const char *errPrefix = errPrefixStr.c_str();
+
+	ensureOrError(mxIsNumeric(a), "%s must be numeric", errPrefix);
+	ensureOrError(!mxIsComplex(a), "%s must be real", errPrefix);
+
+	if (classId == getTypeClassId<Ty>()) {
+		// Type of MATLAB array already matches type, so just use it
+		return std::shared_ptr<Ty>(static_cast<Ty *>(mxGetData(a)), [](auto) {});
+	} else {
+		// Type does *not* match so alloc an array to hold and cast all values
+		size_t n = mxGetNumberOfElements(a);
+		std::shared_ptr<Ty> out(static_cast<Ty *>(mxMalloc(n * sizeof(Ty))), [](auto ptr) { mxFree(ptr); });
+		Ty *outPtr = out.get();
+		// Need to use uint8_t pointer so we can do pointer arithmetic
+		const std::uint8_t *data = static_cast<const std::uint8_t *>(mxGetData(a));
+		for (int i = 0; i < n; ++i) {
+			outPtr[i] = derefAndCast<Ty>(data, classId);
+			data += mxGetElementSize(a);
+		}
+		return out;
+	}
+}
+
+int getMaxCompThreads();
+
+#endif // MATLAB_UTIL_H__

surfseg/matlab/mex_gauss_curvature.cpp

+#include <vector>
+#include <cstring>
+
+#include "mex.h"
+#include "matrix.h"
+
+#include <GEL/CGLA/Vec3f.h>
+
+#include "volume.h"
+#include "manifold_mesh.h"
+#include "matlab_util.h"
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+    ensureOrError(nrhs == 2, "Must supply 2 inputs");
+    ensureOrError(isSize(prhs[0], { -1, 3 }), "Faces must be an N x 3 array");
+    ensureOrError(isSize(prhs[1], { -1, 3 }), "Vertices must be an M x 3 array");
+    Volume<int> faceVol = getVolumeChecked<int>(prhs[0], "Faces");
+    Volume<float> vertVol = getVolumeChecked<float>(prhs[1], "Vertices");
+
+    size_t nverts = vertVol.nx;
+    size_t nfaces = faceVol.nx;
+
+    std::vector<Vec3f> vertices;
+    vertices.reserve(nverts);
+    for (int i = 0; i < nverts; ++i) {
+        vertices.push_back(Vec3f(vertVol.at(i, 0, 0), vertVol.at(i, 1, 0), vertVol.at(i, 2, 0)));
+    }
+    std::vector<std::vector<int>> faces;
+    faces.reserve(nfaces);
+    for (int i = 0; i < nfaces; ++i) {
+        faces.push_back({ faceVol.at(i, 0, 0), faceVol.at(i, 1, 0), faceVol.at(i, 2, 0) });
+    }
+
+    ManifoldMesh mesh;
+    mesh.build(vertices, faces);
+
+    mxArray *mxGaussCurvs = mxCreateNumericMatrix(1, vertVol.nx, mxSINGLE_CLASS, mxREAL);
+    std::vector<float> gaussCurvs = mesh.gaussCurvatures();
+    std::memcpy(mxGetData(mxGaussCurvs), gaussCurvs.data(), nverts * sizeof(float));
+
+    plhs[0] = mxGaussCurvs;
+}
\ No newline at end of file

surfseg/matlab/mex_hausdorff.cpp

+#include <cmath>
+
+#include "mex.h"
+#include "matrix.h"
+
+#include "compute_error.h"
+#include "geomutils.h"
+#undef max
+#undef min
+
+#include "matlab_util.h"
+#include "volume.h"
+#include "util.h"
+
+inline vertex_t vmax(const vertex_t& v1, const vertex_t& v2)
+{
+	vertex_t out;
+	out.x = std::fmax(v1.x, v2.x);
+	out.y = std::fmax(v1.y, v2.y);
+	out.z = std::fmax(v1.z, v2.z);
+	return out;
+}
+
+inline vertex_t vmin(const vertex_t& v1, const vertex_t& v2)
+{
+	vertex_t out;
+	out.x = std::fmin(v1.x, v2.x);
+	out.y = std::fmin(v1.y, v2.y);
+	out.z = std::fmin(v1.z, v2.z);
+	return out;
+}
+
+model_error initModelError()
+{
+	model_error me;
+	me.min_error = 0;
+	me.max_error = 0;
+	me.mean_error = 0;
+	me.mesh = nullptr;
+	me.n_samples = 0;
+	me.fe = nullptr;
+	me.verror = nullptr;
+	me.info = nullptr;
+
+	return me;
+}
+
+dist_surf_surf_stats initDistSurfSurfStats()
+{
+	size3d grid_size;
+	grid_size.x = 0;
+	grid_size.y = 0;
+	grid_size.z = 0;
+
+	dist_surf_surf_stats s;
+	s.st_m1_area = 0;
+	s.m1_area = 0;
+	s.m2_area = 0;
+	s.min_dist = 0;
+	s.max_dist = 0;
+	s.mean_dist = 0;
+	s.rms_dist = 0;
+	s.cell_sz = 0;
+	s.n_t_p_nec = 0;
+	s.m1_samples = 0;
+	s.grid_sz = grid_size;
+	s.n_ne_cells = 0;
+
+	return s;
+}
+
+model buildModel(const Volume<int> faces, const Volume<float> verts)
+{
+	// Init fields and alloc memory
+	int nfaces = faces.nx;
+	int nverts = verts.nx;
+	model out;
+	out.num_faces = nfaces;
+	out.num_vert = nverts;
+	out.faces = static_cast<face_t *>(mxMalloc(nfaces * sizeof(face_t)));
+	out.vertices = static_cast<vertex_t *>(mxMalloc(nverts * sizeof(vertex_t)));
+	out.area = nullptr;
+	out.builtin_normals = 0;
+	out.face_normals = nullptr;
+	out.normals = nullptr;
+	out.total_area = 0.0f;
+	out.tree = nullptr;
+
+	// Fill vertices
+	vertex_t minPt; minPt.x = infOrMax<float>; minPt.y = infOrMax<float>; minPt.z = infOrMax<float>;
+	vertex_t maxPt; maxPt.x = -infOrMax<float>; maxPt.y = -infOrMax<float>; maxPt.z = -infOrMax<float>;
+	for (int i = 0; i < nverts; ++i) {
+		vertex_t pt; pt.x = verts.at(i, 0, 0); pt.y = verts.at(i, 1, 0); pt.z = verts.at(i, 2, 0);
+		minPt = vmin(minPt, pt);
+		maxPt = vmax(maxPt, pt);
+		out.vertices[i] = pt;
+	}
+	out.bBox[0] = minPt;
+	out.bBox[1] = maxPt;
+
+	// Fill faces
+	for (int i = 0; i < nfaces; ++i) {
+		face_t f; f.f0 = faces.at(i, 0, 0) - 1; f.f1 = faces.at(i, 1, 0) - 1; f.f2 = faces.at(i, 2, 0) - 1;
+		out.faces[i] = f;
+	}
+
+	return out;
+}
+
+mxArray *makeFaceErrorStruct(const face_error *fe, int nfaces)
+{
+	static const char *faceErrorNames[7] = {
+		"face_area",
+		"min_error",
+		"max_error",
+		"mean_error",
+		"mean_sqr_error",
+		"serror",
+		"samples_freq"
+	};
+
+	mxArray *mxFaceError = mxCreateStructMatrix(1, nfaces, 7, faceErrorNames);
+	for (int i = 0; i < nfaces; ++i) {
+		mxSetField(mxFaceError, i, "face_area", mxCreateDoubleScalar(fe[i].face_area));
+		mxSetField(mxFaceError, i, "min_error", mxCreateDoubleScalar(fe[i].min_error));
+		mxSetField(mxFaceError, i, "max_error", mxCreateDoubleScalar(fe[i].max_error));
+		mxSetField(mxFaceError, i, "mean_error", mxCreateDoubleScalar(fe[i].mean_error));
+		mxSetField(mxFaceError, i, "mean_sqr_error", mxCreateDoubleScalar(fe[i].mean_sqr_error));
+		mxSetField(mxFaceError, i, "sample_freq", mxCreateDoubleScalar(fe[i].sample_freq));
+		size_t nsamples = fe[i].sample_freq * (fe[i].sample_freq + 1) / 2;
+		mxArray *mxSampleErrors = mxCreateDoubleMatrix(1, nsamples, mxREAL);
+		double *sampleErrorData = mxGetPr(mxSampleErrors);
+		for (int k = 0; k < nsamples; ++k) {
+			sampleErrorData[k] = fe[i].serror[k];
+		}
+		mxSetField(mxFaceError, i, "serror", mxSampleErrors);
+	}
+
+	return mxFaceError;
+}
+
+mxArray *makeModelErrorStruct(const model_error& me)
+{
+	static const char *modelErrorNames[4] = {
+		"min_error",
+		"max_error",
+		"mean_error",
+		"n_samples"/*,
+		"fe",
+		"verror"*/
+	};
+
+	mxArray *mxModelError = mxCreateStructMatrix(1, 1, 4, modelErrorNames);
+	mxSetField(mxModelError, 0, "min_error", mxCreateDoubleScalar(me.min_error));
+	mxSetField(mxModelError, 0, "max_error", mxCreateDoubleScalar(me.max_error));
+	mxSetField(mxModelError, 0, "mean_error", mxCreateDoubleScalar(me.mean_error));
+	mxSetField(mxModelError, 0, "n_samples", mxCreateDoubleScalar(me.n_samples));
+	/*mxSetField(mxModelError, 0, "fe", makeFaceErrorStruct(me.fe, me.mesh->num_faces));
+	mxArray *mxVertError = mxCreateDoubleMatrix(1, me.mesh->num_vert, mxREAL);
+	double *vertErrorData = mxGetPr(mxVertError);
+	for (int i = 0; i < me.mesh->num_vert; ++i) {
+		vertErrorData[i] = me.verror[i];
+	}
+	mxSetField(mxModelError, 0, "verror", mxVertError);*/
+	return mxModelError;
+}
+
+mxArray *makeDistSurfSurfStatsStruct(const dist_surf_surf_stats& stats)
+{
+	static const char *distSurfSurfStatsNames[12] = {
+		"st_m1_area",
+		"m1_area",
+		"m2_area",
+		"min_dist",
+		"max_dist",
+		"mean_dist",
+		"rms_dist",
+		"cell_sz",
+		"n_t_p_nec",
+		"m1_samples",
+		"grid_sz",
+		"n_ne_cells"
+	};
+
+	mxArray *mxDistSurfSurfStats = mxCreateStructMatrix(1, 1, 12, distSurfSurfStatsNames);
+	mxSetField(mxDistSurfSurfStats, 0, "st_m1_area", mxCreateDoubleScalar(stats.st_m1_area));
+	mxSetField(mxDistSurfSurfStats, 0, "m1_area", mxCreateDoubleScalar(stats.m1_area));
+	mxSetField(mxDistSurfSurfStats, 0, "m2_area", mxCreateDoubleScalar(stats.m2_area));
+	mxSetField(mxDistSurfSurfStats, 0, "min_dist", mxCreateDoubleScalar(stats.min_dist));
+	mxSetField(mxDistSurfSurfStats, 0, "max_dist", mxCreateDoubleScalar(stats.max_dist));
+	mxSetField(mxDistSurfSurfStats, 0, "mean_dist", mxCreateDoubleScalar(stats.mean_dist));
+	mxSetField(mxDistSurfSurfStats, 0, "rms_dist", mxCreateDoubleScalar(stats.rms_dist));
+	mxSetField(mxDistSurfSurfStats, 0, "cell_sz", mxCreateDoubleScalar(stats.cell_sz));
+	mxSetField(mxDistSurfSurfStats, 0, "n_t_p_nec", mxCreateDoubleScalar(stats.n_t_p_nec));
+	mxSetField(mxDistSurfSurfStats, 0, "m1_samples", mxCreateDoubleScalar(stats.m1_samples));
+	mxArray *mxGridSz = mxCreateDoubleMatrix(1, 3, mxREAL);
+	double *gridSzData = mxGetPr(mxGridSz);
+	gridSzData[0] = stats.grid_sz.x;
+	gridSzData[1] = stats.grid_sz.y;
+	gridSzData[2] = stats.grid_sz.z;
+	mxSetField(mxDistSurfSurfStats, 0, "grid_sz", mxGridSz);
+	mxSetField(mxDistSurfSurfStats, 0, "n_ne_cells", mxCreateDoubleScalar(stats.n_ne_cells));
+
+	return mxDistSurfSurfStats;
+}
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+	//mexPrintf("Check inputs\n"); mexEvalString("drawnow;");
+	ensureOrError(nrhs == 6, "Must supply 6 inputs");
+	ensureOrError(isSize(prhs[0], { -1, 3 }), "FacesA must be an N x 3 array");
+	ensureOrError(isSize(prhs[1], { -1, 3 }), "VerticesA must be an M x 3 array");
+	ensureOrError(isSize(prhs[2], { -1, 3 }), "FacesB must be an N x 3 array");
+	ensureOrError(isSize(prhs[3], { -1, 3 }), "VerticesB must be an M x 3 array");
+	//mexPrintf("Get inputs\n"); mexEvalString("drawnow;");
+	Volume<int> faceVolA = getVolumeChecked<int>(prhs[0], "FacesA");
+	Volume<float> vertVolA = getVolumeChecked<float>(prhs[1], "VerticesA");
+	Volume<int> faceVolB = getVolumeChecked<int>(prhs[2], "FacesB");
+	Volume<float> vertVolB = getVolumeChecked<float>(prhs[3], "VerticesB");
+	double sampleDensity = getCastScalarChecked<double>(prhs[4], "sampleDensity");
+	int minSampleFreq = getCastScalarChecked<int>(prhs[5], "minSampleFreq");
+
+	// Convert data to MESH data structures
+	//mexPrintf("Init model errors\n"); mexEvalString("drawnow;");
+	model_error me1 = initModelError();
+	model_error me2 = initModelError();
+	//mexPrintf("Build models\n"); mexEvalString("drawnow;");
+	model m1 = buildModel(faceVolA, vertVolA);
+	model m2 = buildModel(faceVolB, vertVolB);
+	me1.mesh = &m1;
+	me2.mesh = &m2;
+	//mexPrintf("Compute distances\n"); mexEvalString("drawnow;");
+	dist_surf_surf_stats stats1 = initDistSurfSurfStats();
+	dist_surf_surf_stats stats2 = initDistSurfSurfStats();
+
+	// Compute distance statistics
+	//mexPrintf("Distance stats\n"); mexEvalString("drawnow;");
+	dist_surf_surf(&me1, &m2, sampleDensity, minSampleFreq, &stats1, 0, nullptr);
+	dist_surf_surf(&me2, &m1, sampleDensity, minSampleFreq, &stats2, 0, nullptr);
+
+	int tmp1, tmp2;
+	//mexPrintf("Vertex errors 1\n"); mexEvalString("drawnow;");
+	calc_vertex_error(&me1, &tmp1, &tmp2);
+	//mexPrintf("Vertex errors 2\n"); mexEvalString("drawnow;");
+	calc_vertex_error(&me2, &tmp1, &tmp2);
+
+	// Place results in MATLAB structs
+	//mexPrintf("Put in plhs\n"); mexEvalString("drawnow;");
+	plhs[0] = makeModelErrorStruct(me1);
+	plhs[1] = makeModelErrorStruct(me2);
+	/*plhs[2] = makeDistSurfSurfStatsStruct(stats1);
+	plhs[3] = makeDistSurfSurfStatsStruct(stats2);*/
+
+	// Free memory
+	//mexPrintf("Free memory\n"); mexEvalString("drawnow;");
+	free(me1.verror);
+	free(me2.verror);
+	free_face_error(me1.fe);
+	free_face_error(me2.fe);
+	mxFree(m1.faces);
+	mxFree(m2.faces);
+	mxFree(m1.vertices);
+	mxFree(m2.vertices);
+
+	//mexPrintf("Returning\n"); mexEvalString("drawnow;");
+}