soc-2008-mxcurioni: merged changes to revision 14747, cosmetic changes for source/blender/freestyle

1 files changed, 647 insertions, 0 deletions

diff --git a/source/blender/freestyle/intern/winged_edge/Curvature.cpp b/source/blender/freestyle/intern/winged_edge/Curvature.cpp
new file mode 100755
index 00000000000..a890fb92c04
--- /dev/null
+++ b/source/blender/freestyle/intern/winged_edge/Curvature.cpp
@@ -0,0 +1,647 @@
+/* GTS - Library for the manipulation of triangulated surfaces
+ * Copyright (C) 1999-2002 Ray Jones, Stéphane Popinet
+ *
+ * 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 library 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.
+ *
+ * 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.
+ */
+
+#include "Curvature.h"
+#include <math.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)
+{
+  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)
+{
+  bool b=false;
+  for (int i=0; i<3; i++)
+    b = b ||
+      ((f->GetEdgeList()[i]->getVec3r() * f->GetEdgeList()[(i+1)%3]->getVec3r()) < 0);
+  return b;
+}
+
+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());
+
+  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);
+}
+
+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;
+
+  Vec3r u (v1->GetVertex()-vo->GetVertex());
+  Vec3r v(v2->GetVertex()-vo->GetVertex());
+
+  udotv = u * v;
+  denom = sqrt(u.squareNorm() * v.squareNorm() - udotv * udotv);
+
+  /* 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)));
+}
+
+/** 
+ * 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.
+ *
+ * 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
+ *
+ * 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)
+{
+	real area = 0.0;
+
+	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 ";
+
+			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.
+ *
+ * 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
+ *
+ * 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)
+{
+	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:
+ * @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.  
+ *
+ * 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().
+ */
+void gts_vertex_principal_curvatures (real Kh, real Kg, 
+				      real * K1, real * K2)
+{
+  real temp = Kh*Kh - Kg;
+
+  if (!K1) return; 
+  if (!K1) return; 
+
+  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;
+
+  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;
+}
+
+/** 
+ * 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. 
+ *
+ * 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)
+{
+  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->getVec3r());
+
+    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 {
+    static real angle(WOEdge* h) {
+        Vec3r e(h->GetbVertex()->GetVertex()-h->GetaVertex()->GetVertex()) ;
+        Vec3r n1 = h->GetbFace()->GetNormal();
+        Vec3r n2 = h->GetaFace()->GetNormal();
+        real sine = (n1 ^ n2) * e / e.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;
+              Vec3r V(h->GetaVertex()->GetVertex()-h->GetbVertex()->GetVertex()) ;
+              if((v == start) || V * (P - O) > 0.0) {
+                    bool isect = sphere_clip_vector(O, radius, P, V) ;
+                    if(h->GetOwner()->GetNumberOfOEdges() == 2) {
+                        real beta = angle(h) ;
+                        nc.accumulate_dihedral_angle(V, beta) ;
+                    }
+                    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();
+          Vec3r hvec(h->GetbVertex()->GetVertex()-h->GetaVertex()->GetVertex());
+          nc.accumulate_dihedral_angle(hvec, angle(h)) ;
+          WOEdge *hprev = h->getPrevOnFace();
+          Vec3r hprevvec(hprev->GetbVertex()->GetVertex()-hprev->GetaVertex()->GetVertex());
+          nc.accumulate_dihedral_angle(hprevvec, angle(hprev)) ;
+        }
+    }
+
+}