From 9f6f514d05fb9ed085fc930521180ae54f036975 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Fri, 27 Jul 2018 16:23:39 -0400
Subject: [PATCH] start on decl data structure

---
 analysis/Minkowski.cpp |   2 -
 analysis/Minkowski.h   |   1 -
 analysis/decl.cpp      | 338 +++++++++++++++++++++++++++++++++++++++++
 analysis/decl.h        |  56 +++++++
 tests/lbpm_uCT_pp.cpp  |   6 +-
 5 files changed, 397 insertions(+), 6 deletions(-)
 create mode 100644 analysis/decl.cpp
 create mode 100644 analysis/decl.h

diff --git a/analysis/Minkowski.cpp b/analysis/Minkowski.cpp
index 7be0ec95..8b9a473d 100644
--- a/analysis/Minkowski.cpp
+++ b/analysis/Minkowski.cpp
@@ -221,7 +221,5 @@ void Minkowski::PrintAll()
 		fprintf(LOGFILE,"%.5g %.5g %.5g %.5g\n",vol_n_global, An_global, Jn_global, euler_global);			// minkowski measures
 		fflush(LOGFILE);
 	}
-
-
 }
 
diff --git a/analysis/Minkowski.h b/analysis/Minkowski.h
index f30d8357..d88747ab 100644
--- a/analysis/Minkowski.h
+++ b/analysis/Minkowski.h
@@ -82,6 +82,5 @@ public:
 	void SortBlobs();
 };
 
-
 #endif
 
diff --git a/analysis/decl.cpp b/analysis/decl.cpp
new file mode 100644
index 00000000..bda8dac0
--- /dev/null
+++ b/analysis/decl.cpp
@@ -0,0 +1,338 @@
+#include "analysis/pmmc.h"
+#include "analysis/decl.h"
+
+/* 
+Double connected edge list (DECL) 
+*/
+
+Vertex::Vertex(){
+	count = 0;
+}
+
+Vertex::~Vertex(){
+	
+}
+void Vertex::add(Point P){
+	vertex_data.push_back(P.x);
+	vertex_data.push_back(P.y);
+	vertex_data.push_back(P.z);
+	count++;
+}
+
+void Vertex::assign(unsigned long int idx, Point P){
+	vertex_data[3*idx] = P.x;
+	vertex_data[3*idx+1] = P.y;
+	vertex_data[3*idx+2] = P.z;
+}
+
+Point Vertex::coords(unsigned long int idx){
+	Point P;
+	P.x = vertex_data[3*idx];
+	P.y = vertex_data[3*idx+1];
+	P.z = vertex_data[3*idx+2];
+	return P;
+}
+
+unsigned long int Halfedge::v1(unsigned long int edge){
+	return HalfEdge(0,edge);
+}
+
+unsigned long int Halfedge::v2(unsigned long int edge){
+	return HalfEdge(1,edge);
+}
+
+unsigned long int Halfedge::face(unsigned long int edge){
+	return HalfEdge(2,edge);
+}
+
+unsigned long int Halfedge::twin(unsigned long int edge){
+	return HalfEdge(3,edge);
+}
+
+unsigned long int Halfedge::prev(unsigned long int edge){
+	return HalfEdge(4,edge);
+}
+
+unsigned long int Halfedge::next(unsigned long int edge){
+	return HalfEdge(5,edge);
+}
+
+Point DECL::TriNormal(int edge)
+{
+	Point P,Q;
+	double ux,uy,uz,vx,vy,vz;
+	double nx,ny,nz,len;
+	if (edge == -1)			P.x = 1.0; P.y = 0.0; P.z = 0.0; // x cube face
+	else if (edge == -2)	P.x = 0.0; P.y = 1.0; P.z = 0.0; // y cube face
+	else if (edge == -3)	P.x = 0.0; P.y = 0.0; P.z = 1.0; // z cube face
+	else{
+		// coordinates for first edge
+		P = vertex.coords(halfedge.v1(edge));
+		Q = vertex.coords(halfedge.v2(edge));
+		ux = Q.x-P.x;
+		uy = Q.y-P.y;
+		uz = Q.z-P.z;
+		// coordinates for second edge
+		P = vertex.coords(halfedge.v1(halfedge.next(edge)));
+		Q = vertex.coords(halfedge.v2(halfedge.next(edge)));
+		vx = Q.x-P.x;
+		vy = Q.y-P.y;
+		vz = Q.z-P.z;
+		// normal vector
+		nx = uy*vz - uz*vy;
+		ny = uz*vx - ux*vz;
+		nz = ux*vy - uy*vx;
+		len = sqrt(nx*nx+ny*ny+nz*nz);
+		P.x = nx/len; P.y = ny/len; P.z = nz/len;
+	}
+	return P;
+}
+
+double DECL::EdgeAngle(int edge)
+{
+	double angle;
+	Point U,V; // triangle normal vectors
+	U = TriNormal(edge);
+	V = TriNormal(halfedge.twin(edge));
+	angle = acos(U.x*V.x + U.y*V.y + U.z*V.z);
+	return angle;
+}
+
+void Isosurface(DoubleArray &A, const double &v)
+{
+	Point P,Q;
+	Point PlaceHolder;
+	double temp;
+	Point C0,C1,C2,C3,C4,C5,C6,C7;
+	
+	int TriangleCount;
+    int NewVertexCount;
+	int CubeIndex;
+	int nTris, nVert;
+	
+	Point VertexList[12];
+	Point NewVertexList[12];
+	int LocalRemap[12];
+	
+	DTMutableList<Point> cellvertices = DTMutableList<Point>(20);
+	IntArray Triangles = IntArray(3,20);
+
+	// Values from array 'A' at the cube corners
+	double CubeValues[8];
+	
+	int Nx = A.size(0);
+	int Ny = A.size(1);
+	int Nz = A.size(2);
+	
+	// Points corresponding to cube corners
+	C0.x = 0.0; C0.y = 0.0; C0.z = 0.0;
+	C1.x = 1.0; C1.y = 0.0; C1.z = 0.0;
+	C2.x = 1.0; C2.y = 1.0; C2.z = 0.0;
+	C3.x = 0.0; C3.y = 1.0; C3.z = 0.0;
+	C4.x = 0.0; C4.y = 0.0; C4.z = 1.0;
+	C5.x = 1.0; C5.y = 0.0; C5.z = 1.0;
+	C6.x = 1.0; C6.y = 1.0; C6.z = 1.0;
+	C7.x = 0.0; C7.y = 1.0; C7.z = 1.0;							
+	
+	for (int k=1; k<Nz-1; k++){		
+		for (int j=1; j<Ny-1; j++){
+			for (int i=1; i<Nx-1; i++){
+				// Set the corner values for this cube
+				CubeValues[0] = A(i,j,k);
+				CubeValues[1] = A(i+1,j,k);
+				CubeValues[2] = A(i+1,j+1,k);
+				CubeValues[3] = A(i,j+1,k);
+				CubeValues[4] = A(i,j,k+1);
+				CubeValues[5] = A(i+1,j,k+1);
+				CubeValues[6] = A(i+1,j+1,k+1);
+				CubeValues[7] = A(i,j+1,k+1);
+
+				//Determine the index into the edge table which
+				//tells us which vertices are inside of the surface
+				CubeIndex = 0;
+				if (CubeValues[0] < 0.0f) CubeIndex |= 1;
+				if (CubeValues[1] < 0.0f) CubeIndex |= 2;
+				if (CubeValues[2] < 0.0f) CubeIndex |= 4;
+				if (CubeValues[3] < 0.0f) CubeIndex |= 8;
+				if (CubeValues[4] < 0.0f) CubeIndex |= 16;
+				if (CubeValues[5] < 0.0f) CubeIndex |= 32;
+				if (CubeValues[6] < 0.0f) CubeIndex |= 64;
+				if (CubeValues[7] < 0.0f) CubeIndex |= 128;
+
+				//Find the vertices where the surface intersects the cube
+				if (edgeTable[CubeIndex] & 1){
+					P = VertexInterp(C0,C1,CubeValues[0],CubeValues[1]);
+					VertexList[0] = P;
+					Q = C0;
+				}
+				if (edgeTable[CubeIndex] & 2){
+					P = VertexInterp(C1,C2,CubeValues[1],CubeValues[2]);
+					VertexList[1] = P;
+					Q = C1;
+				}
+				if (edgeTable[CubeIndex] & 4){
+					P = VertexInterp(C2,C3,CubeValues[2],CubeValues[3]);
+					VertexList[2] =	P;
+					Q = C2;
+				}
+				if (edgeTable[CubeIndex] & 8){
+					P = VertexInterp(C3,C0,CubeValues[3],CubeValues[0]);
+					VertexList[3] =	P;
+					Q = C3;
+				}
+				if (edgeTable[CubeIndex] & 16){
+					P = VertexInterp(C4,C5,CubeValues[4],CubeValues[5]);
+					VertexList[4] =	P;
+					Q = C4;
+				}
+				if (edgeTable[CubeIndex] & 32){
+					P = VertexInterp(C5,C6,CubeValues[5],CubeValues[6]);
+					VertexList[5] = P;
+					Q = C5;
+				}
+				if (edgeTable[CubeIndex] & 64){ 
+					P = VertexInterp(C6,C7,CubeValues[6],CubeValues[7]);
+					VertexList[6] = P;
+					Q = C6;
+				}
+				if (edgeTable[CubeIndex] & 128){
+					P = VertexInterp(C7,C4,CubeValues[7],CubeValues[4]);
+					VertexList[7] =	P;
+					Q = C7;
+				}
+				if (edgeTable[CubeIndex] & 256){
+					P = VertexInterp(C0,C4,CubeValues[0],CubeValues[4]);
+					VertexList[8] =	P;
+					Q = C0;
+				}
+				if (edgeTable[CubeIndex] & 512){
+					P = VertexInterp(C1,C5,CubeValues[1],CubeValues[5]);
+					VertexList[9] =	P;
+					Q = C1;
+				}
+				if (edgeTable[CubeIndex] & 1024){
+					P = VertexInterp(C2,C6,CubeValues[2],CubeValues[6]);
+					VertexList[10] = P;
+					Q = C2;
+				}
+				if (edgeTable[CubeIndex] & 2048){
+					P = VertexInterp(C3,C7,CubeValues[3],CubeValues[7]);
+					VertexList[11] = P;
+					Q = C3;
+				}
+
+				NewVertexCount=0;
+				for (int idx=0;idx<12;idx++)
+					LocalRemap[idx] = -1;
+
+				for (int idx=0;triTable[CubeIndex][idx]!=-1;idx++)
+				{
+					if(LocalRemap[triTable[CubeIndex][idx]] == -1)
+					{
+						NewVertexList[NewVertexCount] = VertexList[triTable[CubeIndex][idx]];
+						LocalRemap[triTable[CubeIndex][idx]] = NewVertexCount;
+						NewVertexCount++;
+					}
+				}
+
+				for (int idx=0;idx<NewVertexCount;idx++) {
+					P = NewVertexList[idx];
+					//P.x  += i;
+					//P.y  += j;
+					//P.z  += k;
+					cellvertices(idx) = P;
+				}
+				nVert = NewVertexCount;
+
+				TriangleCount = 0;
+				for (int idx=0;triTable[CubeIndex][idx]!=-1;idx+=3) {
+					Triangles(0,TriangleCount) = LocalRemap[triTable[CubeIndex][idx+0]];
+					Triangles(1,TriangleCount) = LocalRemap[triTable[CubeIndex][idx+1]];
+					Triangles(2,TriangleCount) = LocalRemap[triTable[CubeIndex][idx+2]];
+					TriangleCount++;
+				}
+				nTris = TriangleCount;
+
+				// Now add the local values to the DECL data structure
+				IntArray HalfEdge(6,nTris*3);
+				DoubleArray EdgeAngles(nTris*3);
+				int idx_edge=0;
+				for (int idx=0; idx<TriangleCount; idx++){
+					int V1 = Triangles(0,idx);
+					int V2 = Triangles(1,idx);
+					int V3 = Triangles(2,idx);
+					// first edge: V1->V2 
+					HalfEdge(0,idx_edge) = V1;  // first vertex
+					HalfEdge(1,idx_edge) = V2;  // second vertex
+					HalfEdge(2,idx_edge) = idx; // triangle
+					HalfEdge(3,idx_edge) = -1;  // twin
+					HalfEdge(4,idx_edge) = idx_edge+2;  // previous edge
+					HalfEdge(5,idx_edge) = idx_edge+1;  // next edge
+					idx_edge++;
+					// second edge: V2->V3
+					HalfEdge(0,idx_edge) = V2;  // first vertex
+					HalfEdge(1,idx_edge) = V3;  // second vertex
+					HalfEdge(2,idx_edge) = idx; // triangle
+					HalfEdge(3,idx_edge) = -1;  // twin
+					HalfEdge(4,idx_edge) = idx_edge-1;  // previous edge
+					HalfEdge(5,idx_edge) = idx_edge+1;  // next edge
+					idx_edge++;
+					// third edge: V3->V1
+					HalfEdge(0,idx_edge) = V3;  // first vertex
+					HalfEdge(1,idx_edge) = V1;  // second vertex
+					HalfEdge(2,idx_edge) = idx; // triangle
+					HalfEdge(3,idx_edge) = -1;  // twin
+					HalfEdge(4,idx_edge) = idx_edge-1;  // previous edge
+					HalfEdge(5,idx_edge) = idx_edge-2;  // next edge
+					idx_edge++;
+				}
+				int EdgeCount=idx_edge;
+				for (int idx=0; idx<EdgeCount; idx++){
+					int V1=HalfEdge(0,idx);
+					int V2=HalfEdge(1,idx);
+					// Find all the twins within the cube
+					for (int jdx=0; idx<EdgeCount; jdx++){
+						if (HalfEdge(1,jdx) == V1 && HalfEdge(0,jdx) == V2){
+							// this is the pair
+							HalfEdge(3,idx) = jdx;
+							HalfEdge(3,jdx) = idx;
+						}
+						if (HalfEdge(1,jdx) == V2 && HalfEdge(0,jdx) == V1 && !(idx==jdx)){
+							std::printf("WARNING: half edges with identical orientation! \n");
+						}
+					}	
+					// Use "ghost" twins if edge is on a cube face
+					P = cellvertices(V1);
+					Q = cellvertices(V2);
+					if (P.x == 0.0 && Q.x == 0.0) HalfEdge(3,idx_edge) = -1;  // ghost twin for x=0 face
+					if (P.x == 1.0 && Q.x == 1.0) HalfEdge(3,idx_edge) = -2;  // ghost twin for x=1 face
+					if (P.y == 0.0 && Q.y == 0.0) HalfEdge(3,idx_edge) = -3;  // ghost twin for y=0 face
+					if (P.y == 1.0 && Q.y == 1.0) HalfEdge(3,idx_edge) = -4;  // ghost twin for y=1 face
+					if (P.z == 0.0 && Q.z == 0.0) HalfEdge(3,idx_edge) = -5;  // ghost twin for z=0 face
+					if (P.z == 1.0 && Q.z == 1.0) HalfEdge(3,idx_edge) = -6;  // ghost twin for z=1 face
+				}
+				// Find all the angles
+				for (int idx=0; idx<EdgeCount; idx++){
+					int V1=HalfEdge(0,idx);
+					int V2=HalfEdge(1,idx);
+					int T1= HalfEdge(2,idx_edge);
+					int twin=HalfEdge(3,idx_edge); 
+					if (twin == -1) ;
+				}
+				
+				// Map vertices to global coordinates
+				for (int idx=0;idx<NewVertexCount;idx++) {
+					P = cellvertices(idx);
+					P.x  += i;
+					P.y  += j;
+					P.z  += k;
+					cellvertices(idx) = P;
+				}
+			}
+		}
+	}
+
+}
+
+
+
diff --git a/analysis/decl.h b/analysis/decl.h
new file mode 100644
index 00000000..9b9df5d6
--- /dev/null
+++ b/analysis/decl.h
@@ -0,0 +1,56 @@
+#include <vector>
+
+/* 
+Doubly-connected edge list (DECL) 
+*/
+
+// Vertex structure
+class Vertex{
+public:
+	Vertex();
+	~Vertex();
+	void add(Point P);
+	void assign(unsigned long int idx, Point P);
+	Point coords(unsigned long int idx);
+	unsigned long int IncidentEdge();
+	unsigned long int count;
+private:
+	std::vector<double> vertex_data;
+};
+
+// Halfedge structure
+// Face
+class Halfedge{
+public:
+	Halfedge();
+	~Halfedge();
+
+	unsigned long int v1(unsigned long int edge);
+	unsigned long int v2(unsigned long int edge);
+	unsigned long int twin(unsigned long int edge);
+	unsigned long int face(unsigned long int edge);
+	unsigned long int next(unsigned long int edge);
+	unsigned long int prev(unsigned long int edge);
+	
+private:
+	Array<unsigned long int> HalfEdge;
+};
+
+// DECL
+class DECL{
+public:
+	DECL();
+	~DECL();
+	
+	unsigned long int face();
+	Vertex vertex;
+	Halfedge halfedge;
+	void AddCube(); // need a function to add new faces based on marching cubes surface
+	
+	double origin(int edge);
+	double EdgeAngle(int edge);
+	Point TriNormal(int edge);
+	
+private:
+	unsigned long int *face_data;
+};
diff --git a/tests/lbpm_uCT_pp.cpp b/tests/lbpm_uCT_pp.cpp
index 34d4825d..391d0f95 100644
--- a/tests/lbpm_uCT_pp.cpp
+++ b/tests/lbpm_uCT_pp.cpp
@@ -105,7 +105,7 @@ int main(int argc, char **argv)
 		//std::vector<size_t> ratio = {4,4,4};
 		// need to set up databases for each level of the mesh
 		std:vector<Database> multidomain_db;
-		
+
 		std::vector<int> Nx(1,nx), Ny(1,ny), Nz(1,nz);
 		while ( Nx.back()%ratio[0]==0 && Nx.back()>8 &&
 				Ny.back()%ratio[1]==0 && Ny.back()>8 &&
@@ -231,12 +231,12 @@ int main(int argc, char **argv)
 				for (int i=1;i<Nx[0]+1;i++) {
 					if (MASK(i,j,k) > 0.f ){
 						auto tmp = LOCVOL[0](i,j,k);
-/*						if ((tmp-background)*(tmp-target) > 0){
+						/*						if ((tmp-background)*(tmp-target) > 0){
 							// direction to background / target is the same
 							if (fabs(tmp-target) > fabs(tmp-background)) tmp=background; // tmp closer to background
 							else                             			 tmp=target;     // tmp closer to target
 						}
-						*/
+						 */
 						if ( tmp > THRESHOLD ) {
 							mean_plus += tmp;
 							count_plus++;