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