ClangFormat: apply to source, most of intern
[blender.git] / source / blender / freestyle / intern / view_map / ViewMapBuilder.cpp
1 /*
2  * This program is free software; you can redistribute it and/or
3  * modify it under the terms of the GNU General Public License
4  * as published by the Free Software Foundation; either version 2
5  * of the License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program; if not, write to the Free Software Foundation,
14  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
15  */
16
17 /** \file
18  * \ingroup freestyle
19  * \brief Class to build silhouette edges from a Winged-Edge structure
20  */
21
22 #include <algorithm>
23 #include <memory>
24 #include <stdexcept>
25 #include <sstream>
26
27 #include "FRS_freestyle.h"
28
29 #include "BoxGrid.h"
30 #include "CulledOccluderSource.h"
31 #include "HeuristicGridDensityProviderFactory.h"
32 #include "OccluderSource.h"
33 #include "SphericalGrid.h"
34 #include "ViewMapBuilder.h"
35
36 #include "../geometry/GridHelpers.h"
37 #include "../geometry/GeomUtils.h"
38
39 #include "../winged_edge/WFillGrid.h"
40
41 #include "BKE_global.h"
42
43 namespace Freestyle {
44
45 // XXX Grmll... G is used as template's typename parameter :/
46 static const Global &_global = G;
47
48 #define LOGGING 0
49
50 using namespace std;
51
52 template<typename G, typename I>
53 static void findOccludee(FEdge *fe,
54                          G & /*grid*/,
55                          I &occluders,
56                          real epsilon,
57                          WFace **oaWFace,
58                          Vec3r &u,
59                          Vec3r &A,
60                          Vec3r &origin,
61                          Vec3r &edgeDir,
62                          vector<WVertex *> &faceVertices)
63 {
64   WFace *face = NULL;
65   if (fe->isSmooth()) {
66     FEdgeSmooth *fes = dynamic_cast<FEdgeSmooth *>(fe);
67     face = (WFace *)fes->face();
68   }
69   WFace *oface;
70   bool skipFace;
71
72   WVertex::incoming_edge_iterator ie;
73
74   *oaWFace = NULL;
75   if (((fe)->getNature() & Nature::SILHOUETTE) || ((fe)->getNature() & Nature::BORDER)) {
76     // we cast a ray from A in the same direction but looking behind
77     Vec3r v(-u[0], -u[1], -u[2]);
78     bool noIntersection = true;
79     real mint = FLT_MAX;
80
81     for (occluders.initAfterTarget(); occluders.validAfterTarget(); occluders.nextOccludee()) {
82 #if LOGGING
83       if (_global.debug & G_DEBUG_FREESTYLE) {
84         cout << "\t\tEvaluating intersection for occludee " << occluders.getWFace() << " and ray "
85              << A << " * " << u << endl;
86       }
87 #endif
88       oface = occluders.getWFace();
89       Polygon3r *p = occluders.getCameraSpacePolygon();
90       real d = -((p->getVertices())[0] * p->getNormal());
91       real t, t_u, t_v;
92
93       if (0 != face) {
94         skipFace = false;
95
96         if (face == oface)
97           continue;
98
99         if (faceVertices.empty())
100           continue;
101
102         for (vector<WVertex *>::iterator fv = faceVertices.begin(), fvend = faceVertices.end();
103              fv != fvend;
104              ++fv) {
105           if ((*fv)->isBoundary())
106             continue;
107           WVertex::incoming_edge_iterator iebegin = (*fv)->incoming_edges_begin();
108           WVertex::incoming_edge_iterator ieend = (*fv)->incoming_edges_end();
109           for (ie = iebegin; ie != ieend; ++ie) {
110             if ((*ie) == 0)
111               continue;
112
113             WFace *sface = (*ie)->GetbFace();
114             if (sface == oface) {
115               skipFace = true;
116               break;
117             }
118           }
119           if (skipFace)
120             break;
121         }
122         if (skipFace)
123           continue;
124       }
125       else {
126         // check whether the edge and the polygon plane are coincident:
127         //-------------------------------------------------------------
128         //first let us compute the plane equation.
129         if (GeomUtils::COINCIDENT ==
130             GeomUtils::intersectRayPlane(origin, edgeDir, p->getNormal(), d, t, epsilon)) {
131 #if LOGGING
132           if (_global.debug & G_DEBUG_FREESTYLE) {
133             cout << "\t\tRejecting occluder for target coincidence." << endl;
134           }
135 #endif
136           continue;
137         }
138       }
139
140       if (p->rayIntersect(A, v, t, t_u, t_v)) {
141 #if LOGGING
142         if (_global.debug & G_DEBUG_FREESTYLE) {
143           cout << "\t\tRay " << A << " * " << v << " intersects at time " << t << endl;
144           cout << "\t\t(v * normal) == " << (v * p->getNormal()) << " for normal "
145                << p->getNormal() << endl;
146         }
147 #endif
148         if (fabs(v * p->getNormal()) > 0.0001) {
149           if ((t > 0.0)) {  // && (t<1.0))
150             if (t < mint) {
151               *oaWFace = oface;
152               mint = t;
153               noIntersection = false;
154               fe->setOccludeeIntersection(Vec3r(A + t * v));
155 #if LOGGING
156               if (_global.debug & G_DEBUG_FREESTYLE) {
157                 cout << "\t\tIs occludee" << endl;
158               }
159 #endif
160             }
161           }
162         }
163
164         occluders.reportDepth(A, v, t);
165       }
166     }
167
168     if (noIntersection)
169       *oaWFace = NULL;
170   }
171 }
172
173 template<typename G, typename I>
174 static void findOccludee(FEdge *fe, G &grid, real epsilon, ViewEdge * /*ve*/, WFace **oaFace)
175 {
176   Vec3r A;
177   Vec3r edgeDir;
178   Vec3r origin;
179   A = Vec3r(((fe)->vertexA()->point3D() + (fe)->vertexB()->point3D()) / 2.0);
180   edgeDir = Vec3r((fe)->vertexB()->point3D() - (fe)->vertexA()->point3D());
181   edgeDir.normalize();
182   origin = Vec3r((fe)->vertexA()->point3D());
183   Vec3r u;
184   if (grid.orthographicProjection()) {
185     u = Vec3r(0.0, 0.0, grid.viewpoint().z() - A.z());
186   }
187   else {
188     u = Vec3r(grid.viewpoint() - A);
189   }
190   u.normalize();
191
192   vector<WVertex *> faceVertices;
193
194   WFace *face = NULL;
195   if (fe->isSmooth()) {
196     FEdgeSmooth *fes = dynamic_cast<FEdgeSmooth *>(fe);
197     face = (WFace *)fes->face();
198   }
199
200   if (face) {
201     face->RetrieveVertexList(faceVertices);
202   }
203
204   I occluders(grid, A, epsilon);
205   findOccludee<G, I>(fe, grid, occluders, epsilon, oaFace, u, A, origin, edgeDir, faceVertices);
206 }
207
208 // computeVisibility takes a pointer to foundOccluders, instead of using a reference,
209 // so that computeVeryFastVisibility can skip the AddOccluders step with minimal overhead.
210 template<typename G, typename I>
211 static int computeVisibility(ViewMap *viewMap,
212                              FEdge *fe,
213                              G &grid,
214                              real epsilon,
215                              ViewEdge * /*ve*/,
216                              WFace **oaWFace,
217                              set<ViewShape *> *foundOccluders)
218 {
219   int qi = 0;
220
221   Vec3r center;
222   Vec3r edgeDir;
223   Vec3r origin;
224
225   center = fe->center3d();
226   edgeDir = Vec3r(fe->vertexB()->point3D() - fe->vertexA()->point3D());
227   edgeDir.normalize();
228   origin = Vec3r(fe->vertexA()->point3D());
229
230   Vec3r vp;
231   if (grid.orthographicProjection()) {
232     vp = Vec3r(center.x(), center.y(), grid.viewpoint().z());
233   }
234   else {
235     vp = Vec3r(grid.viewpoint());
236   }
237   Vec3r u(vp - center);
238   real raylength = u.norm();
239   u.normalize();
240
241   WFace *face = NULL;
242   if (fe->isSmooth()) {
243     FEdgeSmooth *fes = dynamic_cast<FEdgeSmooth *>(fe);
244     face = (WFace *)fes->face();
245   }
246   vector<WVertex *> faceVertices;
247   WVertex::incoming_edge_iterator ie;
248
249   WFace *oface;
250   bool skipFace;
251
252   if (face)
253     face->RetrieveVertexList(faceVertices);
254
255   I occluders(grid, center, epsilon);
256
257   for (occluders.initBeforeTarget(); occluders.validBeforeTarget(); occluders.nextOccluder()) {
258     // If we're dealing with an exact silhouette, check whether we must take care of this occluder of not.
259     // (Indeed, we don't consider the occluders that share at least one vertex with the face containing this edge).
260     //-----------
261     oface = occluders.getWFace();
262     Polygon3r *p = occluders.getCameraSpacePolygon();
263     real t, t_u, t_v;
264 #if LOGGING
265     if (_global.debug & G_DEBUG_FREESTYLE) {
266       cout << "\t\tEvaluating intersection for occluder " << (p->getVertices())[0]
267            << (p->getVertices())[1] << (p->getVertices())[2] << endl
268            << "\t\t\tand ray " << vp << " * " << u << " (center " << center << ")" << endl;
269     }
270 #endif
271
272 #if LOGGING
273     Vec3r v(vp - center);
274     real rl = v.norm();
275     v.normalize();
276     vector<Vec3r> points;
277     // Iterate over vertices, storing projections in points
278     for (vector<WOEdge *>::const_iterator woe = oface->getEdgeList().begin(),
279                                           woend = oface->getEdgeList().end();
280          woe != woend;
281          woe++) {
282       points.push_back(Vec3r((*woe)->GetaVertex()->GetVertex()));
283     }
284     Polygon3r p1(points, oface->GetNormal());
285     Vec3r v1((p1.getVertices())[0]);
286     real d = -(v1 * p->getNormal());
287     if (_global.debug & G_DEBUG_FREESTYLE) {
288       cout << "\t\tp:  " << (p->getVertices())[0] << (p->getVertices())[1] << (p->getVertices())[2]
289            << ", norm: " << p->getNormal() << endl;
290       cout << "\t\tp1: " << (p1.getVertices())[0] << (p1.getVertices())[1] << (p1.getVertices())[2]
291            << ", norm: " << p1.getNormal() << endl;
292     }
293 #else
294     real d = -((p->getVertices())[0] * p->getNormal());
295 #endif
296
297     if (face) {
298 #if LOGGING
299       if (_global.debug & G_DEBUG_FREESTYLE) {
300         cout << "\t\tDetermining face adjacency...";
301       }
302 #endif
303       skipFace = false;
304
305       if (face == oface) {
306 #if LOGGING
307         if (_global.debug & G_DEBUG_FREESTYLE) {
308           cout << "  Rejecting occluder for face concurrency." << endl;
309         }
310 #endif
311         continue;
312       }
313
314       for (vector<WVertex *>::iterator fv = faceVertices.begin(), fvend = faceVertices.end();
315            fv != fvend;
316            ++fv) {
317         if ((*fv)->isBoundary())
318           continue;
319
320         WVertex::incoming_edge_iterator iebegin = (*fv)->incoming_edges_begin();
321         WVertex::incoming_edge_iterator ieend = (*fv)->incoming_edges_end();
322         for (ie = iebegin; ie != ieend; ++ie) {
323           if ((*ie) == 0)
324             continue;
325
326           WFace *sface = (*ie)->GetbFace();
327           //WFace *sfacea = (*ie)->GetaFace();
328           //if ((sface == oface) || (sfacea == oface))
329           if (sface == oface) {
330             skipFace = true;
331             break;
332           }
333         }
334         if (skipFace)
335           break;
336       }
337       if (skipFace) {
338 #if LOGGING
339         if (_global.debug & G_DEBUG_FREESTYLE) {
340           cout << "  Rejecting occluder for face adjacency." << endl;
341         }
342 #endif
343         continue;
344       }
345     }
346     else {
347       // check whether the edge and the polygon plane are coincident:
348       //-------------------------------------------------------------
349       //first let us compute the plane equation.
350       if (GeomUtils::COINCIDENT ==
351           GeomUtils::intersectRayPlane(origin, edgeDir, p->getNormal(), d, t, epsilon)) {
352 #if LOGGING
353         if (_global.debug & G_DEBUG_FREESTYLE) {
354           cout << "\t\tRejecting occluder for target coincidence." << endl;
355         }
356 #endif
357         continue;
358       }
359     }
360
361 #if LOGGING
362     if (_global.debug & G_DEBUG_FREESTYLE) {
363       real x;
364       if (p1.rayIntersect(center, v, x, t_u, t_v)) {
365         cout << "\t\tRay should intersect at time " << (rl - x) << ". Center: " << center
366              << ", V: " << v << ", RL: " << rl << ", T:" << x << endl;
367       }
368       else {
369         cout << "\t\tRay should not intersect. Center: " << center << ", V: " << v
370              << ", RL: " << rl << endl;
371       }
372     }
373 #endif
374
375     if (p->rayIntersect(center, u, t, t_u, t_v)) {
376 #if LOGGING
377       if (_global.debug & G_DEBUG_FREESTYLE) {
378         cout << "\t\tRay " << center << " * " << u << " intersects at time " << t
379              << " (raylength is " << raylength << ")" << endl;
380         cout << "\t\t(u * normal) == " << (u * p->getNormal()) << " for normal " << p->getNormal()
381              << endl;
382       }
383 #endif
384       if (fabs(u * p->getNormal()) > 0.0001) {
385         if ((t > 0.0) && (t < raylength)) {
386 #if LOGGING
387           if (_global.debug & G_DEBUG_FREESTYLE) {
388             cout << "\t\tIs occluder" << endl;
389           }
390 #endif
391           if (foundOccluders != NULL) {
392             ViewShape *vshape = viewMap->viewShape(oface->GetVertex(0)->shape()->GetId());
393             foundOccluders->insert(vshape);
394           }
395           ++qi;
396
397           if (!grid.enableQI())
398             break;
399         }
400
401         occluders.reportDepth(center, u, t);
402       }
403     }
404   }
405
406   // Find occludee
407   findOccludee<G, I>(
408       fe, grid, occluders, epsilon, oaWFace, u, center, origin, edgeDir, faceVertices);
409
410   return qi;
411 }
412
413 // computeCumulativeVisibility returns the lowest x such that the majority of FEdges have QI <= x
414 //
415 // This was probably the original intention of the "normal" algorithm on which computeDetailedVisibility is based.
416 // But because the "normal" algorithm chooses the most popular QI, without considering any other values, a ViewEdge
417 // with FEdges having QIs of 0, 21, 22, 23, 24 and 25 will end up having a total QI of 0, even though most of the
418 // FEdges are heavily occluded. computeCumulativeVisibility will treat this case as a QI of 22 because 3 out of
419 // 6 occluders have QI <= 22.
420
421 template<typename G, typename I>
422 static void computeCumulativeVisibility(ViewMap *ioViewMap,
423                                         G &grid,
424                                         real epsilon,
425                                         RenderMonitor *iRenderMonitor)
426 {
427   vector<ViewEdge *> &vedges = ioViewMap->ViewEdges();
428
429   FEdge *fe, *festart;
430   int nSamples = 0;
431   vector<WFace *> wFaces;
432   WFace *wFace = NULL;
433   unsigned cnt = 0;
434   unsigned cntStep = (unsigned)ceil(0.01f * vedges.size());
435   unsigned tmpQI = 0;
436   unsigned qiClasses[256];
437   unsigned maxIndex, maxCard;
438   unsigned qiMajority;
439   for (vector<ViewEdge *>::iterator ve = vedges.begin(), veend = vedges.end(); ve != veend; ve++) {
440     if (iRenderMonitor) {
441       if (iRenderMonitor->testBreak())
442         break;
443       if (cnt % cntStep == 0) {
444         stringstream ss;
445         ss << "Freestyle: Visibility computations " << (100 * cnt / vedges.size()) << "%";
446         iRenderMonitor->setInfo(ss.str());
447         iRenderMonitor->progress((float)cnt / vedges.size());
448       }
449       cnt++;
450     }
451 #if LOGGING
452     if (_global.debug & G_DEBUG_FREESTYLE) {
453       cout << "Processing ViewEdge " << (*ve)->getId() << endl;
454     }
455 #endif
456     // Find an edge to test
457     if (!(*ve)->isInImage()) {
458       // This view edge has been proscenium culled
459       (*ve)->setQI(255);
460       (*ve)->setaShape(0);
461 #if LOGGING
462       if (_global.debug & G_DEBUG_FREESTYLE) {
463         cout << "\tCulled." << endl;
464       }
465 #endif
466       continue;
467     }
468
469     // Test edge
470     festart = (*ve)->fedgeA();
471     fe = (*ve)->fedgeA();
472     qiMajority = 0;
473     do {
474       if (fe != NULL && fe->isInImage()) {
475         qiMajority++;
476       }
477       fe = fe->nextEdge();
478     } while (fe && fe != festart);
479
480     if (qiMajority == 0) {
481       // There are no occludable FEdges on this ViewEdge
482       // This should be impossible.
483       if (_global.debug & G_DEBUG_FREESTYLE) {
484         cout << "View Edge in viewport without occludable FEdges: " << (*ve)->getId() << endl;
485       }
486       // We can recover from this error:
487       // Treat this edge as fully visible with no occludee
488       (*ve)->setQI(0);
489       (*ve)->setaShape(0);
490       continue;
491     }
492     else {
493       ++qiMajority;
494       qiMajority >>= 1;
495     }
496 #if LOGGING
497     if (_global.debug & G_DEBUG_FREESTYLE) {
498       cout << "\tqiMajority: " << qiMajority << endl;
499     }
500 #endif
501
502     tmpQI = 0;
503     maxIndex = 0;
504     maxCard = 0;
505     nSamples = 0;
506     memset(qiClasses, 0, 256 * sizeof(*qiClasses));
507     set<ViewShape *> foundOccluders;
508
509     fe = (*ve)->fedgeA();
510     do {
511       if (!fe || !fe->isInImage()) {
512         fe = fe->nextEdge();
513         continue;
514       }
515       if ((maxCard < qiMajority)) {
516         //ARB: change &wFace to wFace and use reference in called function
517         tmpQI = computeVisibility<G, I>(
518             ioViewMap, fe, grid, epsilon, *ve, &wFace, &foundOccluders);
519 #if LOGGING
520         if (_global.debug & G_DEBUG_FREESTYLE) {
521           cout << "\tFEdge: visibility " << tmpQI << endl;
522         }
523 #endif
524
525         //ARB: This is an error condition, not an alert condition.
526         // Some sort of recovery or abort is necessary.
527         if (tmpQI >= 256) {
528           cerr << "Warning: too many occluding levels" << endl;
529           //ARB: Wild guess: instead of aborting or corrupting memory, treat as tmpQI == 255
530           tmpQI = 255;
531         }
532
533         if (++qiClasses[tmpQI] > maxCard) {
534           maxCard = qiClasses[tmpQI];
535           maxIndex = tmpQI;
536         }
537       }
538       else {
539         //ARB: FindOccludee is redundant if ComputeRayCastingVisibility has been called
540         //ARB: change &wFace to wFace and use reference in called function
541         findOccludee<G, I>(fe, grid, epsilon, *ve, &wFace);
542 #if LOGGING
543         if (_global.debug & G_DEBUG_FREESTYLE) {
544           cout << "\tFEdge: occludee only (" << (wFace != NULL ? "found" : "not found") << ")"
545                << endl;
546         }
547 #endif
548       }
549
550       // Store test results
551       if (wFace) {
552         vector<Vec3r> vertices;
553         for (int i = 0, numEdges = wFace->numberOfEdges(); i < numEdges; ++i) {
554           vertices.push_back(Vec3r(wFace->GetVertex(i)->GetVertex()));
555         }
556         Polygon3r poly(vertices, wFace->GetNormal());
557         poly.userdata = (void *)wFace;
558         fe->setaFace(poly);
559         wFaces.push_back(wFace);
560         fe->setOccludeeEmpty(false);
561 #if LOGGING
562         if (_global.debug & G_DEBUG_FREESTYLE) {
563           cout << "\tFound occludee" << endl;
564         }
565 #endif
566       }
567       else {
568         fe->setOccludeeEmpty(true);
569       }
570
571       ++nSamples;
572       fe = fe->nextEdge();
573     } while ((maxCard < qiMajority) && (fe) && (fe != festart));
574
575 #if LOGGING
576     if (_global.debug & G_DEBUG_FREESTYLE) {
577       cout << "\tFinished with " << nSamples << " samples, maxCard = " << maxCard << endl;
578     }
579 #endif
580
581     // ViewEdge
582     // qi --
583     // Find the minimum value that is >= the majority of the QI
584     for (unsigned count = 0, i = 0; i < 256; ++i) {
585       count += qiClasses[i];
586       if (count >= qiMajority) {
587         (*ve)->setQI(i);
588         break;
589       }
590     }
591     // occluders --
592     // I would rather not have to go through the effort of creating this set and then copying out its contents.
593     // Is there a reason why ViewEdge::_Occluders cannot be converted to a set<>?
594     for (set<ViewShape *>::iterator o = foundOccluders.begin(), oend = foundOccluders.end();
595          o != oend;
596          ++o) {
597       (*ve)->AddOccluder((*o));
598     }
599 #if LOGGING
600     if (_global.debug & G_DEBUG_FREESTYLE) {
601       cout << "\tConclusion: QI = " << maxIndex << ", " << (*ve)->occluders_size() << " occluders."
602            << endl;
603     }
604 #else
605     (void)maxIndex;
606 #endif
607     // occludee --
608     if (!wFaces.empty()) {
609       if (wFaces.size() <= (float)nSamples / 2.0f) {
610         (*ve)->setaShape(0);
611       }
612       else {
613         ViewShape *vshape = ioViewMap->viewShape(
614             (*wFaces.begin())->GetVertex(0)->shape()->GetId());
615         (*ve)->setaShape(vshape);
616       }
617     }
618
619     wFaces.clear();
620   }
621   if (iRenderMonitor && vedges.size()) {
622     stringstream ss;
623     ss << "Freestyle: Visibility computations " << (100 * cnt / vedges.size()) << "%";
624     iRenderMonitor->setInfo(ss.str());
625     iRenderMonitor->progress((float)cnt / vedges.size());
626   }
627 }
628
629 template<typename G, typename I>
630 static void computeDetailedVisibility(ViewMap *ioViewMap,
631                                       G &grid,
632                                       real epsilon,
633                                       RenderMonitor *iRenderMonitor)
634 {
635   vector<ViewEdge *> &vedges = ioViewMap->ViewEdges();
636
637   FEdge *fe, *festart;
638   int nSamples = 0;
639   vector<WFace *> wFaces;
640   WFace *wFace = NULL;
641   unsigned tmpQI = 0;
642   unsigned qiClasses[256];
643   unsigned maxIndex, maxCard;
644   unsigned qiMajority;
645   for (vector<ViewEdge *>::iterator ve = vedges.begin(), veend = vedges.end(); ve != veend; ve++) {
646     if (iRenderMonitor && iRenderMonitor->testBreak())
647       break;
648 #if LOGGING
649     if (_global.debug & G_DEBUG_FREESTYLE) {
650       cout << "Processing ViewEdge " << (*ve)->getId() << endl;
651     }
652 #endif
653     // Find an edge to test
654     if (!(*ve)->isInImage()) {
655       // This view edge has been proscenium culled
656       (*ve)->setQI(255);
657       (*ve)->setaShape(0);
658 #if LOGGING
659       if (_global.debug & G_DEBUG_FREESTYLE) {
660         cout << "\tCulled." << endl;
661       }
662 #endif
663       continue;
664     }
665
666     // Test edge
667     festart = (*ve)->fedgeA();
668     fe = (*ve)->fedgeA();
669     qiMajority = 0;
670     do {
671       if (fe != NULL && fe->isInImage()) {
672         qiMajority++;
673       }
674       fe = fe->nextEdge();
675     } while (fe && fe != festart);
676
677     if (qiMajority == 0) {
678       // There are no occludable FEdges on this ViewEdge
679       // This should be impossible.
680       if (_global.debug & G_DEBUG_FREESTYLE) {
681         cout << "View Edge in viewport without occludable FEdges: " << (*ve)->getId() << endl;
682       }
683       // We can recover from this error:
684       // Treat this edge as fully visible with no occludee
685       (*ve)->setQI(0);
686       (*ve)->setaShape(0);
687       continue;
688     }
689     else {
690       ++qiMajority;
691       qiMajority >>= 1;
692     }
693 #if LOGGING
694     if (_global.debug & G_DEBUG_FREESTYLE) {
695       cout << "\tqiMajority: " << qiMajority << endl;
696     }
697 #endif
698
699     tmpQI = 0;
700     maxIndex = 0;
701     maxCard = 0;
702     nSamples = 0;
703     memset(qiClasses, 0, 256 * sizeof(*qiClasses));
704     set<ViewShape *> foundOccluders;
705
706     fe = (*ve)->fedgeA();
707     do {
708       if (fe == NULL || !fe->isInImage()) {
709         fe = fe->nextEdge();
710         continue;
711       }
712       if ((maxCard < qiMajority)) {
713         //ARB: change &wFace to wFace and use reference in called function
714         tmpQI = computeVisibility<G, I>(
715             ioViewMap, fe, grid, epsilon, *ve, &wFace, &foundOccluders);
716 #if LOGGING
717         if (_global.debug & G_DEBUG_FREESTYLE) {
718           cout << "\tFEdge: visibility " << tmpQI << endl;
719         }
720 #endif
721
722         //ARB: This is an error condition, not an alert condition.
723         // Some sort of recovery or abort is necessary.
724         if (tmpQI >= 256) {
725           cerr << "Warning: too many occluding levels" << endl;
726           //ARB: Wild guess: instead of aborting or corrupting memory, treat as tmpQI == 255
727           tmpQI = 255;
728         }
729
730         if (++qiClasses[tmpQI] > maxCard) {
731           maxCard = qiClasses[tmpQI];
732           maxIndex = tmpQI;
733         }
734       }
735       else {
736         //ARB: FindOccludee is redundant if ComputeRayCastingVisibility has been called
737         //ARB: change &wFace to wFace and use reference in called function
738         findOccludee<G, I>(fe, grid, epsilon, *ve, &wFace);
739 #if LOGGING
740         if (_global.debug & G_DEBUG_FREESTYLE) {
741           cout << "\tFEdge: occludee only (" << (wFace != NULL ? "found" : "not found") << ")"
742                << endl;
743         }
744 #endif
745       }
746
747       // Store test results
748       if (wFace) {
749         vector<Vec3r> vertices;
750         for (int i = 0, numEdges = wFace->numberOfEdges(); i < numEdges; ++i) {
751           vertices.push_back(Vec3r(wFace->GetVertex(i)->GetVertex()));
752         }
753         Polygon3r poly(vertices, wFace->GetNormal());
754         poly.userdata = (void *)wFace;
755         fe->setaFace(poly);
756         wFaces.push_back(wFace);
757         fe->setOccludeeEmpty(false);
758 #if LOGGING
759         if (_global.debug & G_DEBUG_FREESTYLE) {
760           cout << "\tFound occludee" << endl;
761         }
762 #endif
763       }
764       else {
765         fe->setOccludeeEmpty(true);
766       }
767
768       ++nSamples;
769       fe = fe->nextEdge();
770     } while ((maxCard < qiMajority) && (fe) && (fe != festart));
771
772 #if LOGGING
773     if (_global.debug & G_DEBUG_FREESTYLE) {
774       cout << "\tFinished with " << nSamples << " samples, maxCard = " << maxCard << endl;
775     }
776 #endif
777
778     // ViewEdge
779     // qi --
780     (*ve)->setQI(maxIndex);
781     // occluders --
782     // I would rather not have to go through the effort of creating this this set and then copying out its contents.
783     // Is there a reason why ViewEdge::_Occluders cannot be converted to a set<>?
784     for (set<ViewShape *>::iterator o = foundOccluders.begin(), oend = foundOccluders.end();
785          o != oend;
786          ++o) {
787       (*ve)->AddOccluder((*o));
788     }
789 #if LOGGING
790     if (_global.debug & G_DEBUG_FREESTYLE) {
791       cout << "\tConclusion: QI = " << maxIndex << ", " << (*ve)->occluders_size() << " occluders."
792            << endl;
793     }
794 #endif
795     // occludee --
796     if (!wFaces.empty()) {
797       if (wFaces.size() <= (float)nSamples / 2.0f) {
798         (*ve)->setaShape(0);
799       }
800       else {
801         ViewShape *vshape = ioViewMap->viewShape(
802             (*wFaces.begin())->GetVertex(0)->shape()->GetId());
803         (*ve)->setaShape(vshape);
804       }
805     }
806
807     wFaces.clear();
808   }
809 }
810
811 template<typename G, typename I>
812 static void computeFastVisibility(ViewMap *ioViewMap, G &grid, real epsilon)
813 {
814   vector<ViewEdge *> &vedges = ioViewMap->ViewEdges();
815
816   FEdge *fe, *festart;
817   unsigned nSamples = 0;
818   vector<WFace *> wFaces;
819   WFace *wFace = NULL;
820   unsigned tmpQI = 0;
821   unsigned qiClasses[256];
822   unsigned maxIndex, maxCard;
823   unsigned qiMajority;
824   bool even_test;
825   for (vector<ViewEdge *>::iterator ve = vedges.begin(), veend = vedges.end(); ve != veend; ve++) {
826     // Find an edge to test
827     if (!(*ve)->isInImage()) {
828       // This view edge has been proscenium culled
829       (*ve)->setQI(255);
830       (*ve)->setaShape(0);
831       continue;
832     }
833
834     // Test edge
835     festart = (*ve)->fedgeA();
836     fe = (*ve)->fedgeA();
837
838     even_test = true;
839     qiMajority = 0;
840     do {
841       if (even_test && fe && fe->isInImage()) {
842         qiMajority++;
843         even_test = !even_test;
844       }
845       fe = fe->nextEdge();
846     } while (fe && fe != festart);
847
848     if (qiMajority == 0) {
849       // There are no occludable FEdges on this ViewEdge
850       // This should be impossible.
851       if (_global.debug & G_DEBUG_FREESTYLE) {
852         cout << "View Edge in viewport without occludable FEdges: " << (*ve)->getId() << endl;
853       }
854       // We can recover from this error:
855       // Treat this edge as fully visible with no occludee
856       (*ve)->setQI(0);
857       (*ve)->setaShape(0);
858       continue;
859     }
860     else {
861       ++qiMajority;
862       qiMajority >>= 1;
863     }
864
865     even_test = true;
866     maxIndex = 0;
867     maxCard = 0;
868     nSamples = 0;
869     memset(qiClasses, 0, 256 * sizeof(*qiClasses));
870     set<ViewShape *> foundOccluders;
871
872     fe = (*ve)->fedgeA();
873     do {
874       if (!fe || !fe->isInImage()) {
875         fe = fe->nextEdge();
876         continue;
877       }
878       if (even_test) {
879         if ((maxCard < qiMajority)) {
880           //ARB: change &wFace to wFace and use reference in called function
881           tmpQI = computeVisibility<G, I>(
882               ioViewMap, fe, grid, epsilon, *ve, &wFace, &foundOccluders);
883
884           //ARB: This is an error condition, not an alert condition.
885           // Some sort of recovery or abort is necessary.
886           if (tmpQI >= 256) {
887             cerr << "Warning: too many occluding levels" << endl;
888             //ARB: Wild guess: instead of aborting or corrupting memory, treat as tmpQI == 255
889             tmpQI = 255;
890           }
891
892           if (++qiClasses[tmpQI] > maxCard) {
893             maxCard = qiClasses[tmpQI];
894             maxIndex = tmpQI;
895           }
896         }
897         else {
898           //ARB: FindOccludee is redundant if ComputeRayCastingVisibility has been called
899           //ARB: change &wFace to wFace and use reference in called function
900           findOccludee<G, I>(fe, grid, epsilon, *ve, &wFace);
901         }
902
903         if (wFace) {
904           vector<Vec3r> vertices;
905           for (int i = 0, numEdges = wFace->numberOfEdges(); i < numEdges; ++i) {
906             vertices.push_back(Vec3r(wFace->GetVertex(i)->GetVertex()));
907           }
908           Polygon3r poly(vertices, wFace->GetNormal());
909           poly.userdata = (void *)wFace;
910           fe->setaFace(poly);
911           wFaces.push_back(wFace);
912         }
913         ++nSamples;
914       }
915
916       even_test = !even_test;
917       fe = fe->nextEdge();
918     } while ((maxCard < qiMajority) && (fe) && (fe != festart));
919
920     // qi --
921     (*ve)->setQI(maxIndex);
922
923     // occluders --
924     for (set<ViewShape *>::iterator o = foundOccluders.begin(), oend = foundOccluders.end();
925          o != oend;
926          ++o) {
927       (*ve)->AddOccluder((*o));
928     }
929
930     // occludee --
931     if (!wFaces.empty()) {
932       if (wFaces.size() < nSamples / 2) {
933         (*ve)->setaShape(0);
934       }
935       else {
936         ViewShape *vshape = ioViewMap->viewShape(
937             (*wFaces.begin())->GetVertex(0)->shape()->GetId());
938         (*ve)->setaShape(vshape);
939       }
940     }
941
942     wFaces.clear();
943   }
944 }
945
946 template<typename G, typename I>
947 static void computeVeryFastVisibility(ViewMap *ioViewMap, G &grid, real epsilon)
948 {
949   vector<ViewEdge *> &vedges = ioViewMap->ViewEdges();
950
951   FEdge *fe;
952   unsigned qi = 0;
953   WFace *wFace = 0;
954
955   for (vector<ViewEdge *>::iterator ve = vedges.begin(), veend = vedges.end(); ve != veend; ve++) {
956     // Find an edge to test
957     if (!(*ve)->isInImage()) {
958       // This view edge has been proscenium culled
959       (*ve)->setQI(255);
960       (*ve)->setaShape(0);
961       continue;
962     }
963     fe = (*ve)->fedgeA();
964     // Find a FEdge inside the occluder proscenium to test for visibility
965     FEdge *festart = fe;
966     while (fe && !fe->isInImage() && fe != festart) {
967       fe = fe->nextEdge();
968     }
969
970     // Test edge
971     if (!fe || !fe->isInImage()) {
972       // There are no occludable FEdges on this ViewEdge
973       // This should be impossible.
974       if (_global.debug & G_DEBUG_FREESTYLE) {
975         cout << "View Edge in viewport without occludable FEdges: " << (*ve)->getId() << endl;
976       }
977       // We can recover from this error:
978       // Treat this edge as fully visible with no occludee
979       qi = 0;
980       wFace = NULL;
981     }
982     else {
983       qi = computeVisibility<G, I>(ioViewMap, fe, grid, epsilon, *ve, &wFace, NULL);
984     }
985
986     // Store test results
987     if (wFace) {
988       vector<Vec3r> vertices;
989       for (int i = 0, numEdges = wFace->numberOfEdges(); i < numEdges; ++i) {
990         vertices.push_back(Vec3r(wFace->GetVertex(i)->GetVertex()));
991       }
992       Polygon3r poly(vertices, wFace->GetNormal());
993       poly.userdata = (void *)wFace;
994       fe->setaFace(poly);  // This works because setaFace *copies* the polygon
995       ViewShape *vshape = ioViewMap->viewShape(wFace->GetVertex(0)->shape()->GetId());
996       (*ve)->setaShape(vshape);
997     }
998     else {
999       (*ve)->setaShape(0);
1000     }
1001     (*ve)->setQI(qi);
1002   }
1003 }
1004
1005 void ViewMapBuilder::BuildGrid(WingedEdge &we, const BBox<Vec3r> &bbox, unsigned int sceneNumFaces)
1006 {
1007   _Grid->clear();
1008   Vec3r size;
1009   for (unsigned int i = 0; i < 3; i++) {
1010     size[i] = fabs(bbox.getMax()[i] - bbox.getMin()[i]);
1011     // let make the grid 1/10 bigger to avoid numerical errors while computing triangles/cells intersections.
1012     size[i] += size[i] / 10.0;
1013     if (size[i] == 0) {
1014       if (_global.debug & G_DEBUG_FREESTYLE) {
1015         cout << "Warning: the bbox size is 0 in dimension " << i << endl;
1016       }
1017     }
1018   }
1019   _Grid->configure(Vec3r(bbox.getMin() - size / 20.0), size, sceneNumFaces);
1020
1021   // Fill in the grid:
1022   WFillGrid fillGridRenderer(_Grid, &we);
1023   fillGridRenderer.fillGrid();
1024
1025   // DEBUG
1026   _Grid->displayDebug();
1027 }
1028
1029 ViewMap *ViewMapBuilder::BuildViewMap(WingedEdge &we,
1030                                       visibility_algo iAlgo,
1031                                       real epsilon,
1032                                       const BBox<Vec3r> &bbox,
1033                                       unsigned int sceneNumFaces)
1034 {
1035   _ViewMap = new ViewMap;
1036   _currentId = 1;
1037   _currentFId = 0;
1038   _currentSVertexId = 0;
1039
1040   // Builds initial view edges
1041   computeInitialViewEdges(we);
1042
1043   // Detects cusps
1044   computeCusps(_ViewMap);
1045
1046   // Compute intersections
1047   ComputeIntersections(_ViewMap, sweep_line, epsilon);
1048
1049   // Compute visibility
1050   ComputeEdgesVisibility(_ViewMap, we, bbox, sceneNumFaces, iAlgo, epsilon);
1051
1052   return _ViewMap;
1053 }
1054
1055 static inline real distance2D(const Vec3r &point, const real origin[2])
1056 {
1057   return ::hypot((point[0] - origin[0]), (point[1] - origin[1]));
1058 }
1059
1060 static inline bool crossesProscenium(real proscenium[4], FEdge *fe)
1061 {
1062   Vec2r min(proscenium[0], proscenium[2]);
1063   Vec2r max(proscenium[1], proscenium[3]);
1064   Vec2r A(fe->vertexA()->getProjectedX(), fe->vertexA()->getProjectedY());
1065   Vec2r B(fe->vertexB()->getProjectedX(), fe->vertexB()->getProjectedY());
1066
1067   return GeomUtils::intersect2dSeg2dArea(min, max, A, B);
1068 }
1069
1070 static inline bool insideProscenium(real proscenium[4], const Vec3r &point)
1071 {
1072   return !(point[0] < proscenium[0] || point[0] > proscenium[1] || point[1] < proscenium[2] ||
1073            point[1] > proscenium[3]);
1074 }
1075
1076 void ViewMapBuilder::CullViewEdges(ViewMap *ioViewMap,
1077                                    real viewProscenium[4],
1078                                    real occluderProscenium[4],
1079                                    bool extensiveFEdgeSearch)
1080 {
1081   // Cull view edges by marking them as non-displayable.
1082   // This avoids the complications of trying to delete edges from the ViewMap.
1083
1084   // Non-displayable view edges will be skipped over during visibility calculation.
1085
1086   // View edges will be culled according to their position w.r.t. the viewport proscenium (viewport + 5% border,
1087   // or some such).
1088
1089   // Get proscenium boundary for culling
1090   GridHelpers::getDefaultViewProscenium(viewProscenium);
1091   real prosceniumOrigin[2];
1092   prosceniumOrigin[0] = (viewProscenium[1] - viewProscenium[0]) / 2.0;
1093   prosceniumOrigin[1] = (viewProscenium[3] - viewProscenium[2]) / 2.0;
1094   if (_global.debug & G_DEBUG_FREESTYLE) {
1095     cout << "Proscenium culling:" << endl;
1096     cout << "Proscenium: [" << viewProscenium[0] << ", " << viewProscenium[1] << ", "
1097          << viewProscenium[2] << ", " << viewProscenium[3] << "]" << endl;
1098     cout << "Origin: [" << prosceniumOrigin[0] << ", " << prosceniumOrigin[1] << "]" << endl;
1099   }
1100
1101   // A separate occluder proscenium will also be maintained, starting out the same as the viewport proscenium, and
1102   // expanding as necessary so that it encompasses the center point of at least one feature edge in each retained view
1103   // edge.
1104   // The occluder proscenium will be used later to cull occluding triangles before they are inserted into the Grid.
1105   // The occluder proscenium starts out the same size as the view proscenium
1106   GridHelpers::getDefaultViewProscenium(occluderProscenium);
1107
1108   // N.B. Freestyle is inconsistent in its use of ViewMap::viewedges_container and vector<ViewEdge*>::iterator.
1109   // Probably all occurences of vector<ViewEdge*>::iterator should be replaced ViewMap::viewedges_container
1110   // throughout the code.
1111   // For each view edge
1112   ViewMap::viewedges_container::iterator ve, veend;
1113
1114   for (ve = ioViewMap->ViewEdges().begin(), veend = ioViewMap->ViewEdges().end(); ve != veend;
1115        ve++) {
1116     // Overview:
1117     //    Search for a visible feature edge
1118     //    If none: mark view edge as non-displayable
1119     //    Otherwise:
1120     //        Find a feature edge with center point inside occluder proscenium.
1121     //        If none exists, find the feature edge with center point closest to viewport origin.
1122     //            Expand occluder proscenium to enclose center point.
1123
1124     // For each feature edge, while bestOccluderTarget not found and view edge not visibile
1125     bool bestOccluderTargetFound = false;
1126     FEdge *bestOccluderTarget = NULL;
1127     real bestOccluderDistance = 0.0;
1128     FEdge *festart = (*ve)->fedgeA();
1129     FEdge *fe = festart;
1130     // All ViewEdges start culled
1131     (*ve)->setIsInImage(false);
1132
1133     // For simple visibility calculation: mark a feature edge that is known to have a center point inside the
1134     // occluder proscenium. Cull all other feature edges.
1135     do {
1136       // All FEdges start culled
1137       fe->setIsInImage(false);
1138
1139       // Look for the visible edge that can most easily be included in the occluder proscenium.
1140       if (!bestOccluderTargetFound) {
1141         // If center point is inside occluder proscenium,
1142         if (insideProscenium(occluderProscenium, fe->center2d())) {
1143           // Use this feature edge for visibility deterimination
1144           fe->setIsInImage(true);
1145           // Mark bestOccluderTarget as found
1146           bestOccluderTargetFound = true;
1147           bestOccluderTarget = fe;
1148         }
1149         else {
1150           real d = distance2D(fe->center2d(), prosceniumOrigin);
1151           // If center point is closer to viewport origin than current target
1152           if (bestOccluderTarget == NULL || d < bestOccluderDistance) {
1153             // Then store as bestOccluderTarget
1154             bestOccluderDistance = d;
1155             bestOccluderTarget = fe;
1156           }
1157         }
1158       }
1159
1160       // If feature edge crosses the view proscenium
1161       if (!(*ve)->isInImage() && crossesProscenium(viewProscenium, fe)) {
1162         // Then the view edge will be included in the image
1163         (*ve)->setIsInImage(true);
1164       }
1165       fe = fe->nextEdge();
1166     } while (fe && fe != festart && !(bestOccluderTargetFound && (*ve)->isInImage()));
1167
1168     // Either we have run out of FEdges, or we already have the one edge we need to determine visibility
1169     // Cull all remaining edges.
1170     while (fe && fe != festart) {
1171       fe->setIsInImage(false);
1172       fe = fe->nextEdge();
1173     }
1174
1175     // If bestOccluderTarget was not found inside the occluder proscenium, we need to expand the occluder
1176     // proscenium to include it.
1177     if ((*ve)->isInImage() && bestOccluderTarget != NULL && !bestOccluderTargetFound) {
1178       // Expand occluder proscenium to enclose bestOccluderTarget
1179       Vec3r point = bestOccluderTarget->center2d();
1180       if (point[0] < occluderProscenium[0]) {
1181         occluderProscenium[0] = point[0];
1182       }
1183       else if (point[0] > occluderProscenium[1]) {
1184         occluderProscenium[1] = point[0];
1185       }
1186       if (point[1] < occluderProscenium[2]) {
1187         occluderProscenium[2] = point[1];
1188       }
1189       else if (point[1] > occluderProscenium[3]) {
1190         occluderProscenium[3] = point[1];
1191       }
1192       // Use bestOccluderTarget for visibility determination
1193       bestOccluderTarget->setIsInImage(true);
1194     }
1195   }
1196
1197   // We are done calculating the occluder proscenium.
1198   // Expand the occluder proscenium by an epsilon to avoid rounding errors.
1199   const real epsilon = 1.0e-6;
1200   occluderProscenium[0] -= epsilon;
1201   occluderProscenium[1] += epsilon;
1202   occluderProscenium[2] -= epsilon;
1203   occluderProscenium[3] += epsilon;
1204
1205   // For "Normal" or "Fast" style visibility computation only:
1206
1207   // For more detailed visibility calculation, make a second pass through the view map, marking all feature edges
1208   // with center points inside the final occluder proscenium. All of these feature edges can be considered during
1209   // visibility calculation.
1210
1211   // So far we have only found one FEdge per ViewEdge.  The "Normal" and "Fast" styles of visibility computation
1212   // want to consider many FEdges for each ViewEdge.
1213   // Here we re-scan the view map to find any usable FEdges that we skipped on the first pass, or that have become
1214   // usable because the occluder proscenium has been expanded since the edge was visited on the first pass.
1215   if (extensiveFEdgeSearch) {
1216     // For each view edge,
1217     for (ve = ioViewMap->ViewEdges().begin(), veend = ioViewMap->ViewEdges().end(); ve != veend;
1218          ve++) {
1219       if (!(*ve)->isInImage()) {
1220         continue;
1221       }
1222       // For each feature edge,
1223       FEdge *festart = (*ve)->fedgeA();
1224       FEdge *fe = festart;
1225       do {
1226         // If not (already) visible and center point inside occluder proscenium,
1227         if (!fe->isInImage() && insideProscenium(occluderProscenium, fe->center2d())) {
1228           // Use the feature edge for visibility determination
1229           fe->setIsInImage(true);
1230         }
1231         fe = fe->nextEdge();
1232       } while (fe && fe != festart);
1233     }
1234   }
1235 }
1236
1237 void ViewMapBuilder::computeInitialViewEdges(WingedEdge &we)
1238 {
1239   vector<WShape *> wshapes = we.getWShapes();
1240   SShape *psShape;
1241
1242   for (vector<WShape *>::const_iterator it = wshapes.begin(); it != wshapes.end(); it++) {
1243     if (_pRenderMonitor && _pRenderMonitor->testBreak())
1244       break;
1245
1246     // create the embedding
1247     psShape = new SShape;
1248     psShape->setId((*it)->GetId());
1249     psShape->setName((*it)->getName());
1250     psShape->setLibraryPath((*it)->getLibraryPath());
1251     psShape->setFrsMaterials((*it)->frs_materials());  // FIXME
1252
1253     // create the view shape
1254     ViewShape *vshape = new ViewShape(psShape);
1255     // add this view shape to the view map:
1256     _ViewMap->AddViewShape(vshape);
1257
1258     // we want to number the view edges in a unique way for the while scene.
1259     _pViewEdgeBuilder->setCurrentViewId(_currentId);
1260     // we want to number the feature edges in a unique way for the while scene.
1261     _pViewEdgeBuilder->setCurrentFId(_currentFId);
1262     // we want to number the SVertex in a unique way for the while scene.
1263     _pViewEdgeBuilder->setCurrentSVertexId(_currentFId);
1264     _pViewEdgeBuilder->BuildViewEdges(dynamic_cast<WXShape *>(*it),
1265                                       vshape,
1266                                       _ViewMap->ViewEdges(),
1267                                       _ViewMap->ViewVertices(),
1268                                       _ViewMap->FEdges(),
1269                                       _ViewMap->SVertices());
1270
1271     _currentId = _pViewEdgeBuilder->currentViewId() + 1;
1272     _currentFId = _pViewEdgeBuilder->currentFId() + 1;
1273     _currentSVertexId = _pViewEdgeBuilder->currentSVertexId() + 1;
1274
1275     psShape->ComputeBBox();
1276   }
1277 }
1278
1279 void ViewMapBuilder::computeCusps(ViewMap *ioViewMap)
1280 {
1281   vector<ViewVertex *> newVVertices;
1282   vector<ViewEdge *> newVEdges;
1283   ViewMap::viewedges_container &vedges = ioViewMap->ViewEdges();
1284   ViewMap::viewedges_container::iterator ve = vedges.begin(), veend = vedges.end();
1285   for (; ve != veend; ++ve) {
1286     if (_pRenderMonitor && _pRenderMonitor->testBreak())
1287       break;
1288     if ((!((*ve)->getNature() & Nature::SILHOUETTE)) || (!((*ve)->fedgeA()->isSmooth())))
1289       continue;
1290     FEdge *fe = (*ve)->fedgeA();
1291     FEdge *fefirst = fe;
1292     bool first = true;
1293     bool positive = true;
1294     do {
1295       FEdgeSmooth *fes = dynamic_cast<FEdgeSmooth *>(fe);
1296       Vec3r A((fes)->vertexA()->point3d());
1297       Vec3r B((fes)->vertexB()->point3d());
1298       Vec3r AB(B - A);
1299       AB.normalize();
1300       Vec3r m((A + B) / 2.0);
1301       Vec3r crossP(AB ^ (fes)->normal());
1302       crossP.normalize();
1303       Vec3r viewvector;
1304       if (_orthographicProjection) {
1305         viewvector = Vec3r(0.0, 0.0, m.z() - _viewpoint.z());
1306       }
1307       else {
1308         viewvector = Vec3r(m - _viewpoint);
1309       }
1310       viewvector.normalize();
1311       if (first) {
1312         if (((crossP) * (viewvector)) > 0)
1313           positive = true;
1314         else
1315           positive = false;
1316         first = false;
1317       }
1318       // If we're in a positive part, we need a stronger negative value to change
1319       NonTVertex *cusp = NULL;
1320       if (positive) {
1321         if (((crossP) * (viewvector)) < -0.1) {
1322           // state changes
1323           positive = false;
1324           // creates and insert cusp
1325           cusp = dynamic_cast<NonTVertex *>(
1326               ioViewMap->InsertViewVertex(fes->vertexA(), newVEdges));
1327           if (cusp)
1328             cusp->setNature(cusp->getNature() | Nature::CUSP);
1329         }
1330       }
1331       else {
1332         // If we're in a negative part, we need a stronger negative value to change
1333         if (((crossP) * (viewvector)) > 0.1) {
1334           positive = true;
1335           cusp = dynamic_cast<NonTVertex *>(
1336               ioViewMap->InsertViewVertex(fes->vertexA(), newVEdges));
1337           if (cusp)
1338             cusp->setNature(cusp->getNature() | Nature::CUSP);
1339         }
1340       }
1341       fe = fe->nextEdge();
1342     } while (fe && fe != fefirst);
1343   }
1344   for (ve = newVEdges.begin(), veend = newVEdges.end(); ve != veend; ++ve) {
1345     (*ve)->viewShape()->AddEdge(*ve);
1346     vedges.push_back(*ve);
1347   }
1348 }
1349
1350 void ViewMapBuilder::ComputeCumulativeVisibility(ViewMap *ioViewMap,
1351                                                  WingedEdge &we,
1352                                                  const BBox<Vec3r> &bbox,
1353                                                  real epsilon,
1354                                                  bool cull,
1355                                                  GridDensityProviderFactory &factory)
1356 {
1357   AutoPtr<GridHelpers::Transform> transform;
1358   AutoPtr<OccluderSource> source;
1359
1360   if (_orthographicProjection) {
1361     transform.reset(new BoxGrid::Transform);
1362   }
1363   else {
1364     transform.reset(new SphericalGrid::Transform);
1365   }
1366
1367   if (cull) {
1368     source.reset(new CulledOccluderSource(*transform, we, *ioViewMap, true));
1369   }
1370   else {
1371     source.reset(new OccluderSource(*transform, we));
1372   }
1373
1374   AutoPtr<GridDensityProvider> density(factory.newGridDensityProvider(*source, bbox, *transform));
1375
1376   if (_orthographicProjection) {
1377     BoxGrid grid(*source, *density, ioViewMap, _viewpoint, _EnableQI);
1378     computeCumulativeVisibility<BoxGrid, BoxGrid::Iterator>(
1379         ioViewMap, grid, epsilon, _pRenderMonitor);
1380   }
1381   else {
1382     SphericalGrid grid(*source, *density, ioViewMap, _viewpoint, _EnableQI);
1383     computeCumulativeVisibility<SphericalGrid, SphericalGrid::Iterator>(
1384         ioViewMap, grid, epsilon, _pRenderMonitor);
1385   }
1386 }
1387
1388 void ViewMapBuilder::ComputeDetailedVisibility(ViewMap *ioViewMap,
1389                                                WingedEdge &we,
1390                                                const BBox<Vec3r> &bbox,
1391                                                real epsilon,
1392                                                bool cull,
1393                                                GridDensityProviderFactory &factory)
1394 {
1395   AutoPtr<GridHelpers::Transform> transform;
1396   AutoPtr<OccluderSource> source;
1397
1398   if (_orthographicProjection) {
1399     transform.reset(new BoxGrid::Transform);
1400   }
1401   else {
1402     transform.reset(new SphericalGrid::Transform);
1403   }
1404
1405   if (cull) {
1406     source.reset(new CulledOccluderSource(*transform, we, *ioViewMap, true));
1407   }
1408   else {
1409     source.reset(new OccluderSource(*transform, we));
1410   }
1411
1412   AutoPtr<GridDensityProvider> density(factory.newGridDensityProvider(*source, bbox, *transform));
1413
1414   if (_orthographicProjection) {
1415     BoxGrid grid(*source, *density, ioViewMap, _viewpoint, _EnableQI);
1416     computeDetailedVisibility<BoxGrid, BoxGrid::Iterator>(
1417         ioViewMap, grid, epsilon, _pRenderMonitor);
1418   }
1419   else {
1420     SphericalGrid grid(*source, *density, ioViewMap, _viewpoint, _EnableQI);
1421     computeDetailedVisibility<SphericalGrid, SphericalGrid::Iterator>(
1422         ioViewMap, grid, epsilon, _pRenderMonitor);
1423   }
1424 }
1425
1426 void ViewMapBuilder::ComputeEdgesVisibility(ViewMap *ioViewMap,
1427                                             WingedEdge &we,
1428                                             const BBox<Vec3r> &bbox,
1429                                             unsigned int sceneNumFaces,
1430                                             visibility_algo iAlgo,
1431                                             real epsilon)
1432 {
1433 #if 0
1434   iAlgo = ray_casting; // for testing algorithms equivalence
1435 #endif
1436   switch (iAlgo) {
1437     case ray_casting:
1438       if (_global.debug & G_DEBUG_FREESTYLE) {
1439         cout << "Using ordinary ray casting" << endl;
1440       }
1441       BuildGrid(we, bbox, sceneNumFaces);
1442       ComputeRayCastingVisibility(ioViewMap, epsilon);
1443       break;
1444     case ray_casting_fast:
1445       if (_global.debug & G_DEBUG_FREESTYLE) {
1446         cout << "Using fast ray casting" << endl;
1447       }
1448       BuildGrid(we, bbox, sceneNumFaces);
1449       ComputeFastRayCastingVisibility(ioViewMap, epsilon);
1450       break;
1451     case ray_casting_very_fast:
1452       if (_global.debug & G_DEBUG_FREESTYLE) {
1453         cout << "Using very fast ray casting" << endl;
1454       }
1455       BuildGrid(we, bbox, sceneNumFaces);
1456       ComputeVeryFastRayCastingVisibility(ioViewMap, epsilon);
1457       break;
1458     case ray_casting_culled_adaptive_traditional:
1459       if (_global.debug & G_DEBUG_FREESTYLE) {
1460         cout << "Using culled adaptive grid with heuristic density and traditional QI calculation"
1461              << endl;
1462       }
1463       try {
1464         HeuristicGridDensityProviderFactory factory(0.5f, sceneNumFaces);
1465         ComputeDetailedVisibility(ioViewMap, we, bbox, epsilon, true, factory);
1466       }
1467       catch (...) {
1468         // Last resort catch to make sure RAII semantics hold for OptimizedGrid. Can be replaced with
1469         // try...catch block around main() if the program as a whole is converted to RAII
1470
1471         // This is the little-mentioned caveat of RAII: RAII does not work unless destructors are always
1472         // called, but destructors are only called if all exceptions are caught (or std::terminate() is
1473         // replaced).
1474
1475         // We don't actually handle the exception here, so re-throw it now that our destructors have had a
1476         // chance to run.
1477         throw;
1478       }
1479       break;
1480     case ray_casting_adaptive_traditional:
1481       if (_global.debug & G_DEBUG_FREESTYLE) {
1482         cout
1483             << "Using unculled adaptive grid with heuristic density and traditional QI calculation"
1484             << endl;
1485       }
1486       try {
1487         HeuristicGridDensityProviderFactory factory(0.5f, sceneNumFaces);
1488         ComputeDetailedVisibility(ioViewMap, we, bbox, epsilon, false, factory);
1489       }
1490       catch (...) {
1491         throw;
1492       }
1493       break;
1494     case ray_casting_culled_adaptive_cumulative:
1495       if (_global.debug & G_DEBUG_FREESTYLE) {
1496         cout << "Using culled adaptive grid with heuristic density and cumulative QI calculation"
1497              << endl;
1498       }
1499       try {
1500         HeuristicGridDensityProviderFactory factory(0.5f, sceneNumFaces);
1501         ComputeCumulativeVisibility(ioViewMap, we, bbox, epsilon, true, factory);
1502       }
1503       catch (...) {
1504         throw;
1505       }
1506       break;
1507     case ray_casting_adaptive_cumulative:
1508       if (_global.debug & G_DEBUG_FREESTYLE) {
1509         cout << "Using unculled adaptive grid with heuristic density and cumulative QI calculation"
1510              << endl;
1511       }
1512       try {
1513         HeuristicGridDensityProviderFactory factory(0.5f, sceneNumFaces);
1514         ComputeCumulativeVisibility(ioViewMap, we, bbox, epsilon, false, factory);
1515       }
1516       catch (...) {
1517         throw;
1518       }
1519       break;
1520     default:
1521       break;
1522   }
1523 }
1524
1525 static const unsigned gProgressBarMaxSteps = 10;
1526 static const unsigned gProgressBarMinSize = 2000;
1527
1528 void ViewMapBuilder::ComputeRayCastingVisibility(ViewMap *ioViewMap, real epsilon)
1529 {
1530   vector<ViewEdge *> &vedges = ioViewMap->ViewEdges();
1531   bool progressBarDisplay = false;
1532   unsigned progressBarStep = 0;
1533   unsigned vEdgesSize = vedges.size();
1534   unsigned fEdgesSize = ioViewMap->FEdges().size();
1535
1536   if (_pProgressBar != NULL && fEdgesSize > gProgressBarMinSize) {
1537     unsigned progressBarSteps = min(gProgressBarMaxSteps, vEdgesSize);
1538     progressBarStep = vEdgesSize / progressBarSteps;
1539     _pProgressBar->reset();
1540     _pProgressBar->setLabelText("Computing Ray casting Visibility");
1541     _pProgressBar->setTotalSteps(progressBarSteps);
1542     _pProgressBar->setProgress(0);
1543     progressBarDisplay = true;
1544   }
1545
1546   unsigned counter = progressBarStep;
1547   FEdge *fe, *festart;
1548   int nSamples = 0;
1549   vector<Polygon3r *> aFaces;
1550   Polygon3r *aFace = NULL;
1551   unsigned tmpQI = 0;
1552   unsigned qiClasses[256];
1553   unsigned maxIndex, maxCard;
1554   unsigned qiMajority;
1555   static unsigned timestamp = 1;
1556   for (vector<ViewEdge *>::iterator ve = vedges.begin(), veend = vedges.end(); ve != veend; ve++) {
1557     if (_pRenderMonitor && _pRenderMonitor->testBreak())
1558       break;
1559 #if LOGGING
1560     if (_global.debug & G_DEBUG_FREESTYLE) {
1561       cout << "Processing ViewEdge " << (*ve)->getId() << endl;
1562     }
1563 #endif
1564     festart = (*ve)->fedgeA();
1565     fe = (*ve)->fedgeA();
1566     qiMajority = 1;
1567     do {
1568       qiMajority++;
1569       fe = fe->nextEdge();
1570     } while (fe && fe != festart);
1571     qiMajority >>= 1;
1572 #if LOGGING
1573     if (_global.debug & G_DEBUG_FREESTYLE) {
1574       cout << "\tqiMajority: " << qiMajority << endl;
1575     }
1576 #endif
1577
1578     tmpQI = 0;
1579     maxIndex = 0;
1580     maxCard = 0;
1581     nSamples = 0;
1582     fe = (*ve)->fedgeA();
1583     memset(qiClasses, 0, 256 * sizeof(*qiClasses));
1584     set<ViewShape *> occluders;
1585     do {
1586       if ((maxCard < qiMajority)) {
1587         tmpQI = ComputeRayCastingVisibility(fe, _Grid, epsilon, occluders, &aFace, timestamp++);
1588
1589 #if LOGGING
1590         if (_global.debug & G_DEBUG_FREESTYLE) {
1591           cout << "\tFEdge: visibility " << tmpQI << endl;
1592         }
1593 #endif
1594         //ARB: This is an error condition, not an alert condition.
1595         // Some sort of recovery or abort is necessary.
1596         if (tmpQI >= 256) {
1597           cerr << "Warning: too many occluding levels" << endl;
1598           //ARB: Wild guess: instead of aborting or corrupting memory, treat as tmpQI == 255
1599           tmpQI = 255;
1600         }
1601
1602         if (++qiClasses[tmpQI] > maxCard) {
1603           maxCard = qiClasses[tmpQI];
1604           maxIndex = tmpQI;
1605         }
1606       }
1607       else {
1608         //ARB: FindOccludee is redundant if ComputeRayCastingVisibility has been called
1609         FindOccludee(fe, _Grid, epsilon, &aFace, timestamp++);
1610 #if LOGGING
1611         if (_global.debug & G_DEBUG_FREESTYLE) {
1612           cout << "\tFEdge: occludee only (" << (aFace != NULL ? "found" : "not found") << ")"
1613                << endl;
1614         }
1615 #endif
1616       }
1617
1618       if (aFace) {
1619         fe->setaFace(*aFace);
1620         aFaces.push_back(aFace);
1621         fe->setOccludeeEmpty(false);
1622 #if LOGGING
1623         if (_global.debug & G_DEBUG_FREESTYLE) {
1624           cout << "\tFound occludee" << endl;
1625         }
1626 #endif
1627       }
1628       else {
1629         //ARB: We are arbitrarily using the last observed value for occludee (almost always the value observed
1630         //     for the edge before festart). Is that meaningful?
1631         // ...in fact, _occludeeEmpty seems to be unused.
1632         fe->setOccludeeEmpty(true);
1633       }
1634
1635       ++nSamples;
1636       fe = fe->nextEdge();
1637     } while ((maxCard < qiMajority) && (fe) && (fe != festart));
1638 #if LOGGING
1639     if (_global.debug & G_DEBUG_FREESTYLE) {
1640       cout << "\tFinished with " << nSamples << " samples, maxCard = " << maxCard << endl;
1641     }
1642 #endif
1643
1644     // ViewEdge
1645     // qi --
1646     (*ve)->setQI(maxIndex);
1647     // occluders --
1648     for (set<ViewShape *>::iterator o = occluders.begin(), oend = occluders.end(); o != oend; ++o)
1649       (*ve)->AddOccluder((*o));
1650 #if LOGGING
1651     if (_global.debug & G_DEBUG_FREESTYLE) {
1652       cout << "\tConclusion: QI = " << maxIndex << ", " << (*ve)->occluders_size() << " occluders."
1653            << endl;
1654     }
1655 #endif
1656     // occludee --
1657     if (!aFaces.empty()) {
1658       if (aFaces.size() <= (float)nSamples / 2.0f) {
1659         (*ve)->setaShape(0);
1660       }
1661       else {
1662         vector<Polygon3r *>::iterator p = aFaces.begin();
1663         WFace *wface = (WFace *)((*p)->userdata);
1664         ViewShape *vshape = ioViewMap->viewShape(wface->GetVertex(0)->shape()->GetId());
1665         ++p;
1666         (*ve)->setaShape(vshape);
1667       }
1668     }
1669
1670     if (progressBarDisplay) {
1671       counter--;
1672       if (counter <= 0) {
1673         counter = progressBarStep;
1674         _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
1675       }
1676     }
1677     aFaces.clear();
1678   }
1679 }
1680
1681 void ViewMapBuilder::ComputeFastRayCastingVisibility(ViewMap *ioViewMap, real epsilon)
1682 {
1683   vector<ViewEdge *> &vedges = ioViewMap->ViewEdges();
1684   bool progressBarDisplay = false;
1685   unsigned progressBarStep = 0;
1686   unsigned vEdgesSize = vedges.size();
1687   unsigned fEdgesSize = ioViewMap->FEdges().size();
1688
1689   if (_pProgressBar != NULL && fEdgesSize > gProgressBarMinSize) {
1690     unsigned progressBarSteps = min(gProgressBarMaxSteps, vEdgesSize);
1691     progressBarStep = vEdgesSize / progressBarSteps;
1692     _pProgressBar->reset();
1693     _pProgressBar->setLabelText("Computing Ray casting Visibility");
1694     _pProgressBar->setTotalSteps(progressBarSteps);
1695     _pProgressBar->setProgress(0);
1696     progressBarDisplay = true;
1697   }
1698
1699   unsigned counter = progressBarStep;
1700   FEdge *fe, *festart;
1701   unsigned nSamples = 0;
1702   vector<Polygon3r *> aFaces;
1703   Polygon3r *aFace = NULL;
1704   unsigned tmpQI = 0;
1705   unsigned qiClasses[256];
1706   unsigned maxIndex, maxCard;
1707   unsigned qiMajority;
1708   static unsigned timestamp = 1;
1709   bool even_test;
1710   for (vector<ViewEdge *>::iterator ve = vedges.begin(), veend = vedges.end(); ve != veend; ve++) {
1711     if (_pRenderMonitor && _pRenderMonitor->testBreak())
1712       break;
1713
1714     festart = (*ve)->fedgeA();
1715     fe = (*ve)->fedgeA();
1716     qiMajority = 1;
1717     do {
1718       qiMajority++;
1719       fe = fe->nextEdge();
1720     } while (fe && fe != festart);
1721     if (qiMajority >= 4)
1722       qiMajority >>= 2;
1723     else
1724       qiMajority = 1;
1725
1726     set<ViewShape *> occluders;
1727
1728     even_test = true;
1729     maxIndex = 0;
1730     maxCard = 0;
1731     nSamples = 0;
1732     memset(qiClasses, 0, 256 * sizeof(*qiClasses));
1733     fe = (*ve)->fedgeA();
1734     do {
1735       if (even_test) {
1736         if ((maxCard < qiMajority)) {
1737           tmpQI = ComputeRayCastingVisibility(fe, _Grid, epsilon, occluders, &aFace, timestamp++);
1738
1739           //ARB: This is an error condition, not an alert condition.
1740           // Some sort of recovery or abort is necessary.
1741           if (tmpQI >= 256) {
1742             cerr << "Warning: too many occluding levels" << endl;
1743             //ARB: Wild guess: instead of aborting or corrupting memory, treat as tmpQI == 255
1744             tmpQI = 255;
1745           }
1746
1747           if (++qiClasses[tmpQI] > maxCard) {
1748             maxCard = qiClasses[tmpQI];
1749             maxIndex = tmpQI;
1750           }
1751         }
1752         else {
1753           //ARB: FindOccludee is redundant if ComputeRayCastingVisibility has been called
1754           FindOccludee(fe, _Grid, epsilon, &aFace, timestamp++);
1755         }
1756
1757         if (aFace) {
1758           fe->setaFace(*aFace);
1759           aFaces.push_back(aFace);
1760         }
1761         ++nSamples;
1762         even_test = false;
1763       }
1764       else {
1765         even_test = true;
1766       }
1767       fe = fe->nextEdge();
1768     } while ((maxCard < qiMajority) && (fe) && (fe != festart));
1769
1770     (*ve)->setQI(maxIndex);
1771
1772     if (!aFaces.empty()) {
1773       if (aFaces.size() < nSamples / 2) {
1774         (*ve)->setaShape(0);
1775       }
1776       else {
1777         vector<Polygon3r *>::iterator p = aFaces.begin();
1778         WFace *wface = (WFace *)((*p)->userdata);
1779         ViewShape *vshape = ioViewMap->viewShape(wface->GetVertex(0)->shape()->GetId());
1780         ++p;
1781 #if 0
1782         for (; p != pend; ++p) {
1783           WFace *f = (WFace*)((*p)->userdata);
1784           ViewShape *vs = ioViewMap->viewShape(f->GetVertex(0)->shape()->GetId());
1785           if (vs != vshape) {
1786             sameShape = false;
1787             break;
1788           }
1789         }
1790         if (sameShape)
1791 #endif
1792         (*ve)->setaShape(vshape);
1793       }
1794     }
1795
1796     //(*ve)->setaFace(aFace);
1797
1798     if (progressBarDisplay) {
1799       counter--;
1800       if (counter <= 0) {
1801         counter = progressBarStep;
1802         _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
1803       }
1804     }
1805     aFaces.clear();
1806   }
1807 }
1808
1809 void ViewMapBuilder::ComputeVeryFastRayCastingVisibility(ViewMap *ioViewMap, real epsilon)
1810 {
1811   vector<ViewEdge *> &vedges = ioViewMap->ViewEdges();
1812   bool progressBarDisplay = false;
1813   unsigned progressBarStep = 0;
1814   unsigned vEdgesSize = vedges.size();
1815   unsigned fEdgesSize = ioViewMap->FEdges().size();
1816
1817   if (_pProgressBar != NULL && fEdgesSize > gProgressBarMinSize) {
1818     unsigned progressBarSteps = min(gProgressBarMaxSteps, vEdgesSize);
1819     progressBarStep = vEdgesSize / progressBarSteps;
1820     _pProgressBar->reset();
1821     _pProgressBar->setLabelText("Computing Ray casting Visibility");
1822     _pProgressBar->setTotalSteps(progressBarSteps);
1823     _pProgressBar->setProgress(0);
1824     progressBarDisplay = true;
1825   }
1826
1827   unsigned counter = progressBarStep;
1828   FEdge *fe;
1829   unsigned qi = 0;
1830   Polygon3r *aFace = NULL;
1831   static unsigned timestamp = 1;
1832   for (vector<ViewEdge *>::iterator ve = vedges.begin(), veend = vedges.end(); ve != veend; ve++) {
1833     if (_pRenderMonitor && _pRenderMonitor->testBreak())
1834       break;
1835
1836     set<ViewShape *> occluders;
1837
1838     fe = (*ve)->fedgeA();
1839     qi = ComputeRayCastingVisibility(fe, _Grid, epsilon, occluders, &aFace, timestamp++);
1840     if (aFace) {
1841       fe->setaFace(*aFace);
1842       WFace *wface = (WFace *)(aFace->userdata);
1843       ViewShape *vshape = ioViewMap->viewShape(wface->GetVertex(0)->shape()->GetId());
1844       (*ve)->setaShape(vshape);
1845     }
1846     else {
1847       (*ve)->setaShape(0);
1848     }
1849
1850     (*ve)->setQI(qi);
1851
1852     if (progressBarDisplay) {
1853       counter--;
1854       if (counter <= 0) {
1855         counter = progressBarStep;
1856         _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
1857       }
1858     }
1859   }
1860 }
1861
1862 void ViewMapBuilder::FindOccludee(FEdge *fe,
1863                                   Grid *iGrid,
1864                                   real epsilon,
1865                                   Polygon3r **oaPolygon,
1866                                   unsigned timestamp,
1867                                   Vec3r &u,
1868                                   Vec3r &A,
1869                                   Vec3r &origin,
1870                                   Vec3r &edgeDir,
1871                                   vector<WVertex *> &faceVertices)
1872 {
1873   WFace *face = NULL;
1874   if (fe->isSmooth()) {
1875     FEdgeSmooth *fes = dynamic_cast<FEdgeSmooth *>(fe);
1876     face = (WFace *)fes->face();
1877   }
1878   OccludersSet occluders;
1879   WFace *oface;
1880   bool skipFace;
1881
1882   WVertex::incoming_edge_iterator ie;
1883   OccludersSet::iterator p, pend;
1884
1885   *oaPolygon = NULL;
1886   if (((fe)->getNature() & Nature::SILHOUETTE) || ((fe)->getNature() & Nature::BORDER)) {
1887     occluders.clear();
1888     // we cast a ray from A in the same direction but looking behind
1889     Vec3r v(-u[0], -u[1], -u[2]);
1890     iGrid->castInfiniteRay(A, v, occluders, timestamp);
1891
1892     bool noIntersection = true;
1893     real mint = FLT_MAX;
1894     // we met some occluders, let us fill the aShape field with the first intersected occluder
1895     for (p = occluders.begin(), pend = occluders.end(); p != pend; p++) {
1896       // check whether the edge and the polygon plane are coincident:
1897       //-------------------------------------------------------------
1898       //first let us compute the plane equation.
1899       oface = (WFace *)(*p)->userdata;
1900       Vec3r v1(((*p)->getVertices())[0]);
1901       Vec3r normal((*p)->getNormal());
1902       real d = -(v1 * normal);
1903       real t, t_u, t_v;
1904
1905       if (face) {
1906         skipFace = false;
1907
1908         if (face == oface)
1909           continue;
1910
1911         if (faceVertices.empty())
1912           continue;
1913
1914         for (vector<WVertex *>::iterator fv = faceVertices.begin(), fvend = faceVertices.end();
1915              fv != fvend;
1916              ++fv) {
1917           if ((*fv)->isBoundary())
1918             continue;
1919           WVertex::incoming_edge_iterator iebegin = (*fv)->incoming_edges_begin();
1920           WVertex::incoming_edge_iterator ieend = (*fv)->incoming_edges_end();
1921           for (ie = iebegin; ie != ieend; ++ie) {
1922             if ((*ie) == 0)
1923               continue;
1924
1925             WFace *sface = (*ie)->GetbFace();
1926             if (sface == oface) {
1927               skipFace = true;
1928               break;
1929             }
1930           }
1931           if (skipFace)
1932             break;
1933         }
1934         if (skipFace)
1935           continue;
1936       }
1937       else {
1938         if (GeomUtils::COINCIDENT ==
1939             GeomUtils::intersectRayPlane(origin, edgeDir, normal, d, t, epsilon))
1940           continue;
1941       }
1942       if ((*p)->rayIntersect(A, v, t, t_u, t_v)) {
1943         if (fabs(v * normal) > 0.0001) {
1944           if (t > 0.0) {  // && t < 1.0) {
1945             if (t < mint) {
1946               *oaPolygon = (*p);
1947               mint = t;
1948               noIntersection = false;
1949               fe->setOccludeeIntersection(Vec3r(A + t * v));
1950             }
1951           }
1952         }
1953       }
1954     }
1955
1956     if (noIntersection)
1957       *oaPolygon = NULL;
1958   }
1959 }
1960
1961 void ViewMapBuilder::FindOccludee(
1962     FEdge *fe, Grid *iGrid, real epsilon, Polygon3r **oaPolygon, unsigned timestamp)
1963 {
1964   OccludersSet occluders;
1965
1966   Vec3r A;
1967   Vec3r edgeDir;
1968   Vec3r origin;
1969   A = Vec3r(((fe)->vertexA()->point3D() + (fe)->vertexB()->point3D()) / 2.0);
1970   edgeDir = Vec3r((fe)->vertexB()->point3D() - (fe)->vertexA()->point3D());
1971   edgeDir.normalize();
1972   origin = Vec3r((fe)->vertexA()->point3D());
1973   Vec3r u;
1974   if (_orthographicProjection) {
1975     u = Vec3r(0.0, 0.0, _viewpoint.z() - A.z());
1976   }
1977   else {
1978     u = Vec3r(_viewpoint - A);
1979   }
1980   u.normalize();
1981   if (A < iGrid->getOrigin())
1982     cerr << "Warning: point is out of the grid for fedge " << fe->getId().getFirst() << "-"
1983          << fe->getId().getSecond() << endl;
1984
1985   vector<WVertex *> faceVertices;
1986
1987   WFace *face = NULL;
1988   if (fe->isSmooth()) {
1989     FEdgeSmooth *fes = dynamic_cast<FEdgeSmooth *>(fe);
1990     face = (WFace *)fes->face();
1991   }
1992   if (face)
1993     face->RetrieveVertexList(faceVertices);
1994
1995   return FindOccludee(
1996       fe, iGrid, epsilon, oaPolygon, timestamp, u, A, origin, edgeDir, faceVertices);
1997 }
1998
1999 int ViewMapBuilder::ComputeRayCastingVisibility(FEdge *fe,
2000                                                 Grid *iGrid,
2001                                                 real epsilon,
2002                                                 set<ViewShape *> &oOccluders,
2003                                                 Polygon3r **oaPolygon,
2004                                                 unsigned timestamp)
2005 {
2006   OccludersSet occluders;
2007   int qi = 0;
2008
2009   Vec3r center;
2010   Vec3r edgeDir;
2011   Vec3r origin;
2012
2013   center = fe->center3d();
2014   edgeDir = Vec3r(fe->vertexB()->point3D() - fe->vertexA()->point3D());
2015   edgeDir.normalize();
2016   origin = Vec3r(fe->vertexA()->point3D());
2017   // Is the edge outside the view frustum ?
2018   Vec3r gridOrigin(iGrid->getOrigin());
2019   Vec3r gridExtremity(iGrid->getOrigin() + iGrid->gridSize());
2020
2021   if ((center.x() < gridOrigin.x()) || (center.y() < gridOrigin.y()) ||
2022       (center.z() < gridOrigin.z()) || (center.x() > gridExtremity.x()) ||
2023       (center.y() > gridExtremity.y()) || (center.z() > gridExtremity.z())) {
2024     cerr << "Warning: point is out of the grid for fedge " << fe->getId() << endl;
2025     //return 0;
2026   }
2027
2028 #if 0
2029   Vec3r A(fe->vertexA()->point2d());
2030   Vec3r B(fe->vertexB()->point2d());
2031   int viewport[4];
2032   SilhouetteGeomEngine::retrieveViewport(viewport);
2033   if ((A.x() < viewport[0]) || (A.x() > viewport[2]) || (A.y() < viewport[1]) || (A.y() > viewport[3]) ||
2034       (B.x() < viewport[0]) || (B.x() > viewport[2]) || (B.y() < viewport[1]) || (B.y() > viewport[3]))
2035   {
2036     cerr << "Warning: point is out of the grid for fedge " << fe->getId() << endl;
2037     //return 0;
2038   }
2039 #endif
2040
2041   Vec3r vp;
2042   if (_orthographicProjection) {
2043     vp = Vec3r(center.x(), center.y(), _viewpoint.z());
2044   }
2045   else {
2046     vp = Vec3r(_viewpoint);
2047   }
2048   Vec3r u(vp - center);
2049   real raylength = u.norm();
2050   u.normalize();
2051 #if 0
2052   if (_global.debug & G_DEBUG_FREESTYLE) {
2053     cout << "grid origin " << iGrid->getOrigin().x() << "," << iGrid->getOrigin().y() << ","
2054          << iGrid->getOrigin().z() << endl;
2055     cout << "center " << center.x() << "," << center.y() << "," << center.z() << endl;
2056   }
2057 #endif
2058
2059   iGrid->castRay(center, vp, occluders, timestamp);
2060
2061   WFace *face = NULL;
2062   if (fe->isSmooth()) {
2063     FEdgeSmooth *fes = dynamic_cast<FEdgeSmooth *>(fe);
2064     face = (WFace *)fes->face();
2065   }
2066   vector<WVertex *> faceVertices;
2067   WVertex::incoming_edge_iterator ie;
2068
2069   WFace *oface;
2070   bool skipFace;
2071   OccludersSet::iterator p, pend;
2072   if (face)
2073     face->RetrieveVertexList(faceVertices);
2074
2075   for (p = occluders.begin(), pend = occluders.end(); p != pend; p++) {
2076     // If we're dealing with an exact silhouette, check whether we must take care of this occluder of not.
2077     // (Indeed, we don't consider the occluders that share at least one vertex with the face containing this edge).
2078     //-----------
2079     oface = (WFace *)(*p)->userdata;
2080 #if LOGGING
2081     if (_global.debug & G_DEBUG_FREESTYLE) {
2082       cout << "\t\tEvaluating intersection for occluder " << ((*p)->getVertices())[0]
2083            << ((*p)->getVertices())[1] << ((*p)->getVertices())[2] << endl
2084            << "\t\t\tand ray " << vp << " * " << u << " (center " << center << ")" << endl;
2085     }
2086 #endif
2087     Vec3r v1(((*p)->getVertices())[0]);
2088     Vec3r normal((*p)->getNormal());
2089     real d = -(v1 * normal);
2090     real t, t_u, t_v;
2091
2092 #if LOGGING
2093     if (_global.debug & G_DEBUG_FREESTYLE) {
2094       cout << "\t\tp:  " << ((*p)->getVertices())[0] << ((*p)->getVertices())[1]
2095            << ((*p)->getVertices())[2] << ", norm: " << (*p)->getNormal() << endl;
2096     }
2097 #endif
2098
2099     if (face) {
2100 #if LOGGING
2101       if (_global.debug & G_DEBUG_FREESTYLE) {
2102         cout << "\t\tDetermining face adjacency...";
2103       }
2104 #endif
2105       skipFace = false;
2106
2107       if (face == oface) {
2108 #if LOGGING
2109         if (_global.debug & G_DEBUG_FREESTYLE) {
2110           cout << "  Rejecting occluder for face concurrency." << endl;
2111         }
2112 #endif
2113         continue;
2114       }
2115
2116       for (vector<WVertex *>::iterator fv = faceVertices.begin(), fvend = faceVertices.end();
2117            fv != fvend;
2118            ++fv) {
2119         if ((*fv)->isBoundary())
2120           continue;
2121
2122         WVertex::incoming_edge_iterator iebegin = (*fv)->incoming_edges_begin();
2123         WVertex::incoming_edge_iterator ieend = (*fv)->incoming_edges_end();
2124         for (ie = iebegin; ie != ieend; ++ie) {
2125           if ((*ie) == 0)
2126             continue;
2127
2128           WFace *sface = (*ie)->GetbFace();
2129           //WFace *sfacea = (*ie)->GetaFace();
2130           //if ((sface == oface) || (sfacea == oface)) {
2131           if (sface == oface) {
2132             skipFace = true;
2133             break;
2134           }
2135         }
2136         if (skipFace)
2137           break;
2138       }
2139       if (skipFace) {
2140 #if LOGGING
2141         if (_global.debug & G_DEBUG_FREESTYLE) {
2142           cout << "  Rejecting occluder for face adjacency." << endl;
2143         }
2144 #endif
2145         continue;
2146       }
2147     }
2148     else {
2149       // check whether the edge and the polygon plane are coincident:
2150       //-------------------------------------------------------------
2151       //first let us compute the plane equation.
2152
2153       if (GeomUtils::COINCIDENT ==
2154           GeomUtils::intersectRayPlane(origin, edgeDir, normal, d, t, epsilon)) {
2155 #if LOGGING
2156         if (_global.debug & G_DEBUG_FREESTYLE) {
2157           cout << "\t\tRejecting occluder for target coincidence." << endl;
2158         }
2159 #endif
2160         continue;
2161       }
2162     }
2163
2164     if ((*p)->rayIntersect(center, u, t, t_u, t_v)) {
2165 #if LOGGING
2166       if (_global.debug & G_DEBUG_FREESTYLE) {
2167         cout << "\t\tRay " << vp << " * " << u << " intersects at time " << t << " (raylength is "
2168              << raylength << ")" << endl;
2169         cout << "\t\t(u * normal) == " << (u * normal) << " for normal " << normal << endl;
2170       }
2171 #endif
2172       if (fabs(u * normal) > 0.0001) {
2173         if ((t > 0.0) && (t < raylength)) {
2174 #if LOGGING
2175           if (_global.debug & G_DEBUG_FREESTYLE) {
2176             cout << "\t\tIs occluder" << endl;
2177           }
2178 #endif
2179           WFace *f = (WFace *)((*p)->userdata);
2180           ViewShape *vshape = _ViewMap->viewShape(f->GetVertex(0)->shape()->GetId());
2181           oOccluders.insert(vshape);
2182           ++qi;
2183           if (!_EnableQI)
2184             break;
2185         }
2186       }
2187     }
2188   }
2189
2190   // Find occludee
2191   FindOccludee(fe, iGrid, epsilon, oaPolygon, timestamp, u, center, origin, edgeDir, faceVertices);
2192
2193   return qi;
2194 }
2195
2196 void ViewMapBuilder::ComputeIntersections(ViewMap *ioViewMap,
2197                                           intersection_algo iAlgo,
2198                                           real epsilon)
2199 {
2200   switch (iAlgo) {
2201     case sweep_line:
2202       ComputeSweepLineIntersections(ioViewMap, epsilon);
2203       break;
2204     default:
2205       break;
2206   }
2207 #if 0
2208   if (_global.debug & G_DEBUG_FREESTYLE) {
2209     ViewMap::viewvertices_container& vvertices = ioViewMap->ViewVertices();
2210     for (ViewMap::viewvertices_container::iterator vv = vvertices.begin(), vvend = vvertices.end();
2211          vv != vvend; ++vv)
2212     {
2213       if ((*vv)->getNature() == Nature::T_VERTEX) {
2214         TVertex *tvertex = (TVertex *)(*vv);
2215         cout << "TVertex " << tvertex->getId() << " has :" << endl;
2216         cout << "FrontEdgeA: " << tvertex->frontEdgeA().first << endl;
2217         cout << "FrontEdgeB: " << tvertex->frontEdgeB().first << endl;
2218         cout << "BackEdgeA: " << tvertex->backEdgeA().first << endl;
2219         cout << "BackEdgeB: " << tvertex->backEdgeB().first << endl << endl;
2220       }
2221     }
2222   }
2223 #endif
2224 }
2225
2226 struct less_SVertex2D : public binary_function<SVertex *, SVertex *, bool> {
2227   real epsilon;
2228
2229   less_SVertex2D(real eps) : binary_function<SVertex *, SVertex *, bool>()
2230   {
2231     epsilon = eps;
2232   }
2233
2234   bool operator()(SVertex *x, SVertex *y)
2235   {
2236     Vec3r A = x->point2D();
2237     Vec3r B = y->point2D();
2238     for (unsigned int i = 0; i < 3; i++) {
2239       if ((fabs(A[i] - B[i])) < epsilon)
2240         continue;
2241       if (A[i] < B[i])
2242         return true;
2243       if (A[i] > B[i])
2244         return false;
2245     }
2246     return false;
2247   }
2248 };
2249
2250 typedef Segment<FEdge *, Vec3r> segment;
2251 typedef Intersection<segment> intersection;
2252
2253 struct less_Intersection : public binary_function<intersection *, intersection *, bool> {
2254   segment *edge;
2255
2256   less_Intersection(segment *iEdge) : binary_function<intersection *, intersection *, bool>()
2257   {
2258     edge = iEdge;
2259   }
2260
2261   bool operator()(intersection *x, intersection *y)
2262   {
2263     real tx = x->getParameter(edge);
2264     real ty = y->getParameter(edge);
2265     if (tx > ty)
2266       return true;
2267     return false;
2268   }
2269 };
2270
2271 struct silhouette_binary_rule : public binary_rule<segment, segment> {
2272   silhouette_binary_rule() : binary_rule<segment, segment>()
2273   {
2274   }
2275
2276   virtual bool operator()(segment &s1, segment &s2)
2277   {
2278     FEdge *f1 = s1.edge();
2279     FEdge *f2 = s2.edge();
2280
2281     if ((!(((f1)->getNature() & Nature::SILHOUETTE) || ((f1)->getNature() & Nature::BORDER))) &&
2282         (!(((f2)->getNature() & Nature::SILHOUETTE) || ((f2)->getNature() & Nature::BORDER)))) {
2283       return false;
2284     }
2285
2286     return true;
2287   }
2288 };
2289
2290 void ViewMapBuilder::ComputeSweepLineIntersections(ViewMap *ioViewMap, real epsilon)
2291 {
2292   vector<SVertex *> &svertices = ioViewMap->SVertices();
2293   bool progressBarDisplay = false;
2294   unsigned sVerticesSize = svertices.size();
2295   unsigned fEdgesSize = ioViewMap->FEdges().size();
2296 #if 0
2297   if (_global.debug & G_DEBUG_FREESTYLE) {
2298     ViewMap::fedges_container& fedges = ioViewMap->FEdges();
2299     for (ViewMap::fedges_container::const_iterator f = fedges.begin(), end = fedges.end(); f != end; ++f) {
2300       cout << (*f)->aMaterialIndex() << "-" << (*f)->bMaterialIndex() << endl;
2301     }
2302   }
2303 #endif
2304   unsigned progressBarStep = 0;
2305
2306   if (_pProgressBar != NULL && fEdgesSize > gProgressBarMinSize) {
2307     unsigned progressBarSteps = min(gProgressBarMaxSteps, sVerticesSize);
2308     progressBarStep = sVerticesSize / progressBarSteps;
2309     _pProgressBar->reset();
2310     _pProgressBar->setLabelText("Computing Sweep Line Intersections");
2311     _pProgressBar->setTotalSteps(progressBarSteps);
2312     _pProgressBar->setProgress(0);
2313     progressBarDisplay = true;
2314   }
2315
2316   unsigned counter = progressBarStep;
2317
2318   sort(svertices.begin(), svertices.end(), less_SVertex2D(epsilon));
2319
2320   SweepLine<FEdge *, Vec3r> SL;
2321
2322   vector<FEdge *> &ioEdges = ioViewMap->FEdges();
2323
2324   vector<segment *> segments;
2325
2326   vector<FEdge *>::iterator fe, fend;
2327
2328   for (fe = ioEdges.begin(), fend = ioEdges.end(); fe != fend; fe++) {
2329     segment *s = new segment((*fe), (*fe)->vertexA()->point2D(), (*fe)->vertexB()->point2D());
2330     (*fe)->userdata = s;
2331     segments.push_back(s);
2332   }
2333
2334   vector<segment *> vsegments;
2335   for (vector<SVertex *>::iterator sv = svertices.begin(), svend = svertices.end(); sv != svend;
2336        sv++) {
2337     if (_pRenderMonitor && _pRenderMonitor->testBreak())
2338       break;
2339
2340     const vector<FEdge *> &vedges = (*sv)->fedges();
2341
2342     for (vector<FEdge *>::const_iterator sve = vedges.begin(), sveend = vedges.end();
2343          sve != sveend;
2344          sve++) {
2345       vsegments.push_back((segment *)((*sve)->userdata));
2346     }
2347
2348     Vec3r evt((*sv)->point2D());
2349     silhouette_binary_rule sbr;
2350     SL.process(evt, vsegments, sbr, epsilon);
2351
2352     if (progressBarDisplay) {
2353       counter--;
2354       if (counter <= 0) {
2355         counter = progressBarStep;
2356         _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
2357       }
2358     }
2359     vsegments.clear();
2360   }
2361
2362   if (_pRenderMonitor && _pRenderMonitor->testBreak()) {
2363     // delete segments
2364     if (!segments.empty()) {
2365       vector<segment *>::iterator s, send;
2366       for (s = segments.begin(), send = segments.end(); s != send; s++) {
2367         delete *s;
2368       }
2369     }
2370     return;
2371   }
2372
2373   // reset userdata:
2374   for (fe = ioEdges.begin(), fend = ioEdges.end(); fe != fend; fe++)
2375     (*fe)->userdata = NULL;
2376
2377   // list containing the new edges resulting from splitting operations.
2378   vector<FEdge *> newEdges;
2379
2380   // retrieve the intersected edges:
2381   vector<segment *> &iedges = SL.intersectedEdges();
2382   // retrieve the intersections:
2383   vector<intersection *> &intersections = SL.intersections();
2384
2385   int id = 0;
2386   // create a view vertex for each intersection and linked this one with the intersection object
2387   vector<intersection *>::iterator i, iend;
2388   for (i = intersections.begin(), iend = intersections.end(); i != iend; i++) {
2389     FEdge *fA = (*i)->EdgeA->edge();
2390     FEdge *fB = (*i)->EdgeB->edge();
2391
2392     Vec3r A1 = fA->vertexA()->point3D();
2393     Vec3r A2 = fA->vertexB()->point3D();
2394     Vec3r B1 = fB->vertexA()->point3D();
2395     Vec3r B2 = fB->vertexB()->point3D();
2396
2397     Vec3r a1 = fA->vertexA()->point2D();
2398     Vec3r a2 = fA->vertexB()->point2D();
2399     Vec3r b1 = fB->vertexA()->point2D();
2400     Vec3r b2 = fB->vertexB()->point2D();
2401
2402     real ta = (*i)->tA;
2403     real tb = (*i)->tB;
2404
2405     if ((ta < -epsilon) || (ta > 1 + epsilon))
2406       cerr << "Warning: 2D intersection out of range for edge " << fA->vertexA()->getId() << " - "
2407            << fA->vertexB()->getId() << endl;
2408
2409     if ((tb < -epsilon) || (tb > 1 + epsilon))
2410       cerr << "Warning: 2D intersection out of range for edge " << fB->vertexA()->getId() << " - "
2411            << fB->vertexB()->getId() << endl;
2412
2413     real Ta = SilhouetteGeomEngine::ImageToWorldParameter(fA, ta);
2414     real Tb = SilhouetteGeomEngine::ImageToWorldParameter(fB, tb);
2415
2416     if ((Ta < -epsilon) || (Ta > 1 + epsilon))
2417       cerr << "Warning: 3D intersection out of range for edge " << fA->vertexA()->getId() << " - "
2418            << fA->vertexB()->getId() << endl;
2419
2420     if ((Tb < -epsilon) || (Tb > 1 + epsilon))
2421       cerr << "Warning: 3D intersection out of range for edge " << fB->vertexA()->getId() << " - "
2422            << fB->vertexB()->getId() << endl;
2423
2424 #if 0
2425     if (_global.debug & G_DEBUG_FREESTYLE) {
2426       if ((Ta < -epsilon) || (Ta > 1 + epsilon) || (Tb < -epsilon) || (Tb > 1 + epsilon)) {
2427         printf("ta %.12e\n", ta);
2428         printf("tb %.12e\n", tb);
2429         printf("a1 %e, %e -- a2 %e, %e\n", a1[0], a1[1], a2[0], a2[1]);
2430         printf("b1 %e, %e -- b2 %e, %e\n", b1[0], b1[1], b2[0], b2[1]);
2431         //printf("line([%e, %e], [%e, %e]);\n", a1[0], a2[0], a1[1], a2[1]);
2432         //printf("line([%e, %e], [%e, %e]);\n", b1[0], b2[0], b1[1], b2[1]);
2433         if ((Ta < -epsilon) || (Ta > 1 + epsilon))
2434           printf("Ta %.12e\n", Ta);
2435         if ((Tb < -epsilon) || (Tb > 1 + epsilon))
2436           printf("Tb %.12e\n", Tb);
2437         printf("A1 %e, %e, %e -- A2 %e, %e, %e\n", A1[0], A1[1], A1[2], A2[0], A2[1], A2[2]);
2438         printf("B1 %e, %e, %e -- B2 %e, %e, %e\n", B1[0], B1[1], B1[2], B2[0], B2[1], B2[2]);
2439       }
2440     }
2441 #endif
2442
2443     TVertex *tvertex = ioViewMap->CreateTVertex(Vec3r(A1 + Ta * (A2 - A1)),
2444                                                 Vec3r(a1 + ta * (a2 - a1)),
2445                                                 fA,
2446                                                 Vec3r(B1 + Tb * (B2 - B1)),
2447                                                 Vec3r(b1 + tb * (b2 - b1)),
2448                                                 fB,
2449                                                 id);
2450
2451     (*i)->userdata = tvertex;
2452     ++id;
2453   }
2454
2455   progressBarStep = 0;
2456
2457   if (progressBarDisplay) {
2458     unsigned iEdgesSize = iedges.size();
2459     unsigned progressBarSteps = min(gProgressBarMaxSteps, iEdgesSize);
2460     progressBarStep = iEdgesSize / progressBarSteps;
2461     _pProgressBar->reset();
2462     _pProgressBar->setLabelText("Splitting intersected edges");
2463     _pProgressBar->setTotalSteps(progressBarSteps);
2464     _pProgressBar->setProgress(0);
2465   }
2466
2467   counter = progressBarStep;
2468
2469   vector<TVertex *> edgeVVertices;
2470   vector<ViewEdge *> newVEdges;
2471   vector<segment *>::iterator s, send;
2472   for (s = iedges.begin(), send = iedges.end(); s != send; s++) {
2473     edgeVVertices.clear();
2474     newEdges.clear();
2475     newVEdges.clear();
2476
2477     FEdge *fedge = (*s)->edge();
2478     ViewEdge *vEdge = fedge->viewedge();
2479     ViewShape *shape = vEdge->viewShape();
2480
2481     vector<intersection *> &eIntersections = (*s)->intersections();
2482     // we first need to sort these intersections from farther to closer to A
2483     sort(eIntersections.begin(), eIntersections.end(), less_Intersection(*s));
2484     for (i = eIntersections.begin(), iend = eIntersections.end(); i != iend; i++)
2485       edgeVVertices.push_back((TVertex *)(*i)->userdata);
2486
2487     shape->SplitEdge(fedge, edgeVVertices, ioViewMap->FEdges(), ioViewMap->ViewEdges());
2488
2489     if (progressBarDisplay) {
2490       counter--;
2491       if (counter <= 0) {
2492         counter = progressBarStep;
2493         _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
2494       }
2495     }
2496   }
2497
2498   // reset userdata:
2499   for (fe = ioEdges.begin(), fend = ioEdges.end(); fe != fend; fe++)
2500     (*fe)->userdata = NULL;
2501
2502   // delete segments
2503   if (!segments.empty()) {
2504     for (s = segments.begin(), send = segments.end(); s != send; s++) {
2505       delete *s;
2506     }
2507   }
2508 }
2509
2510 } /* namespace Freestyle */