cleanup: style/indentation
[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
218         for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++) {
219                 WOEdge *e = (*itE)->getPrevOnFace();
220                 WVertex *v1 = e->GetaVertex();
221                 WVertex *v2 = e->GetbVertex();
222                 angle_sum += angle_from_cotan(v, v1, v2);
223         }
224
225         *Kg = (2.0 * M_PI - angle_sum) / area;
226
227         return true;
228 }
229
230 /*! gts_vertex_principal_curvatures:
231  *  @Kh: mean curvature.
232  *  @Kg: Gaussian curvature.
233  *  @K1: first principal curvature.
234  *  @K2: second principal curvature.
235  *
236  *  Computes the principal curvatures at a point given the mean and Gaussian curvatures at that point.
237  *
238  *  The mean curvature can be computed as one-half the magnitude of the vector computed by
239  *  gts_vertex_mean_curvature_normal().
240  *
241  *  The Gaussian curvature can be computed with gts_vertex_gaussian_curvature().
242  */
243 void gts_vertex_principal_curvatures (real Kh, real Kg, real *K1, real *K2)
244 {
245         real temp = Kh * Kh - Kg;
246
247         if (!K1 || !K2)
248                 return;
249
250         if (temp < 0.0)
251                 temp = 0.0;
252         temp = sqrt (temp);
253         *K1 = Kh + temp;
254         *K2 = Kh - temp;
255 }
256
257 /* from Maple */
258 static void linsolve(real m11, real m12, real b1, real m21, real m22, real b2, real *x1, real *x2)
259 {
260         real temp;
261
262         temp = 1.0 / (m21 * m12 - m11 * m22);
263         *x1 = (m12 * b2 - m22 * b1) * temp;
264         *x2 = (m11 * b2 - m21 * b1) * temp;
265 }
266
267 /* from Maple - largest eigenvector of [a b; b c] */
268 static void eigenvector(real a, real b, real c, Vec3r e)
269 {
270         if (b == 0.0) {
271                 e[0] = 0.0;
272         }
273         else {
274                 e[0] = -(c - a - sqrt(c * c - 2 * a * c + a * a + 4 * b * b)) / (2 * b);
275         }
276         e[1] = 1.0;
277         e[2] = 0.0;
278 }
279
280 /*! gts_vertex_principal_directions:
281  *  @v: a #WVertex.
282  *  @s: a #GtsSurface.
283  *  @Kh: mean curvature normal (a #Vec3r).
284  *  @Kg: Gaussian curvature (a real).
285  *  @e1: first principal curvature direction (direction of largest curvature).
286  *  @e2: second principal curvature direction.
287  *
288  *  Computes the principal curvature directions at a point given @Kh and @Kg, the mean curvature normal and
289  *  Gaussian curvatures at that point, computed with gts_vertex_mean_curvature_normal() and
290  *  gts_vertex_gaussian_curvature(), respectively.
291  *
292  *  Note that this computation is very approximate and tends to be unstable. Smoothing of the surface or the principal
293  *  directions may be necessary to achieve reasonable results.
294  */
295 void gts_vertex_principal_directions(WVertex *v, Vec3r Kh, real Kg, Vec3r &e1, Vec3r &e2)
296 {
297         Vec3r N;
298         real normKh;
299
300         Vec3r basis1, basis2, d, eig;
301         real ve2, vdotN;
302         real aterm_da, bterm_da, cterm_da, const_da;
303         real aterm_db, bterm_db, cterm_db, const_db;
304         real a, b, c;
305         real K1, K2;
306         real *weights, *kappas, *d1s, *d2s;
307         int edge_count;
308         real err_e1, err_e2;
309         int e;
310         WVertex::incoming_edge_iterator itE;
311
312         /* compute unit normal */
313         normKh = Kh.norm();
314
315         if (normKh > 0.0) {
316                 Kh.normalize();
317         }
318         else {
319                 /* This vertex is a point of zero mean curvature (flat or saddle point). Compute a normal by averaging
320                  * the adjacent triangles
321                  */
322                 N[0] = N[1] = N[2] = 0.0;
323
324                 for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++)
325                         N = Vec3r(N + (*itE)->GetaFace()->GetNormal());
326                 real normN = N.norm();
327                 if (normN <= 0.0)
328                         return;
329                 N.normalize();
330         }
331
332         /* construct a basis from N: */
333         /* set basis1 to any component not the largest of N */
334         basis1[0] =  basis1[1] =  basis1[2] = 0.0;
335         if (fabs (N[0]) > fabs (N[1]))
336                 basis1[1] = 1.0;
337         else
338                 basis1[0] = 1.0;
339
340         /* make basis2 orthogonal to N */
341         basis2 = (N ^ basis1);
342         basis2.normalize();
343
344         /* make basis1 orthogonal to N and basis2 */
345         basis1 = (N ^ basis2);
346         basis1.normalize();
347
348         aterm_da = bterm_da = cterm_da = const_da = 0.0;
349         aterm_db = bterm_db = cterm_db = const_db = 0.0;
350         int nb_edges = v->GetEdges().size();
351
352         weights = (real *)malloc(sizeof(real) * nb_edges);
353         kappas = (real *)malloc(sizeof(real) * nb_edges);
354         d1s = (real *)malloc(sizeof(real) * nb_edges);
355         d2s = (real *)malloc(sizeof(real) * nb_edges);
356         edge_count = 0;
357
358         for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++) {
359                 WOEdge *e;
360                 WFace *f1, *f2;
361                 real weight, kappa, d1, d2;
362                 Vec3r vec_edge;
363                 if (!*itE)
364                         continue;
365                 e = *itE;
366
367                 /* since this vertex passed the tests in gts_vertex_mean_curvature_normal(), this should be true. */
368                 //g_assert(gts_edge_face_number (e, s) == 2);
369
370                 /* identify the two triangles bordering e in s */
371                 f1 = e->GetaFace();
372                 f2 = e->GetbFace();
373
374                 /* We are solving for the values of the curvature tensor
375                 *     B = [ a b ; b c ].
376                 *  The computations here are from section 5 of [Meyer et al 2002].
377                 *
378                 *  The first step is to calculate the linear equations governing the values of (a,b,c). These can be computed
379                 *  by setting the derivatives of the error E to zero (section 5.3).
380                 *
381                 *  Since a + c = norm(Kh), we only compute the linear equations for dE/da and dE/db. (NB: [Meyer et al 2002]
382                 *  has the equation a + b = norm(Kh), but I'm almost positive this is incorrect).
383                 *
384                 *  Note that the w_ij (defined in section 5.2) are all scaled by (1/8*A_mixed). We drop this uniform scale
385                 *  factor because the solution of the linear equations doesn't rely on it.
386                 *
387                 *  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
388                 *  terms that are the constant factors in the equations.
389                 */
390
391                 /* find the vector from v along edge e */
392                 vec_edge = Vec3r(-1 * e->GetVec());
393
394                 ve2 = vec_edge.squareNorm();
395                 vdotN = vec_edge * N;
396
397                 /* section 5.2 - There is a typo in the computation of kappa. The edges should be x_j-x_i. */
398                 kappa = 2.0 * vdotN / ve2;
399
400                 /* section 5.2 */
401
402                 /* I don't like performing a minimization where some of the weights can be negative (as can be the case
403                 *  if f1 or f2 are obtuse). To ensure all-positive weights, we check for obtuseness. */
404                 weight = 0.0;
405                 if (!triangle_obtuse(v, f1)) {
406                         weight += ve2 * cotan(f1->GetNextOEdge(e->twin())->GetbVertex(), e->GetaVertex(), e->GetbVertex()) / 8.0;
407                 }
408                 else {
409                         if (angle_obtuse(v, f1)) {
410                                 weight += ve2 * f1->getArea() / 4.0;
411                         }
412                         else {
413                                 weight += ve2 * f1->getArea() / 8.0;
414                         }
415                 }
416
417                 if (!triangle_obtuse(v, f2)) {
418                         weight += ve2 * cotan (f2->GetNextOEdge(e)->GetbVertex(), e->GetaVertex(), e->GetbVertex()) / 8.0;
419                 }
420                 else {
421                         if (angle_obtuse(v, f2)) {
422                                 weight += ve2 * f1->getArea() / 4.0;
423                         }
424                         else {
425                                 weight += ve2 * f1->getArea() / 8.0;
426                         }
427                 }
428
429                 /* projection of edge perpendicular to N (section 5.3) */
430                 d[0] = vec_edge[0] - vdotN * N[0];
431                 d[1] = vec_edge[1] - vdotN * N[1];
432                 d[2] = vec_edge[2] - vdotN * N[2];
433                 d.normalize();
434
435                 /* not explicit in the paper, but necessary. Move d to 2D basis. */
436                 d1 = d * basis1;
437                 d2 = d * basis2;
438
439                 /* store off the curvature, direction of edge, and weights for later use */
440                 weights[edge_count] = weight;
441                 kappas[edge_count] = kappa;
442                 d1s[edge_count] = d1;
443                 d2s[edge_count] = d2;
444                 edge_count++;
445
446                 /* Finally, update the linear equations */
447                 aterm_da += weight * d1 * d1 * d1 * d1;
448                 bterm_da += weight * d1 * d1 * 2 * d1 * d2;
449                 cterm_da += weight * d1 * d1 * d2 * d2;
450                 const_da += weight * d1 * d1 * (-kappa);
451
452                 aterm_db += weight * d1 * d2 * d1 * d1;
453                 bterm_db += weight * d1 * d2 * 2 * d1 * d2;
454                 cterm_db += weight * d1 * d2 * d2 * d2;
455                 const_db += weight * d1 * d2 * (-kappa);
456         }
457
458         /* now use the identity (Section 5.3) a + c = |Kh| = 2 * kappa_h */
459         aterm_da -= cterm_da;
460         const_da += cterm_da * normKh;
461
462         aterm_db -= cterm_db;
463         const_db += cterm_db * normKh;
464
465         /* check for solvability of the linear system */
466         if (((aterm_da * bterm_db - aterm_db * bterm_da) != 0.0) && ((const_da != 0.0) || (const_db != 0.0))) {
467                 linsolve(aterm_da, bterm_da, -const_da, aterm_db, bterm_db, -const_db, &a, &b);
468
469                 c = normKh - a;
470
471                 eigenvector(a, b, c, eig);
472         }
473         else {
474                 /* region of v is planar */
475                 eig[0] = 1.0;
476                 eig[1] = 0.0;
477         }
478
479         /* Although the eigenvectors of B are good estimates of the principal directions, it seems that which one is
480          * attached to which curvature direction is a bit arbitrary. This may be a bug in my implementation, or just
481          * a side-effect of the inaccuracy of B due to the discrete nature of the sampling.
482          *
483          * To overcome this behavior, we'll evaluate which assignment best matches the given eigenvectors by comparing
484          * the curvature estimates computed above and the curvatures calculated from the discrete differential operators.
485          */
486
487         gts_vertex_principal_curvatures(0.5 * normKh, Kg, &K1, &K2);
488
489         err_e1 = err_e2 = 0.0;
490         /* loop through the values previously saved */
491         for (e = 0; e < edge_count; e++) {
492                 real weight, kappa, d1, d2;
493                 real temp1, temp2;
494                 real delta;
495
496                 weight = weights[e];
497                 kappa = kappas[e];
498                 d1 = d1s[e];
499                 d2 = d2s[e];
500
501                 temp1 = fabs (eig[0] * d1 + eig[1] * d2);
502                 temp1 = temp1 * temp1;
503                 temp2 = fabs (eig[1] * d1 - eig[0] * d2);
504                 temp2 = temp2 * temp2;
505
506                 /* err_e1 is for K1 associated with e1 */
507                 delta = K1 * temp1 + K2 * temp2 - kappa;
508                 err_e1 += weight * delta * delta;
509
510                 /* err_e2 is for K1 associated with e2 */
511                 delta = K2 * temp1 + K1 * temp2 - kappa;
512                 err_e2 += weight * delta * delta;
513         }
514         free (weights);
515         free (kappas);
516         free (d1s);
517         free (d2s);
518
519         /* rotate eig by a right angle if that would decrease the error */
520         if (err_e2 < err_e1) {
521                 real temp = eig[0];
522
523                 eig[0] = eig[1];
524                 eig[1] = -temp;
525         }
526
527         e1[0] = eig[0] * basis1[0] + eig[1] * basis2[0];
528         e1[1] = eig[0] * basis1[1] + eig[1] * basis2[1];
529         e1[2] = eig[0] * basis1[2] + eig[1] * basis2[2];
530         e1.normalize();
531
532         /* make N,e1,e2 a right handed coordinate sytem */
533         e2 =  N ^ e1;
534         e2.normalize();
535 }
536
537 namespace OGF {
538
539 #if 0
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->GetVec();
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 #endif
555
556 // precondition1: P is inside the sphere
557 // precondition2: P,V points to the outside of the sphere (i.e. OP.V > 0)
558 static bool sphere_clip_vector(const Vec3r& O, real r, const Vec3r& P, Vec3r& V)
559 {
560         Vec3r W = P - O;
561         real a = V.squareNorm();
562         real b = 2.0 * V * W;
563         real c = W.squareNorm() - r * r;
564         real delta = b * b - 4 * a * c;
565         if (delta < 0) {
566                 // Should not happen, but happens sometimes (numerical precision)
567                 return true;
568         }
569         real t = - b + ::sqrt(delta) / (2.0 * a);
570         if (t < 0.0) {
571                 // Should not happen, but happens sometimes (numerical precision)
572                 return true;
573         }
574         if (t >= 1.0) {
575                 // Inside the sphere
576                 return false;
577         }
578
579         V[0] = (t * V.x());
580         V[1] = (t * V.y());
581         V[2] = (t * V.z());
582
583         return true;
584 }
585
586 // TODO: check optimizations:
587 // use marking ? (measure *timings* ...)
588 void compute_curvature_tensor(WVertex *start, real radius, NormalCycle& nc)
589 {
590         // in case we have a non-manifold vertex, skip it...
591         if (start->isBoundary())
592                 return;
593
594         std::set<WVertex*> vertices;
595         const Vec3r& O = start->GetVertex();
596         std::stack<WVertex*> S;
597         S.push(start);
598         vertices.insert(start);
599         while (!S.empty()) {
600                 WVertex *v = S.top();
601                 S.pop();
602                 if (v->isBoundary())
603                         continue;
604                 const Vec3r& P = v->GetVertex();
605                 WVertex::incoming_edge_iterator woeit = v->incoming_edges_begin();
606                 WVertex::incoming_edge_iterator woeitend = v->incoming_edges_end();
607                 for (; woeit != woeitend; ++woeit) {
608                         WOEdge *h = *woeit;
609                         if ((v == start) || h->GetVec() * (O - P) > 0.0) {
610                                 Vec3r V(-1 * h->GetVec());
611                                 bool isect = sphere_clip_vector(O, radius, P, V);
612                                 assert (h->GetOwner()->GetNumberOfOEdges() == 2); // Because otherwise v->isBoundary() would be true
613                                 nc.accumulate_dihedral_angle(V, h->GetAngle());
614
615                                 if (!isect) {
616                                         WVertex *w = h->GetaVertex();
617                                         if (vertices.find(w) == vertices.end()) {
618                                                 vertices.insert(w);
619                                                 S.push(w);
620                                         }
621                                 }
622                         }
623                 }
624         }
625 }
626
627 void compute_curvature_tensor_one_ring(WVertex *start, NormalCycle& nc)
628 {
629         // in case we have a non-manifold vertex, skip it...
630         if (start->isBoundary())
631                 return;
632
633         WVertex::incoming_edge_iterator woeit = start->incoming_edges_begin();
634         WVertex::incoming_edge_iterator woeitend = start->incoming_edges_end();
635         for (; woeit != woeitend; ++woeit) {
636                 WOEdge *h = (*woeit)->twin();
637                 nc.accumulate_dihedral_angle(h->GetVec(), h->GetAngle());
638                 WOEdge *hprev = h->getPrevOnFace();
639                 nc.accumulate_dihedral_angle(hprev->GetVec(), hprev->GetAngle());
640         }
641 }
642
643 }  // OGF namespace
644
645 } /* namespace Freestyle */