From 1963fbfd3125949d570d6a4c8c32b5bba8d228ba Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Sat, 15 Sep 2018 06:24:00 -0400
Subject: [PATCH] decl

---
 analysis/Minkowski.cpp |  6 +++---
 analysis/decl.cpp      | 21 ++++++++++++---------
 2 files changed, 15 insertions(+), 12 deletions(-)

diff --git a/analysis/Minkowski.cpp b/analysis/Minkowski.cpp
index 61a9463c..352a7464 100644
--- a/analysis/Minkowski.cpp
+++ b/analysis/Minkowski.cpp
@@ -210,9 +210,9 @@ void Minkowski::ComputeScalar(const DoubleArray Field, const double isovalue)
 					//if (0.08333333333333*(a1*s1+a2*s2+a3*s3) < 0.f){
 					//double intcurv=0.08333333333333*(a1*s1+a2*s2+a3*s3);
 					//double surfarea=sqrt(s*(s-s1)*(s-s2)*(s-s3));
-					  //printf("   (%i,%i,%i) PQ(%i,%i)={%f,%f,%f} {%f,%f,%f} a=%f l=%f \n",i,j,k,e1,object.halfedge.twin(e1),P1.x,P1.y,P1.z,P2.x,P2.y,P2.z,a1,s1);
-					  //printf("   (%i,%i,%i) QR(%i,%i)={%f,%f,%f} {%f,%f,%f} a=%f l=%f \n",i,j,k,e2,object.halfedge.twin(e2),P2.x,P2.y,P2.z,P3.x,P3.y,P3.z,a2,s2);
-					  //printf("   (%i,%i,%i) RP(%i,%i)={%f,%f,%f} {%f,%f,%f} a=%f l=%f \n",i,j,k,e3,object.halfedge.twin(e3),P3.x,P3.y,P3.z,P1.x,P1.y,P1.z,a3,s3);
+					//printf("   (%i,%i,%i) PQ(%i,%i)={%f,%f,%f} {%f,%f,%f} a=%f l=%f \n",i,j,k,e1,object.halfedge.twin(e1),P1.x,P1.y,P1.z,P2.x,P2.y,P2.z,a1,s1);
+					// printf("   (%i,%i,%i) QR(%i,%i)={%f,%f,%f} {%f,%f,%f} a=%f l=%f \n",i,j,k,e2,object.halfedge.twin(e2),P2.x,P2.y,P2.z,P3.x,P3.y,P3.z,a2,s2);
+					// printf("   (%i,%i,%i) RP(%i,%i)={%f,%f,%f} {%f,%f,%f} a=%f l=%f \n",i,j,k,e3,object.halfedge.twin(e3),P3.x,P3.y,P3.z,P1.x,P1.y,P1.z,a3,s3);
 					  //}
 					// Euler characteristic (half edge rule: one face - 0.5*(three edges))
 					Xi -= 0.5;
diff --git a/analysis/decl.cpp b/analysis/decl.cpp
index 1c84a62d..2c41d125 100644
--- a/analysis/decl.cpp
+++ b/analysis/decl.cpp
@@ -359,6 +359,7 @@ Point DECL::TriNormal(int edge)
 double DECL::EdgeAngle(int edge)
 {
 	double angle;
+	double nx,ny,nz;
 	double dotprod,length,hypotenuse;
 	Point P,Q,R; // triangle vertices
 	Point U,V,W; // normal vectors
@@ -377,14 +378,15 @@ double DECL::EdgeAngle(int edge)
 		W.y /= length;
 		W.z /= length;
 		// edge normal within the plane of the cube face
-		double nx = W.y*V.z - W.z*V.y;
-		double ny = W.z*V.x - W.x*V.z;
-		double nz = W.x*V.y - W.y*V.x;
-		double len = sqrt(nx*nx+ny*ny+nz*nz);
+		nx = W.y*V.z - W.z*V.y;
+		ny = W.z*V.x - W.x*V.z;
+		nz = W.x*V.y - W.y*V.x;
+		length = sqrt(nx*nx+ny*ny+nz*nz);
 		// new value for V is this normal vector
-		V.x = nx/len; V.y = ny/len; V.z = nz/len;
+		V.x = nx/length; V.y = ny/length; V.z = nz/length;
 		dotprod = U.x*V.x + U.y*V.y + U.z*V.z;
 		if (dotprod < 0.f){
+		  //printf("negative dot product on face\n");
 			dotprod=-dotprod;
 			V.x = -V.x; V.y = -V.y; V.z = -V.z;
 		}
@@ -406,10 +408,11 @@ double DECL::EdgeAngle(int edge)
 		angle = 0.5*acos(dotprod);
 	}
 	// determine if angle is concave or convex based on edge normal
-	W = 0.5*(P+Q)-R; // vector that lies in plane of triangle
-	//hypotenuse = sqrt(W.x*W.x+W.y*W.y+W.z*W.z + V.x*V.x+V.y*V.y+V.z*V.z); // hypotenuse of right triangle
-	//length = sqrt((W.x+V.x)*(W.x+V.x) + (W.y+V.y)*(W.y+V.y) + (W.z+V.z)*(W.z+V.z));
-	//if (length > hypotenuse){
+	W.x = (P.y-Q.y)*U.z - (P.z-Q.z)*U.y;
+	W.y = (P.z-Q.z)*U.x - (P.x-Q.x)*U.z;
+	W.z = (P.x-Q.x)*U.y - (P.y-Q.y)*U.x;
+        //length = sqrt(nx*nx+ny*ny+nz*nz);
+
 	if (W.x*V.x + W.y*V.y + W.z*V.z < 0.f){
 		// concave
 		angle = -angle;