65ebdeac3a8cf957cef1c498c068762fe4685ddf
[blender.git] / source / blender / freestyle / intern / winged_edge / Curvature.cpp
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * This Code is Copyright (C) 2010 Blender Foundation.
19  * All rights reserved.
20  *
21  * The Original Code is:
22  *     GTS - Library for the manipulation of triangulated surfaces
23  *     Copyright (C) 1999 Stephane Popinet
24  * and:
25  *     OGF/Graphite: Geometry and Graphics Programming Library + Utilities
26  *     Copyright (C) 2000-2003 Bruno Levy
27  *     Contact: Bruno Levy levy@loria.fr
28  *         ISA Project
29  *         LORIA, INRIA Lorraine,
30  *         Campus Scientifique, BP 239
31  *         54506 VANDOEUVRE LES NANCY CEDEX
32  *         FRANCE
33  *
34  * Contributor(s): none yet.
35  *
36  * ***** END GPL LICENSE BLOCK *****
37  */
38
39 /** \file blender/freestyle/intern/winged_edge/Curvature.cpp
40  *  \ingroup freestyle
41  *  \brief GTS - Library for the manipulation of triangulated surfaces
42  *  \author Stephane Popinet
43  *  \date 1999
44  *  \brief OGF/Graphite: Geometry and Graphics Programming Library + Utilities
45  *  \author Bruno Levy
46  *  \date 2000-2003
47  */
48
49 #include <assert.h>
50 #include <cstdlib> // for malloc and free
51 #include <math.h>
52 #include <set>
53 #include <stack>
54
55 #include "Curvature.h"
56 #include "WEdge.h"
57
58 #include "../geometry/normal_cycle.h"
59
60 #include "../system/FreestyleConfig.h"
61
62 static bool angle_obtuse(WVertex *v, WFace *f)
63 {
64         WOEdge *e;
65         f->getOppositeEdge(v, e);
66
67         Vec3r vec1(e->GetaVertex()->GetVertex()-v->GetVertex());
68         Vec3r vec2(e->GetbVertex()->GetVertex()-v->GetVertex());
69         return ((vec1 * vec2) < 0);
70 }
71
72 // FIXME
73 // WVvertex is useless but kept for history reasons
74 static bool triangle_obtuse(WVertex *, WFace *f)
75 {
76         bool b = false;
77         for (int i = 0; i < 3; i++)
78                 b = b || ((f->getEdgeList()[i]->GetVec() * f->getEdgeList()[(i + 1) % 3]->GetVec()) < 0);
79         return b;
80 }
81
82 static real cotan(WVertex *vo, WVertex *v1, WVertex *v2)
83 {
84         /* cf. Appendix B of [Meyer et al 2002] */
85         real udotv, denom;
86
87         Vec3r u(v1->GetVertex() - vo->GetVertex());
88         Vec3r v(v2->GetVertex() - vo->GetVertex());
89
90         udotv = u * v;
91         denom = sqrt(u.squareNorm() * v.squareNorm() - udotv * udotv);
92
93         /* denom can be zero if u==v.  Returning 0 is acceptable, based on the callers of this function below. */
94         if (denom == 0.0)
95                 return 0.0;
96         return (udotv / denom);
97 }
98
99 static real angle_from_cotan(WVertex *vo, WVertex *v1, WVertex *v2)
100 {
101         /* cf. Appendix B and the caption of Table 1 from [Meyer et al 2002] */
102         real udotv, denom;
103
104         Vec3r u (v1->GetVertex() - vo->GetVertex());
105         Vec3r v(v2->GetVertex() - vo->GetVertex());
106
107         udotv = u * v;
108         denom = sqrt(u.squareNorm() * v.squareNorm() - udotv * udotv);
109
110         /* Note: I assume this is what they mean by using atan2(). -Ray Jones */
111
112         /* tan = denom/udotv = y/x (see man page for atan2) */
113         return (fabs(atan2(denom, udotv)));
114 }
115
116 /*! gts_vertex_mean_curvature_normal:
117  *  @v: a #WVertex.
118  *  @s: a #GtsSurface.
119  *  @Kh: the Mean Curvature Normal at @v.
120  *
121  *  Computes the Discrete Mean Curvature Normal approximation at @v.
122  *  The mean curvature at @v is half the magnitude of the vector @Kh.
123  *
124  *  Note: the normal computed is not unit length, and may point either into or out of the surface, depending on
125  *  the curvature at @v. It is the responsibility of the caller of the function to use the mean curvature normal
126  *  appropriately.
127  *
128  *  This approximation is from the paper:
129  *      Discrete Differential-Geometry Operators for Triangulated 2-Manifolds
130  *      Mark Meyer, Mathieu Desbrun, Peter Schroder, Alan H. Barr
131  *      VisMath '02, Berlin (Germany)
132  *      http://www-grail.usc.edu/pubs.html
133  *
134  *  Returns: %TRUE if the operator could be evaluated, %FALSE if the evaluation failed for some reason (@v is
135  *  boundary or is the endpoint of a non-manifold edge.)
136  */
137 bool gts_vertex_mean_curvature_normal(WVertex *v, Vec3r &Kh)
138 {
139         real area = 0.0;
140
141         if (!v)
142                 return false;
143
144         /* this operator is not defined for boundary edges */
145         if (v->isBoundary())
146                 return false;
147
148         WVertex::incoming_edge_iterator itE;
149
150         for (itE=v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++)
151                 area += (*itE)->GetaFace()->getArea();
152
153         Kh = Vec3r(0.0, 0.0, 0.0);
154
155         for (itE=v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++) {
156                 WOEdge *e = (*itE)->getPrevOnFace();
157 #if 0
158                 if ((e->GetaVertex() == v) || (e->GetbVertex() == v))
159                         cerr<< "BUG ";
160 #endif
161                 WVertex *v1 = e->GetaVertex();
162                 WVertex *v2 = e->GetbVertex();
163                 real temp;
164
165                 temp = cotan(v1, v, v2);
166                 Kh = Vec3r(Kh + temp * (v2->GetVertex() - v->GetVertex()));
167
168                 temp = cotan(v2, v, v1);
169                 Kh = Vec3r(Kh + temp * (v1->GetVertex() - v->GetVertex()));
170         }
171         if (area > 0.0) {
172                 Kh[0] /= 2 * area;
173                 Kh[1] /= 2 * area;
174                 Kh[2] /= 2 * area;
175         }
176         else {
177                 return false;
178         }
179
180         return true;
181 }
182
183 /*! gts_vertex_gaussian_curvature:
184  *  @v: a #WVertex.
185  *  @s: a #GtsSurface.
186  *  @Kg: the Discrete Gaussian Curvature approximation at @v.
187  *
188  *  Computes the Discrete Gaussian Curvature approximation at @v.
189  *
190  *  This approximation is from the paper:
191  *      Discrete Differential-Geometry Operators for Triangulated 2-Manifolds
192  *      Mark Meyer, Mathieu Desbrun, Peter Schroder, Alan H. Barr
193  *      VisMath '02, Berlin (Germany)
194  *      http://www-grail.usc.edu/pubs.html
195  *
196  *  Returns: %TRUE if the operator could be evaluated, %FALSE if the evaluation failed for some reason (@v is
197  * boundary or is the endpoint of a non-manifold edge.)
198  */
199 bool gts_vertex_gaussian_curvature(WVertex *v, real *Kg)
200 {
201         real area = 0.0;
202         real angle_sum = 0.0;
203
204         if (!v)
205                 return false;
206         if (!Kg)
207                 return false;
208
209         /* this operator is not defined for boundary edges */
210         if (v->isBoundary()) {
211                 *Kg = 0.0;
212                 return false;
213         }
214
215         WVertex::incoming_edge_iterator itE;
216         for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++)
217                  area += (*itE)->GetaFace()->getArea();
218
219         for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++) {
220                 WOEdge *e = (*itE)->getPrevOnFace();
221                 WVertex *v1 = e->GetaVertex();
222                 WVertex *v2 = e->GetbVertex();
223                 angle_sum += angle_from_cotan(v, v1, v2);
224         }
225
226         *Kg = (2.0 * M_PI - angle_sum) / area;
227
228         return true;
229 }
230
231 /*! gts_vertex_principal_curvatures:
232  *  @Kh: mean curvature.
233  *  @Kg: Gaussian curvature.
234  *  @K1: first principal curvature.
235  *  @K2: second principal curvature.
236  *
237  *  Computes the principal curvatures at a point given the mean and Gaussian curvatures at that point.
238  *
239  *  The mean curvature can be computed as one-half the magnitude of the vector computed by
240  *  gts_vertex_mean_curvature_normal().
241  *
242  *  The Gaussian curvature can be computed with gts_vertex_gaussian_curvature().
243  */
244 void gts_vertex_principal_curvatures (real Kh, real Kg, real *K1, real *K2)
245 {
246         real temp = Kh * Kh - Kg;
247
248         if (!K1 || !K2)
249                 return;
250
251         if (temp < 0.0)
252                 temp = 0.0;
253         temp = sqrt (temp);
254         *K1 = Kh + temp;
255         *K2 = Kh - temp;
256 }
257
258 /* from Maple */
259 static void linsolve(real m11, real m12, real b1, real m21, real m22, real b2, real *x1, real *x2)
260 {
261         real temp;
262
263         temp = 1.0 / (m21 * m12 - m11 * m22);
264         *x1 = (m12 * b2 - m22 * b1) * temp;
265         *x2 = (m11 * b2 - m21 * b1) * temp;
266 }
267
268 /* from Maple - largest eigenvector of [a b; b c] */
269 static void eigenvector(real a, real b, real c, Vec3r e)
270 {
271         if (b == 0.0) {
272                 e[0] = 0.0;
273         }
274         else {
275                 e[0] = -(c - a - sqrt(c * c - 2 * a * c + a * a + 4 * b * b)) / (2 * b);
276         }
277         e[1] = 1.0;
278         e[2] = 0.0;
279 }
280
281 /*! gts_vertex_principal_directions:
282  *  @v: a #WVertex.
283  *  @s: a #GtsSurface.
284  *  @Kh: mean curvature normal (a #Vec3r).
285  *  @Kg: Gaussian curvature (a real).
286  *  @e1: first principal curvature direction (direction of largest curvature).
287  *  @e2: second principal curvature direction.
288  *
289  *  Computes the principal curvature directions at a point given @Kh and @Kg, the mean curvature normal and
290  *  Gaussian curvatures at that point, computed with gts_vertex_mean_curvature_normal() and
291  *  gts_vertex_gaussian_curvature(), respectively.
292  *
293  *  Note that this computation is very approximate and tends to be unstable. Smoothing of the surface or the principal
294  *  directions may be necessary to achieve reasonable results.
295  */
296 void gts_vertex_principal_directions(WVertex *v, Vec3r Kh, real Kg, Vec3r &e1, Vec3r &e2)
297 {
298         Vec3r N;
299         real normKh;
300
301         Vec3r basis1, basis2, d, eig;
302         real ve2, vdotN;
303         real aterm_da, bterm_da, cterm_da, const_da;
304         real aterm_db, bterm_db, cterm_db, const_db;
305         real a, b, c;
306         real K1, K2;
307         real *weights, *kappas, *d1s, *d2s;
308         int edge_count;
309         real err_e1, err_e2;
310         int e;
311         WVertex::incoming_edge_iterator itE;
312
313         /* compute unit normal */
314         normKh = Kh.norm();
315
316         if (normKh > 0.0) {
317                 Kh.normalize();
318         }
319         else {
320                 /* This vertex is a point of zero mean curvature (flat or saddle point). Compute a normal by averaging
321                  * the adjacent triangles
322                  */
323                 N[0] = N[1] = N[2] = 0.0;
324
325                 for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++)
326                         N = Vec3r(N + (*itE)->GetaFace()->GetNormal());
327                 real normN = N.norm();
328                 if (normN <= 0.0)
329                         return;
330                 N.normalize();
331         }
332
333         /* construct a basis from N: */
334         /* set basis1 to any component not the largest of N */
335         basis1[0] =  basis1[1] =  basis1[2] = 0.0;
336         if (fabs (N[0]) > fabs (N[1]))
337                 basis1[1] = 1.0;
338         else
339                 basis1[0] = 1.0;
340
341         /* make basis2 orthogonal to N */
342         basis2 = (N ^ basis1);
343         basis2.normalize();
344
345         /* make basis1 orthogonal to N and basis2 */
346         basis1 = (N ^ basis2);
347         basis1.normalize();
348
349         aterm_da = bterm_da = cterm_da = const_da = 0.0;
350         aterm_db = bterm_db = cterm_db = const_db = 0.0;
351         int nb_edges=v->GetEdges().size();
352
353         weights = (real *)malloc(sizeof (real) * nb_edges);
354         kappas = (real *)malloc(sizeof (real) * nb_edges);
355         d1s = (real *)malloc(sizeof (real) * nb_edges);
356         d2s = (real *)malloc(sizeof (real) * nb_edges);
357         edge_count = 0;
358
359         for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++) {
360                 WOEdge *e;
361                 WFace *f1, *f2;
362                 real weight, kappa, d1, d2;
363                 Vec3r vec_edge;
364                 if (!*itE)
365                         continue;
366                 e = *itE;
367
368                 /* since this vertex passed the tests in gts_vertex_mean_curvature_normal(), this should be true. */
369                 //g_assert(gts_edge_face_number (e, s) == 2);
370
371                 /* identify the two triangles bordering e in s */
372                 f1 = e->GetaFace();
373                 f2 = e->GetbFace();
374
375                 /* We are solving for the values of the curvature tensor
376                 *     B = [ a b ; b c ].
377                 *  The computations here are from section 5 of [Meyer et al 2002].
378                 *
379                 *  The first step is to calculate the linear equations governing the values of (a,b,c). These can be computed
380                 *  by setting the derivatives of the error E to zero (section 5.3).
381                 *
382                 *  Since a + c = norm(Kh), we only compute the linear equations for dE/da and dE/db. (NB: [Meyer et al 2002]
383                 *  has the equation a + b = norm(Kh), but I'm almost positive this is incorrect).
384                 *
385                 *  Note that the w_ij (defined in section 5.2) are all scaled by (1/8*A_mixed). We drop this uniform scale
386                 *  factor because the solution of the linear equations doesn't rely on it.
387                 *
388                 *  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
389                 *  terms that are the constant factors in the equations.
390                 */
391
392                 /* find the vector from v along edge e */
393                 vec_edge = Vec3r(-1 * e->GetVec());
394
395                 ve2 = vec_edge.squareNorm();
396                 vdotN = vec_edge * N;
397
398                 /* section 5.2 - There is a typo in the computation of kappa. The edges should be x_j-x_i. */
399                 kappa = 2.0 * vdotN / ve2;
400
401                 /* section 5.2 */
402
403                 /* I don't like performing a minimization where some of the weights can be negative (as can be the case
404                 *  if f1 or f2 are obtuse). To ensure all-positive weights, we check for obtuseness. */
405                 weight = 0.0;
406                 if (!triangle_obtuse(v, f1)) {
407                         weight += ve2 * cotan(f1->GetNextOEdge(e->twin())->GetbVertex(), e->GetaVertex(), e->GetbVertex()) / 8.0;
408                 }
409                 else {
410                         if (angle_obtuse(v, f1)) {
411                                 weight += ve2 * f1->getArea() / 4.0;
412                         }
413                         else {
414                                 weight += ve2 * f1->getArea() / 8.0;
415                         }
416                 }
417
418                 if (!triangle_obtuse(v, f2)) {
419                         weight += ve2 * cotan (f2->GetNextOEdge(e)->GetbVertex(), e->GetaVertex(), e->GetbVertex()) / 8.0;
420                 }
421                 else {
422                         if (angle_obtuse(v, f2)) {
423                                 weight += ve2 * f1->getArea() / 4.0;
424                         }
425                         else {
426                                 weight += ve2 * f1->getArea() / 8.0;
427                         }
428                 }
429
430                 /* projection of edge perpendicular to N (section 5.3) */
431                 d[0] = vec_edge[0] - vdotN * N[0];
432                 d[1] = vec_edge[1] - vdotN * N[1];
433                 d[2] = vec_edge[2] - vdotN * N[2];
434                 d.normalize();
435
436                 /* not explicit in the paper, but necessary. Move d to 2D basis. */
437                 d1 = d * basis1;
438                 d2 = d * basis2;
439
440                 /* store off the curvature, direction of edge, and weights for later use */
441                 weights[edge_count] = weight;
442                 kappas[edge_count] = kappa;
443                 d1s[edge_count] = d1;
444                 d2s[edge_count] = d2;
445                 edge_count++;
446
447                 /* Finally, update the linear equations */
448                 aterm_da += weight * d1 * d1 * d1 * d1;
449                 bterm_da += weight * d1 * d1 * 2 * d1 * d2;
450                 cterm_da += weight * d1 * d1 * d2 * d2;
451                 const_da += weight * d1 * d1 * (-kappa);
452
453                 aterm_db += weight * d1 * d2 * d1 * d1;
454                 bterm_db += weight * d1 * d2 * 2 * d1 * d2;
455                 cterm_db += weight * d1 * d2 * d2 * d2;
456                 const_db += weight * d1 * d2 * (-kappa);
457         }
458
459         /* now use the identity (Section 5.3) a + c = |Kh| = 2 * kappa_h */
460         aterm_da -= cterm_da;
461         const_da += cterm_da * normKh;
462
463         aterm_db -= cterm_db;
464         const_db += cterm_db * normKh;
465
466         /* check for solvability of the linear system */
467         if (((aterm_da * bterm_db - aterm_db * bterm_da) != 0.0) && ((const_da != 0.0) || (const_db != 0.0))) {
468                 linsolve(aterm_da, bterm_da, -const_da, aterm_db, bterm_db, -const_db, &a, &b);
469
470                 c = normKh - a;
471
472                 eigenvector(a, b, c, eig);
473         }
474         else {
475                 /* region of v is planar */
476                 eig[0] = 1.0;
477                 eig[1] = 0.0;
478         }
479
480         /* Although the eigenvectors of B are good estimates of the principal directions, it seems that which one is
481          * attached to which curvature direction is a bit arbitrary. This may be a bug in my implementation, or just
482          * a side-effect of the inaccuracy of B due to the discrete nature of the sampling.
483          *
484          * To overcome this behavior, we'll evaluate which assignment best matches the given eigenvectors by comparing
485          * the curvature estimates computed above and the curvatures calculated from the discrete differential operators.
486          */
487
488         gts_vertex_principal_curvatures(0.5 * normKh, Kg, &K1, &K2);
489
490         err_e1 = err_e2 = 0.0;
491         /* loop through the values previously saved */
492         for (e = 0; e < edge_count; e++) {
493                 real weight, kappa, d1, d2;
494                 real temp1, temp2;
495                 real delta;
496
497                 weight = weights[e];
498                 kappa = kappas[e];
499                 d1 = d1s[e];
500                 d2 = d2s[e];
501
502                 temp1 = fabs (eig[0] * d1 + eig[1] * d2);
503                 temp1 = temp1 * temp1;
504                 temp2 = fabs (eig[1] * d1 - eig[0] * d2);
505                 temp2 = temp2 * temp2;
506
507                 /* err_e1 is for K1 associated with e1 */
508                 delta = K1 * temp1 + K2 * temp2 - kappa;
509                 err_e1 += weight * delta * delta;
510
511                 /* err_e2 is for K1 associated with e2 */
512                 delta = K2 * temp1 + K1 * temp2 - kappa;
513                 err_e2 += weight * delta * delta;
514         }
515         free (weights);
516         free (kappas);
517         free (d1s);
518         free (d2s);
519
520         /* rotate eig by a right angle if that would decrease the error */
521         if (err_e2 < err_e1) {
522                 real temp = eig[0];
523
524                 eig[0] = eig[1];
525                 eig[1] = -temp;
526         }
527
528         e1[0] = eig[0] * basis1[0] + eig[1] * basis2[0];
529         e1[1] = eig[0] * basis1[1] + eig[1] * basis2[1];
530         e1[2] = eig[0] * basis1[2] + eig[1] * basis2[2];
531         e1.normalize();
532
533         /* make N,e1,e2 a right handed coordinate sytem */
534         e2 =  N ^ e1;
535         e2.normalize();
536 }
537
538 namespace OGF {
539
540 inline static real angle(WOEdge *h)
541 {
542         const Vec3r& n1 = h->GetbFace()->GetNormal();
543         const Vec3r& n2 = h->GetaFace()->GetNormal();
544         const Vec3r v = h->getVec3r();
545         real sine = (n1 ^ n2) * v / v.norm();
546         if (sine >= 1.0) {
547                 return M_PI / 2.0;
548         }
549         if (sine <= -1.0) {
550                 return -M_PI / 2.0;
551         }
552         return ::asin(sine);
553 }
554
555 // precondition1: P is inside the sphere
556 // precondition2: P,V points to the outside of the sphere (i.e. OP.V > 0)
557 static bool sphere_clip_vector(const Vec3r& O, real r, const Vec3r& P, Vec3r& V)
558 {
559         Vec3r W = P - O;
560         real a = V.squareNorm();
561         real b = 2.0 * V * W;
562         real c = W.squareNorm() - r * r;
563         real delta = b * b - 4 * a * c;
564         if (delta < 0) {
565                 // Should not happen, but happens sometimes (numerical precision)
566                 return true;
567         }
568         real t = - b + ::sqrt(delta) / (2.0 * a);
569         if (t < 0.0) {
570                 // Should not happen, but happens sometimes (numerical precision)
571                 return true;
572         }
573         if (t >= 1.0) {
574                 // Inside the sphere
575                 return false;
576         }
577
578         V[0] = (t * V.x());
579         V[1] = (t * V.y());
580         V[2] = (t * V.z());
581
582         return true;
583 }
584
585 // TODO: check optimizations:
586 // use marking ? (measure *timings* ...)
587 void compute_curvature_tensor(WVertex *start, real radius, NormalCycle& nc)
588 {
589         // TODO: for some reason, the WVertex 'start' may have no associated edges
590         // (i.e., WVertex::_EdgeList is empty), which causes a crash due to a call
591         // of WVertex::_EdgeList.front().
592         if (start->GetEdges().empty())
593                 return;
594
595         // in case we have a non-manifold vertex, skip it...
596         if (start->isBoundary())
597                 return;
598
599         std::set<WVertex*> vertices;
600         const Vec3r& O = start->GetVertex();
601         std::stack<WVertex*> S;
602         S.push(start);
603         vertices.insert(start);
604         while (!S.empty()) {
605                 WVertex *v = S.top();
606                 S.pop();
607                 if (v->isBoundary())
608                         continue;
609                 const Vec3r& P = v->GetVertex();
610                 WVertex::incoming_edge_iterator woeit = v->incoming_edges_begin();
611                 WVertex::incoming_edge_iterator woeitend = v->incoming_edges_end();
612                 for (; woeit != woeitend; ++woeit) {
613                         WOEdge *h = *woeit;
614                         if ((v == start) || h->GetVec() * (O - P) > 0.0) {
615                                 Vec3r V(-1 * h->GetVec());
616                                 bool isect = sphere_clip_vector(O, radius, P, V);
617                                 assert (h->GetOwner()->GetNumberOfOEdges() == 2); // Because otherwise v->isBoundary() would be true
618                                 nc.accumulate_dihedral_angle(V, h->GetAngle());
619
620                                 if (!isect) {
621                                         WVertex *w = h->GetaVertex();
622                                         if (vertices.find(w) == vertices.end()) {
623                                                 vertices.insert(w);
624                                                 S.push(w);
625                                         }
626                                 }
627                         }
628                 }
629         }
630 }
631
632 void compute_curvature_tensor_one_ring(WVertex *start, NormalCycle& nc)
633 {
634         // in case we have a non-manifold vertex, skip it...
635         if (start->isBoundary())
636                 return;
637
638         WVertex::incoming_edge_iterator woeit = start->incoming_edges_begin();
639         WVertex::incoming_edge_iterator woeitend = start->incoming_edges_end();
640         for (; woeit != woeitend; ++woeit) {
641                 WOEdge *h = (*woeit)->twin();
642                 nc.accumulate_dihedral_angle(h->GetVec(), h->GetAngle());
643                 WOEdge *hprev = h->getPrevOnFace();
644                 nc.accumulate_dihedral_angle(hprev->GetVec(), hprev->GetAngle());
645         }
646 }
647
648 }  // OGF namespace