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