diff options
author | Tamito Kajiyama <rd6t-kjym@asahi-net.or.jp> | 2013-01-02 05:55:30 +0400 |
---|---|---|
committer | Tamito Kajiyama <rd6t-kjym@asahi-net.or.jp> | 2013-01-02 05:55:30 +0400 |
commit | 699da2fb0d9012cef5e45cc1b547a01fd92dbc1c (patch) | |
tree | eb6e3de1f90387af1a2628bbc79d51327e3ac248 /source/blender/freestyle/intern/winged_edge/Curvature.cpp | |
parent | 520ab93465d8056bf6d4a4115743daaee4c2fd66 (diff) |
Yet another big style clean-up patch by Bastien Montagne, thanks a lot!
Now the code style is acceptable for the merge now, according to Bastien.
Thanks again Bastien for having this done! :)
Diffstat (limited to 'source/blender/freestyle/intern/winged_edge/Curvature.cpp')
-rw-r--r-- | source/blender/freestyle/intern/winged_edge/Curvature.cpp | 1112 |
1 files changed, 554 insertions, 558 deletions
diff --git a/source/blender/freestyle/intern/winged_edge/Curvature.cpp b/source/blender/freestyle/intern/winged_edge/Curvature.cpp index ce8fc5b98fa..acefe1aa5fc 100644 --- a/source/blender/freestyle/intern/winged_edge/Curvature.cpp +++ b/source/blender/freestyle/intern/winged_edge/Curvature.cpp @@ -1,646 +1,642 @@ -/* GTS - Library for the manipulation of triangulated surfaces - * Copyright (C) 1999-2002 Ray Jones, Stéphane Popinet +/* + * ***** BEGIN GPL LICENSE BLOCK ***** * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Library General Public - * License as published by the Free Software Foundation; either - * version 2 of the License, or (at your option) any later version. + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. * - * This library is distributed in the hope that it will be useful, + * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Library General Public License for more details. + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. * - * You should have received a copy of the GNU Library General Public - * License along with this library; if not, write to the - * Free Software Foundation, Inc., 59 Temple Place - Suite 330, - * Boston, MA 02111-1307, USA. + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + * + * This Code is Copyright (C) 2010 Blender Foundation. + * All rights reserved. + * + * The Original Code is: + * GTS - Library for the manipulation of triangulated surfaces + * Copyright (C) 1999 Stephane Popinet + * and: + * OGF/Graphite: Geometry and Graphics Programming Library + Utilities + * Copyright (C) 2000-2003 Bruno Levy + * Contact: Bruno Levy levy@loria.fr + * ISA Project + * LORIA, INRIA Lorraine, + * Campus Scientifique, BP 239 + * 54506 VANDOEUVRE LES NANCY CEDEX + * FRANCE + * + * Contributor(s): none yet. + * + * ***** END GPL LICENSE BLOCK ***** */ +/** \file blender/freestyle/intern/winged_edge/Curvature.cpp + * \ingroup freestyle + * \brief GTS - Library for the manipulation of triangulated surfaces + * \author Stephane Popinet + * \date 1999 + * \brief OGF/Graphite: Geometry and Graphics Programming Library + Utilities + * \author Bruno Levy + * \date 2000-2003 + */ + +#include <assert.h> #include <cstdlib> // for malloc and free -#include "Curvature.h" #include <math.h> -#include <assert.h> -#include "WEdge.h" -#include "../system/FreestyleConfig.h" -#include "../geometry/normal_cycle.h" #include <set> #include <stack> -static bool angle_obtuse (WVertex * v, WFace * f) +#include "Curvature.h" +#include "WEdge.h" + +#include "../geometry/normal_cycle.h" + +#include "../system/FreestyleConfig.h" + +static bool angle_obtuse(WVertex *v, WFace *f) { - WOEdge * e; - f->getOppositeEdge (v, e); - - Vec3r vec1(e->GetaVertex()->GetVertex()-v->GetVertex()); - Vec3r vec2(e->GetbVertex()->GetVertex()-v->GetVertex()); - return ((vec1 * vec2) < 0); + WOEdge *e; + f->getOppositeEdge(v, e); + + 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) +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) +static real cotan(WVertex *vo, WVertex *v1, WVertex *v2) { - /* cf. Appendix B of [Meyer et al 2002] */ - real udotv, denom; - - Vec3r u(v1->GetVertex()- vo->GetVertex()); - Vec3r v(v2->GetVertex()- vo->GetVertex()); + /* cf. Appendix B of [Meyer et al 2002] */ + real udotv, denom; - udotv = u * v; - denom = sqrt(u.squareNorm() * v.squareNorm() - udotv * udotv); + Vec3r u(v1->GetVertex() - vo->GetVertex()); + Vec3r v(v2->GetVertex() - vo->GetVertex()); - /* 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); + udotv = u * v; + denom = sqrt(u.squareNorm() * v.squareNorm() - udotv * udotv); - 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) +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: - * @v: a #WVertex. - * @s: a #GtsSurface. - * @Kh: the Mean Curvature Normal at @v. +/*! gts_vertex_mean_curvature_normal: + * @v: a #WVertex. + * @s: a #GtsSurface. + * @Kh: the Mean Curvature Normal at @v. * - * Computes the Discrete Mean Curvature Normal approximation at @v. - * The mean curvature at @v is half the magnitude of the vector @Kh. + * Computes the Discrete Mean Curvature Normal approximation at @v. + * The mean curvature at @v is half the magnitude of the vector @Kh. * - * Note: the normal computed is not unit length, and may point either - * into or out of the surface, depending on the curvature at @v. It - * is the responsibility of the caller of the function to use the mean - * curvature normal appropriately. + * Note: the normal computed is not unit length, and may point either into or out of the surface, depending on + * the curvature at @v. It is the responsibility of the caller of the function to use the mean curvature normal + * appropriately. * - * This approximation is from the paper: - * Discrete Differential-Geometry Operators for Triangulated 2-Manifolds - * Mark Meyer, Mathieu Desbrun, Peter Schroder, Alan H. Barr - * VisMath '02, Berlin (Germany) - * http://www-grail.usc.edu/pubs.html + * This approximation is from the paper: + * Discrete Differential-Geometry Operators for Triangulated 2-Manifolds + * Mark Meyer, Mathieu Desbrun, Peter Schroder, Alan H. Barr + * VisMath '02, Berlin (Germany) + * http://www-grail.usc.edu/pubs.html * - * Returns: %TRUE if the operator could be evaluated, %FALSE if the - * evaluation failed for some reason (@v is boundary or is the - * endpoint of a non-manifold edge.) + * Returns: %TRUE if the operator could be evaluated, %FALSE if the evaluation failed for some reason (@v is + * boundary or is the endpoint of a non-manifold edge.) */ -bool gts_vertex_mean_curvature_normal (WVertex * v, Vec3r &Kh) +bool gts_vertex_mean_curvature_normal(WVertex *v, Vec3r &Kh) { 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; - - WVertex::incoming_edge_iterator itE; - - for (itE=v->incoming_edges_begin(); - itE!=v->incoming_edges_end(); itE++) - area+=(*itE)->GetaFace()->getArea(); - - Kh=Vec3r(0.0, 0.0, 0.0); - - for (itE=v->incoming_edges_begin(); - itE!=v->incoming_edges_end(); itE++) - { - WOEdge * e = (*itE)->getPrevOnFace(); - //if ((e->GetaVertex()==v) || (e->GetbVertex()==v)) cerr<< "BUG "; + if (v->isBoundary()) + return false; - WVertex * v1 = e->GetaVertex(); - WVertex * v2 = e->GetbVertex(); - real temp; - - temp = cotan (v1, v, v2); - Kh = Vec3r(Kh+temp*(v2->GetVertex()-v->GetVertex())); + WVertex::incoming_edge_iterator itE; - 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; + for (itE=v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++) + area += (*itE)->GetaFace()->getArea(); + + Kh = Vec3r(0.0, 0.0, 0.0); + + 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 "; +#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; } -/** - * gts_vertex_gaussian_curvature: - * @v: a #WVertex. - * @s: a #GtsSurface. - * @Kg: the Discrete Gaussian Curvature approximation at @v. +/*! gts_vertex_gaussian_curvature: + * @v: a #WVertex. + * @s: a #GtsSurface. + * @Kg: the Discrete Gaussian Curvature approximation at @v. * - * Computes the Discrete Gaussian Curvature approximation at @v. + * Computes the Discrete Gaussian Curvature approximation at @v. * - * This approximation is from the paper: - * Discrete Differential-Geometry Operators for Triangulated 2-Manifolds - * Mark Meyer, Mathieu Desbrun, Peter Schroder, Alan H. Barr - * VisMath '02, Berlin (Germany) - * http://www-grail.usc.edu/pubs.html + * This approximation is from the paper: + * Discrete Differential-Geometry Operators for Triangulated 2-Manifolds + * Mark Meyer, Mathieu Desbrun, Peter Schroder, Alan H. Barr + * VisMath '02, Berlin (Germany) + * http://www-grail.usc.edu/pubs.html * - * Returns: %TRUE if the operator could be evaluated, %FALSE if the - * evaluation failed for some reason (@v is boundary or is the - * endpoint of a non-manifold edge.) + * Returns: %TRUE if the operator could be evaluated, %FALSE if the evaluation failed for some reason (@v is + * boundary or is the endpoint of a non-manifold edge.) */ -bool gts_vertex_gaussian_curvature (WVertex * v, real * Kg) +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; + if (!v) + return false; + if (!Kg) + return false; + + /* this operator is not defined for boundary edges */ + if (v->isBoundary()) { + *Kg = 0.0; + 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; + 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: - * @Kh: mean curvature. - * @Kg: Gaussian curvature. - * @K1: first principal curvature. - * @K2: second principal curvature. +/*! gts_vertex_principal_curvatures: + * @Kh: mean curvature. + * @Kg: Gaussian curvature. + * @K1: first principal curvature. + * @K2: second principal curvature. * - * Computes the principal curvatures at a point given the mean and - * Gaussian curvatures at that point. + * Computes the principal curvatures at a point given the mean and Gaussian curvatures at that point. * - * The mean curvature can be computed as one-half the magnitude of the - * vector computed by gts_vertex_mean_curvature_normal(). + * The mean curvature can be computed as one-half the magnitude of the vector computed by + * gts_vertex_mean_curvature_normal(). * - * The Gaussian curvature can be computed with - * gts_vertex_gaussian_curvature(). + * 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) return; - if (!K1) 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) +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) +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: - * @v: a #WVertex. - * @s: a #GtsSurface. - * @Kh: mean curvature normal (a #Vec3r). - * @Kg: Gaussian curvature (a real). - * @e1: first principal curvature direction (direction of largest curvature). - * @e2: second principal curvature direction. +/*! gts_vertex_principal_directions: + * @v: a #WVertex. + * @s: a #GtsSurface. + * @Kh: mean curvature normal (a #Vec3r). + * @Kg: Gaussian curvature (a real). + * @e1: first principal curvature direction (direction of largest curvature). + * @e2: second principal curvature direction. * - * Computes the principal curvature directions at a point given @Kh - * and @Kg, the mean curvature normal and Gaussian curvatures at that - * point, computed with gts_vertex_mean_curvature_normal() and - * gts_vertex_gaussian_curvature(), respectively. + * Computes the principal curvature directions at a point given @Kh and @Kg, the mean curvature normal and + * Gaussian curvatures at that point, computed with gts_vertex_mean_curvature_normal() and + * gts_vertex_gaussian_curvature(), respectively. * - * Note that this computation is very approximate and tends to be - * unstable. Smoothing of the surface or the principal directions may - * be necessary to achieve reasonable results. + * Note that this computation is very approximate and tends to be unstable. Smoothing of the surface or the principal + * directions may be necessary to achieve reasonable results. */ -void gts_vertex_principal_directions (WVertex * v, - Vec3r Kh, real Kg, - Vec3r &e1, Vec3r &e2) +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 { - inline static real angle(WOEdge * h) { - const Vec3r& n1 = h->GetbFace()->GetNormal(); - const Vec3r& n2 = h->GetaFace()->GetNormal(); - const Vec3r v = h->getVec3r(); - 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) ; - } - - // 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 - ) { - - 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 - ) { - // 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 - ) { - // 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()) ; - } - } +inline static real angle(WOEdge *h) +{ + const Vec3r& n1 = h->GetbFace()->GetNormal(); + const Vec3r& n2 = h->GetaFace()->GetNormal(); + const Vec3r v = h->getVec3r(); + 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); +} + +// 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) +{ + 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) +{ + // 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) +{ + // 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 |