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

convex test for faces

Browse Source
This commit is contained in:
James E McClure 2018-09-14 18:21:36 -04:00
parent 7935d2ccb7
commit f29df4b8aa
1 changed files with 15 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

35
analysis/decl.cpp
View File

@ -383,33 +383,28 @@ double DECL::EdgeAngle(int edge)
double len = sqrt(nx*nx+ny*ny+nz*nz); double len = sqrt(nx*nx+ny*ny+nz*nz);
// new value for V is this normal vector // new value for V is this normal vector
V.x = nx/len; V.y = ny/len; V.z = nz/len; V.x = nx/len; V.y = ny/len; V.z = nz/len;
dotprod = U.x*V.x + U.y*V.y + U.z*V.z; /* project onto plane of cube face also works
if (dotprod > 1.f) dotprod=1.f; W = U - dotprod*V;
if (dotprod < 0.f) dotprod=-dotprod;
angle = acos(dotprod);
/*W = U - dotprod*V;
length = sqrt(W.x*W.x+W.y*W.y+W.z*W.z); // for normalization 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; 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;
if (dotprod < -1.f) dotprod=-1.f; if (dotprod < -1.f) dotprod=-1.f;
angle = acos(dotprod); angle = acos(dotprod);
*/ */
} }
else{ dotprod=U.x*V.x + U.y*V.y + U.z*V.z;
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; if (dotprod < -1.f) dotprod=-1.f;
if (dotprod < -1.f) dotprod=-1.f; angle = acos(dotprod);
angle = acos(dotprod); // determine if angle is concave or convex based on edge normal
// determine if angle is concave or convex based on edge normal W = 0.5*(P+Q)-R; // vector that lies in plane of triangle
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
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));
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){
if (length > hypotenuse){ // concave
// concave angle = -angle;
angle = -angle;
}
angle *= 0.5; // half edge value
} }
angle *= 0.5; // half edge value
//printf(" %f, %f (Edge=%i, twin=%i)\n U={%f, %f, %f}, V={%f, %f, %f}\n",angle,dotprod,edge,halfedge.twin(edge),U.x,U.y,U.z,V.x,V.y,V.z); //printf(" %f, %f (Edge=%i, twin=%i)\n U={%f, %f, %f}, V={%f, %f, %f}\n",angle,dotprod,edge,halfedge.twin(edge),U.x,U.y,U.z,V.x,V.y,V.z);
return angle; return angle;
} }