More code cleanup in intern/dualcon.

Removed a lot of unused code, added comments and some clearer
naming. Minor code shuffles and style cleanup too.

4 files changed, 1389 insertions, 1864 deletions

diff --git a/intern/dualcon/intern/Projections.cpp b/intern/dualcon/intern/Projections.cpp
index e50065c004d..7e7d5e0081c 100644
--- a/intern/dualcon/intern/Projections.cpp
+++ b/intern/dualcon/intern/Projections.cpp
@@ -71,3 +71,308 @@ const int facemap[6][4] = {
 	{0, 2, 4, 6},
 	{1, 3, 5, 7}
 };
+
+/**
+ * Method to perform cross-product
+ */
+static void crossProduct(int64_t res[3], const int64_t a[3], const int64_t b[3])
+{
+	res[0] = a[1] * b[2] - a[2] * b[1];
+	res[1] = a[2] * b[0] - a[0] * b[2];
+	res[2] = a[0] * b[1] - a[1] * b[0];
+}
+
+static void crossProduct(double res[3], const double a[3], const double b[3])
+{
+	res[0] = a[1] * b[2] - a[2] * b[1];
+	res[1] = a[2] * b[0] - a[0] * b[2];
+	res[2] = a[0] * b[1] - a[1] * b[0];
+}
+
+/**
+ * Method to perform dot product
+ */
+int64_t dotProduct(const int64_t a[3], const int64_t b[3])
+{
+	return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
+}
+
+void normalize(double a[3])
+{
+	double mag = a[0] * a[0] + a[1] * a[1] + a[2] * a[2];
+	if (mag > 0) {
+		mag = sqrt(mag);
+		a[0] /= mag;
+		a[1] /= mag;
+		a[2] /= mag;
+	}
+}
+
+/* Create projection axes for cube+triangle intersection testing.
+ *    0, 1, 2: cube face normals
+ *    
+ *          3: triangle normal
+ *          
+ *    4, 5, 6,
+ *    7, 8, 9,
+ * 10, 11, 12: cross of each triangle edge vector with each cube
+ *             face normal
+ */
+static void create_projection_axes(int64_t axes[NUM_AXES][3], const int64_t tri[3][3])
+{
+	/* Cube face normals */
+	axes[0][0] = 1;
+	axes[0][1] = 0;
+	axes[0][2] = 0;
+	axes[1][0] = 0;
+	axes[1][1] = 1;
+	axes[1][2] = 0;
+	axes[2][0] = 0;
+	axes[2][1] = 0;
+	axes[2][2] = 1;
+
+	/* Get triangle edge vectors */
+	int64_t tri_edges[3][3];
+	for (int i = 0; i < 3; i++) {
+		for (int j = 0; j < 3; j++)
+			tri_edges[i][j] = tri[(i + 1) % 3][j] - tri[i][j];
+	}
+
+	/* Triangle normal */
+	crossProduct(axes[3], tri_edges[0], tri_edges[1]);
+
+	// Face edges and triangle edges
+	int ct = 4;
+	for (int i = 0; i < 3; i++) {
+		for (int j = 0; j < 3; j++) {
+			crossProduct(axes[ct], axes[j], tri_edges[i]);
+			ct++;
+		}
+	}
+}
+
+/**
+ * Construction from a cube (axes aligned) and triangle
+ */
+CubeTriangleIsect::CubeTriangleIsect(int64_t cube[2][3], int64_t tri[3][3], int64_t error, int triind)
+{
+	int i;
+	inherit = new TriangleProjection;
+	inherit->index = triind;
+
+	int64_t axes[NUM_AXES][3];
+	create_projection_axes(axes, tri);
+
+	/* Normalize face normal and store */
+	double dedge1[] = {(double)tri[1][0] - (double)tri[0][0],
+					   (double)tri[1][1] - (double)tri[0][1],
+					   (double)tri[1][2] - (double)tri[0][2]};
+	double dedge2[] = {(double)tri[2][0] - (double)tri[1][0],
+					   (double)tri[2][1] - (double)tri[1][1],
+					   (double)tri[2][2] - (double)tri[1][2]};
+	crossProduct(inherit->norm, dedge1, dedge2);
+	normalize(inherit->norm);
+
+	int64_t cubeedge[3][3];
+	for (i = 0; i < 3; i++) {
+		for (int j = 0; j < 3; j++) {
+			cubeedge[i][j] = 0;
+		}
+		cubeedge[i][i] = cube[1][i] - cube[0][i];
+	}
+
+	/* Project the cube on to each axis */
+	for (int axis = 0; axis < NUM_AXES; axis++) {
+		CubeProjection &cube_proj = cubeProj[axis];
+
+		/* Origin */
+		cube_proj.origin = dotProduct(axes[axis], cube[0]);
+
+		/* 3 direction vectors */
+		for (i = 0; i < 3; i++)
+			cube_proj.edges[i] = dotProduct(axes[axis], cubeedge[i]);
+
+		/* Offsets of 2 ends of cube projection */
+		int64_t max = 0;
+		int64_t min = 0;
+		for (i = 1; i < 8; i++) {
+			int64_t proj = (vertmap[i][0] * cube_proj.edges[0] +
+							vertmap[i][1] * cube_proj.edges[1] +
+							vertmap[i][2] * cube_proj.edges[2]);
+			if (proj > max) {
+				max = proj;
+			}
+			if (proj < min) {
+				min = proj;
+			}
+		}
+		cube_proj.min = min;
+		cube_proj.max = max;
+
+	}
+
+	/* Project the triangle on to each axis */
+	for (int axis = 0; axis < NUM_AXES; axis++) {
+		const int64_t vts[3] = {dotProduct(axes[axis], tri[0]),
+								dotProduct(axes[axis], tri[1]),
+								dotProduct(axes[axis], tri[2])};
+
+		// Triangle
+		inherit->tri_proj[axis][0] = vts[0];
+		inherit->tri_proj[axis][1] = vts[0];
+		for (i = 1; i < 3; i++) {
+			if (vts[i] < inherit->tri_proj[axis][0])
+				inherit->tri_proj[axis][0] = vts[i];
+			
+			if (vts[i] > inherit->tri_proj[axis][1])
+				inherit->tri_proj[axis][1] = vts[i];
+		}
+	}
+}
+
+/**
+ * Construction
+ * from a parent CubeTriangleIsect object and the index of the children
+ */
+CubeTriangleIsect::CubeTriangleIsect(CubeTriangleIsect *parent)
+{
+	// Copy inheritable projections
+	this->inherit = parent->inherit;
+
+	// Shrink cube projections
+	for (int i = 0; i < NUM_AXES; i++) {
+		cubeProj[i].origin = parent->cubeProj[i].origin;
+
+		for (int j = 0; j < 3; j++)
+			cubeProj[i].edges[j] = parent->cubeProj[i].edges[j] >> 1;
+		
+		cubeProj[i].min = parent->cubeProj[i].min >> 1;
+		cubeProj[i].max = parent->cubeProj[i].max >> 1;
+	}
+}
+
+unsigned char CubeTriangleIsect::getBoxMask( )
+{
+	int i, j, k;
+	int bmask[3][2] = {{0, 0}, {0, 0}, {0, 0}};
+	unsigned char boxmask = 0;
+	int64_t child_len = cubeProj[0].edges[0] >> 1;
+
+	for (i = 0; i < 3; i++) {
+		int64_t mid = cubeProj[i].origin + child_len;
+
+		// Check bounding box
+		if (mid >= inherit->tri_proj[i][0]) {
+			bmask[i][0] = 1;
+		}
+		if (mid <= inherit->tri_proj[i][1]) {
+			bmask[i][1] = 1;
+		}
+
+	}
+
+	// Fill in masks
+	int ct = 0;
+	for (i = 0; i < 2; i++) {
+		for (j = 0; j < 2; j++) {
+			for (k = 0; k < 2; k++) {
+				boxmask |= ( (bmask[0][i] & bmask[1][j] & bmask[2][k]) << ct);
+				ct++;
+			}
+		}
+	}
+
+	// Return bounding box masks
+	return boxmask;
+}
+
+
+/**
+ * Shifting a cube to a new origin
+ */
+void CubeTriangleIsect::shift(int off[3])
+{
+	for (int i = 0; i < NUM_AXES; i++) {
+		cubeProj[i].origin += (off[0] * cubeProj[i].edges[0] +
+							   off[1] * cubeProj[i].edges[1] +
+							   off[2] * cubeProj[i].edges[2]);
+	}
+}
+
+/**
+ * Method to test intersection of the triangle and the cube
+ */
+int CubeTriangleIsect::isIntersecting() const
+{
+	for (int i = 0; i < NUM_AXES; i++) {
+		/*
+		  int64_t proj0 = cubeProj[i][0] +
+		  vertmap[inherit->cubeEnds[i][0]][0] * cubeProj[i][1] +
+		  vertmap[inherit->cubeEnds[i][0]][1] * cubeProj[i][2] +
+		  vertmap[inherit->cubeEnds[i][0]][2] * cubeProj[i][3] ;
+		  int64_t proj1 = cubeProj[i][0] +
+		  vertmap[inherit->cubeEnds[i][1]][0] * cubeProj[i][1] +
+		  vertmap[inherit->cubeEnds[i][1]][1] * cubeProj[i][2] +
+		  vertmap[inherit->cubeEnds[i][1]][2] * cubeProj[i][3] ;
+		*/
+
+		int64_t proj0 = cubeProj[i].origin + cubeProj[i].min;
+		int64_t proj1 = cubeProj[i].origin + cubeProj[i].max;
+
+		if (proj0 > inherit->tri_proj[i][1] ||
+			proj1 < inherit->tri_proj[i][0]) {
+			return 0;
+		}
+	}
+
+	return 1;
+}
+
+int CubeTriangleIsect::isIntersectingPrimary(int edgeInd) const
+{
+	for (int i = 0; i < NUM_AXES; i++) {
+
+		int64_t proj0 = cubeProj[i].origin;
+		int64_t proj1 = cubeProj[i].origin + cubeProj[i].edges[edgeInd];
+
+		if (proj0 < proj1) {
+			if (proj0 > inherit->tri_proj[i][1] ||
+				proj1 < inherit->tri_proj[i][0]) {
+				return 0;
+			}
+		}
+		else {
+			if (proj1 > inherit->tri_proj[i][1] ||
+				proj0 < inherit->tri_proj[i][0]) {
+				return 0;
+			}
+		}
+
+	}
+
+	// printf( "Intersecting: %d %d\n", edgemap[edgeInd][0], edgemap[edgeInd][1] )  ;
+	return 1;
+}
+
+float CubeTriangleIsect::getIntersectionPrimary(int edgeInd) const
+{
+	int i = 3;
+
+
+	int64_t proj0 = cubeProj[i].origin;
+	int64_t proj1 = cubeProj[i].origin + cubeProj[i].edges[edgeInd];
+	int64_t proj2 = inherit->tri_proj[i][1];
+	int64_t d = proj1 - proj0;
+	double alpha;
+
+	if (d == 0)
+		alpha = 0.5;
+	else {
+		alpha = (double)((proj2 - proj0)) / (double)d;
+
+		if (alpha < 0 || alpha > 1)
+			alpha = 0.5;
+	}
+
+	return (float)alpha;
+}
diff --git a/intern/dualcon/intern/Projections.h b/intern/dualcon/intern/Projections.h
index 7740b0d1634..2cc3320f8e4 100644
--- a/intern/dualcon/intern/Projections.h
+++ b/intern/dualcon/intern/Projections.h
@@ -32,808 +32,99 @@
 #if defined(_WIN32) && !defined(__MINGW32__)
 #define isnan(n) _isnan(n)
 #define LONG __int64
+#define int64_t __int64
 #else
 #include <stdint.h>
-#define LONG int64_t
 #endif
-#define UCHAR unsigned char
 
 /**
- * Structures and classes for computing projections of triangles
- * onto separating axes during scan conversion
- *
- * @author Tao Ju
- */
-
+* Structures and classes for computing projections of triangles onto
+* separating axes during scan conversion
+*
+* @author Tao Ju
+*/
 
 extern const int vertmap[8][3];
 extern const int centmap[3][3][3][2];
 extern const int edgemap[12][2];
 extern const int facemap[6][4];
 
+/* Axes:
+ *  0,  1,  2: cube face normals
+ *    
+ *          3: triangle normal
+ *          
+ *  4,  5,  6,
+ *  7,  8,  9,
+ * 10, 11, 12: cross of each triangle edge vector with each cube
+ *             face normal
+ */
+#define NUM_AXES 13
+
 /**
  * Structure for the projections inheritable from parent
  */
-struct InheritableProjections {
-	/// Projections of triangle
-	LONG trigProj[13][2];
-
-	/// Projections of triangle vertices on primary axes
-	LONG trigVertProj[13][3];
-
-	/// Projections of triangle edges
-	LONG trigEdgeProj[13][3][2];
+struct TriangleProjection {
+	/// Projections of triangle (min and max)
+	int64_t tri_proj[NUM_AXES][2];
 
 	/// Normal of the triangle
 	double norm[3];
-	double normA, normB;
-
-	/// End points along each axis
-	//int cubeEnds[13][2] ;
 
-	/// Error range on each axis
-	/// LONG errorProj[13];
-
-#ifdef CONTAINS_INDEX
 	/// Index of polygon
 	int index;
-#endif
+};
+
+/* This is a projection for the cube against a single projection
+   axis, see CubeTriangleIsect.cubeProj */
+struct CubeProjection {
+	int64_t origin;
+	int64_t edges[3];
+	int64_t min, max;
 };
 
 
 /**
  * Class for projections of cube / triangle vertices on the separating axes
  */
-class Projections
+class CubeTriangleIsect
 {
 public:
-/// Inheritable portion
-InheritableProjections *inherit;
+	/// Inheritable portion
+	TriangleProjection *inherit;
 
-/// Projections of the cube vertices
-LONG cubeProj[13][6];
+	/// Projections of the cube vertices
+	CubeProjection cubeProj[NUM_AXES];
 
 public:
+	CubeTriangleIsect() {}
 
-Projections( )
-{
-}
-
-/**
- * Construction
- * from a cube (axes aligned) and triangle
- */
-Projections(LONG cube[2][3], LONG trig[3][3], LONG error, int triind)
-{
-	int i, j;
-	inherit = new InheritableProjections;
-#ifdef CONTAINS_INDEX
-	inherit->index = triind;
-#endif
-	/// Create axes
-	LONG axes[13][3];
-
-	// Cube faces
-	axes[0][0] = 1;
-	axes[0][1] = 0;
-	axes[0][2] = 0;
-
-	axes[1][0] = 0;
-	axes[1][1] = 1;
-	axes[1][2] = 0;
-
-	axes[2][0] = 0;
-	axes[2][1] = 0;
-	axes[2][2] = 1;
-
-	// Triangle face
-	LONG trigedge[3][3];
-	for (i = 0; i < 3; i++)
-	{
-		for (j = 0; j < 3; j++)
-		{
-			trigedge[i][j] = trig[(i + 1) % 3][j] - trig[i][j];
-		}
-	}
-	crossProduct(trigedge[0], trigedge[1], axes[3]);
-
-	/// Normalize face normal and store
-	double dedge1[] = { (double) trig[1][0] - (double) trig[0][0],
-		                (double) trig[1][1] - (double) trig[0][1],
-		                (double) trig[1][2] - (double) trig[0][2] };
-	double dedge2[] = { (double) trig[2][0] - (double) trig[1][0],
-		                (double) trig[2][1] - (double) trig[1][1],
-		                (double) trig[2][2] - (double) trig[1][2] };
-	crossProduct(dedge1, dedge2, inherit->norm);
-	normalize(inherit->norm);
-//		inherit->normA = norm[ 0 ] ;
-//		inherit->normB = norm[ 2 ] > 0 ? norm[ 1 ] : 2 + norm[ 1 ] ;
-
-	// Face edges and triangle edges
-	int ct = 4;
-	for (i = 0; i < 3; i++)
-		for (j = 0; j < 3; j++)
-		{
-			crossProduct(axes[j], trigedge[i], axes[ct]);
-			ct++;
-		}
-
-	/// Generate projections
-	LONG cubeedge[3][3];
-	for (i = 0; i < 3; i++)
-	{
-		for (j = 0; j < 3; j++)
-		{
-			cubeedge[i][j] = 0;
-		}
-		cubeedge[i][i] = cube[1][i] - cube[0][i];
-	}
-
-	for (j = 0; j < 13; j++)
-	{
-		// Origin
-		cubeProj[j][0] = dotProduct(axes[j], cube[0]);
-
-		// 3 direction vectors
-		for (i = 1; i < 4; i++)
-		{
-			cubeProj[j][i] = dotProduct(axes[j], cubeedge[i - 1]);
-		}
-
-		// Offsets of 2 ends of cube projection
-		LONG max = 0;
-		LONG min = 0;
-		for (i = 1; i < 8; i++)
-		{
-			LONG proj = vertmap[i][0] * cubeProj[j][1] + vertmap[i][1] * cubeProj[j][2] + vertmap[i][2] * cubeProj[j][3];
-			if (proj > max)
-			{
-				max = proj;
-			}
-			if (proj < min)
-			{
-				min = proj;
-			}
-		}
-		cubeProj[j][4] = min;
-		cubeProj[j][5] = max;
-
-	}
-
-	for (j = 0; j < 13; j++)
-	{
-		LONG vts[3] = { dotProduct(axes[j], trig[0]),
-			            dotProduct(axes[j], trig[1]),
-			            dotProduct(axes[j], trig[2])  };
-
-		// Vertex
-		inherit->trigVertProj[j][0] = vts[0];
-		inherit->trigVertProj[j][1] = vts[1];
-		inherit->trigVertProj[j][2] = vts[2];
-
-		// Edge
-		for (i = 0; i < 3; i++)
-		{
-			if (vts[i] < vts[(i + 1) % 3])
-			{
-				inherit->trigEdgeProj[j][i][0] = vts[i];
-				inherit->trigEdgeProj[j][i][1] = vts[(i + 1) % 3];
-			}
-			else {
-				inherit->trigEdgeProj[j][i][1] = vts[i];
-				inherit->trigEdgeProj[j][i][0] = vts[(i + 1) % 3];
-			}
-		}
-
-		// Triangle
-		inherit->trigProj[j][0] = vts[0];
-		inherit->trigProj[j][1] = vts[0];
-		for (i = 1; i < 3; i++)
-		{
-			if (vts[i] < inherit->trigProj[j][0])
-			{
-				inherit->trigProj[j][0] = vts[i];
-			}
-			if (vts[i] > inherit->trigProj[j][1])
-			{
-				inherit->trigProj[j][1] = vts[i];
-			}
-		}
-	}
-
-}
-
-/**
- * Construction
- * from a parent Projections object and the index of the children
- */
-Projections (Projections *parent)
-{
-	// Copy inheritable projections
-	this->inherit = parent->inherit;
-
-	// Shrink cube projections
-	for (int i = 0; i < 13; i++)
-	{
-		cubeProj[i][0] = parent->cubeProj[i][0];
-		for (int j = 1; j < 6; j++)
-		{
-			cubeProj[i][j] = parent->cubeProj[i][j] >> 1;
-		}
-	}
-};
-
-Projections (Projections *parent, int box[3], int depth)
-{
-	int mask =  (1 << depth) - 1;
-	int nbox[3] = { box[0] & mask, box[1] & mask, box[2] & mask };
-
-	// Copy inheritable projections
-	this->inherit = parent->inherit;
-
-	// Shrink cube projections
-	for (int i = 0; i < 13; i++)
-	{
-		for (int j = 1; j < 6; j++)
-		{
-			cubeProj[i][j] = parent->cubeProj[i][j] >> depth;
-		}
-
-		cubeProj[i][0] = parent->cubeProj[i][0] + nbox[0] * cubeProj[i][1] + nbox[1] * cubeProj[i][2] + nbox[2] * cubeProj[i][3];
-	}
-};
-
-/**
- * Testing intersection based on vertex/edge masks
- */
-int getIntersectionMasks(UCHAR cedgemask, UCHAR& edgemask)
-{
-	int i, j;
-	edgemask = cedgemask;
-
-	// Pre-processing
-	/*
-	   if ( cvertmask & 1 )
-	   {
-	    edgemask |= 5 ;
-	   }
-	   if ( cvertmask & 2 )
-	   {
-	    edgemask |= 3 ;
-	   }
-	   if ( cvertmask & 4 )
-	   {
-	    edgemask |= 6 ;
-	   }
-
+	/**
+	 * Construction from a cube (axes aligned) and triangle
 	 */
-
-	// Test axes for edge intersection
-	UCHAR bit = 1;
-	for (j = 0; j < 3; j++)
-	{
-		if (edgemask & bit)
-		{
-			for (i = 0; i < 13; i++)
-			{
-				LONG proj0 = cubeProj[i][0] + cubeProj[i][4];
-				LONG proj1 = cubeProj[i][0] + cubeProj[i][5];
-
-				if (proj0 > inherit->trigEdgeProj[i][j][1] ||
-				    proj1 < inherit->trigEdgeProj[i][j][0])
-				{
-					edgemask &= (~bit);
-					break;
-				}
-			}
-		}
-		bit <<= 1;
-	}
-
-	/*
-	   if ( edgemask != 0 )
-	   {
-	    printf("%d %d\n", cedgemask, edgemask) ;
-	   }
+	CubeTriangleIsect(int64_t cube[2][3], int64_t trig[3][3], int64_t error, int triind);
+	
+	/**
+	 * Construction from a parent CubeTriangleIsect object and the index of
+	 * the children
 	 */
+	CubeTriangleIsect(CubeTriangleIsect *parent);
+	
+	unsigned char getBoxMask( );
 
-	// Test axes for triangle intersection
-	if (edgemask)
-	{
-		return 1;
-	}
-
-	for (i = 3; i < 13; i++)
-	{
-		LONG proj0 = cubeProj[i][0] + cubeProj[i][4];
-		LONG proj1 = cubeProj[i][0] + cubeProj[i][5];
-
-		if (proj0 > inherit->trigProj[i][1] ||
-		    proj1 < inherit->trigProj[i][0])
-		{
-			return 0;
-		}
-	}
-
-	return 1;
-}
-
-/**
- * Retrieving children masks using PRIMARY AXES
- */
-UCHAR getChildrenMasks(UCHAR cvertmask, UCHAR vertmask[8])
-{
-	int i, j, k;
-	int bmask[3][2] = {{0, 0}, {0, 0}, {0, 0}};
-	int vmask[3][3][2] = {{{0, 0}, {0, 0}, {0, 0}}, {{0, 0}, {0, 0}, {0, 0}}, {{0, 0}, {0, 0}, {0, 0}}};
-	UCHAR boxmask = 0;
-	LONG len = cubeProj[0][1] >> 1;
-
-	for (i = 0; i < 3; i++)
-	{
-		LONG mid = cubeProj[i][0] + len;
-
-		// Check bounding box
-		if (mid >= inherit->trigProj[i][0])
-		{
-			bmask[i][0] = 1;
-		}
-		if (mid <= inherit->trigProj[i][1])
-		{
-			bmask[i][1] = 1;
-		}
-
-		// Check vertex mask
-		if (cvertmask)
-		{
-			for (j = 0; j < 3; j++)
-			{
-				if (cvertmask & (1 << j) )
-				{
-					// Only check if it's contained this node
-					if (mid >= inherit->trigVertProj[i][j])
-					{
-						vmask[i][j][0] = 1;
-					}
-					if (mid <= inherit->trigVertProj[i][j])
-					{
-						vmask[i][j][1] = 1;
-					}
-				}
-			}
-		}
-
-		/*
-		   // Check edge mask
-		   if ( cedgemask )
-		   {
-		    for ( j = 0 ; j < 3 ; j ++ )
-		    {
-		        if ( cedgemask & ( 1 << j ) )
-		        {
-		            // Only check if it's contained this node
-		            if ( mid >= inherit->trigEdgeProj[i][j][0] )
-		            {
-		                emask[i][j][0] = 1 ;
-		            }
-		            if ( mid <= inherit->trigEdgeProj[i][j][1] )
-		            {
-		                emask[i][j][1] = 1 ;
-		            }
-		        }
-		    }
-		   }
-		 */
-
-	}
-
-	// Fill in masks
-	int ct = 0;
-	for (i = 0; i < 2; i++)
-		for (j = 0; j < 2; j++)
-			for (k = 0; k < 2; k++)
-			{
-				boxmask |= ( (bmask[0][i] & bmask[1][j] & bmask[2][k]) << ct);
-				vertmask[ct] = ((vmask[0][0][i] & vmask[1][0][j] & vmask[2][0][k]) |
-				                ((vmask[0][1][i] & vmask[1][1][j] & vmask[2][1][k]) << 1) |
-				                ((vmask[0][2][i] & vmask[1][2][j] & vmask[2][2][k]) << 2) );
-				/*
-				   edgemask[ct] = (( emask[0][0][i] & emask[1][0][j] & emask[2][0][k] ) |
-				               (( emask[0][1][i] & emask[1][1][j] & emask[2][1][k] ) << 1 ) |
-				               (( emask[0][2][i] & emask[1][2][j] & emask[2][2][k] ) << 2 ) ) ;
-				   edgemask[ct] = cedgemask ;
-				 */
-				ct++;
-			}
-
-	// Return bounding box masks
-	return boxmask;
-}
-
-UCHAR getBoxMask( )
-{
-	int i, j, k;
-	int bmask[3][2] = {{0, 0}, {0, 0}, {0, 0}};
-	UCHAR boxmask = 0;
-	LONG len = cubeProj[0][1] >> 1;
-
-	for (i = 0; i < 3; i++)
-	{
-		LONG mid = cubeProj[i][0] + len;
-
-		// Check bounding box
-		if (mid >= inherit->trigProj[i][0])
-		{
-			bmask[i][0] = 1;
-		}
-		if (mid <= inherit->trigProj[i][1])
-		{
-			bmask[i][1] = 1;
-		}
-
-	}
-
-	// Fill in masks
-	int ct = 0;
-	for (i = 0; i < 2; i++)
-		for (j = 0; j < 2; j++)
-			for (k = 0; k < 2; k++)
-			{
-				boxmask |= ( (bmask[0][i] & bmask[1][j] & bmask[2][k]) << ct);
-				ct++;
-			}
-
-	// Return bounding box masks
-	return boxmask;
-}
-
-
-/**
- * Get projections for sub-cubes (simple axes)
- */
-void getSubProjectionsSimple(Projections *p[8])
-{
-	// Process the axes cooresponding to the triangle's normal
-	int ind = 3;
-	LONG len = cubeProj[0][1] >> 1;
-	LONG trigproj[3] = { cubeProj[ind][1] >> 1, cubeProj[ind][2] >> 1, cubeProj[ind][3] >> 1 };
-
-	int ct = 0;
-	for (int i = 0; i < 2; i++)
-		for (int j = 0; j < 2; j++)
-			for (int k = 0; k < 2; k++)
-			{
-				p[ct] = new Projections( );
-				p[ct]->inherit = inherit;
-
-				p[ct]->cubeProj[0][0] = cubeProj[0][0] + i * len;
-				p[ct]->cubeProj[1][0] = cubeProj[1][0] + j * len;
-				p[ct]->cubeProj[2][0] = cubeProj[2][0] + k * len;
-				p[ct]->cubeProj[0][1] = len;
-
-				for (int m = 1; m < 4; m++)
-				{
-					p[ct]->cubeProj[ind][m] = trigproj[m - 1];
-				}
-				p[ct]->cubeProj[ind][0] = cubeProj[ind][0] + i * trigproj[0] + j * trigproj[1] + k * trigproj[2];
-
-				ct++;
-			}
-}
-
-/**
- * Shifting a cube to a new origin
- */
-void shift(int off[3])
-{
-	for (int i = 0; i < 13; i++)
-	{
-		cubeProj[i][0] += off[0] * cubeProj[i][1] + off[1] * cubeProj[i][2] + off[2] * cubeProj[i][3];
-	}
-}
-
-void shiftNoPrimary(int off[3])
-{
-	for (int i = 3; i < 13; i++)
-	{
-		cubeProj[i][0] += off[0] * cubeProj[i][1] + off[1] * cubeProj[i][2] + off[2] * cubeProj[i][3];
-	}
-}
-
-/**
- * Method to test intersection of the triangle and the cube
- */
-int isIntersecting( )
-{
-	for (int i = 0; i < 13; i++)
-	{
-		/*
-		   LONG proj0 = cubeProj[i][0] +
-		    vertmap[inherit->cubeEnds[i][0]][0] * cubeProj[i][1] +
-		    vertmap[inherit->cubeEnds[i][0]][1] * cubeProj[i][2] +
-		    vertmap[inherit->cubeEnds[i][0]][2] * cubeProj[i][3] ;
-		   LONG proj1 = cubeProj[i][0] +
-		    vertmap[inherit->cubeEnds[i][1]][0] * cubeProj[i][1] +
-		    vertmap[inherit->cubeEnds[i][1]][1] * cubeProj[i][2] +
-		    vertmap[inherit->cubeEnds[i][1]][2] * cubeProj[i][3] ;
-		 */
-
-		LONG proj0 = cubeProj[i][0] + cubeProj[i][4];
-		LONG proj1 = cubeProj[i][0] + cubeProj[i][5];
-
-		if (proj0 > inherit->trigProj[i][1] ||
-		    proj1 < inherit->trigProj[i][0])
-		{
-			return 0;
-		}
-	}
-
-	return 1;
-};
-
-int isIntersectingNoPrimary( )
-{
-	for (int i = 3; i < 13; i++)
-	{
-		/*
-		   LONG proj0 = cubeProj[i][0] +
-		    vertmap[inherit->cubeEnds[i][0]][0] * cubeProj[i][1] +
-		    vertmap[inherit->cubeEnds[i][0]][1] * cubeProj[i][2] +
-		    vertmap[inherit->cubeEnds[i][0]][2] * cubeProj[i][3] ;
-		   LONG proj1 = cubeProj[i][0] +
-		    vertmap[inherit->cubeEnds[i][1]][0] * cubeProj[i][1] +
-		    vertmap[inherit->cubeEnds[i][1]][1] * cubeProj[i][2] +
-		    vertmap[inherit->cubeEnds[i][1]][2] * cubeProj[i][3] ;
-		 */
-
-		LONG proj0 = cubeProj[i][0] + cubeProj[i][4];
-		LONG proj1 = cubeProj[i][0] + cubeProj[i][5];
-
-		if (proj0 > inherit->trigProj[i][1] ||
-		    proj1 < inherit->trigProj[i][0])
-		{
-			return 0;
-		}
-	}
-
-	return 1;
-};
-
-/**
- * Method to test intersection of the triangle and one edge
- */
-int isIntersecting(int edgeInd)
-{
-	for (int i = 0; i < 13; i++)
-	{
-
-		LONG proj0 = cubeProj[i][0] +
-		             vertmap[edgemap[edgeInd][0]][0] * cubeProj[i][1] +
-		             vertmap[edgemap[edgeInd][0]][1] * cubeProj[i][2] +
-		             vertmap[edgemap[edgeInd][0]][2] * cubeProj[i][3];
-		LONG proj1 = cubeProj[i][0] +
-		             vertmap[edgemap[edgeInd][1]][0] * cubeProj[i][1] +
-		             vertmap[edgemap[edgeInd][1]][1] * cubeProj[i][2] +
-		             vertmap[edgemap[edgeInd][1]][2] * cubeProj[i][3];
-
-
-		if (proj0 < proj1)
-		{
-			if (proj0 > inherit->trigProj[i][1] ||
-			    proj1 < inherit->trigProj[i][0])
-			{
-				return 0;
-			}
-		}
-		else {
-			if (proj1 > inherit->trigProj[i][1] ||
-			    proj0 < inherit->trigProj[i][0])
-			{
-				return 0;
-			}
-		}
-	}
-
-	// printf( "Intersecting: %d %d\n", edgemap[edgeInd][0], edgemap[edgeInd][1] )  ;
-	return 1;
-};
-
-/**
- * Method to test intersection of one triangle edge and one cube face
- */
-int isIntersecting(int edgeInd, int faceInd)
-{
-	for (int i = 0; i < 13; i++)
-	{
-		LONG trigproj0 = inherit->trigVertProj[i][edgeInd];
-		LONG trigproj1 = inherit->trigVertProj[i][(edgeInd + 1) % 3];
-
-		if (trigproj0 < trigproj1)
-		{
-			int t1 = 1, t2 = 1;
-			for (int j = 0; j < 4; j++)
-			{
-				LONG proj = cubeProj[i][0] +
-				            vertmap[facemap[faceInd][j]][0] * cubeProj[i][1] +
-				            vertmap[facemap[faceInd][j]][1] * cubeProj[i][2] +
-				            vertmap[facemap[faceInd][j]][2] * cubeProj[i][3];
-				if (proj >= trigproj0)
-				{
-					t1 = 0;
-				}
-				if (proj <= trigproj1)
-				{
-					t2 = 0;
-				}
-			}
-			if (t1 || t2)
-			{
-				return 0;
-			}
-		}
-		else {
-			int t1 = 1, t2 = 1;
-			for (int j = 0; j < 4; j++)
-			{
-				LONG proj = cubeProj[i][0] +
-				            vertmap[facemap[faceInd][j]][0] * cubeProj[i][1] +
-				            vertmap[facemap[faceInd][j]][1] * cubeProj[i][2] +
-				            vertmap[facemap[faceInd][j]][2] * cubeProj[i][3];
-				if (proj >= trigproj1)
-				{
-					t1 = 0;
-				}
-				if (proj <= trigproj0)
-				{
-					t2 = 0;
-				}
-			}
-			if (t1 || t2)
-			{
-				return 0;
-			}
-		}
-	}
-
-	return 1;
-};
-
-
-int isIntersectingPrimary(int edgeInd)
-{
-	for (int i = 0; i < 13; i++)
-	{
-
-		LONG proj0 = cubeProj[i][0];
-		LONG proj1 = cubeProj[i][0] + cubeProj[i][edgeInd + 1];
-
-		if (proj0 < proj1)
-		{
-			if (proj0 > inherit->trigProj[i][1] ||
-			    proj1 < inherit->trigProj[i][0])
-			{
-				return 0;
-			}
-		}
-		else {
-			if (proj1 > inherit->trigProj[i][1] ||
-			    proj0 < inherit->trigProj[i][0])
-			{
-				return 0;
-			}
-		}
-
-	}
-
-	// printf( "Intersecting: %d %d\n", edgemap[edgeInd][0], edgemap[edgeInd][1] )  ;
-	return 1;
-};
-
-double getIntersection(int edgeInd)
-{
-	int i = 3;
-
-	LONG proj0 = cubeProj[i][0] +
-	             vertmap[edgemap[edgeInd][0]][0] * cubeProj[i][1] +
-	             vertmap[edgemap[edgeInd][0]][1] * cubeProj[i][2] +
-	             vertmap[edgemap[edgeInd][0]][2] * cubeProj[i][3];
-	LONG proj1 = cubeProj[i][0] +
-	             vertmap[edgemap[edgeInd][1]][0] * cubeProj[i][1] +
-	             vertmap[edgemap[edgeInd][1]][1] * cubeProj[i][2] +
-	             vertmap[edgemap[edgeInd][1]][2] * cubeProj[i][3];
-	LONG proj2 = inherit->trigProj[i][1];
-
-	/*
-	   if ( proj0 < proj1 )
-	   {
-	    if ( proj2 < proj0 || proj2 > proj1 )
-	    {
-	        return -1 ;
-	    }
-	   }
-	   else
-	   {
-	    if ( proj2 < proj1 || proj2 > proj0 )
-	    {
-	        return -1 ;
-	    }
-	   }
+	/**
+	 * Shifting a cube to a new origin
 	 */
+	void shift(int off[3]);
 
-	double alpha = (double)(proj2 - proj0) / (double)(proj1 - proj0);
-	/*
-	   if ( alpha < 0 )
-	   {
-	    alpha = 0.5 ;
-	   }
-	   else if ( alpha > 1 )
-	   {
-	    alpha = 0.5 ;
-	   }
+	/**
+	 * Method to test intersection of the triangle and the cube
 	 */
+	int isIntersecting() const;
 
-	return alpha;
-};
-
-float getIntersectionPrimary(int edgeInd)
-{
-	int i = 3;
-
-
-	LONG proj0 = cubeProj[i][0];
-	LONG proj1 = cubeProj[i][0] + cubeProj[i][edgeInd + 1];
-	LONG proj2 = inherit->trigProj[i][1];
-	LONG d = proj1 - proj0;
-	double alpha;
-
-	if (d == 0)
-		alpha = 0.5;
-	else {
-		alpha = (double)((proj2 - proj0)) / (double)d;
-
-		if (alpha < 0 || alpha > 1)
-			alpha = 0.5;
-	}
-
-	return (float)alpha;
-};
-
-/**
- * Method to perform cross-product
- */
-void crossProduct(LONG a[3], LONG b[3], LONG res[3])
-{
-	res[0] = a[1] * b[2] - a[2] * b[1];
-	res[1] = a[2] * b[0] - a[0] * b[2];
-	res[2] = a[0] * b[1] - a[1] * b[0];
-}
-void crossProduct(double a[3], double b[3], double res[3])
-{
-	res[0] = a[1] * b[2] - a[2] * b[1];
-	res[1] = a[2] * b[0] - a[0] * b[2];
-	res[2] = a[0] * b[1] - a[1] * b[0];
-}
-
-/**
- * Method to perform dot product
- */
-LONG dotProduct(LONG a[3], LONG b[3])
-{
-	return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
-}
-
-void normalize(double a[3])
-{
-	double mag = a[0] * a[0] + a[1] * a[1] + a[2] * a[2];
-	if (mag > 0)
-	{
-		mag = sqrt(mag);
-		a[0] /= mag;
-		a[1] /= mag;
-		a[2] /= mag;
-	}
-}
+	int isIntersectingPrimary(int edgeInd) const;
 
+	float getIntersectionPrimary(int edgeInd) const;
 };
 
 #endif
diff --git a/intern/dualcon/intern/octree.cpp b/intern/dualcon/intern/octree.cpp
index c5f50d8af4e..1ad502b018d 100644
--- a/intern/dualcon/intern/octree.cpp
+++ b/intern/dualcon/intern/octree.cpp
@@ -113,7 +113,7 @@ void Octree::scanConvert()
 	start = clock();
 #endif
 
-	addTrian();
+	addAllTriangles();
 	resetMinimalEdges();
 	preparePrimalEdgesMask(&root->internal);
 
@@ -257,7 +257,7 @@ void Octree::resetMinimalEdges()
 	cellProcParity(root, 0, maxDepth);
 }
 
-void Octree::addTrian()
+void Octree::addAllTriangles()
 {
 	Triangle *trian;
 	int count = 0;
@@ -273,7 +273,7 @@ void Octree::addTrian()
 	while ((trian = reader->getNextTriangle()) != NULL) {
 		// Drop triangles
 		{
-			addTrian(trian, count);
+			addTriangle(trian, count);
 		}
 		delete trian;
 
@@ -316,48 +316,60 @@ void Octree::addTrian()
 	putchar(13);
 }
 
-void Octree::addTrian(Triangle *trian, int triind)
+/* Prepare a triangle for insertion into the octree; call the other
+   addTriangle() to (recursively) build the octree */
+void Octree::addTriangle(Triangle *trian, int triind)
 {
 	int i, j;
 
-	// Blowing up the triangle to the grid
-	float mid[3] = {0, 0, 0};
-	for (i = 0; i < 3; i++)
-		for (j = 0; j < 3; j++) {
+	/* Project the triangle's coordinates into the grid */
+	for (i = 0; i < 3; i++) {
+		for (j = 0; j < 3; j++)
 			trian->vt[i][j] = dimen * (trian->vt[i][j] - origin[j]) / range;
-			mid[j] += trian->vt[i][j] / 3;
-		}
-
-	// Generate projections
-	LONG cube[2][3] = {{0, 0, 0}, {dimen, dimen, dimen}};
-	LONG trig[3][3];
+	}
 
-	for (i = 0; i < 3; i++)
-		for (j = 0; j < 3; j++) {
-			trig[i][j] = (LONG)(trian->vt[i][j]);
-			// Perturb end points, if set so
-		}
+	/* Generate projections */
+	int64_t cube[2][3] = {{0, 0, 0}, {dimen, dimen, dimen}};
+	int64_t trig[3][3];
+	for (i = 0; i < 3; i++) {
+		for (j = 0; j < 3; j++)
+			trig[i][j] = (int64_t)(trian->vt[i][j]);
+	}
 
-	// Add to the octree
-	// int start[3] = {0, 0, 0};
-	LONG errorvec = (LONG)(0);
-	Projections *proj = new Projections(cube, trig, errorvec, triind);
-	root = (Node *)addTrian(&root->internal, proj, maxDepth);
+	/* Add triangle to the octree */
+	int64_t errorvec = (int64_t)(0);
+	CubeTriangleIsect *proj = new CubeTriangleIsect(cube, trig, errorvec, triind);
+	root = (Node *)addTriangle(&root->internal, proj, maxDepth);
 
 	delete proj->inherit;
 	delete proj;
 }
 
+void print_depth(int height, int maxDepth)
+{
+	for (int i = 0; i < maxDepth - height; i++)
+		printf("  ");
+}
 
-InternalNode *Octree::addTrian(InternalNode *node, Projections *p, int height)
+InternalNode *Octree::addTriangle(InternalNode *node, CubeTriangleIsect *p, int height)
 {
 	int i;
-	int vertdiff[8][3] = {{0, 0, 0}, {0, 0, 1}, {0, 1, -1}, {0, 0, 1}, {1, -1, -1}, {0, 0, 1}, {0, 1, -1}, {0, 0, 1}};
-	UCHAR boxmask = p->getBoxMask();
-	Projections *subp = new Projections(p);
-
+	const int vertdiff[8][3] = {
+		{0,  0,  0},
+		{0,  0,  1},
+		{0,  1, -1},
+		{0,  0,  1},
+		{1, -1, -1},
+		{0,  0,  1},
+		{0,  1, -1},
+		{0,  0,  1}};
+	unsigned char boxmask = p->getBoxMask();
+	CubeTriangleIsect *subp = new CubeTriangleIsect(p);
+	
 	int count = 0;
 	int tempdiff[3] = {0, 0, 0};
+
+	/* Check triangle against each of the input node's children */
 	for (i = 0; i < 8; i++) {
 		tempdiff[0] += vertdiff[i][0];
 		tempdiff[1] += vertdiff[i][1];
@@ -370,30 +382,23 @@ InternalNode *Octree::addTrian(InternalNode *node, Projections *p, int height)
 
 			/* Pruning using intersection test */
 			if (subp->isIntersecting()) {
-			// if(subp->getIntersectionMasks(cedgemask, edgemask))
 				if (!hasChild(node, i)) {
-					if (height == 1) {
+					if (height == 1)
 						node = addLeafChild(node, i, count, createLeaf(0));
-					}
-					else {
+					else
 						node = addInternalChild(node, i, count, createInternal(0));
-					}
 				}
 				Node *chd = getChild(node, count);
 
-				if (!isLeaf(node, i)) {
-					// setChild(node, count, addTrian(chd, subp, height - 1, vertmask[i], edgemask));
-					setChild(node, count, (Node *)addTrian(&chd->internal, subp, height - 1));
-				}
-				else {
+				if (node->is_child_leaf(i))
 					setChild(node, count, (Node *)updateCell(&chd->leaf, subp));
-				}
+				else
+					setChild(node, count, (Node *)addTriangle(&chd->internal, subp, height - 1));
 			}
 		}
 
-		if (hasChild(node, i)) {
+		if (hasChild(node, i))
 			count++;
-		}
 	}
 
 	delete subp;
@@ -401,7 +406,7 @@ InternalNode *Octree::addTrian(InternalNode *node, Projections *p, int height)
 	return node;
 }
 
-LeafNode *Octree::updateCell(LeafNode *node, Projections *p)
+LeafNode *Octree::updateCell(LeafNode *node, CubeTriangleIsect *p)
 {
 	int i;
 
@@ -426,13 +431,6 @@ LeafNode *Octree::updateCell(LeafNode *node, Projections *p)
 		else {
 			offs[newc] = getEdgeOffsetNormal(node, oldc, a[newc], b[newc], c[newc]);
 
-//			if(p->isIntersectingPrimary(i))
-			{
-				// dc_printf("Multiple intersections!\n");
-
-//				setPatchEdge(node, i);
-			}
-
 			oldc++;
 			newc++;
 		}
@@ -451,7 +449,7 @@ void Octree::preparePrimalEdgesMask(InternalNode *node)
 	int count = 0;
 	for (int i = 0; i < 8; i++) {
 		if (hasChild(node, i)) {
-			if (isLeaf(node, i))
+			if (node->is_child_leaf(i))
 				createPrimalEdgesMask(&getChild(node, count)->leaf);
 			else
 				preparePrimalEdgesMask(&getChild(node, count)->internal);
@@ -495,7 +493,7 @@ Node *Octree::trace(Node *newnode, int *st, int len, int depth, PathList *& path
 			nst[i][j] = st[j] + len * vertmap[i][j];
 		}
 
-		if (chd[i] == NULL || isLeaf(&newnode->internal, i)) {
+		if (chd[i] == NULL || newnode->internal.is_child_leaf(i)) {
 			chdpaths[i] = NULL;
 		}
 		else {
@@ -1411,7 +1409,7 @@ Node *Octree::locateCell(InternalNode *node, int st[3], int len, int ori[3], int
 	if (hasChild(node, ind)) {
 		int count = getChildCount(node, ind);
 		Node *chd = getChild(node, count);
-		if (isLeaf(node, ind)) {
+		if (node->is_child_leaf(ind)) {
 			rleaf = chd;
 			rlen = len;
 		}
@@ -2367,7 +2365,7 @@ void Octree::edgeProcContour(Node *node[4], int leaf[4], int depth[4], int maxde
 					de[j] = depth[j];
 				}
 				else {
-					le[j] = isLeaf(&node[j]->internal, c[j]);
+					le[j] = node[j]->internal.is_child_leaf(c[j]);
 					ne[j] = chd[j][c[j]];
 					de[j] = depth[j] - 1;
 				}
@@ -2410,7 +2408,7 @@ void Octree::faceProcContour(Node *node[2], int leaf[2], int depth[2], int maxde
 					df[j] = depth[j];
 				}
 				else {
-					lf[j] = isLeaf(&node[j]->internal, c[j]);
+					lf[j] = node[j]->internal.is_child_leaf(c[j]);
 					nf[j] = chd[j][c[j]];
 					df[j] = depth[j] - 1;
 				}
@@ -2436,7 +2434,7 @@ void Octree::faceProcContour(Node *node[2], int leaf[2], int depth[2], int maxde
 					de[j] = depth[order[j]];
 				}
 				else {
-					le[j] = isLeaf(&node[order[j]]->internal, c[j]);
+					le[j] = node[order[j]]->internal.is_child_leaf(c[j]);
 					ne[j] = chd[order[j]][c[j]];
 					de[j] = depth[order[j]] - 1;
 				}
@@ -2467,7 +2465,7 @@ void Octree::cellProcContour(Node *node, int leaf, int depth)
 
 		// 8 Cell calls
 		for (i = 0; i < 8; i++) {
-			cellProcContour(chd[i], isLeaf(&node->internal, i), depth - 1);
+			cellProcContour(chd[i], node->internal.is_child_leaf(i), depth - 1);
 		}
 
 		// 12 face calls
@@ -2477,8 +2475,8 @@ void Octree::cellProcContour(Node *node, int leaf, int depth)
 		for (i = 0; i < 12; i++) {
 			int c[2] = {cellProcFaceMask[i][0], cellProcFaceMask[i][1]};
 
-			lf[0] = isLeaf(&node->internal, c[0]);
-			lf[1] = isLeaf(&node->internal, c[1]);
+			lf[0] = node->internal.is_child_leaf(c[0]);
+			lf[1] = node->internal.is_child_leaf(c[1]);
 
 			nf[0] = chd[c[0]];
 			nf[1] = chd[c[1]];
@@ -2494,7 +2492,7 @@ void Octree::cellProcContour(Node *node, int leaf, int depth)
 			int c[4] = {cellProcEdgeMask[i][0], cellProcEdgeMask[i][1], cellProcEdgeMask[i][2], cellProcEdgeMask[i][3]};
 
 			for (int j = 0; j < 4; j++) {
-				le[j] = isLeaf(&node->internal, c[j]);
+				le[j] = node->internal.is_child_leaf(c[j]);
 				ne[j] = chd[c[j]];
 			}
 
@@ -2563,7 +2561,7 @@ void Octree::edgeProcParity(Node *node[4], int leaf[4], int depth[4], int maxdep
 					de[j] = depth[j];
 				}
 				else {
-					le[j] = isLeaf(&node[j]->internal, c[j]);
+					le[j] = node[j]->internal.is_child_leaf(c[j]);
 					ne[j] = chd[j][c[j]];
 					de[j] = depth[j] - 1;
 
@@ -2608,7 +2606,7 @@ void Octree::faceProcParity(Node *node[2], int leaf[2], int depth[2], int maxdep
 					df[j] = depth[j];
 				}
 				else {
-					lf[j] = isLeaf(&node[j]->internal, c[j]);
+					lf[j] = node[j]->internal.is_child_leaf(c[j]);
 					nf[j] = chd[j][c[j]];
 					df[j] = depth[j] - 1;
 				}
@@ -2634,7 +2632,7 @@ void Octree::faceProcParity(Node *node[2], int leaf[2], int depth[2], int maxdep
 					de[j] = depth[order[j]];
 				}
 				else {
-					le[j] = isLeaf((InternalNode *)(node[order[j]]), c[j]);
+					le[j] = node[order[j]]->internal.is_child_leaf(c[j]);
 					ne[j] = chd[order[j]][c[j]];
 					de[j] = depth[order[j]] - 1;
 				}
@@ -2665,7 +2663,7 @@ void Octree::cellProcParity(Node *node, int leaf, int depth)
 
 		// 8 Cell calls
 		for (i = 0; i < 8; i++) {
-			cellProcParity(chd[i], isLeaf((InternalNode *)node, i), depth - 1);
+			cellProcParity(chd[i], node->internal.is_child_leaf(i), depth - 1);
 		}
 
 		// 12 face calls
@@ -2675,8 +2673,8 @@ void Octree::cellProcParity(Node *node, int leaf, int depth)
 		for (i = 0; i < 12; i++) {
 			int c[2] = {cellProcFaceMask[i][0], cellProcFaceMask[i][1]};
 
-			lf[0] = isLeaf((InternalNode *)node, c[0]);
-			lf[1] = isLeaf((InternalNode *)node, c[1]);
+			lf[0] = node->internal.is_child_leaf(c[0]);
+			lf[1] = node->internal.is_child_leaf(c[1]);
 
 			nf[0] = chd[c[0]];
 			nf[1] = chd[c[1]];
@@ -2692,7 +2690,7 @@ void Octree::cellProcParity(Node *node, int leaf, int depth)
 			int c[4] = {cellProcEdgeMask[i][0], cellProcEdgeMask[i][1], cellProcEdgeMask[i][2], cellProcEdgeMask[i][3]};
 
 			for (int j = 0; j < 4; j++) {
-				le[j] = isLeaf((InternalNode *)node, c[j]);
+				le[j] = node->internal.is_child_leaf(c[j]);
 				ne[j] = chd[c[j]];
 			}
 
diff --git a/intern/dualcon/intern/octree.h b/intern/dualcon/intern/octree.h
index 35d24a074bb..550d584baa7 100644
--- a/intern/dualcon/intern/octree.h
+++ b/intern/dualcon/intern/octree.h
@@ -65,6 +65,12 @@ struct InternalNode {
 
 	/* Can have up to eight children */
 	Node *children[0];
+
+	/// Test if child is leaf
+	int is_child_leaf(int index) const
+	{
+		return (child_is_leaf >> index) & 1;
+	}
 };
 
 
@@ -145,1277 +151,1202 @@ struct PathList {
  */
 class Octree
 {
-public:
-/* Public members */
-
-/// Memory allocators
-VirtualMemoryAllocator *alloc[9];
-VirtualMemoryAllocator *leafalloc[4];
+ public:
+	/* Public members */
 
-/// Root node
-Node *root;
+	/// Memory allocators
+	VirtualMemoryAllocator *alloc[9];
+	VirtualMemoryAllocator *leafalloc[4];
 
-/// Model reader
-ModelReader *reader;
+	/// Root node
+	Node *root;
 
-/// Marching cubes table
-Cubes *cubes;
+	/// Model reader
+	ModelReader *reader;
 
-/// Length of grid
-int dimen;
-int mindimen, minshift;
+	/// Marching cubes table
+	Cubes *cubes;
 
-/// Maximum depth
-int maxDepth;
+	/// Length of grid
+	int dimen;
+	int mindimen, minshift;
 
-/// The lower corner of the bounding box and the size
-float origin[3];
-float range;
+	/// Maximum depth
+	int maxDepth;
 
-/// Counting information
-int nodeCount;
-int nodeSpace;
-int nodeCounts[9];
+	/// The lower corner of the bounding box and the size
+	float origin[3];
+	float range;
 
-int actualQuads, actualVerts;
+	/// Counting information
+	int nodeCount;
+	int nodeSpace;
+	int nodeCounts[9];
 
-PathList *ringList;
+	int actualQuads, actualVerts;
 
-int maxTrianglePerCell;
-int outType;     // 0 for OFF, 1 for PLY, 2 for VOL
+	PathList *ringList;
 
-// For flood filling
-int use_flood_fill;
-float thresh;
+	int maxTrianglePerCell;
+	int outType;     // 0 for OFF, 1 for PLY, 2 for VOL
 
-int use_manifold;
+	// For flood filling
+	int use_flood_fill;
+	float thresh;
 
-float hermite_num;
+	int use_manifold;
 
-DualConMode mode;
-
-public:
-/**
- * Construtor
- */
-Octree(ModelReader *mr,
-       DualConAllocOutput alloc_output_func,
-       DualConAddVert add_vert_func,
-       DualConAddQuad add_quad_func,
-       DualConFlags flags, DualConMode mode, int depth,
-       float threshold, float hermite_num);
-
-/**
- * Destructor
- */
-~Octree();
-
-/**
- * Scan convert
- */
-void scanConvert();
+	float hermite_num;
 
-void *getOutputMesh() {
-	return output_mesh;
-}
+	DualConMode mode;
 
-private:
-/* Helper functions */
-
-/**
- * Initialize memory allocators
- */
-void initMemory();
-
-/**
- * Release memory
- */
-void freeMemory();
-
-/**
- * Print memory usage
- */
-void printMemUsage();
+ public:
+	/**
+	 * Construtor
+	 */
+	Octree(ModelReader *mr,
+		   DualConAllocOutput alloc_output_func,
+		   DualConAddVert add_vert_func,
+		   DualConAddQuad add_quad_func,
+		   DualConFlags flags, DualConMode mode, int depth,
+		   float threshold, float hermite_num);
+
+	/**
+	 * Destructor
+	 */
+	~Octree();
 
+	/**
+	 * Scan convert
+	 */
+	void scanConvert();
 
-/**
- * Methods to set / restore minimum edges
- */
-void resetMinimalEdges();
+	void *getOutputMesh() {
+		return output_mesh;
+	}
 
-void cellProcParity(Node *node, int leaf, int depth);
-void faceProcParity(Node * node[2], int leaf[2], int depth[2], int maxdep, int dir);
-void edgeProcParity(Node * node[4], int leaf[4], int depth[4], int maxdep, int dir);
+ private:
+	/* Helper functions */
 
-void processEdgeParity(LeafNode * node[4], int depths[4], int maxdep, int dir);
+	/**
+	 * Initialize memory allocators
+	 */
+	void initMemory();
 
-/**
- * Add triangles to the tree
- */
-void addTrian();
-void addTrian(Triangle *trian, int triind);
-InternalNode *addTrian(InternalNode *node, Projections *p, int height);
+	/**
+	 * Release memory
+	 */
+	void freeMemory();
 
-/**
- * Method to update minimizer in a cell: update edge intersections instead
- */
-LeafNode *updateCell(LeafNode *node, Projections *p);
+	/**
+	 * Print memory usage
+	 */
+	void printMemUsage();
 
-/* Routines to detect and patch holes */
-int numRings;
-int totRingLengths;
-int maxRingLength;
 
-/**
- * Entry routine.
- */
-void trace();
-/**
- * Trace the given node, find patches and fill them in
- */
-Node *trace(Node *node, int *st, int len, int depth, PathList *& paths);
-/**
- * Look for path on the face and add to paths
- */
-void findPaths(Node * node[2], int leaf[2], int depth[2], int *st[2], int maxdep, int dir, PathList * &paths);
-/**
- * Combine two list1 and list2 into list1 using connecting paths list3,
- * while closed paths are appended to rings
- */
-void combinePaths(PathList *& list1, PathList *list2, PathList *paths, PathList *& rings);
-/**
- * Helper function: combine current paths in list1 and list2 to a single path and append to list3
- */
-PathList *combineSinglePath(PathList *& head1, PathList *pre1, PathList *& list1, PathList *& head2, PathList *pre2, PathList *& list2);
+	/**
+	 * Methods to set / restore minimum edges
+	 */
+	void resetMinimalEdges();
 
-/**
- * Functions to patch rings in a node
- */
-Node *patch(Node * node, int st[3], int len, PathList * rings);
-Node *patchSplit(Node * node, int st[3], int len, PathList * rings, int dir, PathList * &nrings1, PathList * &nrings2);
-Node *patchSplitSingle(Node * node, int st[3], int len, PathElement * head, int dir, PathList * &nrings1, PathList * &nrings2);
-Node *connectFace(Node * node, int st[3], int len, int dir, PathElement * f1, PathElement * f2);
-Node *locateCell(InternalNode * node, int st[3], int len, int ori[3], int dir, int side, Node * &rleaf, int rst[3], int& rlen);
-void compressRing(PathElement *& ring);
-void getFacePoint(PathElement *leaf, int dir, int& x, int& y, float& p, float& q);
-LeafNode *patchAdjacent(InternalNode * node, int len, int st1[3], LeafNode * leaf1, int st2[3], LeafNode * leaf2, int walkdir, int inc, int dir, int side, float alpha);
-int findPair(PathElement *head, int pos, int dir, PathElement *& pre1, PathElement *& pre2);
-int getSide(PathElement *e, int pos, int dir);
-int isEqual(PathElement *e1, PathElement *e2);
-void preparePrimalEdgesMask(InternalNode *node);
-void testFacePoint(PathElement *e1, PathElement *e2);
+	void cellProcParity(Node *node, int leaf, int depth);
+	void faceProcParity(Node * node[2], int leaf[2], int depth[2], int maxdep, int dir);
+	void edgeProcParity(Node * node[4], int leaf[4], int depth[4], int maxdep, int dir);
 
-/**
- * Path-related functions
- */
-void deletePath(PathList *& head, PathList *pre, PathList *& curr);
-void printPath(PathList *path);
-void printPath(PathElement *path);
-void printElement(PathElement *ele);
-void printPaths(PathList *path);
-void checkElement(PathElement *ele);
-void checkPath(PathElement *path);
+	void processEdgeParity(LeafNode * node[4], int depths[4], int maxdep, int dir);
 
+	/**
+	 * Add triangles to the tree
+	 */
+	void addAllTriangles();
+	void addTriangle(Triangle *trian, int triind);
+	InternalNode *addTriangle(InternalNode *node, CubeTriangleIsect *p, int height);
 
-/**
- * Routines to build signs to create a partitioned volume
- *(after patching rings)
- */
-void buildSigns();
-void buildSigns(unsigned char table[], Node * node, int isLeaf, int sg, int rvalue[8]);
+	/**
+	 * Method to update minimizer in a cell: update edge intersections instead
+	 */
+	LeafNode *updateCell(LeafNode *node, CubeTriangleIsect *p);
 
-/************************************************************************/
-/* To remove disconnected components */
-/************************************************************************/
-void floodFill();
-void clearProcessBits(Node *node, int height);
-int floodFill(LeafNode * leaf, int st[3], int len, int height, int threshold);
-int floodFill(Node * node, int st[3], int len, int height, int threshold);
+	/* Routines to detect and patch holes */
+	int numRings;
+	int totRingLengths;
+	int maxRingLength;
 
-/**
- * Write out polygon file
- */
-void writeOut();
+	/**
+	 * Entry routine.
+	 */
+	void trace();
+	/**
+	 * Trace the given node, find patches and fill them in
+	 */
+	Node *trace(Node *node, int *st, int len, int depth, PathList *& paths);
+	/**
+	 * Look for path on the face and add to paths
+	 */
+	void findPaths(Node * node[2], int leaf[2], int depth[2], int *st[2], int maxdep, int dir, PathList * &paths);
+	/**
+	 * Combine two list1 and list2 into list1 using connecting paths list3,
+	 * while closed paths are appended to rings
+	 */
+	void combinePaths(PathList *& list1, PathList *list2, PathList *paths, PathList *& rings);
+	/**
+	 * Helper function: combine current paths in list1 and list2 to a single path and append to list3
+	 */
+	PathList *combineSinglePath(PathList *& head1, PathList *pre1, PathList *& list1, PathList *& head2, PathList *pre2, PathList *& list2);
 
-void countIntersection(Node *node, int height, int& nedge, int& ncell, int& nface);
-void generateMinimizer(Node * node, int st[3], int len, int height, int& offset);
-void computeMinimizer(LeafNode * leaf, int st[3], int len, float rvalue[3]);
-/**
- * Traversal functions to generate polygon model
- * op: 0 for counting, 1 for writing OBJ, 2 for writing OFF, 3 for writing PLY
- */
-void cellProcContour(Node *node, int leaf, int depth);
-void faceProcContour(Node * node[2], int leaf[2], int depth[2], int maxdep, int dir);
-void edgeProcContour(Node * node[4], int leaf[4], int depth[4], int maxdep, int dir);
-void processEdgeWrite(Node * node[4], int depths[4], int maxdep, int dir);
-
-/* output callbacks/data */
-DualConAllocOutput alloc_output;
-DualConAddVert add_vert;
-DualConAddQuad add_quad;
-void *output_mesh;
-
-private:
-/************ Operators for all nodes ************/
-
-/// Lookup table
-int numChildrenTable[256];
-int childrenCountTable[256][8];
-int childrenIndexTable[256][8];
-int numEdgeTable[8];
-int edgeCountTable[8][3];
-
-/// Build up lookup table
-void buildTable()
-{
-	for (int i = 0; i < 256; i++)
+	/**
+	 * Functions to patch rings in a node
+	 */
+	Node *patch(Node * node, int st[3], int len, PathList * rings);
+	Node *patchSplit(Node * node, int st[3], int len, PathList * rings, int dir, PathList * &nrings1, PathList * &nrings2);
+	Node *patchSplitSingle(Node * node, int st[3], int len, PathElement * head, int dir, PathList * &nrings1, PathList * &nrings2);
+	Node *connectFace(Node * node, int st[3], int len, int dir, PathElement * f1, PathElement * f2);
+	Node *locateCell(InternalNode * node, int st[3], int len, int ori[3], int dir, int side, Node * &rleaf, int rst[3], int& rlen);
+	void compressRing(PathElement *& ring);
+	void getFacePoint(PathElement *leaf, int dir, int& x, int& y, float& p, float& q);
+	LeafNode *patchAdjacent(InternalNode * node, int len, int st1[3], LeafNode * leaf1, int st2[3], LeafNode * leaf2, int walkdir, int inc, int dir, int side, float alpha);
+	int findPair(PathElement *head, int pos, int dir, PathElement *& pre1, PathElement *& pre2);
+	int getSide(PathElement *e, int pos, int dir);
+	int isEqual(PathElement *e1, PathElement *e2);
+	void preparePrimalEdgesMask(InternalNode *node);
+	void testFacePoint(PathElement *e1, PathElement *e2);
+
+	/**
+	 * Path-related functions
+	 */
+	void deletePath(PathList *& head, PathList *pre, PathList *& curr);
+	void printPath(PathList *path);
+	void printPath(PathElement *path);
+	void printElement(PathElement *ele);
+	void printPaths(PathList *path);
+	void checkElement(PathElement *ele);
+	void checkPath(PathElement *path);
+
+
+	/**
+	 * Routines to build signs to create a partitioned volume
+	 *(after patching rings)
+	 */
+	void buildSigns();
+	void buildSigns(unsigned char table[], Node * node, int isLeaf, int sg, int rvalue[8]);
+
+	/************************************************************************/
+	/* To remove disconnected components */
+	/************************************************************************/
+	void floodFill();
+	void clearProcessBits(Node *node, int height);
+	int floodFill(LeafNode * leaf, int st[3], int len, int height, int threshold);
+	int floodFill(Node * node, int st[3], int len, int height, int threshold);
+
+	/**
+	 * Write out polygon file
+	 */
+	void writeOut();
+
+	void countIntersection(Node *node, int height, int& nedge, int& ncell, int& nface);
+	void generateMinimizer(Node * node, int st[3], int len, int height, int& offset);
+	void computeMinimizer(LeafNode * leaf, int st[3], int len, float rvalue[3]);
+	/**
+	 * Traversal functions to generate polygon model
+	 * op: 0 for counting, 1 for writing OBJ, 2 for writing OFF, 3 for writing PLY
+	 */
+	void cellProcContour(Node *node, int leaf, int depth);
+	void faceProcContour(Node * node[2], int leaf[2], int depth[2], int maxdep, int dir);
+	void edgeProcContour(Node * node[4], int leaf[4], int depth[4], int maxdep, int dir);
+	void processEdgeWrite(Node * node[4], int depths[4], int maxdep, int dir);
+
+	/* output callbacks/data */
+	DualConAllocOutput alloc_output;
+	DualConAddVert add_vert;
+	DualConAddQuad add_quad;
+	void *output_mesh;
+
+ private:
+	/************ Operators for all nodes ************/
+
+	/// Lookup table
+	int numChildrenTable[256];
+	int childrenCountTable[256][8];
+	int childrenIndexTable[256][8];
+	int numEdgeTable[8];
+	int edgeCountTable[8][3];
+
+	/// Build up lookup table
+	void buildTable()
 	{
-		numChildrenTable[i] = 0;
-		int count = 0;
-		for (int j = 0; j < 8; j++)
-		{
-			numChildrenTable[i] += ((i >> j) & 1);
-			childrenCountTable[i][j] = count;
-			childrenIndexTable[i][count] = j;
-			count += ((i >> j) & 1);
+		for (int i = 0; i < 256; i++) {
+			numChildrenTable[i] = 0;
+			int count = 0;
+			for (int j = 0; j < 8; j++) {
+				numChildrenTable[i] += ((i >> j) & 1);
+				childrenCountTable[i][j] = count;
+				childrenIndexTable[i][count] = j;
+				count += ((i >> j) & 1);
+			}
 		}
-	}
 
-	for (int i = 0; i < 8; i++)
-	{
-		numEdgeTable[i] = 0;
-		int count = 0;
-		for (int j = 0; j < 3; j++)
-		{
-			numEdgeTable[i] += ((i >> j) & 1);
-			edgeCountTable[i][j] = count;
-			count += ((i >> j) & 1);
+		for (int i = 0; i < 8; i++) {
+			numEdgeTable[i] = 0;
+			int count = 0;
+			for (int j = 0; j < 3; j++) {
+				numEdgeTable[i] += ((i >> j) & 1);
+				edgeCountTable[i][j] = count;
+				count += ((i >> j) & 1);
+			}
 		}
 	}
-}
 
-int getSign(Node *node, int height, int index)
-{
-	if (height == 0)
+	int getSign(Node *node, int height, int index)
 	{
-		return getSign(&node->leaf, index);
-	}
-	else {
-		if (hasChild(&node->internal, index))
-		{
-			return getSign(getChild(&node->internal, getChildCount(&node->internal, index)),
-			               height - 1,
-			               index);
+		if (height == 0) {
+			return getSign(&node->leaf, index);
 		}
 		else {
-			return getSign(getChild(&node->internal, 0),
-			               height - 1,
-			               7 - getChildIndex(&node->internal, 0));
+			if (hasChild(&node->internal, index)) {
+				return getSign(getChild(&node->internal, getChildCount(&node->internal, index)),
+							   height - 1,
+							   index);
+			}
+			else {
+				return getSign(getChild(&node->internal, 0),
+							   height - 1,
+							   7 - getChildIndex(&node->internal, 0));
+			}
 		}
 	}
-}
 
-/************ Operators for leaf nodes ************/
+	/************ Operators for leaf nodes ************/
 
-void printInfo(int st[3])
-{
-	printf("INFO AT: %d %d %d\n", st[0] >> minshift, st[1] >> minshift, st[2] >> minshift);
-	LeafNode *leaf = (LeafNode *)locateLeafCheck(st);
-	if (leaf)
-		printInfo(leaf);
-	else
-		printf("Leaf not exists!\n");
-}
-
-void printInfo(const LeafNode *leaf)
-{
-	/*
-	   printf("Edge mask: ");
-	   for(int i = 0; i < 12; i ++)
-	   {
-	    printf("%d ", getEdgeParity(leaf, i));
-	   }
-	   printf("\n");
-	   printf("Stored edge mask: ");
-	   for(i = 0; i < 3; i ++)
-	   {
-	    printf("%d ", getStoredEdgesParity(leaf, i));
-	   }
-	   printf("\n");
-	 */
-	printf("Sign mask: ");
-	for (int i = 0; i < 8; i++)
+	void printInfo(int st[3])
 	{
-		printf("%d ", getSign(leaf, i));
+		printf("INFO AT: %d %d %d\n", st[0] >> minshift, st[1] >> minshift, st[2] >> minshift);
+		LeafNode *leaf = (LeafNode *)locateLeafCheck(st);
+		if (leaf)
+			printInfo(leaf);
+		else
+			printf("Leaf not exists!\n");
 	}
-	printf("\n");
-
-}
 
-/// Retrieve signs
-int getSign(const LeafNode *leaf, int index)
-{
-	return ((leaf->signs >> index) & 1);
-}
-
-/// Set sign
-void setSign(LeafNode *leaf, int index)
-{
-	leaf->signs |= (1 << index);
-}
+	void printInfo(const LeafNode *leaf)
+	{
+		/*
+		  printf("Edge mask: ");
+		  for(int i = 0; i < 12; i ++)
+		  {
+		  printf("%d ", getEdgeParity(leaf, i));
+		  }
+		  printf("\n");
+		  printf("Stored edge mask: ");
+		  for(i = 0; i < 3; i ++)
+		  {
+		  printf("%d ", getStoredEdgesParity(leaf, i));
+		  }
+		  printf("\n");
+		*/
+		printf("Sign mask: ");
+		for (int i = 0; i < 8; i++) {
+			printf("%d ", getSign(leaf, i));
+		}
+		printf("\n");
 
-void setSign(LeafNode *leaf, int index, int sign)
-{
-	leaf->signs &= (~(1 << index));
-	leaf->signs |= ((sign & 1) << index);
-}
+	}
 
-int getSignMask(const LeafNode *leaf)
-{
-	return leaf->signs;
-}
+	/// Retrieve signs
+	int getSign(const LeafNode *leaf, int index)
+	{
+		return ((leaf->signs >> index) & 1);
+	}
 
-void setInProcessAll(int st[3], int dir)
-{
-	int nst[3], eind;
-	for (int i = 0; i < 4; i++)
+	/// Set sign
+	void setSign(LeafNode *leaf, int index)
 	{
-		nst[0] = st[0] + dirCell[dir][i][0] * mindimen;
-		nst[1] = st[1] + dirCell[dir][i][1] * mindimen;
-		nst[2] = st[2] + dirCell[dir][i][2] * mindimen;
-		eind = dirEdge[dir][i];
+		leaf->signs |= (1 << index);
+	}
 
-		LeafNode *cell = locateLeafCheck(nst);
-		assert(cell);
+	void setSign(LeafNode *leaf, int index, int sign)
+	{
+		leaf->signs &= (~(1 << index));
+		leaf->signs |= ((sign & 1) << index);
+	}
 
-		setInProcess(cell, eind);
+	int getSignMask(const LeafNode *leaf)
+	{
+		return leaf->signs;
 	}
-}
 
-void flipParityAll(int st[3], int dir)
-{
-	int nst[3], eind;
-	for (int i = 0; i < 4; i++)
+	void setInProcessAll(int st[3], int dir)
 	{
-		nst[0] = st[0] + dirCell[dir][i][0] * mindimen;
-		nst[1] = st[1] + dirCell[dir][i][1] * mindimen;
-		nst[2] = st[2] + dirCell[dir][i][2] * mindimen;
-		eind = dirEdge[dir][i];
+		int nst[3], eind;
+		for (int i = 0; i < 4; i++) {
+			nst[0] = st[0] + dirCell[dir][i][0] * mindimen;
+			nst[1] = st[1] + dirCell[dir][i][1] * mindimen;
+			nst[2] = st[2] + dirCell[dir][i][2] * mindimen;
+			eind = dirEdge[dir][i];
+
+			LeafNode *cell = locateLeafCheck(nst);
+			assert(cell);
 
-		LeafNode *cell = locateLeaf(nst);
-		flipEdge(cell, eind);
+			setInProcess(cell, eind);
+		}
 	}
-}
 
-void setInProcess(LeafNode *leaf, int eind)
-{
-	assert(eind >= 0 && eind <= 11);
+	void flipParityAll(int st[3], int dir)
+	{
+		int nst[3], eind;
+		for (int i = 0; i < 4; i++) {
+			nst[0] = st[0] + dirCell[dir][i][0] * mindimen;
+			nst[1] = st[1] + dirCell[dir][i][1] * mindimen;
+			nst[2] = st[2] + dirCell[dir][i][2] * mindimen;
+			eind = dirEdge[dir][i];
+
+			LeafNode *cell = locateLeaf(nst);
+			flipEdge(cell, eind);
+		}
+	}
 
-	leaf->flood_fill |= (1 << eind);
-}
+	void setInProcess(LeafNode *leaf, int eind)
+	{
+		assert(eind >= 0 && eind <= 11);
 
-void setOutProcess(LeafNode *leaf, int eind)
-{
-	assert(eind >= 0 && eind <= 11);
+		leaf->flood_fill |= (1 << eind);
+	}
 
-	leaf->flood_fill &= ~(1 << eind);
-}
+	void setOutProcess(LeafNode *leaf, int eind)
+	{
+		assert(eind >= 0 && eind <= 11);
 
-int isInProcess(LeafNode *leaf, int eind)
-{
-	assert(eind >= 0 && eind <= 11);
+		leaf->flood_fill &= ~(1 << eind);
+	}
 
-	return (leaf->flood_fill >> eind) & 1;
-}
+	int isInProcess(LeafNode *leaf, int eind)
+	{
+		assert(eind >= 0 && eind <= 11);
 
-/// Generate signs at the corners from the edge parity
-void generateSigns(LeafNode *leaf, unsigned char table[], int start)
-{
-	leaf->signs = table[leaf->edge_parity];
+		return (leaf->flood_fill >> eind) & 1;
+	}
 
-	if ((start ^ leaf->signs) & 1)
+	/// Generate signs at the corners from the edge parity
+	void generateSigns(LeafNode *leaf, unsigned char table[], int start)
 	{
-		leaf->signs = ~(leaf->signs);
+		leaf->signs = table[leaf->edge_parity];
+
+		if ((start ^ leaf->signs) & 1) {
+			leaf->signs = ~(leaf->signs);
+		}
 	}
-}
 
-/// Get edge parity
-int getEdgeParity(LeafNode *leaf, int index)
-{
-	assert(index >= 0 && index <= 11);
+	/// Get edge parity
+	int getEdgeParity(LeafNode *leaf, int index)
+	{
+		assert(index >= 0 && index <= 11);
 
-	return (leaf->edge_parity >> index) & 1;
-}
+		return (leaf->edge_parity >> index) & 1;
+	}
 
-/// Get edge parity on a face
-int getFaceParity(LeafNode *leaf, int index)
-{
-	int a = getEdgeParity(leaf, faceMap[index][0]) +
+	/// Get edge parity on a face
+	int getFaceParity(LeafNode *leaf, int index)
+	{
+		int a = getEdgeParity(leaf, faceMap[index][0]) +
 	        getEdgeParity(leaf, faceMap[index][1]) +
 	        getEdgeParity(leaf, faceMap[index][2]) +
 	        getEdgeParity(leaf, faceMap[index][3]);
-	return (a & 1);
-}
-int getFaceEdgeNum(LeafNode *leaf, int index)
-{
-	int a = getEdgeParity(leaf, faceMap[index][0]) +
+		return (a & 1);
+	}
+	int getFaceEdgeNum(LeafNode *leaf, int index)
+	{
+		int a = getEdgeParity(leaf, faceMap[index][0]) +
 	        getEdgeParity(leaf, faceMap[index][1]) +
 	        getEdgeParity(leaf, faceMap[index][2]) +
 	        getEdgeParity(leaf, faceMap[index][3]);
-	return a;
-}
-
-/// Set edge parity
-void flipEdge(LeafNode *leaf, int index)
-{
-	assert(index >= 0 && index <= 11);
-
-	leaf->edge_parity ^= (1 << index);
-}
-
-/// Set 1
-void setEdge(LeafNode *leaf, int index)
-{
-	assert(index >= 0 && index <= 11);
+		return a;
+	}
 
-	leaf->edge_parity |= (1 << index);
-}
+	/// Set edge parity
+	void flipEdge(LeafNode *leaf, int index)
+	{
+		assert(index >= 0 && index <= 11);
 
-/// Set 0
-void resetEdge(LeafNode *leaf, int index)
-{
-	assert(index >= 0 && index <= 11);
+		leaf->edge_parity ^= (1 << index);
+	}
 
-	leaf->edge_parity &= ~(1 << index);
-}
+	/// Set 1
+	void setEdge(LeafNode *leaf, int index)
+	{
+		assert(index >= 0 && index <= 11);
 
-/// Flipping with a new intersection offset
-void createPrimalEdgesMask(LeafNode *leaf)
-{
-	leaf->primary_edge_intersections = getPrimalEdgesMask2(leaf);
-}
+		leaf->edge_parity |= (1 << index);
+	}
 
-void setStoredEdgesParity(LeafNode *leaf, int pindex)
-{
-	assert(pindex <= 2 && pindex >= 0);
+	/// Set 0
+	void resetEdge(LeafNode *leaf, int index)
+	{
+		assert(index >= 0 && index <= 11);
 
-	leaf->primary_edge_intersections |= (1 << pindex);
-}
-int getStoredEdgesParity(LeafNode *leaf, int pindex)
-{
-	assert(pindex <= 2 && pindex >= 0);
+		leaf->edge_parity &= ~(1 << index);
+	}
 
-	return (leaf->primary_edge_intersections >> pindex) & 1;
-}
+	/// Flipping with a new intersection offset
+	void createPrimalEdgesMask(LeafNode *leaf)
+	{
+		leaf->primary_edge_intersections = getPrimalEdgesMask2(leaf);
+	}
 
-LeafNode *flipEdge(LeafNode *leaf, int index, float alpha)
-{
-	flipEdge(leaf, index);
+	void setStoredEdgesParity(LeafNode *leaf, int pindex)
+	{
+		assert(pindex <= 2 && pindex >= 0);
 
-	if ((index & 3) == 0)
+		leaf->primary_edge_intersections |= (1 << pindex);
+	}
+	int getStoredEdgesParity(LeafNode *leaf, int pindex)
 	{
-		int ind = index / 4;
-		if (getEdgeParity(leaf, index) && !getStoredEdgesParity(leaf, ind))
-		{
-			// Create a new node
-			int num = getNumEdges(leaf) + 1;
-			setStoredEdgesParity(leaf, ind);
-			int count = getEdgeCount(leaf, ind);
-			LeafNode *nleaf = createLeaf(num);
-			*nleaf = *leaf;
+		assert(pindex <= 2 && pindex >= 0);
 
-			setEdgeOffset(nleaf, alpha, count);
+		return (leaf->primary_edge_intersections >> pindex) & 1;
+	}
 
-			if (num > 1)
-			{
-				float *pts = leaf->edge_intersections;
-				float *npts = nleaf->edge_intersections;
-				for (int i = 0; i < count; i++)
-				{
-					for (int j = 0; j < EDGE_FLOATS; j++)
-					{
-						npts[i * EDGE_FLOATS + j] = pts[i * EDGE_FLOATS + j];
+	LeafNode *flipEdge(LeafNode *leaf, int index, float alpha)
+	{
+		flipEdge(leaf, index);
+
+		if ((index & 3) == 0) {
+			int ind = index / 4;
+			if (getEdgeParity(leaf, index) && !getStoredEdgesParity(leaf, ind)) {
+				// Create a new node
+				int num = getNumEdges(leaf) + 1;
+				setStoredEdgesParity(leaf, ind);
+				int count = getEdgeCount(leaf, ind);
+				LeafNode *nleaf = createLeaf(num);
+				*nleaf = *leaf;
+
+				setEdgeOffset(nleaf, alpha, count);
+
+				if (num > 1) {
+					float *pts = leaf->edge_intersections;
+					float *npts = nleaf->edge_intersections;
+					for (int i = 0; i < count; i++) {
+						for (int j = 0; j < EDGE_FLOATS; j++) {
+							npts[i * EDGE_FLOATS + j] = pts[i * EDGE_FLOATS + j];
+						}
 					}
-				}
-				for (int i = count + 1; i < num; i++)
-				{
-					for (int j = 0; j < EDGE_FLOATS; j++)
-					{
-						npts[i * EDGE_FLOATS + j] = pts[(i - 1) * EDGE_FLOATS + j];
+					for (int i = count + 1; i < num; i++) {
+						for (int j = 0; j < EDGE_FLOATS; j++) {
+							npts[i * EDGE_FLOATS + j] = pts[(i - 1) * EDGE_FLOATS + j];
+						}
 					}
 				}
-			}
 
 
-			removeLeaf(num - 1, (LeafNode *)leaf);
-			leaf = nleaf;
+				removeLeaf(num - 1, (LeafNode *)leaf);
+				leaf = nleaf;
+			}
 		}
-	}
 
-	return leaf;
-}
+		return leaf;
+	}
 
-/// Update parent link
-void updateParent(InternalNode *node, int len, int st[3], LeafNode *leaf)
-{
-	// First, locate the parent
-	int count;
-	InternalNode *parent = locateParent(node, len, st, count);
+	/// Update parent link
+	void updateParent(InternalNode *node, int len, int st[3], LeafNode *leaf)
+	{
+		// First, locate the parent
+		int count;
+		InternalNode *parent = locateParent(node, len, st, count);
 
-	// Update
-	setChild(parent, count, (Node *)leaf);
-}
+		// Update
+		setChild(parent, count, (Node *)leaf);
+	}
 
-void updateParent(InternalNode *node, int len, int st[3])
-{
-	if (len == dimen)
+	void updateParent(InternalNode *node, int len, int st[3])
 	{
-		root = (Node *)node;
-		return;
-	}
+		if (len == dimen) {
+			root = (Node *)node;
+			return;
+		}
 
-	// First, locate the parent
-	int count;
-	InternalNode *parent = locateParent(len, st, count);
+		// First, locate the parent
+		int count;
+		InternalNode *parent = locateParent(len, st, count);
 
-	// UPdate
-	setChild(parent, count, (Node *)node);
-}
+		// UPdate
+		setChild(parent, count, (Node *)node);
+	}
 
-/// Find edge intersection on a given edge
-int getEdgeIntersectionByIndex(int st[3], int index, float pt[3], int check)
-{
-	// First, locat the leaf
-	LeafNode *leaf;
-	if (check)
+	/// Find edge intersection on a given edge
+	int getEdgeIntersectionByIndex(int st[3], int index, float pt[3], int check)
 	{
-		leaf = locateLeafCheck(st);
+		// First, locat the leaf
+		LeafNode *leaf;
+		if (check) {
+			leaf = locateLeafCheck(st);
+		}
+		else {
+			leaf = locateLeaf(st);
+		}
+
+		if (leaf && getStoredEdgesParity(leaf, index)) {
+			float off = getEdgeOffset(leaf, getEdgeCount(leaf, index));
+			pt[0] = (float) st[0];
+			pt[1] = (float) st[1];
+			pt[2] = (float) st[2];
+			pt[index] += off * mindimen;
+
+			return 1;
+		}
+		else {
+			return 0;
+		}
 	}
-	else {
-		leaf = locateLeaf(st);
+
+	/// Retrieve number of edges intersected
+	int getPrimalEdgesMask(LeafNode *leaf)
+	{
+		return leaf->primary_edge_intersections;
 	}
 
-	if (leaf && getStoredEdgesParity(leaf, index))
+	int getPrimalEdgesMask2(LeafNode *leaf)
 	{
-		float off = getEdgeOffset(leaf, getEdgeCount(leaf, index));
-		pt[0] = (float) st[0];
-		pt[1] = (float) st[1];
-		pt[2] = (float) st[2];
-		pt[index] += off * mindimen;
+		return (((leaf->edge_parity &   0x1) >> 0) |
+				((leaf->edge_parity &  0x10) >> 3) |
+				((leaf->edge_parity & 0x100) >> 6));
+	}
 
-		return 1;
+	/// Get the count for a primary edge
+	int getEdgeCount(LeafNode *leaf, int index)
+	{
+		return edgeCountTable[getPrimalEdgesMask(leaf)][index];
 	}
-	else {
-		return 0;
+	int getNumEdges(LeafNode *leaf)
+	{
+		return numEdgeTable[getPrimalEdgesMask(leaf)];
 	}
-}
-
-/// Retrieve number of edges intersected
-int getPrimalEdgesMask(LeafNode *leaf)
-{
-	return leaf->primary_edge_intersections;
-}
-
-int getPrimalEdgesMask2(LeafNode *leaf)
-{
-	return (((leaf->edge_parity &   0x1) >> 0) |
-	        ((leaf->edge_parity &  0x10) >> 3) |
-	        ((leaf->edge_parity & 0x100) >> 6));
-}
-
-/// Get the count for a primary edge
-int getEdgeCount(LeafNode *leaf, int index)
-{
-	return edgeCountTable[getPrimalEdgesMask(leaf)][index];
-}
-int getNumEdges(LeafNode *leaf)
-{
-	return numEdgeTable[getPrimalEdgesMask(leaf)];
-}
-
-int getNumEdges2(LeafNode *leaf)
-{
-	return numEdgeTable[getPrimalEdgesMask2(leaf)];
-}
 
-/// Set edge intersection
-void setEdgeOffset(LeafNode *leaf, float pt, int count)
-{
-	float *pts = leaf->edge_intersections;
-	pts[EDGE_FLOATS * count] = pt;
-	pts[EDGE_FLOATS * count + 1] = 0;
-	pts[EDGE_FLOATS * count + 2] = 0;
-	pts[EDGE_FLOATS * count + 3] = 0;
-}
-
-/// Set multiple edge intersections
-void setEdgeOffsets(LeafNode *leaf, float pt[3], int len)
-{
-	float *pts = leaf->edge_intersections;
-	for (int i = 0; i < len; i++)
+	int getNumEdges2(LeafNode *leaf)
 	{
-		pts[i] = pt[i];
+		return numEdgeTable[getPrimalEdgesMask2(leaf)];
 	}
-}
 
-/// Retrieve edge intersection
-float getEdgeOffset(LeafNode *leaf, int count)
-{
-	return leaf->edge_intersections[4 * count];
-}
-
-/// Update method
-LeafNode *updateEdgeOffsets(LeafNode *leaf, int oldlen, int newlen, float offs[3])
-{
-	// First, create a new leaf node
-	LeafNode *nleaf = createLeaf(newlen);
-	*nleaf = *leaf;
+	/// Set edge intersection
+	void setEdgeOffset(LeafNode *leaf, float pt, int count)
+	{
+		float *pts = leaf->edge_intersections;
+		pts[EDGE_FLOATS * count] = pt;
+		pts[EDGE_FLOATS * count + 1] = 0;
+		pts[EDGE_FLOATS * count + 2] = 0;
+		pts[EDGE_FLOATS * count + 3] = 0;
+	}
 
-	// Next, fill in the offsets
-	setEdgeOffsets(nleaf, offs, newlen);
+	/// Set multiple edge intersections
+	void setEdgeOffsets(LeafNode *leaf, float pt[3], int len)
+	{
+		float *pts = leaf->edge_intersections;
+		for (int i = 0; i < len; i++) {
+			pts[i] = pt[i];
+		}
+	}
 
-	// Finally, delete the old leaf
-	removeLeaf(oldlen, leaf);
+	/// Retrieve edge intersection
+	float getEdgeOffset(LeafNode *leaf, int count)
+	{
+		return leaf->edge_intersections[4 * count];
+	}
 
-	return nleaf;
-}
+	/// Update method
+	LeafNode *updateEdgeOffsets(LeafNode *leaf, int oldlen, int newlen, float offs[3])
+	{
+		// First, create a new leaf node
+		LeafNode *nleaf = createLeaf(newlen);
+		*nleaf = *leaf;
 
-/// Set minimizer index
-void setMinimizerIndex(LeafNode *leaf, int index)
-{
-	leaf->minimizer_index = index;
-}
+		// Next, fill in the offsets
+		setEdgeOffsets(nleaf, offs, newlen);
 
-/// Get minimizer index
-int getMinimizerIndex(LeafNode *leaf)
-{
-	return leaf->minimizer_index;
-}
+		// Finally, delete the old leaf
+		removeLeaf(oldlen, leaf);
 
-int getMinimizerIndex(LeafNode *leaf, int eind)
-{
-	int add = manifold_table[getSignMask(leaf)].pairs[eind][0] - 1;
-	assert(add >= 0);
-	return leaf->minimizer_index + add;
-}
+		return nleaf;
+	}
 
-void getMinimizerIndices(LeafNode *leaf, int eind, int inds[2])
-{
-	const int *add = manifold_table[getSignMask(leaf)].pairs[eind];
-	inds[0] = leaf->minimizer_index + add[0] - 1;
-	if (add[0] == add[1])
+	/// Set minimizer index
+	void setMinimizerIndex(LeafNode *leaf, int index)
 	{
-		inds[1] = -1;
+		leaf->minimizer_index = index;
 	}
-	else {
-		inds[1] = leaf->minimizer_index + add[1] - 1;
+
+	/// Get minimizer index
+	int getMinimizerIndex(LeafNode *leaf)
+	{
+		return leaf->minimizer_index;
 	}
-}
 
+	int getMinimizerIndex(LeafNode *leaf, int eind)
+	{
+		int add = manifold_table[getSignMask(leaf)].pairs[eind][0] - 1;
+		assert(add >= 0);
+		return leaf->minimizer_index + add;
+	}
 
-/// Set edge intersection
-void setEdgeOffsetNormal(LeafNode *leaf, float pt, float a, float b, float c, int count)
-{
-	float *pts = leaf->edge_intersections;
-	pts[4 * count] = pt;
-	pts[4 * count + 1] = a;
-	pts[4 * count + 2] = b;
-	pts[4 * count + 3] = c;
-}
-
-float getEdgeOffsetNormal(LeafNode *leaf, int count, float& a, float& b, float& c)
-{
-	float *pts = leaf->edge_intersections;
-	a = pts[4 * count + 1];
-	b = pts[4 * count + 2];
-	c = pts[4 * count + 3];
-	return pts[4 * count];
-}
-
-/// Set multiple edge intersections
-void setEdgeOffsetsNormals(LeafNode *leaf, float pt[], float a[], float b[], float c[], int len)
-{
-	float *pts = leaf->edge_intersections;
-	for (int i = 0; i < len; i++)
+	void getMinimizerIndices(LeafNode *leaf, int eind, int inds[2])
 	{
-		if (pt[i] > 1 || pt[i] < 0)
-		{
-			printf("\noffset: %f\n", pt[i]);
+		const int *add = manifold_table[getSignMask(leaf)].pairs[eind];
+		inds[0] = leaf->minimizer_index + add[0] - 1;
+		if (add[0] == add[1]) {
+			inds[1] = -1;
+		}
+		else {
+			inds[1] = leaf->minimizer_index + add[1] - 1;
 		}
-		pts[i * 4] = pt[i];
-		pts[i * 4 + 1] = a[i];
-		pts[i * 4 + 2] = b[i];
-		pts[i * 4 + 3] = c[i];
 	}
-}
-
-/// Retrieve complete edge intersection
-void getEdgeIntersectionByIndex(LeafNode *leaf, int index, int st[3], int len, float pt[3], float nm[3])
-{
-	int count = getEdgeCount(leaf, index);
-	float *pts = leaf->edge_intersections;
-
-	float off = pts[4 * count];
-
-	pt[0] = (float) st[0];
-	pt[1] = (float) st[1];
-	pt[2] = (float) st[2];
-	pt[index] += (off * len);
-
-	nm[0] = pts[4 * count + 1];
-	nm[1] = pts[4 * count + 2];
-	nm[2] = pts[4 * count + 3];
-}
-
-float getEdgeOffsetNormalByIndex(LeafNode *leaf, int index, float nm[3])
-{
-	int count = getEdgeCount(leaf, index);
-	float *pts = leaf->edge_intersections;
-
-	float off = pts[4 * count];
-
-	nm[0] = pts[4 * count + 1];
-	nm[1] = pts[4 * count + 2];
-	nm[2] = pts[4 * count + 3];
 
-	return off;
-}
 
-void fillEdgeIntersections(LeafNode *leaf, int st[3], int len, float pts[12][3], float norms[12][3])
-{
-	int i;
-	// int stt[3] = {0, 0, 0};
-
-	// The three primal edges are easy
-	int pmask[3] = {0, 4, 8};
-	for (i = 0; i < 3; i++)
+	/// Set edge intersection
+	void setEdgeOffsetNormal(LeafNode *leaf, float pt, float a, float b, float c, int count)
 	{
-		if (getEdgeParity(leaf, pmask[i]))
-		{
-			// getEdgeIntersectionByIndex(leaf, i, stt, 1, pts[pmask[i]], norms[pmask[i]]);
-			getEdgeIntersectionByIndex(leaf, i, st, len, pts[pmask[i]], norms[pmask[i]]);
-		}
+		float *pts = leaf->edge_intersections;
+		pts[4 * count] = pt;
+		pts[4 * count + 1] = a;
+		pts[4 * count + 2] = b;
+		pts[4 * count + 3] = c;
 	}
 
-	// 3 face adjacent cubes
-	int fmask[3][2] = {{6, 10}, {2, 9}, {1, 5}};
-	int femask[3][2] = {{1, 2}, {0, 2}, {0, 1}};
-	for (i = 0; i < 3; i++)
+	float getEdgeOffsetNormal(LeafNode *leaf, int count, float& a, float& b, float& c)
 	{
-		int e1 = getEdgeParity(leaf, fmask[i][0]);
-		int e2 = getEdgeParity(leaf, fmask[i][1]);
-		if (e1 || e2)
-		{
-			int nst[3] = {st[0], st[1], st[2]};
-			nst[i] += len;
-			// int nstt[3] = {0, 0, 0};
-			// nstt[i] += 1;
-			LeafNode *node = locateLeaf(nst);
+		float *pts = leaf->edge_intersections;
+		a = pts[4 * count + 1];
+		b = pts[4 * count + 2];
+		c = pts[4 * count + 3];
+		return pts[4 * count];
+	}
 
-			if (e1)
-			{
-				// getEdgeIntersectionByIndex(node, femask[i][0], nstt, 1, pts[fmask[i][0]], norms[fmask[i][0]]);
-				getEdgeIntersectionByIndex(node, femask[i][0], nst, len, pts[fmask[i][0]], norms[fmask[i][0]]);
-			}
-			if (e2)
-			{
-				// getEdgeIntersectionByIndex(node, femask[i][1], nstt, 1, pts[fmask[i][1]], norms[fmask[i][1]]);
-				getEdgeIntersectionByIndex(node, femask[i][1], nst, len, pts[fmask[i][1]], norms[fmask[i][1]]);
+	/// Set multiple edge intersections
+	void setEdgeOffsetsNormals(LeafNode *leaf, const float pt[],
+							   const float a[], const float b[],
+							   const float c[], int len)
+	{
+		float *pts = leaf->edge_intersections;
+		for (int i = 0; i < len; i++) {
+			if (pt[i] > 1 || pt[i] < 0) {
+				printf("\noffset: %f\n", pt[i]);
 			}
+			pts[i * 4] = pt[i];
+			pts[i * 4 + 1] = a[i];
+			pts[i * 4 + 2] = b[i];
+			pts[i * 4 + 3] = c[i];
 		}
 	}
 
-	// 3 edge adjacent cubes
-	int emask[3] = {3, 7, 11};
-	int eemask[3] = {0, 1, 2};
-	for (i = 0; i < 3; i++)
+	/// Retrieve complete edge intersection
+	void getEdgeIntersectionByIndex(LeafNode *leaf, int index, int st[3], int len, float pt[3], float nm[3])
 	{
-		if (getEdgeParity(leaf, emask[i]))
-		{
-			int nst[3] = {st[0] + len, st[1] + len, st[2] + len};
-			nst[i] -= len;
-			// int nstt[3] = {1, 1, 1};
-			// nstt[i] -= 1;
-			LeafNode *node = locateLeaf(nst);
+		int count = getEdgeCount(leaf, index);
+		float *pts = leaf->edge_intersections;
 
-			// getEdgeIntersectionByIndex(node, eemask[i], nstt, 1, pts[emask[i]], norms[emask[i]]);
-			getEdgeIntersectionByIndex(node, eemask[i], nst, len, pts[emask[i]], norms[emask[i]]);
-		}
-	}
-}
+		float off = pts[4 * count];
 
+		pt[0] = (float) st[0];
+		pt[1] = (float) st[1];
+		pt[2] = (float) st[2];
+		pt[index] += (off * len);
 
-void fillEdgeIntersections(LeafNode *leaf, int st[3], int len, float pts[12][3], float norms[12][3], int parity[12])
-{
-	int i;
-	for (i = 0; i < 12; i++)
-	{
-		parity[i] = 0;
+		nm[0] = pts[4 * count + 1];
+		nm[1] = pts[4 * count + 2];
+		nm[2] = pts[4 * count + 3];
 	}
-	// int stt[3] = {0, 0, 0};
 
-	// The three primal edges are easy
-	int pmask[3] = {0, 4, 8};
-	for (i = 0; i < 3; i++)
+	float getEdgeOffsetNormalByIndex(LeafNode *leaf, int index, float nm[3])
 	{
-		if (getStoredEdgesParity(leaf, i))
-		{
-			// getEdgeIntersectionByIndex(leaf, i, stt, 1, pts[pmask[i]], norms[pmask[i]]);
-			getEdgeIntersectionByIndex(leaf, i, st, len, pts[pmask[i]], norms[pmask[i]]);
-			parity[pmask[i]] = 1;
-		}
+		int count = getEdgeCount(leaf, index);
+		float *pts = leaf->edge_intersections;
+
+		float off = pts[4 * count];
+
+		nm[0] = pts[4 * count + 1];
+		nm[1] = pts[4 * count + 2];
+		nm[2] = pts[4 * count + 3];
+
+		return off;
 	}
 
-	// 3 face adjacent cubes
-	int fmask[3][2] = {{6, 10}, {2, 9}, {1, 5}};
-	int femask[3][2] = {{1, 2}, {0, 2}, {0, 1}};
-	for (i = 0; i < 3; i++)
+	void fillEdgeIntersections(LeafNode *leaf, int st[3], int len, float pts[12][3], float norms[12][3])
 	{
-		{
-			int nst[3] = {st[0], st[1], st[2]};
-			nst[i] += len;
-			// int nstt[3] = {0, 0, 0};
-			// nstt[i] += 1;
-			LeafNode *node = locateLeafCheck(nst);
-			if (node == NULL)
-			{
-				continue;
+		int i;
+		// int stt[3] = {0, 0, 0};
+
+		// The three primal edges are easy
+		int pmask[3] = {0, 4, 8};
+		for (i = 0; i < 3; i++) {
+			if (getEdgeParity(leaf, pmask[i])) {
+				// getEdgeIntersectionByIndex(leaf, i, stt, 1, pts[pmask[i]], norms[pmask[i]]);
+				getEdgeIntersectionByIndex(leaf, i, st, len, pts[pmask[i]], norms[pmask[i]]);
 			}
+		}
 
-			int e1 = getStoredEdgesParity(node, femask[i][0]);
-			int e2 = getStoredEdgesParity(node, femask[i][1]);
-
-			if (e1)
-			{
-				// getEdgeIntersectionByIndex(node, femask[i][0], nstt, 1, pts[fmask[i][0]], norms[fmask[i][0]]);
-				getEdgeIntersectionByIndex(node, femask[i][0], nst, len, pts[fmask[i][0]], norms[fmask[i][0]]);
-				parity[fmask[i][0]] = 1;
-			}
-			if (e2)
-			{
-				// getEdgeIntersectionByIndex(node, femask[i][1], nstt, 1, pts[fmask[i][1]], norms[fmask[i][1]]);
-				getEdgeIntersectionByIndex(node, femask[i][1], nst, len, pts[fmask[i][1]], norms[fmask[i][1]]);
-				parity[fmask[i][1]] = 1;
+		// 3 face adjacent cubes
+		int fmask[3][2] = {{6, 10}, {2, 9}, {1, 5}};
+		int femask[3][2] = {{1, 2}, {0, 2}, {0, 1}};
+		for (i = 0; i < 3; i++) {
+			int e1 = getEdgeParity(leaf, fmask[i][0]);
+			int e2 = getEdgeParity(leaf, fmask[i][1]);
+			if (e1 || e2) {
+				int nst[3] = {st[0], st[1], st[2]};
+				nst[i] += len;
+				// int nstt[3] = {0, 0, 0};
+				// nstt[i] += 1;
+				LeafNode *node = locateLeaf(nst);
+
+				if (e1) {
+					// getEdgeIntersectionByIndex(node, femask[i][0], nstt, 1, pts[fmask[i][0]], norms[fmask[i][0]]);
+					getEdgeIntersectionByIndex(node, femask[i][0], nst, len, pts[fmask[i][0]], norms[fmask[i][0]]);
+				}
+				if (e2) {
+					// getEdgeIntersectionByIndex(node, femask[i][1], nstt, 1, pts[fmask[i][1]], norms[fmask[i][1]]);
+					getEdgeIntersectionByIndex(node, femask[i][1], nst, len, pts[fmask[i][1]], norms[fmask[i][1]]);
+				}
 			}
 		}
-	}
 
-	// 3 edge adjacent cubes
-	int emask[3] = {3, 7, 11};
-	int eemask[3] = {0, 1, 2};
-	for (i = 0; i < 3; i++)
-	{
-//			if(getEdgeParity(leaf, emask[i]))
-		{
-			int nst[3] = {st[0] + len, st[1] + len, st[2] + len};
-			nst[i] -= len;
-			// int nstt[3] = {1, 1, 1};
-			// nstt[i] -= 1;
-			LeafNode *node = locateLeafCheck(nst);
-			if (node == NULL)
-			{
-				continue;
-			}
+		// 3 edge adjacent cubes
+		int emask[3] = {3, 7, 11};
+		int eemask[3] = {0, 1, 2};
+		for (i = 0; i < 3; i++) {
+			if (getEdgeParity(leaf, emask[i])) {
+				int nst[3] = {st[0] + len, st[1] + len, st[2] + len};
+				nst[i] -= len;
+				// int nstt[3] = {1, 1, 1};
+				// nstt[i] -= 1;
+				LeafNode *node = locateLeaf(nst);
 
-			if (getStoredEdgesParity(node, eemask[i]))
-			{
 				// getEdgeIntersectionByIndex(node, eemask[i], nstt, 1, pts[emask[i]], norms[emask[i]]);
 				getEdgeIntersectionByIndex(node, eemask[i], nst, len, pts[emask[i]], norms[emask[i]]);
-				parity[emask[i]] = 1;
 			}
 		}
 	}
-}
 
-void fillEdgeOffsetsNormals(LeafNode *leaf, int st[3], int len, float pts[12], float norms[12][3], int parity[12])
-{
-	int i;
-	for (i = 0; i < 12; i++)
-	{
-		parity[i] = 0;
-	}
-	// int stt[3] = {0, 0, 0};
 
-	// The three primal edges are easy
-	int pmask[3] = {0, 4, 8};
-	for (i = 0; i < 3; i++)
+	void fillEdgeIntersections(LeafNode *leaf, int st[3], int len, float pts[12][3], float norms[12][3], int parity[12])
 	{
-		if (getStoredEdgesParity(leaf, i))
-		{
-			pts[pmask[i]] = getEdgeOffsetNormalByIndex(leaf, i, norms[pmask[i]]);
-			parity[pmask[i]] = 1;
+		int i;
+		for (i = 0; i < 12; i++) {
+			parity[i] = 0;
+		}
+		// int stt[3] = {0, 0, 0};
+
+		// The three primal edges are easy
+		int pmask[3] = {0, 4, 8};
+		for (i = 0; i < 3; i++) {
+			if (getStoredEdgesParity(leaf, i)) {
+				// getEdgeIntersectionByIndex(leaf, i, stt, 1, pts[pmask[i]], norms[pmask[i]]);
+				getEdgeIntersectionByIndex(leaf, i, st, len, pts[pmask[i]], norms[pmask[i]]);
+				parity[pmask[i]] = 1;
+			}
 		}
-	}
 
-	// 3 face adjacent cubes
-	int fmask[3][2] = {{6, 10}, {2, 9}, {1, 5}};
-	int femask[3][2] = {{1, 2}, {0, 2}, {0, 1}};
-	for (i = 0; i < 3; i++)
-	{
-		{
-			int nst[3] = {st[0], st[1], st[2]};
-			nst[i] += len;
-			// int nstt[3] = {0, 0, 0};
-			// nstt[i] += 1;
-			LeafNode *node = locateLeafCheck(nst);
-			if (node == NULL)
+		// 3 face adjacent cubes
+		int fmask[3][2] = {{6, 10}, {2, 9}, {1, 5}};
+		int femask[3][2] = {{1, 2}, {0, 2}, {0, 1}};
+		for (i = 0; i < 3; i++) {
 			{
-				continue;
-			}
+				int nst[3] = {st[0], st[1], st[2]};
+				nst[i] += len;
+				// int nstt[3] = {0, 0, 0};
+				// nstt[i] += 1;
+				LeafNode *node = locateLeafCheck(nst);
+				if (node == NULL) {
+					continue;
+				}
 
-			int e1 = getStoredEdgesParity(node, femask[i][0]);
-			int e2 = getStoredEdgesParity(node, femask[i][1]);
+				int e1 = getStoredEdgesParity(node, femask[i][0]);
+				int e2 = getStoredEdgesParity(node, femask[i][1]);
 
-			if (e1)
-			{
-				pts[fmask[i][0]] = getEdgeOffsetNormalByIndex(node, femask[i][0], norms[fmask[i][0]]);
-				parity[fmask[i][0]] = 1;
+				if (e1) {
+					// getEdgeIntersectionByIndex(node, femask[i][0], nstt, 1, pts[fmask[i][0]], norms[fmask[i][0]]);
+					getEdgeIntersectionByIndex(node, femask[i][0], nst, len, pts[fmask[i][0]], norms[fmask[i][0]]);
+					parity[fmask[i][0]] = 1;
+				}
+				if (e2) {
+					// getEdgeIntersectionByIndex(node, femask[i][1], nstt, 1, pts[fmask[i][1]], norms[fmask[i][1]]);
+					getEdgeIntersectionByIndex(node, femask[i][1], nst, len, pts[fmask[i][1]], norms[fmask[i][1]]);
+					parity[fmask[i][1]] = 1;
+				}
 			}
-			if (e2)
+		}
+
+		// 3 edge adjacent cubes
+		int emask[3] = {3, 7, 11};
+		int eemask[3] = {0, 1, 2};
+		for (i = 0; i < 3; i++) {
+			//			if(getEdgeParity(leaf, emask[i]))
 			{
-				pts[fmask[i][1]] = getEdgeOffsetNormalByIndex(node, femask[i][1], norms[fmask[i][1]]);
-				parity[fmask[i][1]] = 1;
+				int nst[3] = {st[0] + len, st[1] + len, st[2] + len};
+				nst[i] -= len;
+				// int nstt[3] = {1, 1, 1};
+				// nstt[i] -= 1;
+				LeafNode *node = locateLeafCheck(nst);
+				if (node == NULL) {
+					continue;
+				}
+
+				if (getStoredEdgesParity(node, eemask[i])) {
+					// getEdgeIntersectionByIndex(node, eemask[i], nstt, 1, pts[emask[i]], norms[emask[i]]);
+					getEdgeIntersectionByIndex(node, eemask[i], nst, len, pts[emask[i]], norms[emask[i]]);
+					parity[emask[i]] = 1;
+				}
 			}
 		}
 	}
 
-	// 3 edge adjacent cubes
-	int emask[3] = {3, 7, 11};
-	int eemask[3] = {0, 1, 2};
-	for (i = 0; i < 3; i++)
+	void fillEdgeOffsetsNormals(LeafNode *leaf, int st[3], int len, float pts[12], float norms[12][3], int parity[12])
 	{
-//			if(getEdgeParity(leaf, emask[i]))
-		{
-			int nst[3] = {st[0] + len, st[1] + len, st[2] + len};
-			nst[i] -= len;
-			// int nstt[3] = {1, 1, 1};
-			// nstt[i] -= 1;
-			LeafNode *node = locateLeafCheck(nst);
-			if (node == NULL)
-			{
-				continue;
+		int i;
+		for (i = 0; i < 12; i++) {
+			parity[i] = 0;
+		}
+		// int stt[3] = {0, 0, 0};
+
+		// The three primal edges are easy
+		int pmask[3] = {0, 4, 8};
+		for (i = 0; i < 3; i++) {
+			if (getStoredEdgesParity(leaf, i)) {
+				pts[pmask[i]] = getEdgeOffsetNormalByIndex(leaf, i, norms[pmask[i]]);
+				parity[pmask[i]] = 1;
 			}
+		}
 
-			if (getStoredEdgesParity(node, eemask[i]))
+		// 3 face adjacent cubes
+		int fmask[3][2] = {{6, 10}, {2, 9}, {1, 5}};
+		int femask[3][2] = {{1, 2}, {0, 2}, {0, 1}};
+		for (i = 0; i < 3; i++) {
 			{
-				pts[emask[i]] = getEdgeOffsetNormalByIndex(node, eemask[i], norms[emask[i]]);
-				parity[emask[i]] = 1;
+				int nst[3] = {st[0], st[1], st[2]};
+				nst[i] += len;
+				// int nstt[3] = {0, 0, 0};
+				// nstt[i] += 1;
+				LeafNode *node = locateLeafCheck(nst);
+				if (node == NULL) {
+					continue;
+				}
+
+				int e1 = getStoredEdgesParity(node, femask[i][0]);
+				int e2 = getStoredEdgesParity(node, femask[i][1]);
+
+				if (e1) {
+					pts[fmask[i][0]] = getEdgeOffsetNormalByIndex(node, femask[i][0], norms[fmask[i][0]]);
+					parity[fmask[i][0]] = 1;
+				}
+				if (e2) {
+					pts[fmask[i][1]] = getEdgeOffsetNormalByIndex(node, femask[i][1], norms[fmask[i][1]]);
+					parity[fmask[i][1]] = 1;
+				}
 			}
 		}
-	}
-}
 
+		// 3 edge adjacent cubes
+		int emask[3] = {3, 7, 11};
+		int eemask[3] = {0, 1, 2};
+		for (i = 0; i < 3; i++) {
+			//			if(getEdgeParity(leaf, emask[i]))
+			{
+				int nst[3] = {st[0] + len, st[1] + len, st[2] + len};
+				nst[i] -= len;
+				// int nstt[3] = {1, 1, 1};
+				// nstt[i] -= 1;
+				LeafNode *node = locateLeafCheck(nst);
+				if (node == NULL) {
+					continue;
+				}
 
-/// Update method
-LeafNode *updateEdgeOffsetsNormals(LeafNode *leaf, int oldlen, int newlen, float offs[3], float a[3], float b[3], float c[3])
-{
-	// First, create a new leaf node
-	LeafNode *nleaf = createLeaf(newlen);
-	*nleaf = *leaf;
+				if (getStoredEdgesParity(node, eemask[i])) {
+					pts[emask[i]] = getEdgeOffsetNormalByIndex(node, eemask[i], norms[emask[i]]);
+					parity[emask[i]] = 1;
+				}
+			}
+		}
+	}
 
-	// Next, fill in the offsets
-	setEdgeOffsetsNormals(nleaf, offs, a, b, c, newlen);
 
-	// Finally, delete the old leaf
-	removeLeaf(oldlen, leaf);
+	/// Update method
+	LeafNode *updateEdgeOffsetsNormals(LeafNode *leaf, int oldlen, int newlen, float offs[3], float a[3], float b[3], float c[3])
+	{
+		// First, create a new leaf node
+		LeafNode *nleaf = createLeaf(newlen);
+		*nleaf = *leaf;
 
-	return nleaf;
-}
+		// Next, fill in the offsets
+		setEdgeOffsetsNormals(nleaf, offs, a, b, c, newlen);
 
-/// Locate a leaf
-/// WARNING: assuming this leaf already exists!
+		// Finally, delete the old leaf
+		removeLeaf(oldlen, leaf);
 
-LeafNode *locateLeaf(int st[3])
-{
-	Node *node = (Node *)root;
-	for (int i = GRID_DIMENSION - 1; i > GRID_DIMENSION - maxDepth - 1; i--)
-	{
-		int index = (((st[0] >> i) & 1) << 2) |
-		            (((st[1] >> i) & 1) << 1) |
-		            (((st[2] >> i) & 1));
-		node = getChild(&node->internal, getChildCount(&node->internal, index));
+		return nleaf;
 	}
 
-	return &node->leaf;
-}
+	/// Locate a leaf
+	/// WARNING: assuming this leaf already exists!
 
-LeafNode *locateLeaf(InternalNode *parent, int len, int st[3])
-{
-	Node *node = (Node *)parent;
-	int index;
-	for (int i = len / 2; i >= mindimen; i >>= 1)
+	LeafNode *locateLeaf(int st[3])
 	{
-		index = (((st[0] & i) ? 4 : 0) |
-		         ((st[1] & i) ? 2 : 0) |
-		         ((st[2] & i) ? 1 : 0));
-		node = getChild(&node->internal,
-		                getChildCount(&node->internal, index));
-	}
+		Node *node = (Node *)root;
+		for (int i = GRID_DIMENSION - 1; i > GRID_DIMENSION - maxDepth - 1; i--) {
+			int index = (((st[0] >> i) & 1) << 2) |
+				(((st[1] >> i) & 1) << 1) |
+				(((st[2] >> i) & 1));
+			node = getChild(&node->internal, getChildCount(&node->internal, index));
+		}
 
-	return &node->leaf;
-}
+		return &node->leaf;
+	}
 
-LeafNode *locateLeafCheck(int st[3])
-{
-	Node *node = (Node *)root;
-	for (int i = GRID_DIMENSION - 1; i > GRID_DIMENSION - maxDepth - 1; i--)
-	{
-		int index = (((st[0] >> i) & 1) << 2) |
-		            (((st[1] >> i) & 1) << 1) |
-		            (((st[2] >> i) & 1));
-		if (!hasChild(&node->internal, index))
-		{
-			return NULL;
+	LeafNode *locateLeaf(InternalNode *parent, int len, int st[3])
+	{
+		Node *node = (Node *)parent;
+		int index;
+		for (int i = len / 2; i >= mindimen; i >>= 1) {
+			index = (((st[0] & i) ? 4 : 0) |
+					 ((st[1] & i) ? 2 : 0) |
+					 ((st[2] & i) ? 1 : 0));
+			node = getChild(&node->internal,
+							getChildCount(&node->internal, index));
 		}
-		node = getChild(&node->internal, getChildCount(&node->internal, index));
-	}
 
-	return &node->leaf;
-}
+		return &node->leaf;
+	}
 
-InternalNode *locateParent(int len, int st[3], int& count)
-{
-	InternalNode *node = (InternalNode *)root;
-	InternalNode *pre = NULL;
-	int index = 0;
-	for (int i = dimen / 2; i >= len; i >>= 1)
+	LeafNode *locateLeafCheck(int st[3])
 	{
-		index = (((st[0] & i) ? 4 : 0) |
-		         ((st[1] & i) ? 2 : 0) |
-		         ((st[2] & i) ? 1 : 0));
-		pre = node;
-		node = &getChild(node, getChildCount(node, index))->internal;
-	}
+		Node *node = (Node *)root;
+		for (int i = GRID_DIMENSION - 1; i > GRID_DIMENSION - maxDepth - 1; i--) {
+			int index = (((st[0] >> i) & 1) << 2) |
+				(((st[1] >> i) & 1) << 1) |
+				(((st[2] >> i) & 1));
+			if (!hasChild(&node->internal, index)) {
+				return NULL;
+			}
+			node = getChild(&node->internal, getChildCount(&node->internal, index));
+		}
 
-	count = getChildCount(pre, index);
-	return pre;
-}
+		return &node->leaf;
+	}
 
-InternalNode *locateParent(InternalNode *parent, int len, int st[3], int& count)
-{
-	InternalNode *node = parent;
-	InternalNode *pre = NULL;
-	int index = 0;
-	for (int i = len / 2; i >= mindimen; i >>= 1)
+	InternalNode *locateParent(int len, int st[3], int& count)
 	{
-		index = (((st[0] & i) ? 4 : 0) |
-		         ((st[1] & i) ? 2 : 0) |
-		         ((st[2] & i) ? 1 : 0));
-		pre = node;
-		node = (InternalNode *)getChild(node, getChildCount(node, index));
-	}
+		InternalNode *node = (InternalNode *)root;
+		InternalNode *pre = NULL;
+		int index = 0;
+		for (int i = dimen / 2; i >= len; i >>= 1) {
+			index = (((st[0] & i) ? 4 : 0) |
+					 ((st[1] & i) ? 2 : 0) |
+					 ((st[2] & i) ? 1 : 0));
+			pre = node;
+			node = &getChild(node, getChildCount(node, index))->internal;
+		}
 
-	count = getChildCount(pre, index);
-	return pre;
-}
+		count = getChildCount(pre, index);
+		return pre;
+	}
 
-/************ Operators for internal nodes ************/
+	InternalNode *locateParent(InternalNode *parent, int len, int st[3], int& count)
+	{
+		InternalNode *node = parent;
+		InternalNode *pre = NULL;
+		int index = 0;
+		for (int i = len / 2; i >= mindimen; i >>= 1) {
+			index = (((st[0] & i) ? 4 : 0) |
+					 ((st[1] & i) ? 2 : 0) |
+					 ((st[2] & i) ? 1 : 0));
+			pre = node;
+			node = (InternalNode *)getChild(node, getChildCount(node, index));
+		}
 
-/// If child index exists
-int hasChild(InternalNode *node, int index)
-{
-	return (node->has_child >> index) & 1;
-}
+		count = getChildCount(pre, index);
+		return pre;
+	}
 
-/// Test if child is leaf
-int isLeaf(InternalNode *node, int index)
-{
-	return (node->child_is_leaf >> index) & 1;
-}
+	/************ Operators for internal nodes ************/
 
-/// Get the pointer to child index
-Node *getChild(InternalNode *node, int count)
-{
-	return node->children[count];
-};
+	/// If child index exists
+	int hasChild(InternalNode *node, int index)
+	{
+		return (node->has_child >> index) & 1;
+	}
 
-/// Get total number of children
-int getNumChildren(InternalNode *node)
-{
-	return numChildrenTable[node->has_child];
-}
+	/// Get the pointer to child index
+	Node *getChild(InternalNode *node, int count)
+	{
+		return node->children[count];
+	};
 
-/// Get the count of children
-int getChildCount(InternalNode *node, int index)
-{
-	return childrenCountTable[node->has_child][index];
-}
-int getChildIndex(InternalNode *node, int count)
-{
-	return childrenIndexTable[node->has_child][count];
-}
-int *getChildCounts(InternalNode *node)
-{
-	return childrenCountTable[node->has_child];
-}
+	/// Get total number of children
+	int getNumChildren(InternalNode *node)
+	{
+		return numChildrenTable[node->has_child];
+	}
 
-/// Get all children
-void fillChildren(InternalNode *node, Node *children[8], int leaf[8])
-{
-	int count = 0;
-	for (int i = 0; i < 8; i++)
-	{
-		leaf[i] = isLeaf(node, i);
-		if (hasChild(node, i))
-		{
-			children[i] = getChild(node, count);
-			count++;
-		}
-		else {
-			children[i] = NULL;
-			leaf[i] = 0;
-		}
+	/// Get the count of children
+	int getChildCount(InternalNode *node, int index)
+	{
+		return childrenCountTable[node->has_child][index];
+	}
+	int getChildIndex(InternalNode *node, int count)
+	{
+		return childrenIndexTable[node->has_child][count];
+	}
+	int *getChildCounts(InternalNode *node)
+	{
+		return childrenCountTable[node->has_child];
 	}
-}
 
-/// Sets the child pointer
-void setChild(InternalNode *node, int count, Node *chd)
-{
-	node->children[count] = chd;
-}
-void setInternalChild(InternalNode *node, int index, int count, InternalNode *chd)
-{
-	setChild(node, count, (Node *)chd);
-	node->has_child |= (1 << index);
-}
-void setLeafChild(InternalNode *node, int index, int count, LeafNode *chd)
-{
-	setChild(node, count, (Node *)chd);
-	node->has_child |= (1 << index);
-	node->child_is_leaf |= (1 << index);
-}
-
-/// Add a kid to an existing internal node
-/// Fix me: can we do this without wasting memory ?
-/// Fixed: using variable memory
-InternalNode *addChild(InternalNode *node, int index, Node *child, int aLeaf)
-{
-	// Create new internal node
-	int num = getNumChildren(node);
-	InternalNode *rnode = createInternal(num + 1);
-
-	// Establish children
-	int i;
-	int count1 = 0, count2 = 0;
-	for (i = 0; i < 8; i++)
-	{
-		if (i == index)
-		{
-			if (aLeaf)
-			{
-				setLeafChild(rnode, i, count2, &child->leaf);
-			}
-			else {
-				setInternalChild(rnode, i, count2, &child->internal);
-			}
-			count2++;
-		}
-		else if (hasChild(node, i))
-		{
-			if (isLeaf(node, i))
-			{
-				setLeafChild(rnode, i, count2, &getChild(node, count1)->leaf);
+	/// Get all children
+	void fillChildren(InternalNode *node, Node *children[8], int leaf[8])
+	{
+		int count = 0;
+		for (int i = 0; i < 8; i++) {
+			leaf[i] = node->is_child_leaf(i);
+			if (hasChild(node, i)) {
+				children[i] = getChild(node, count);
+				count++;
 			}
 			else {
-				setInternalChild(rnode, i, count2, &getChild(node, count1)->internal);
+				children[i] = NULL;
+				leaf[i] = 0;
 			}
-			count1++;
-			count2++;
 		}
 	}
 
-	removeInternal(num, node);
-	return rnode;
-}
+	/// Sets the child pointer
+	void setChild(InternalNode *node, int count, Node *chd)
+	{
+		node->children[count] = chd;
+	}
+	void setInternalChild(InternalNode *node, int index, int count, InternalNode *chd)
+	{
+		setChild(node, count, (Node *)chd);
+		node->has_child |= (1 << index);
+	}
+	void setLeafChild(InternalNode *node, int index, int count, LeafNode *chd)
+	{
+		setChild(node, count, (Node *)chd);
+		node->has_child |= (1 << index);
+		node->child_is_leaf |= (1 << index);
+	}
 
-/// Allocate a node
-InternalNode *createInternal(int length)
-{
-	InternalNode *inode = (InternalNode *)alloc[length]->allocate();
-	inode->has_child = 0;
-	inode->child_is_leaf = 0;
-	return inode;
-}
+	/// Add a kid to an existing internal node
+	/// Fix me: can we do this without wasting memory ?
+	/// Fixed: using variable memory
+	InternalNode *addChild(InternalNode *node, int index, Node *child, int aLeaf)
+	{
+		// Create new internal node
+		int num = getNumChildren(node);
+		InternalNode *rnode = createInternal(num + 1);
 
-LeafNode *createLeaf(int length)
-{
-	assert(length <= 3);
+		// Establish children
+		int i;
+		int count1 = 0, count2 = 0;
+		for (i = 0; i < 8; i++) {
+			if (i == index) {
+				if (aLeaf) {
+					setLeafChild(rnode, i, count2, &child->leaf);
+				}
+				else {
+					setInternalChild(rnode, i, count2, &child->internal);
+				}
+				count2++;
+			}
+			else if (hasChild(node, i)) {
+				if (node->is_child_leaf(i)) {
+					setLeafChild(rnode, i, count2, &getChild(node, count1)->leaf);
+				}
+				else {
+					setInternalChild(rnode, i, count2, &getChild(node, count1)->internal);
+				}
+				count1++;
+				count2++;
+			}
+		}
 
-	LeafNode *lnode = (LeafNode *)leafalloc[length]->allocate();
-	lnode->edge_parity = 0;
-	lnode->primary_edge_intersections = 0;
-	lnode->signs = 0;
+		removeInternal(num, node);
+		return rnode;
+	}
 
-	return lnode;
-}
+	/// Allocate a node
+	InternalNode *createInternal(int length)
+	{
+		InternalNode *inode = (InternalNode *)alloc[length]->allocate();
+		inode->has_child = 0;
+		inode->child_is_leaf = 0;
+		return inode;
+	}
 
-void removeInternal(int num, InternalNode *node)
-{
-	alloc[num]->deallocate(node);
-}
+	LeafNode *createLeaf(int length)
+	{
+		assert(length <= 3);
 
-void removeLeaf(int num, LeafNode *leaf)
-{
-	assert(num >= 0 && num <= 3);
-	leafalloc[num]->deallocate(leaf);
-}
+		LeafNode *lnode = (LeafNode *)leafalloc[length]->allocate();
+		lnode->edge_parity = 0;
+		lnode->primary_edge_intersections = 0;
+		lnode->signs = 0;
 
-/// Add a leaf (by creating a new par node with the leaf added)
-InternalNode *addLeafChild(InternalNode *par, int index, int count,
-                           LeafNode *leaf)
-{
-	int num = getNumChildren(par) + 1;
-	InternalNode *npar = createInternal(num);
-	*npar = *par;
+		return lnode;
+	}
 
-	if (num == 1)
+	void removeInternal(int num, InternalNode *node)
 	{
-		setLeafChild(npar, index, 0, leaf);
+		alloc[num]->deallocate(node);
 	}
-	else {
-		int i;
-		for (i = 0; i < count; i++)
-		{
-			setChild(npar, i, getChild(par, i));
-		}
-		setLeafChild(npar, index, count, leaf);
-		for (i = count + 1; i < num; i++)
-		{
-			setChild(npar, i, getChild(par, i - 1));
-		}
+
+	void removeLeaf(int num, LeafNode *leaf)
+	{
+		assert(num >= 0 && num <= 3);
+		leafalloc[num]->deallocate(leaf);
 	}
 
-	removeInternal(num - 1, par);
-	return npar;
-}
+	/// Add a leaf (by creating a new par node with the leaf added)
+	InternalNode *addLeafChild(InternalNode *par, int index, int count,
+							   LeafNode *leaf)
+	{
+		int num = getNumChildren(par) + 1;
+		InternalNode *npar = createInternal(num);
+		*npar = *par;
 
-InternalNode *addInternalChild(InternalNode *par, int index, int count,
-                               InternalNode *node)
-{
-	int num = getNumChildren(par) + 1;
-	InternalNode *npar = createInternal(num);
-	*npar = *par;
+		if (num == 1) {
+			setLeafChild(npar, index, 0, leaf);
+		}
+		else {
+			int i;
+			for (i = 0; i < count; i++) {
+				setChild(npar, i, getChild(par, i));
+			}
+			setLeafChild(npar, index, count, leaf);
+			for (i = count + 1; i < num; i++) {
+				setChild(npar, i, getChild(par, i - 1));
+			}
+		}
 
-	if (num == 1)
-	{
-		setInternalChild(npar, index, 0, node);
+		removeInternal(num - 1, par);
+		return npar;
 	}
-	else {
-		int i;
-		for (i = 0; i < count; i++)
-		{
-			setChild(npar, i, getChild(par, i));
+
+	InternalNode *addInternalChild(InternalNode *par, int index, int count,
+								   InternalNode *node)
+	{
+		int num = getNumChildren(par) + 1;
+		InternalNode *npar = createInternal(num);
+		*npar = *par;
+
+		if (num == 1) {
+			setInternalChild(npar, index, 0, node);
 		}
-		setInternalChild(npar, index, count, node);
-		for (i = count + 1; i < num; i++)
-		{
-			setChild(npar, i, getChild(par, i - 1));
+		else {
+			int i;
+			for (i = 0; i < count; i++) {
+				setChild(npar, i, getChild(par, i));
+			}
+			setInternalChild(npar, index, count, node);
+			for (i = count + 1; i < num; i++) {
+				setChild(npar, i, getChild(par, i - 1));
+			}
 		}
-	}
 
-	removeInternal(num - 1, par);
-	return npar;
-}
+		removeInternal(num - 1, par);
+		return npar;
+	}
 };
 
 #endif