diff --git a/analysis/decl.cpp b/analysis/decl.cpp index e817a877..db149b9d 100644 --- a/analysis/decl.cpp +++ b/analysis/decl.cpp @@ -383,33 +383,28 @@ double DECL::EdgeAngle(int edge) 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 = U.x*V.x + U.y*V.y + U.z*V.z; - if (dotprod > 1.f) dotprod=1.f; - if (dotprod < 0.f) dotprod=-dotprod; - angle = acos(dotprod); - /*W = U - dotprod*V; + /* project onto plane of cube face also works + 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); - */ + */ } - 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); - // 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){ - // concave - angle = -angle; - } - angle *= 0.5; // half edge value + 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); + // 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){ + // concave + angle = -angle; } + 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); return angle; }