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