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