diff options
author | Campbell Barton <ideasman42@gmail.com> | 2019-04-17 07:17:24 +0300 |
---|---|---|
committer | Campbell Barton <ideasman42@gmail.com> | 2019-04-17 07:21:24 +0300 |
commit | e12c08e8d170b7ca40f204a5b0423c23a9fbc2c1 (patch) | |
tree | 8cf3453d12edb177a218ef8009357518ec6cab6a /source/blender/freestyle/intern/winged_edge/Curvature.cpp | |
parent | b3dabc200a4b0399ec6b81f2ff2730d07b44fcaa (diff) |
ClangFormat: apply to source, most of intern
Apply clang format as proposed in T53211.
For details on usage and instructions for migrating branches
without conflicts, see:
https://wiki.blender.org/wiki/Tools/ClangFormat
Diffstat (limited to 'source/blender/freestyle/intern/winged_edge/Curvature.cpp')
-rw-r--r-- | source/blender/freestyle/intern/winged_edge/Curvature.cpp | 880 |
1 files changed, 443 insertions, 437 deletions
diff --git a/source/blender/freestyle/intern/winged_edge/Curvature.cpp b/source/blender/freestyle/intern/winged_edge/Curvature.cpp index 28b38ae4854..b42c3f9b0e6 100644 --- a/source/blender/freestyle/intern/winged_edge/Curvature.cpp +++ b/source/blender/freestyle/intern/winged_edge/Curvature.cpp @@ -34,7 +34,7 @@ */ #include <assert.h> -#include <cstdlib> // for malloc and free +#include <cstdlib> // for malloc and free #include <set> #include <stack> @@ -49,56 +49,56 @@ namespace Freestyle { static bool angle_obtuse(WVertex *v, WFace *f) { - WOEdge *e; - f->getOppositeEdge(v, e); + WOEdge *e; + f->getOppositeEdge(v, e); - Vec3r vec1(e->GetaVertex()->GetVertex() - v->GetVertex()); - Vec3r vec2(e->GetbVertex()->GetVertex() - v->GetVertex()); - return ((vec1 * vec2) < 0); + Vec3r vec1(e->GetaVertex()->GetVertex() - v->GetVertex()); + Vec3r vec2(e->GetbVertex()->GetVertex() - v->GetVertex()); + return ((vec1 * vec2) < 0); } // FIXME // WVvertex is useless but kept for history reasons static bool triangle_obtuse(WVertex *, WFace *f) { - bool b = false; - for (int i = 0; i < 3; i++) - b = b || ((f->getEdgeList()[i]->GetVec() * f->getEdgeList()[(i + 1) % 3]->GetVec()) < 0); - return b; + bool b = false; + for (int i = 0; i < 3; i++) + b = b || ((f->getEdgeList()[i]->GetVec() * f->getEdgeList()[(i + 1) % 3]->GetVec()) < 0); + return b; } static real cotan(WVertex *vo, WVertex *v1, WVertex *v2) { - /* cf. Appendix B of [Meyer et al 2002] */ - real udotv, denom; + /* cf. Appendix B of [Meyer et al 2002] */ + real udotv, denom; - Vec3r u(v1->GetVertex() - vo->GetVertex()); - Vec3r v(v2->GetVertex() - vo->GetVertex()); + Vec3r u(v1->GetVertex() - vo->GetVertex()); + Vec3r v(v2->GetVertex() - vo->GetVertex()); - udotv = u * v; - denom = sqrt(u.squareNorm() * v.squareNorm() - udotv * udotv); + udotv = u * v; + denom = sqrt(u.squareNorm() * v.squareNorm() - udotv * udotv); - /* denom can be zero if u==v. Returning 0 is acceptable, based on the callers of this function below. */ - if (denom == 0.0) - return 0.0; - return (udotv / denom); + /* denom can be zero if u==v. Returning 0 is acceptable, based on the callers of this function below. */ + if (denom == 0.0) + return 0.0; + return (udotv / denom); } static real angle_from_cotan(WVertex *vo, WVertex *v1, WVertex *v2) { - /* cf. Appendix B and the caption of Table 1 from [Meyer et al 2002] */ - real udotv, denom; + /* cf. Appendix B and the caption of Table 1 from [Meyer et al 2002] */ + real udotv, denom; - Vec3r u (v1->GetVertex() - vo->GetVertex()); - Vec3r v(v2->GetVertex() - vo->GetVertex()); + Vec3r u(v1->GetVertex() - vo->GetVertex()); + Vec3r v(v2->GetVertex() - vo->GetVertex()); - udotv = u * v; - denom = sqrt(u.squareNorm() * v.squareNorm() - udotv * udotv); + udotv = u * v; + denom = sqrt(u.squareNorm() * v.squareNorm() - udotv * udotv); - /* Note: I assume this is what they mean by using atan2(). -Ray Jones */ + /* Note: I assume this is what they mean by using atan2(). -Ray Jones */ - /* tan = denom/udotv = y/x (see man page for atan2) */ - return (fabs(atan2(denom, udotv))); + /* tan = denom/udotv = y/x (see man page for atan2) */ + return (fabs(atan2(denom, udotv))); } /*! gts_vertex_mean_curvature_normal: @@ -124,48 +124,48 @@ static real angle_from_cotan(WVertex *vo, WVertex *v1, WVertex *v2) */ bool gts_vertex_mean_curvature_normal(WVertex *v, Vec3r &Kh) { - real area = 0.0; + real area = 0.0; - if (!v) - return false; + if (!v) + return false; - /* this operator is not defined for boundary edges */ - if (v->isBoundary()) - return false; + /* this operator is not defined for boundary edges */ + if (v->isBoundary()) + return false; - WVertex::incoming_edge_iterator itE; + WVertex::incoming_edge_iterator itE; - for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++) - area += (*itE)->GetaFace()->getArea(); + for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++) + area += (*itE)->GetaFace()->getArea(); - Kh = Vec3r(0.0, 0.0, 0.0); + Kh = Vec3r(0.0, 0.0, 0.0); - for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++) { - WOEdge *e = (*itE)->getPrevOnFace(); + for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++) { + WOEdge *e = (*itE)->getPrevOnFace(); #if 0 - if ((e->GetaVertex() == v) || (e->GetbVertex() == v)) - cerr<< "BUG "; + if ((e->GetaVertex() == v) || (e->GetbVertex() == v)) + cerr<< "BUG "; #endif - WVertex *v1 = e->GetaVertex(); - WVertex *v2 = e->GetbVertex(); - real temp; - - temp = cotan(v1, v, v2); - Kh = Vec3r(Kh + temp * (v2->GetVertex() - v->GetVertex())); - - temp = cotan(v2, v, v1); - Kh = Vec3r(Kh + temp * (v1->GetVertex() - v->GetVertex())); - } - if (area > 0.0) { - Kh[0] /= 2 * area; - Kh[1] /= 2 * area; - Kh[2] /= 2 * area; - } - else { - return false; - } - - return true; + WVertex *v1 = e->GetaVertex(); + WVertex *v2 = e->GetbVertex(); + real temp; + + temp = cotan(v1, v, v2); + Kh = Vec3r(Kh + temp * (v2->GetVertex() - v->GetVertex())); + + temp = cotan(v2, v, v1); + Kh = Vec3r(Kh + temp * (v1->GetVertex() - v->GetVertex())); + } + if (area > 0.0) { + Kh[0] /= 2 * area; + Kh[1] /= 2 * area; + Kh[2] /= 2 * area; + } + else { + return false; + } + + return true; } /*! gts_vertex_gaussian_curvature: @@ -186,35 +186,35 @@ bool gts_vertex_mean_curvature_normal(WVertex *v, Vec3r &Kh) */ bool gts_vertex_gaussian_curvature(WVertex *v, real *Kg) { - real area = 0.0; - real angle_sum = 0.0; - - if (!v) - return false; - if (!Kg) - return false; - - /* this operator is not defined for boundary edges */ - if (v->isBoundary()) { - *Kg = 0.0; - return false; - } - - WVertex::incoming_edge_iterator itE; - for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++) { - area += (*itE)->GetaFace()->getArea(); - } - - for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++) { - WOEdge *e = (*itE)->getPrevOnFace(); - WVertex *v1 = e->GetaVertex(); - WVertex *v2 = e->GetbVertex(); - angle_sum += angle_from_cotan(v, v1, v2); - } - - *Kg = (2.0 * M_PI - angle_sum) / area; - - return true; + real area = 0.0; + real angle_sum = 0.0; + + if (!v) + return false; + if (!Kg) + return false; + + /* this operator is not defined for boundary edges */ + if (v->isBoundary()) { + *Kg = 0.0; + return false; + } + + WVertex::incoming_edge_iterator itE; + for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++) { + area += (*itE)->GetaFace()->getArea(); + } + + for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++) { + WOEdge *e = (*itE)->getPrevOnFace(); + WVertex *v1 = e->GetaVertex(); + WVertex *v2 = e->GetbVertex(); + angle_sum += angle_from_cotan(v, v1, v2); + } + + *Kg = (2.0 * M_PI - angle_sum) / area; + + return true; } /*! gts_vertex_principal_curvatures: @@ -230,41 +230,41 @@ bool gts_vertex_gaussian_curvature(WVertex *v, real *Kg) * * The Gaussian curvature can be computed with gts_vertex_gaussian_curvature(). */ -void gts_vertex_principal_curvatures (real Kh, real Kg, real *K1, real *K2) +void gts_vertex_principal_curvatures(real Kh, real Kg, real *K1, real *K2) { - real temp = Kh * Kh - Kg; + real temp = Kh * Kh - Kg; - if (!K1 || !K2) - return; + if (!K1 || !K2) + return; - if (temp < 0.0) - temp = 0.0; - temp = sqrt (temp); - *K1 = Kh + temp; - *K2 = Kh - temp; + if (temp < 0.0) + temp = 0.0; + temp = sqrt(temp); + *K1 = Kh + temp; + *K2 = Kh - temp; } /* from Maple */ static void linsolve(real m11, real m12, real b1, real m21, real m22, real b2, real *x1, real *x2) { - real temp; + real temp; - temp = 1.0 / (m21 * m12 - m11 * m22); - *x1 = (m12 * b2 - m22 * b1) * temp; - *x2 = (m11 * b2 - m21 * b1) * temp; + temp = 1.0 / (m21 * m12 - m11 * m22); + *x1 = (m12 * b2 - m22 * b1) * temp; + *x2 = (m11 * b2 - m21 * b1) * temp; } /* from Maple - largest eigenvector of [a b; b c] */ static void eigenvector(real a, real b, real c, Vec3r e) { - if (b == 0.0) { - e[0] = 0.0; - } - else { - e[0] = -(c - a - sqrt(c * c - 2 * a * c + a * a + 4 * b * b)) / (2 * b); - } - e[1] = 1.0; - e[2] = 0.0; + if (b == 0.0) { + e[0] = 0.0; + } + else { + e[0] = -(c - a - sqrt(c * c - 2 * a * c + a * a + 4 * b * b)) / (2 * b); + } + e[1] = 1.0; + e[2] = 0.0; } /*! gts_vertex_principal_directions: @@ -284,244 +284,249 @@ static void eigenvector(real a, real b, real c, Vec3r e) */ void gts_vertex_principal_directions(WVertex *v, Vec3r Kh, real Kg, Vec3r &e1, Vec3r &e2) { - Vec3r N; - real normKh; - - Vec3r basis1, basis2, d, eig; - real ve2, vdotN; - real aterm_da, bterm_da, cterm_da, const_da; - real aterm_db, bterm_db, cterm_db, const_db; - real a, b, c; - real K1, K2; - real *weights, *kappas, *d1s, *d2s; - int edge_count; - real err_e1, err_e2; - int e; - WVertex::incoming_edge_iterator itE; - - /* compute unit normal */ - normKh = Kh.norm(); - - if (normKh > 0.0) { - Kh.normalize(); - } - else { - /* This vertex is a point of zero mean curvature (flat or saddle point). Compute a normal by averaging - * the adjacent triangles - */ - N[0] = N[1] = N[2] = 0.0; - - for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++) - N = Vec3r(N + (*itE)->GetaFace()->GetNormal()); - real normN = N.norm(); - if (normN <= 0.0) - return; - N.normalize(); - } - - /* construct a basis from N: */ - /* set basis1 to any component not the largest of N */ - basis1[0] = basis1[1] = basis1[2] = 0.0; - if (fabs (N[0]) > fabs (N[1])) - basis1[1] = 1.0; - else - basis1[0] = 1.0; - - /* make basis2 orthogonal to N */ - basis2 = (N ^ basis1); - basis2.normalize(); - - /* make basis1 orthogonal to N and basis2 */ - basis1 = (N ^ basis2); - basis1.normalize(); - - aterm_da = bterm_da = cterm_da = const_da = 0.0; - aterm_db = bterm_db = cterm_db = const_db = 0.0; - int nb_edges = v->GetEdges().size(); - - weights = (real *)malloc(sizeof(real) * nb_edges); - kappas = (real *)malloc(sizeof(real) * nb_edges); - d1s = (real *)malloc(sizeof(real) * nb_edges); - d2s = (real *)malloc(sizeof(real) * nb_edges); - edge_count = 0; - - for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++) { - WOEdge *e; - WFace *f1, *f2; - real weight, kappa, d1, d2; - Vec3r vec_edge; - if (!*itE) - continue; - e = *itE; - - /* since this vertex passed the tests in gts_vertex_mean_curvature_normal(), this should be true. */ - //g_assert(gts_edge_face_number (e, s) == 2); - - /* identify the two triangles bordering e in s */ - f1 = e->GetaFace(); - f2 = e->GetbFace(); - - /* We are solving for the values of the curvature tensor - * B = [ a b ; b c ]. - * The computations here are from section 5 of [Meyer et al 2002]. - * - * The first step is to calculate the linear equations governing the values of (a,b,c). These can be computed - * by setting the derivatives of the error E to zero (section 5.3). - * - * Since a + c = norm(Kh), we only compute the linear equations for dE/da and dE/db. (NB: [Meyer et al 2002] - * has the equation a + b = norm(Kh), but I'm almost positive this is incorrect). - * - * Note that the w_ij (defined in section 5.2) are all scaled by (1/8*A_mixed). We drop this uniform scale - * factor because the solution of the linear equations doesn't rely on it. - * - * The terms of the linear equations are xterm_dy with x in {a,b,c} and y in {a,b}. There are also const_dy - * terms that are the constant factors in the equations. - */ - - /* find the vector from v along edge e */ - vec_edge = Vec3r(-1 * e->GetVec()); - - ve2 = vec_edge.squareNorm(); - vdotN = vec_edge * N; - - /* section 5.2 - There is a typo in the computation of kappa. The edges should be x_j-x_i. */ - kappa = 2.0 * vdotN / ve2; - - /* section 5.2 */ - - /* I don't like performing a minimization where some of the weights can be negative (as can be the case - * if f1 or f2 are obtuse). To ensure all-positive weights, we check for obtuseness. */ - weight = 0.0; - if (!triangle_obtuse(v, f1)) { - weight += ve2 * cotan(f1->GetNextOEdge(e->twin())->GetbVertex(), e->GetaVertex(), e->GetbVertex()) / 8.0; - } - else { - if (angle_obtuse(v, f1)) { - weight += ve2 * f1->getArea() / 4.0; - } - else { - weight += ve2 * f1->getArea() / 8.0; - } - } - - if (!triangle_obtuse(v, f2)) { - weight += ve2 * cotan (f2->GetNextOEdge(e)->GetbVertex(), e->GetaVertex(), e->GetbVertex()) / 8.0; - } - else { - if (angle_obtuse(v, f2)) { - weight += ve2 * f1->getArea() / 4.0; - } - else { - weight += ve2 * f1->getArea() / 8.0; - } - } - - /* projection of edge perpendicular to N (section 5.3) */ - d[0] = vec_edge[0] - vdotN * N[0]; - d[1] = vec_edge[1] - vdotN * N[1]; - d[2] = vec_edge[2] - vdotN * N[2]; - d.normalize(); - - /* not explicit in the paper, but necessary. Move d to 2D basis. */ - d1 = d * basis1; - d2 = d * basis2; - - /* store off the curvature, direction of edge, and weights for later use */ - weights[edge_count] = weight; - kappas[edge_count] = kappa; - d1s[edge_count] = d1; - d2s[edge_count] = d2; - edge_count++; - - /* Finally, update the linear equations */ - aterm_da += weight * d1 * d1 * d1 * d1; - bterm_da += weight * d1 * d1 * 2 * d1 * d2; - cterm_da += weight * d1 * d1 * d2 * d2; - const_da += weight * d1 * d1 * (-kappa); - - aterm_db += weight * d1 * d2 * d1 * d1; - bterm_db += weight * d1 * d2 * 2 * d1 * d2; - cterm_db += weight * d1 * d2 * d2 * d2; - const_db += weight * d1 * d2 * (-kappa); - } - - /* now use the identity (Section 5.3) a + c = |Kh| = 2 * kappa_h */ - aterm_da -= cterm_da; - const_da += cterm_da * normKh; - - aterm_db -= cterm_db; - const_db += cterm_db * normKh; - - /* check for solvability of the linear system */ - if (((aterm_da * bterm_db - aterm_db * bterm_da) != 0.0) && ((const_da != 0.0) || (const_db != 0.0))) { - linsolve(aterm_da, bterm_da, -const_da, aterm_db, bterm_db, -const_db, &a, &b); - - c = normKh - a; - - eigenvector(a, b, c, eig); - } - else { - /* region of v is planar */ - eig[0] = 1.0; - eig[1] = 0.0; - } - - /* Although the eigenvectors of B are good estimates of the principal directions, it seems that which one is - * attached to which curvature direction is a bit arbitrary. This may be a bug in my implementation, or just - * a side-effect of the inaccuracy of B due to the discrete nature of the sampling. - * - * To overcome this behavior, we'll evaluate which assignment best matches the given eigenvectors by comparing - * the curvature estimates computed above and the curvatures calculated from the discrete differential operators. - */ - - gts_vertex_principal_curvatures(0.5 * normKh, Kg, &K1, &K2); - - err_e1 = err_e2 = 0.0; - /* loop through the values previously saved */ - for (e = 0; e < edge_count; e++) { - real weight, kappa, d1, d2; - real temp1, temp2; - real delta; - - weight = weights[e]; - kappa = kappas[e]; - d1 = d1s[e]; - d2 = d2s[e]; - - temp1 = fabs (eig[0] * d1 + eig[1] * d2); - temp1 = temp1 * temp1; - temp2 = fabs (eig[1] * d1 - eig[0] * d2); - temp2 = temp2 * temp2; - - /* err_e1 is for K1 associated with e1 */ - delta = K1 * temp1 + K2 * temp2 - kappa; - err_e1 += weight * delta * delta; - - /* err_e2 is for K1 associated with e2 */ - delta = K2 * temp1 + K1 * temp2 - kappa; - err_e2 += weight * delta * delta; - } - free (weights); - free (kappas); - free (d1s); - free (d2s); - - /* rotate eig by a right angle if that would decrease the error */ - if (err_e2 < err_e1) { - real temp = eig[0]; - - eig[0] = eig[1]; - eig[1] = -temp; - } - - e1[0] = eig[0] * basis1[0] + eig[1] * basis2[0]; - e1[1] = eig[0] * basis1[1] + eig[1] * basis2[1]; - e1[2] = eig[0] * basis1[2] + eig[1] * basis2[2]; - e1.normalize(); - - /* make N,e1,e2 a right handed coordinate sytem */ - e2 = N ^ e1; - e2.normalize(); + Vec3r N; + real normKh; + + Vec3r basis1, basis2, d, eig; + real ve2, vdotN; + real aterm_da, bterm_da, cterm_da, const_da; + real aterm_db, bterm_db, cterm_db, const_db; + real a, b, c; + real K1, K2; + real *weights, *kappas, *d1s, *d2s; + int edge_count; + real err_e1, err_e2; + int e; + WVertex::incoming_edge_iterator itE; + + /* compute unit normal */ + normKh = Kh.norm(); + + if (normKh > 0.0) { + Kh.normalize(); + } + else { + /* This vertex is a point of zero mean curvature (flat or saddle point). Compute a normal by averaging + * the adjacent triangles + */ + N[0] = N[1] = N[2] = 0.0; + + for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++) + N = Vec3r(N + (*itE)->GetaFace()->GetNormal()); + real normN = N.norm(); + if (normN <= 0.0) + return; + N.normalize(); + } + + /* construct a basis from N: */ + /* set basis1 to any component not the largest of N */ + basis1[0] = basis1[1] = basis1[2] = 0.0; + if (fabs(N[0]) > fabs(N[1])) + basis1[1] = 1.0; + else + basis1[0] = 1.0; + + /* make basis2 orthogonal to N */ + basis2 = (N ^ basis1); + basis2.normalize(); + + /* make basis1 orthogonal to N and basis2 */ + basis1 = (N ^ basis2); + basis1.normalize(); + + aterm_da = bterm_da = cterm_da = const_da = 0.0; + aterm_db = bterm_db = cterm_db = const_db = 0.0; + int nb_edges = v->GetEdges().size(); + + weights = (real *)malloc(sizeof(real) * nb_edges); + kappas = (real *)malloc(sizeof(real) * nb_edges); + d1s = (real *)malloc(sizeof(real) * nb_edges); + d2s = (real *)malloc(sizeof(real) * nb_edges); + edge_count = 0; + + for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++) { + WOEdge *e; + WFace *f1, *f2; + real weight, kappa, d1, d2; + Vec3r vec_edge; + if (!*itE) + continue; + e = *itE; + + /* since this vertex passed the tests in gts_vertex_mean_curvature_normal(), this should be true. */ + //g_assert(gts_edge_face_number (e, s) == 2); + + /* identify the two triangles bordering e in s */ + f1 = e->GetaFace(); + f2 = e->GetbFace(); + + /* We are solving for the values of the curvature tensor + * B = [ a b ; b c ]. + * The computations here are from section 5 of [Meyer et al 2002]. + * + * The first step is to calculate the linear equations governing the values of (a,b,c). These can be computed + * by setting the derivatives of the error E to zero (section 5.3). + * + * Since a + c = norm(Kh), we only compute the linear equations for dE/da and dE/db. (NB: [Meyer et al 2002] + * has the equation a + b = norm(Kh), but I'm almost positive this is incorrect). + * + * Note that the w_ij (defined in section 5.2) are all scaled by (1/8*A_mixed). We drop this uniform scale + * factor because the solution of the linear equations doesn't rely on it. + * + * The terms of the linear equations are xterm_dy with x in {a,b,c} and y in {a,b}. There are also const_dy + * terms that are the constant factors in the equations. + */ + + /* find the vector from v along edge e */ + vec_edge = Vec3r(-1 * e->GetVec()); + + ve2 = vec_edge.squareNorm(); + vdotN = vec_edge * N; + + /* section 5.2 - There is a typo in the computation of kappa. The edges should be x_j-x_i. */ + kappa = 2.0 * vdotN / ve2; + + /* section 5.2 */ + + /* I don't like performing a minimization where some of the weights can be negative (as can be the case + * if f1 or f2 are obtuse). To ensure all-positive weights, we check for obtuseness. */ + weight = 0.0; + if (!triangle_obtuse(v, f1)) { + weight += ve2 * + cotan( + f1->GetNextOEdge(e->twin())->GetbVertex(), e->GetaVertex(), e->GetbVertex()) / + 8.0; + } + else { + if (angle_obtuse(v, f1)) { + weight += ve2 * f1->getArea() / 4.0; + } + else { + weight += ve2 * f1->getArea() / 8.0; + } + } + + if (!triangle_obtuse(v, f2)) { + weight += ve2 * cotan(f2->GetNextOEdge(e)->GetbVertex(), e->GetaVertex(), e->GetbVertex()) / + 8.0; + } + else { + if (angle_obtuse(v, f2)) { + weight += ve2 * f1->getArea() / 4.0; + } + else { + weight += ve2 * f1->getArea() / 8.0; + } + } + + /* projection of edge perpendicular to N (section 5.3) */ + d[0] = vec_edge[0] - vdotN * N[0]; + d[1] = vec_edge[1] - vdotN * N[1]; + d[2] = vec_edge[2] - vdotN * N[2]; + d.normalize(); + + /* not explicit in the paper, but necessary. Move d to 2D basis. */ + d1 = d * basis1; + d2 = d * basis2; + + /* store off the curvature, direction of edge, and weights for later use */ + weights[edge_count] = weight; + kappas[edge_count] = kappa; + d1s[edge_count] = d1; + d2s[edge_count] = d2; + edge_count++; + + /* Finally, update the linear equations */ + aterm_da += weight * d1 * d1 * d1 * d1; + bterm_da += weight * d1 * d1 * 2 * d1 * d2; + cterm_da += weight * d1 * d1 * d2 * d2; + const_da += weight * d1 * d1 * (-kappa); + + aterm_db += weight * d1 * d2 * d1 * d1; + bterm_db += weight * d1 * d2 * 2 * d1 * d2; + cterm_db += weight * d1 * d2 * d2 * d2; + const_db += weight * d1 * d2 * (-kappa); + } + + /* now use the identity (Section 5.3) a + c = |Kh| = 2 * kappa_h */ + aterm_da -= cterm_da; + const_da += cterm_da * normKh; + + aterm_db -= cterm_db; + const_db += cterm_db * normKh; + + /* check for solvability of the linear system */ + if (((aterm_da * bterm_db - aterm_db * bterm_da) != 0.0) && + ((const_da != 0.0) || (const_db != 0.0))) { + linsolve(aterm_da, bterm_da, -const_da, aterm_db, bterm_db, -const_db, &a, &b); + + c = normKh - a; + + eigenvector(a, b, c, eig); + } + else { + /* region of v is planar */ + eig[0] = 1.0; + eig[1] = 0.0; + } + + /* Although the eigenvectors of B are good estimates of the principal directions, it seems that which one is + * attached to which curvature direction is a bit arbitrary. This may be a bug in my implementation, or just + * a side-effect of the inaccuracy of B due to the discrete nature of the sampling. + * + * To overcome this behavior, we'll evaluate which assignment best matches the given eigenvectors by comparing + * the curvature estimates computed above and the curvatures calculated from the discrete differential operators. + */ + + gts_vertex_principal_curvatures(0.5 * normKh, Kg, &K1, &K2); + + err_e1 = err_e2 = 0.0; + /* loop through the values previously saved */ + for (e = 0; e < edge_count; e++) { + real weight, kappa, d1, d2; + real temp1, temp2; + real delta; + + weight = weights[e]; + kappa = kappas[e]; + d1 = d1s[e]; + d2 = d2s[e]; + + temp1 = fabs(eig[0] * d1 + eig[1] * d2); + temp1 = temp1 * temp1; + temp2 = fabs(eig[1] * d1 - eig[0] * d2); + temp2 = temp2 * temp2; + + /* err_e1 is for K1 associated with e1 */ + delta = K1 * temp1 + K2 * temp2 - kappa; + err_e1 += weight * delta * delta; + + /* err_e2 is for K1 associated with e2 */ + delta = K2 * temp1 + K1 * temp2 - kappa; + err_e2 += weight * delta * delta; + } + free(weights); + free(kappas); + free(d1s); + free(d2s); + + /* rotate eig by a right angle if that would decrease the error */ + if (err_e2 < err_e1) { + real temp = eig[0]; + + eig[0] = eig[1]; + eig[1] = -temp; + } + + e1[0] = eig[0] * basis1[0] + eig[1] * basis2[0]; + e1[1] = eig[0] * basis1[1] + eig[1] * basis2[1]; + e1[2] = eig[0] * basis1[2] + eig[1] * basis2[2]; + e1.normalize(); + + /* make N,e1,e2 a right handed coordinate sytem */ + e2 = N ^ e1; + e2.normalize(); } namespace OGF { @@ -529,107 +534,108 @@ namespace OGF { #if 0 inline static real angle(WOEdge *h) { - const Vec3r& n1 = h->GetbFace()->GetNormal(); - const Vec3r& n2 = h->GetaFace()->GetNormal(); - const Vec3r v = h->GetVec(); - real sine = (n1 ^ n2) * v / v.norm(); - if (sine >= 1.0) { - return M_PI / 2.0; - } - if (sine <= -1.0) { - return -M_PI / 2.0; - } - return ::asin(sine); + const Vec3r& n1 = h->GetbFace()->GetNormal(); + const Vec3r& n2 = h->GetaFace()->GetNormal(); + const Vec3r v = h->GetVec(); + real sine = (n1 ^ n2) * v / v.norm(); + if (sine >= 1.0) { + return M_PI / 2.0; + } + if (sine <= -1.0) { + return -M_PI / 2.0; + } + return ::asin(sine); } #endif // precondition1: P is inside the sphere // precondition2: P,V points to the outside of the sphere (i.e. OP.V > 0) -static bool sphere_clip_vector(const Vec3r& O, real r, const Vec3r& P, Vec3r& V) +static bool sphere_clip_vector(const Vec3r &O, real r, const Vec3r &P, Vec3r &V) { - Vec3r W = P - O; - real a = V.squareNorm(); - real b = 2.0 * V * W; - real c = W.squareNorm() - r * r; - real delta = b * b - 4 * a * c; - if (delta < 0) { - // Should not happen, but happens sometimes (numerical precision) - return true; - } - real t = - b + ::sqrt(delta) / (2.0 * a); - if (t < 0.0) { - // Should not happen, but happens sometimes (numerical precision) - return true; - } - if (t >= 1.0) { - // Inside the sphere - return false; - } - - V[0] = (t * V.x()); - V[1] = (t * V.y()); - V[2] = (t * V.z()); - - return true; + Vec3r W = P - O; + real a = V.squareNorm(); + real b = 2.0 * V * W; + real c = W.squareNorm() - r * r; + real delta = b * b - 4 * a * c; + if (delta < 0) { + // Should not happen, but happens sometimes (numerical precision) + return true; + } + real t = -b + ::sqrt(delta) / (2.0 * a); + if (t < 0.0) { + // Should not happen, but happens sometimes (numerical precision) + return true; + } + if (t >= 1.0) { + // Inside the sphere + return false; + } + + V[0] = (t * V.x()); + V[1] = (t * V.y()); + V[2] = (t * V.z()); + + return true; } // TODO: check optimizations: // use marking ? (measure *timings* ...) -void compute_curvature_tensor(WVertex *start, real radius, NormalCycle& nc) +void compute_curvature_tensor(WVertex *start, real radius, NormalCycle &nc) { - // in case we have a non-manifold vertex, skip it... - if (start->isBoundary()) - return; - - std::set<WVertex*> vertices; - const Vec3r& O = start->GetVertex(); - std::stack<WVertex*> S; - S.push(start); - vertices.insert(start); - while (!S.empty()) { - WVertex *v = S.top(); - S.pop(); - if (v->isBoundary()) - continue; - const Vec3r& P = v->GetVertex(); - WVertex::incoming_edge_iterator woeit = v->incoming_edges_begin(); - WVertex::incoming_edge_iterator woeitend = v->incoming_edges_end(); - for (; woeit != woeitend; ++woeit) { - WOEdge *h = *woeit; - if ((v == start) || h->GetVec() * (O - P) > 0.0) { - Vec3r V(-1 * h->GetVec()); - bool isect = sphere_clip_vector(O, radius, P, V); - assert (h->GetOwner()->GetNumberOfOEdges() == 2); // Because otherwise v->isBoundary() would be true - nc.accumulate_dihedral_angle(V, h->GetAngle()); - - if (!isect) { - WVertex *w = h->GetaVertex(); - if (vertices.find(w) == vertices.end()) { - vertices.insert(w); - S.push(w); - } - } - } - } - } + // in case we have a non-manifold vertex, skip it... + if (start->isBoundary()) + return; + + std::set<WVertex *> vertices; + const Vec3r &O = start->GetVertex(); + std::stack<WVertex *> S; + S.push(start); + vertices.insert(start); + while (!S.empty()) { + WVertex *v = S.top(); + S.pop(); + if (v->isBoundary()) + continue; + const Vec3r &P = v->GetVertex(); + WVertex::incoming_edge_iterator woeit = v->incoming_edges_begin(); + WVertex::incoming_edge_iterator woeitend = v->incoming_edges_end(); + for (; woeit != woeitend; ++woeit) { + WOEdge *h = *woeit; + if ((v == start) || h->GetVec() * (O - P) > 0.0) { + Vec3r V(-1 * h->GetVec()); + bool isect = sphere_clip_vector(O, radius, P, V); + assert(h->GetOwner()->GetNumberOfOEdges() == + 2); // Because otherwise v->isBoundary() would be true + nc.accumulate_dihedral_angle(V, h->GetAngle()); + + if (!isect) { + WVertex *w = h->GetaVertex(); + if (vertices.find(w) == vertices.end()) { + vertices.insert(w); + S.push(w); + } + } + } + } + } } -void compute_curvature_tensor_one_ring(WVertex *start, NormalCycle& nc) +void compute_curvature_tensor_one_ring(WVertex *start, NormalCycle &nc) { - // in case we have a non-manifold vertex, skip it... - if (start->isBoundary()) - return; - - WVertex::incoming_edge_iterator woeit = start->incoming_edges_begin(); - WVertex::incoming_edge_iterator woeitend = start->incoming_edges_end(); - for (; woeit != woeitend; ++woeit) { - WOEdge *h = (*woeit)->twin(); - nc.accumulate_dihedral_angle(h->GetVec(), h->GetAngle()); - WOEdge *hprev = h->getPrevOnFace(); - nc.accumulate_dihedral_angle(hprev->GetVec(), hprev->GetAngle()); - } + // in case we have a non-manifold vertex, skip it... + if (start->isBoundary()) + return; + + WVertex::incoming_edge_iterator woeit = start->incoming_edges_begin(); + WVertex::incoming_edge_iterator woeitend = start->incoming_edges_end(); + for (; woeit != woeitend; ++woeit) { + WOEdge *h = (*woeit)->twin(); + nc.accumulate_dihedral_angle(h->GetVec(), h->GetAngle()); + WOEdge *hprev = h->getPrevOnFace(); + nc.accumulate_dihedral_angle(hprev->GetVec(), hprev->GetAngle()); + } } -} // OGF namespace +} // namespace OGF } /* namespace Freestyle */ |