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

decl

Browse Source
This commit is contained in:
James E McClure 2018-09-15 06:24:00 -04:00
parent 8b04fad860
commit 1963fbfd31
2 changed files with 15 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
analysis/Minkowski.cpp
View File

@ -210,9 +210,9 @@ void Minkowski::ComputeScalar(const DoubleArray Field, const double isovalue)
//if (0.08333333333333*(a1*s1+a2*s2+a3*s3) < 0.f){ //if (0.08333333333333*(a1*s1+a2*s2+a3*s3) < 0.f){
//double intcurv=0.08333333333333*(a1*s1+a2*s2+a3*s3); //double intcurv=0.08333333333333*(a1*s1+a2*s2+a3*s3);
//double surfarea=sqrt(s*(s-s1)*(s-s2)*(s-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) 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) 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) 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)) // Euler characteristic (half edge rule: one face - 0.5*(three edges))
Xi -= 0.5; Xi -= 0.5;

21
analysis/decl.cpp
View File

@ -359,6 +359,7 @@ Point DECL::TriNormal(int edge)
double DECL::EdgeAngle(int edge) double DECL::EdgeAngle(int edge)
{ {
double angle; double angle;
double nx,ny,nz;
double dotprod,length,hypotenuse; double dotprod,length,hypotenuse;
Point P,Q,R; // triangle vertices Point P,Q,R; // triangle vertices
Point U,V,W; // normal vectors Point U,V,W; // normal vectors
@ -377,14 +378,15 @@ double DECL::EdgeAngle(int edge)
W.y /= length; W.y /= length;
W.z /= length; W.z /= length;
// edge normal within the plane of the cube face // edge normal within the plane of the cube face
double nx = W.y*V.z - W.z*V.y; nx = W.y*V.z - W.z*V.y;
double ny = W.z*V.x - W.x*V.z; ny = W.z*V.x - W.x*V.z;
double nz = W.x*V.y - W.y*V.x; nz = W.x*V.y - W.y*V.x;
double len = sqrt(nx*nx+ny*ny+nz*nz); length = 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/length; V.y = ny/length; V.z = nz/length;
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 < 0.f){ if (dotprod < 0.f){
//printf("negative dot product on face\n");
dotprod=-dotprod; dotprod=-dotprod;
V.x = -V.x; V.y = -V.y; V.z = -V.z; 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); angle = 0.5*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.x = (P.y-Q.y)*U.z - (P.z-Q.z)*U.y;
//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 W.y = (P.z-Q.z)*U.x - (P.x-Q.x)*U.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)); W.z = (P.x-Q.x)*U.y - (P.y-Q.y)*U.x;
//if (length > hypotenuse){ //length = sqrt(nx*nx+ny*ny+nz*nz);
if (W.x*V.x + W.y*V.y + W.z*V.z < 0.f){ if (W.x*V.x + W.y*V.y + W.z*V.z < 0.f){
// concave // concave
angle = -angle; angle = -angle;