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