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