From 5b4f72497966e6585aa7fef2af50c117ddd6966d Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Wed, 12 Sep 2018 21:00:34 -0400
Subject: [PATCH] concave / convex test

---
 analysis/decl.cpp | 15 ++++++++++++++-
 1 file changed, 14 insertions(+), 1 deletion(-)

diff --git a/analysis/decl.cpp b/analysis/decl.cpp
index 0acd352e..131f2c53 100644
--- a/analysis/decl.cpp
+++ b/analysis/decl.cpp
@@ -360,12 +360,25 @@ Point DECL::TriNormal(int edge)
 double DECL::EdgeAngle(int edge)
 {
 	double angle;
-	Point U,V; // triangle normal vectors
+	point P,Q,R; // triangle corners
+	Point U,V,W; // triangle normal vectors
 	U = TriNormal(edge);
 	V = TriNormal(halfedge.twin(edge));
 	double dotprod=U.x*V.x + U.y*V.y + U.z*V.z;
 	if (dotprod > 1.f) dotprod=1.f;
 	angle =  acos(dotprod);
+	// triangle corners
+	P=halfedge.v1(edge);
+	Q=halfedge.v2(edge);
+	R=halfedge.v2(halfedge.next(edge));
+	// determine if angle is concave or convex based on edge normal
+	W = 0.5*(P+Q)-R; // vector that lies in plane of triangle
+	double hypotenuse = sqrt(W.x*V.x+W.y*V.y+W.z*V.z); // hypotenuse of right triangle
+	double 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){
+		// concave
+		angle = -angle;
+	}
 	if (!(edge<0)) angle *= 0.5; // half edge value
 	//1.570796326794897
 	//printf("  %f, %f: U={%f, %f, %f}, V={%f, %f, %f}\n",angle,dotprod,U.y,U.z,V.x,V.y,V.z);