OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

using different edge normal

Browse Source
This commit is contained in:
James E McClure
2018-09-14 16:34:32 -04:00
parent 07538d26ad
commit 62e1614812
1 changed files with 28 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files

42
analysis/decl.cpp
View File

@@ -359,33 +359,47 @@ Point DECL::TriNormal(int edge)
double DECL::EdgeAngle(int edge)
{
double angle;
double length,hypotenuse;
Point P,Q,R; // triangle corners
Point U,V,W; // triangle normal vectors
double dotprod,length,hypotenuse;
Point P,Q,R; // triangle vertices
Point U,V,W; // normal vectors
int e2 = halfedge.next(edge);
int e3 = halfedge.next(e2);
P=vertex.coords(halfedge.v1(edge));
Q=vertex.coords(halfedge.v1(e2));
R=vertex.coords(halfedge.v1(e3));
U = TriNormal(edge);
V = TriNormal(halfedge.twin(edge));
double dotprod=U.x*V.x + U.y*V.y + U.z*V.z;
if (halfedge.twin(edge) < 0 ){
// compute projection onto plane
W = U - dotprod*V;
// compute edge normal in plane of cube face
W = Q - P; // edge tangent vector
length = sqrt(W.x*W.x+W.y*W.y+W.z*W.z);
W.x /= length;
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);
// new value for V is this normal vector
V.x = nx/len; V.y = ny/len; V.z = nz/len;
dotprod = dotprod=U.x*V.x + U.y*V.y + U.z*V.z;
if (dotprod > 1.f) dotprod=1.f;
if (dotprod < -1.f) dotprod=-1.f;
angle = acos(dotprod);
/*W = U - dotprod*V;
length = sqrt(W.x*W.x+W.y*W.y+W.z*W.z); // for normalization
dotprod = (U.x*W.x + U.y*W.y + U.z*W.z)/length;
if (dotprod > 1.f) dotprod=1.f;
if (dotprod < -1.f) dotprod=-1.f;
angle = acos(dotprod);
// turn in outward normal at cube face is pi/2 from each side of the cube
//angle-=1.570796326794897;
*/
}
else{
dotprod=U.x*V.x + U.y*V.y + U.z*V.z
if (dotprod > 1.f) dotprod=1.f;
if (dotprod < -1.f) dotprod=-1.f;
angle = acos(dotprod);
// triangle corners
int e2 = halfedge.next(edge);
int e3 = halfedge.next(e2);
P=vertex.coords(halfedge.v1(edge));
Q=vertex.coords(halfedge.v1(e2));
R=vertex.coords(halfedge.v1(e3));
// 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