Freestyle: fix wrong arg order, and cleanup confusing loop (both reported by coverity).
[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->setLibraryPath((*it)->getLibraryPath());
1208                 psShape->setFrsMaterials((*it)->frs_materials()); // FIXME
1209
1210                 // create the view shape
1211                 ViewShape *vshape = new ViewShape(psShape);
1212                 // add this view shape to the view map:
1213                 _ViewMap->AddViewShape(vshape);
1214
1215                 // we want to number the view edges in a unique way for the while scene.
1216                 _pViewEdgeBuilder->setCurrentViewId(_currentId);
1217                 // we want to number the feature edges in a unique way for the while scene.
1218                 _pViewEdgeBuilder->setCurrentFId(_currentFId);
1219                 // we want to number the SVertex in a unique way for the while scene.
1220                 _pViewEdgeBuilder->setCurrentSVertexId(_currentFId);
1221                 _pViewEdgeBuilder->BuildViewEdges(dynamic_cast<WXShape*>(*it), vshape, _ViewMap->ViewEdges(),
1222                                                   _ViewMap->ViewVertices(), _ViewMap->FEdges(), _ViewMap->SVertices());
1223
1224                 _currentId = _pViewEdgeBuilder->currentViewId() + 1;
1225                 _currentFId = _pViewEdgeBuilder->currentFId() + 1;
1226                 _currentSVertexId = _pViewEdgeBuilder->currentSVertexId() + 1;
1227
1228                 psShape->ComputeBBox();
1229         }
1230 }
1231
1232 void ViewMapBuilder::computeCusps(ViewMap *ioViewMap)
1233 {
1234         vector<ViewVertex*> newVVertices;
1235         vector<ViewEdge*> newVEdges;
1236         ViewMap::viewedges_container& vedges = ioViewMap->ViewEdges();
1237         ViewMap::viewedges_container::iterator ve = vedges.begin(), veend = vedges.end();
1238         for (; ve != veend; ++ve) {
1239                 if (_pRenderMonitor && _pRenderMonitor->testBreak())
1240                         break;
1241                 if ((!((*ve)->getNature() & Nature::SILHOUETTE)) || (!((*ve)->fedgeA()->isSmooth())))
1242                         continue;
1243                 FEdge *fe = (*ve)->fedgeA();
1244                 FEdge *fefirst = fe;
1245                 bool first = true;
1246                 bool positive = true;
1247                 do {
1248                         FEdgeSmooth *fes = dynamic_cast<FEdgeSmooth*>(fe);
1249                         Vec3r A((fes)->vertexA()->point3d());
1250                         Vec3r B((fes)->vertexB()->point3d());
1251                         Vec3r AB(B - A);
1252                         AB.normalize();
1253                         Vec3r m((A + B) / 2.0);
1254                         Vec3r crossP(AB ^ (fes)->normal());
1255                         crossP.normalize();
1256                         Vec3r viewvector;
1257                         if (_orthographicProjection) {
1258                                 viewvector = Vec3r(0.0, 0.0, m.z() - _viewpoint.z());
1259                         }
1260                         else {
1261                                 viewvector = Vec3r(m - _viewpoint);
1262                         }
1263                         viewvector.normalize();
1264                         if (first) {
1265                                 if (((crossP) * (viewvector)) > 0)
1266                                         positive = true;
1267                                 else
1268                                         positive = false;
1269                                 first = false;
1270                         }
1271                         // If we're in a positive part, we need a stronger negative value to change
1272                         NonTVertex *cusp = NULL;
1273                         if (positive) {
1274                                 if (((crossP) * (viewvector)) < -0.1) {
1275                                         // state changes
1276                                         positive = false;
1277                                         // creates and insert cusp
1278                                         cusp = dynamic_cast<NonTVertex*>(ioViewMap->InsertViewVertex(fes->vertexA(), newVEdges));
1279                                         if (cusp)
1280                                                 cusp->setNature(cusp->getNature() | Nature::CUSP);
1281                                 }
1282                         }
1283                         else {
1284                                 // If we're in a negative part, we need a stronger negative value to change
1285                                 if (((crossP) * (viewvector)) > 0.1) {
1286                                         positive = true;
1287                                         cusp = dynamic_cast<NonTVertex*>(ioViewMap->InsertViewVertex(fes->vertexA(), newVEdges));
1288                                         if (cusp)
1289                                                 cusp->setNature(cusp->getNature() | Nature::CUSP);
1290                                 }
1291                         }
1292                         fe = fe->nextEdge();
1293                 } while (fe && fe != fefirst);
1294         }
1295         for (ve = newVEdges.begin(), veend = newVEdges.end(); ve != veend; ++ve) {
1296                 (*ve)->viewShape()->AddEdge(*ve);
1297                 vedges.push_back(*ve);
1298         }
1299 }
1300
1301 void ViewMapBuilder::ComputeCumulativeVisibility(ViewMap *ioViewMap, WingedEdge& we, const BBox<Vec3r>& bbox,
1302                                                  real epsilon, bool cull, GridDensityProviderFactory& factory)
1303 {
1304         AutoPtr<GridHelpers::Transform> transform;
1305         AutoPtr<OccluderSource> source;
1306
1307         if (_orthographicProjection) {
1308                 transform.reset(new BoxGrid::Transform);
1309         }
1310         else {
1311                 transform.reset(new SphericalGrid::Transform);
1312         }
1313
1314         if (cull) {
1315                 source.reset(new CulledOccluderSource(*transform, we, *ioViewMap, true));
1316         }
1317         else {
1318                 source.reset(new OccluderSource(*transform, we));
1319         }
1320
1321         AutoPtr<GridDensityProvider> density(factory.newGridDensityProvider(*source, bbox, *transform));
1322
1323         if (_orthographicProjection) {
1324                 BoxGrid grid(*source, *density, ioViewMap, _viewpoint, _EnableQI);
1325                 computeCumulativeVisibility<BoxGrid, BoxGrid::Iterator>(ioViewMap, grid, epsilon, _pRenderMonitor);
1326         }
1327         else {
1328                 SphericalGrid grid(*source, *density, ioViewMap, _viewpoint, _EnableQI);
1329                 computeCumulativeVisibility<SphericalGrid, SphericalGrid::Iterator>(ioViewMap, grid, epsilon, _pRenderMonitor);
1330         }
1331 }
1332
1333 void ViewMapBuilder::ComputeDetailedVisibility(ViewMap *ioViewMap, WingedEdge& we, const BBox<Vec3r>& bbox,
1334                                                real epsilon, bool cull, GridDensityProviderFactory& factory)
1335 {
1336         AutoPtr<GridHelpers::Transform> transform;
1337         AutoPtr<OccluderSource> source;
1338
1339         if (_orthographicProjection) {
1340                 transform.reset(new BoxGrid::Transform);
1341         }
1342         else {
1343                 transform.reset(new SphericalGrid::Transform);
1344         }
1345
1346         if (cull) {
1347                 source.reset(new CulledOccluderSource(*transform, we, *ioViewMap, true));
1348         }
1349         else {
1350                 source.reset(new OccluderSource(*transform, we));
1351         }
1352
1353         AutoPtr<GridDensityProvider> density(factory.newGridDensityProvider(*source, bbox, *transform));
1354
1355         if (_orthographicProjection) {
1356                 BoxGrid grid(*source, *density, ioViewMap, _viewpoint, _EnableQI);
1357                 computeDetailedVisibility<BoxGrid, BoxGrid::Iterator>(ioViewMap, grid, epsilon, _pRenderMonitor);
1358         }
1359         else {
1360                 SphericalGrid grid(*source, *density, ioViewMap, _viewpoint, _EnableQI);
1361                 computeDetailedVisibility<SphericalGrid, SphericalGrid::Iterator>(ioViewMap, grid, epsilon, _pRenderMonitor);
1362         }
1363 }
1364
1365 void ViewMapBuilder::ComputeEdgesVisibility(ViewMap *ioViewMap, WingedEdge& we, const BBox<Vec3r>& bbox,
1366                                             unsigned int sceneNumFaces, visibility_algo iAlgo, real epsilon)
1367 {
1368 #if 0
1369         iAlgo = ray_casting; // for testing algorithms equivalence
1370 #endif
1371         switch (iAlgo) {
1372                 case ray_casting:
1373                         if (_global.debug & G_DEBUG_FREESTYLE) {
1374                                 cout << "Using ordinary ray casting" << endl;
1375                         }
1376                         BuildGrid(we, bbox, sceneNumFaces);
1377                         ComputeRayCastingVisibility(ioViewMap, epsilon);
1378                         break;
1379                 case ray_casting_fast:
1380                         if (_global.debug & G_DEBUG_FREESTYLE) {
1381                                 cout << "Using fast ray casting" << endl;
1382                         }
1383                         BuildGrid(we, bbox, sceneNumFaces);
1384                         ComputeFastRayCastingVisibility(ioViewMap, epsilon);
1385                         break;
1386                 case ray_casting_very_fast:
1387                         if (_global.debug & G_DEBUG_FREESTYLE) {
1388                                 cout << "Using very fast ray casting" << endl;
1389                         }
1390                         BuildGrid(we, bbox, sceneNumFaces);
1391                         ComputeVeryFastRayCastingVisibility(ioViewMap, epsilon);
1392                         break;
1393                 case ray_casting_culled_adaptive_traditional:
1394                         if (_global.debug & G_DEBUG_FREESTYLE) {
1395                                 cout << "Using culled adaptive grid with heuristic density and traditional QI calculation" << endl;
1396                         }
1397                         try {
1398                                 HeuristicGridDensityProviderFactory factory(0.5f, sceneNumFaces);
1399                                 ComputeDetailedVisibility(ioViewMap, we, bbox, epsilon, true, factory);
1400                         }
1401                         catch (...) {
1402                                 // Last resort catch to make sure RAII semantics hold for OptimizedGrid. Can be replaced with
1403                                 // try...catch block around main() if the program as a whole is converted to RAII
1404
1405                                 // This is the little-mentioned caveat of RAII: RAII does not work unless destructors are always
1406                                 // called, but destructors are only called if all exceptions are caught (or std::terminate() is
1407                                 // replaced).
1408
1409                                 // We don't actually handle the exception here, so re-throw it now that our destructors have had a
1410                                 // chance to run.
1411                                 throw;
1412                         }
1413                         break;
1414                 case ray_casting_adaptive_traditional:
1415                         if (_global.debug & G_DEBUG_FREESTYLE) {
1416                                 cout << "Using unculled adaptive grid with heuristic density and traditional QI calculation" << endl;
1417                         }
1418                         try {
1419                                 HeuristicGridDensityProviderFactory factory(0.5f, sceneNumFaces);
1420                                 ComputeDetailedVisibility(ioViewMap, we, bbox, epsilon, false, factory);
1421                         }
1422                         catch (...) {
1423                                 throw;
1424                         }
1425                         break;
1426                 case ray_casting_culled_adaptive_cumulative:
1427                         if (_global.debug & G_DEBUG_FREESTYLE) {
1428                                 cout << "Using culled adaptive grid with heuristic density and cumulative QI calculation" << endl;
1429                         }
1430                         try {
1431                                 HeuristicGridDensityProviderFactory factory(0.5f, sceneNumFaces);
1432                                 ComputeCumulativeVisibility(ioViewMap, we, bbox, epsilon, true, factory);
1433                         }
1434                         catch (...) {
1435                                 throw;
1436                         }
1437                         break;
1438                 case ray_casting_adaptive_cumulative:
1439                         if (_global.debug & G_DEBUG_FREESTYLE) {
1440                                 cout << "Using unculled adaptive grid with heuristic density and cumulative QI calculation" << endl;
1441                         }
1442                         try {
1443                                 HeuristicGridDensityProviderFactory factory(0.5f, sceneNumFaces);
1444                                 ComputeCumulativeVisibility(ioViewMap, we, bbox, epsilon, false, factory);
1445                         }
1446                         catch (...) {
1447                                 throw;
1448                         }
1449                         break;
1450                 default:
1451                         break;
1452         }
1453 }
1454
1455 static const unsigned gProgressBarMaxSteps = 10;
1456 static const unsigned gProgressBarMinSize = 2000;
1457
1458 void ViewMapBuilder::ComputeRayCastingVisibility(ViewMap *ioViewMap, real epsilon)
1459 {
1460         vector<ViewEdge*>& vedges = ioViewMap->ViewEdges();
1461         bool progressBarDisplay = false;
1462         unsigned progressBarStep = 0;
1463         unsigned vEdgesSize = vedges.size();
1464         unsigned fEdgesSize = ioViewMap->FEdges().size();
1465
1466         if (_pProgressBar != NULL && fEdgesSize > gProgressBarMinSize) {
1467                 unsigned progressBarSteps = min(gProgressBarMaxSteps, vEdgesSize);
1468                 progressBarStep = vEdgesSize / progressBarSteps;
1469                 _pProgressBar->reset();
1470                 _pProgressBar->setLabelText("Computing Ray casting Visibility");
1471                 _pProgressBar->setTotalSteps(progressBarSteps);
1472                 _pProgressBar->setProgress(0);
1473                 progressBarDisplay = true;
1474         }
1475
1476         unsigned counter = progressBarStep;
1477         FEdge *fe, *festart;
1478         int nSamples = 0;
1479         vector<Polygon3r*> aFaces;
1480         Polygon3r *aFace = NULL;
1481         unsigned tmpQI = 0;
1482         unsigned qiClasses[256];
1483         unsigned maxIndex, maxCard;
1484         unsigned qiMajority;
1485         static unsigned timestamp = 1;
1486         for (vector<ViewEdge*>::iterator ve = vedges.begin(), veend = vedges.end(); ve != veend; ve++) {
1487                 if (_pRenderMonitor && _pRenderMonitor->testBreak())
1488                         break;
1489 #if LOGGING
1490                 if (_global.debug & G_DEBUG_FREESTYLE) {
1491                         cout << "Processing ViewEdge " << (*ve)->getId() << endl;
1492                 }
1493 #endif
1494                 festart = (*ve)->fedgeA();
1495                 fe = (*ve)->fedgeA();
1496                 qiMajority = 1;
1497                 do {
1498                         qiMajority++;
1499                         fe = fe->nextEdge();
1500                 } while (fe && fe != festart);
1501                 qiMajority >>= 1;
1502 #if LOGGING
1503                 if (_global.debug & G_DEBUG_FREESTYLE) {
1504                         cout << "\tqiMajority: " << qiMajority << endl;
1505                 }
1506 #endif
1507
1508                 tmpQI = 0;
1509                 maxIndex = 0;
1510                 maxCard = 0;
1511                 nSamples = 0;
1512                 fe = (*ve)->fedgeA();
1513                 memset(qiClasses, 0, 256 * sizeof(*qiClasses));
1514                 set<ViewShape*> occluders;
1515                 do {
1516                         if ((maxCard < qiMajority)) {
1517                                 tmpQI = ComputeRayCastingVisibility(fe, _Grid, epsilon, occluders, &aFace, timestamp++);
1518
1519 #if LOGGING
1520                                 if (_global.debug & G_DEBUG_FREESTYLE) {
1521                                         cout << "\tFEdge: visibility " << tmpQI << endl;
1522                                 }
1523 #endif
1524                                 //ARB: This is an error condition, not an alert condition.
1525                                 // Some sort of recovery or abort is necessary.
1526                                 if (tmpQI >= 256) {
1527                                         cerr << "Warning: too many occluding levels" << endl;
1528                                         //ARB: Wild guess: instead of aborting or corrupting memory, treat as tmpQI == 255
1529                                         tmpQI = 255;
1530                                 }
1531
1532                                 if (++qiClasses[tmpQI] > maxCard) {
1533                                         maxCard = qiClasses[tmpQI];
1534                                         maxIndex = tmpQI;
1535                                 }
1536                         }
1537                         else {
1538                                 //ARB: FindOccludee is redundant if ComputeRayCastingVisibility has been called
1539                                 FindOccludee(fe, _Grid, epsilon, &aFace, timestamp++);
1540 #if LOGGING
1541                                 if (_global.debug & G_DEBUG_FREESTYLE) {
1542                                         cout << "\tFEdge: occludee only (" << (aFace != NULL ? "found" : "not found") << ")" << endl;
1543                                 }
1544 #endif
1545                         }
1546
1547                         if (aFace) {
1548                                 fe->setaFace(*aFace);
1549                                 aFaces.push_back(aFace);
1550                                 fe->setOccludeeEmpty(false);
1551 #if LOGGING
1552                                 if (_global.debug & G_DEBUG_FREESTYLE) {
1553                                         cout << "\tFound occludee" << endl;
1554                                 }
1555 #endif
1556                         }
1557                         else {
1558                                 //ARB: We are arbitrarily using the last observed value for occludee (almost always the value observed
1559                                 //     for the edge before festart). Is that meaningful?
1560                                 // ...in fact, _occludeeEmpty seems to be unused.
1561                                 fe->setOccludeeEmpty(true);
1562                         }
1563
1564                         ++nSamples;
1565                         fe = fe->nextEdge();
1566                 } while ((maxCard < qiMajority) && (fe) && (fe != festart));
1567 #if LOGGING
1568                 if (_global.debug & G_DEBUG_FREESTYLE) {
1569                         cout << "\tFinished with " << nSamples << " samples, maxCard = " << maxCard << endl;
1570                 }
1571 #endif
1572
1573                 // ViewEdge
1574                 // qi --
1575                 (*ve)->setQI(maxIndex);
1576                 // occluders --
1577                 for (set<ViewShape*>::iterator o = occluders.begin(), oend = occluders.end(); o != oend; ++o)
1578                         (*ve)->AddOccluder((*o));
1579 #if LOGGING
1580                 if (_global.debug & G_DEBUG_FREESTYLE) {
1581                         cout << "\tConclusion: QI = " << maxIndex << ", " << (*ve)->occluders_size() << " occluders." << endl;
1582                 }
1583 #endif
1584                 // occludee --
1585                 if (!aFaces.empty()) {
1586                         if (aFaces.size() <= (float)nSamples / 2.0f) {
1587                                 (*ve)->setaShape(0);
1588                         }
1589                         else {
1590                                 vector<Polygon3r*>::iterator p = aFaces.begin();
1591                                 WFace *wface = (WFace *)((*p)->userdata);
1592                                 ViewShape *vshape = ioViewMap->viewShape(wface->GetVertex(0)->shape()->GetId());
1593                                 ++p;
1594                                 (*ve)->setaShape(vshape);
1595                         }
1596                 }
1597
1598                 if (progressBarDisplay) {
1599                         counter--;
1600                         if (counter <= 0) {
1601                                 counter = progressBarStep;
1602                                 _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
1603                         }
1604                 }
1605                 aFaces.clear();
1606         }
1607 }
1608
1609 void ViewMapBuilder::ComputeFastRayCastingVisibility(ViewMap *ioViewMap, real epsilon)
1610 {
1611         vector<ViewEdge*>& vedges = ioViewMap->ViewEdges();
1612         bool progressBarDisplay = false;
1613         unsigned progressBarStep = 0;
1614         unsigned vEdgesSize = vedges.size();
1615         unsigned fEdgesSize = ioViewMap->FEdges().size();
1616
1617         if (_pProgressBar != NULL && fEdgesSize > gProgressBarMinSize) {
1618                 unsigned progressBarSteps = min(gProgressBarMaxSteps, vEdgesSize);
1619                 progressBarStep = vEdgesSize / progressBarSteps;
1620                 _pProgressBar->reset();
1621                 _pProgressBar->setLabelText("Computing Ray casting Visibility");
1622                 _pProgressBar->setTotalSteps(progressBarSteps);
1623                 _pProgressBar->setProgress(0);
1624                 progressBarDisplay = true;
1625         }
1626
1627         unsigned counter = progressBarStep;
1628         FEdge *fe, *festart;
1629         unsigned nSamples = 0;
1630         vector<Polygon3r*> aFaces;
1631         Polygon3r *aFace = NULL;
1632         unsigned tmpQI = 0;
1633         unsigned qiClasses[256];
1634         unsigned maxIndex, maxCard;
1635         unsigned qiMajority;
1636         static unsigned timestamp = 1;
1637         bool even_test;
1638         for (vector<ViewEdge*>::iterator ve = vedges.begin(), veend = vedges.end(); ve != veend; ve++) {
1639                 if (_pRenderMonitor && _pRenderMonitor->testBreak())
1640                         break;
1641
1642                 festart = (*ve)->fedgeA();
1643                 fe = (*ve)->fedgeA();
1644                 qiMajority = 1;
1645                 do {
1646                         qiMajority++;
1647                         fe = fe->nextEdge();
1648                 } while (fe && fe != festart);
1649                 if (qiMajority >= 4)
1650                         qiMajority >>= 2;
1651                 else
1652                         qiMajority = 1;
1653
1654                 set<ViewShape*> occluders;
1655
1656                 even_test = true;
1657                 maxIndex = 0;
1658                 maxCard = 0;
1659                 nSamples = 0;
1660                 memset(qiClasses, 0, 256 * sizeof(*qiClasses));
1661                 fe = (*ve)->fedgeA();
1662                 do {
1663                         if (even_test) {
1664                                 if ((maxCard < qiMajority)) {
1665                                         tmpQI = ComputeRayCastingVisibility(fe, _Grid, epsilon, occluders, &aFace, timestamp++);
1666
1667                                         //ARB: This is an error condition, not an alert condition.
1668                                         // Some sort of recovery or abort is necessary.
1669                                         if (tmpQI >= 256) {
1670                                                 cerr << "Warning: too many occluding levels" << endl;
1671                                                 //ARB: Wild guess: instead of aborting or corrupting memory, treat as tmpQI == 255
1672                                                 tmpQI = 255;
1673                                         }
1674
1675                                         if (++qiClasses[tmpQI] > maxCard) {
1676                                                 maxCard = qiClasses[tmpQI];
1677                                                 maxIndex = tmpQI;
1678                                         }
1679                                 }
1680                                 else {
1681                                         //ARB: FindOccludee is redundant if ComputeRayCastingVisibility has been called
1682                                         FindOccludee(fe, _Grid, epsilon, &aFace, timestamp++);
1683                                 }
1684
1685                                 if (aFace) {
1686                                         fe->setaFace(*aFace);
1687                                         aFaces.push_back(aFace);
1688                                 }
1689                                 ++nSamples;
1690                                 even_test = false;
1691                         }
1692                         else {
1693                                 even_test = true;
1694                         }
1695                         fe = fe->nextEdge();
1696                 } while ((maxCard < qiMajority) && (fe) && (fe != festart));
1697
1698                 (*ve)->setQI(maxIndex);
1699
1700                 if (!aFaces.empty()) {
1701                         if (aFaces.size() < nSamples / 2) {
1702                                 (*ve)->setaShape(0);
1703                         }
1704                         else {
1705                                 vector<Polygon3r*>::iterator p = aFaces.begin();
1706                                 WFace *wface = (WFace *)((*p)->userdata);
1707                                 ViewShape *vshape = ioViewMap->viewShape(wface->GetVertex(0)->shape()->GetId());
1708                                 ++p;
1709 #if 0
1710                                 for (; p != pend; ++p) {
1711                                         WFace *f = (WFace*)((*p)->userdata);
1712                                         ViewShape *vs = ioViewMap->viewShape(f->GetVertex(0)->shape()->GetId());
1713                                         if (vs != vshape) {
1714                                                 sameShape = false;
1715                                                 break;
1716                                         }
1717                                 }
1718                                 if (sameShape)
1719 #endif
1720                                 (*ve)->setaShape(vshape);
1721                         }
1722                 }
1723
1724                 //(*ve)->setaFace(aFace);
1725
1726                 if (progressBarDisplay) {
1727                         counter--;
1728                         if (counter <= 0) {
1729                                 counter = progressBarStep;
1730                                 _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
1731                         }
1732                 }
1733                 aFaces.clear();
1734         }
1735 }
1736
1737 void ViewMapBuilder::ComputeVeryFastRayCastingVisibility(ViewMap *ioViewMap, real epsilon)
1738 {
1739         vector<ViewEdge*>& vedges = ioViewMap->ViewEdges();
1740         bool progressBarDisplay = false;
1741         unsigned progressBarStep = 0;
1742         unsigned vEdgesSize = vedges.size();
1743         unsigned fEdgesSize = ioViewMap->FEdges().size();
1744
1745         if (_pProgressBar != NULL && fEdgesSize > gProgressBarMinSize) {
1746                 unsigned progressBarSteps = min(gProgressBarMaxSteps, vEdgesSize);
1747                 progressBarStep = vEdgesSize / progressBarSteps;
1748                 _pProgressBar->reset();
1749                 _pProgressBar->setLabelText("Computing Ray casting Visibility");
1750                 _pProgressBar->setTotalSteps(progressBarSteps);
1751                 _pProgressBar->setProgress(0);
1752                 progressBarDisplay = true;
1753         }
1754
1755         unsigned counter = progressBarStep;
1756         FEdge *fe;
1757         unsigned qi = 0;
1758         Polygon3r *aFace = NULL;
1759         static unsigned timestamp = 1;
1760         for (vector<ViewEdge*>::iterator ve = vedges.begin(), veend = vedges.end(); ve != veend; ve++) {
1761                 if (_pRenderMonitor && _pRenderMonitor->testBreak())
1762                         break;
1763
1764                 set<ViewShape*> occluders;
1765
1766                 fe = (*ve)->fedgeA();
1767                 qi = ComputeRayCastingVisibility(fe, _Grid, epsilon, occluders, &aFace, timestamp++);
1768                 if (aFace) {
1769                         fe->setaFace(*aFace);
1770                         WFace *wface = (WFace *)(aFace->userdata);
1771                         ViewShape *vshape = ioViewMap->viewShape(wface->GetVertex(0)->shape()->GetId());
1772                         (*ve)->setaShape(vshape);
1773                 }
1774                 else {
1775                         (*ve)->setaShape(0);
1776                 }
1777
1778                 (*ve)->setQI(qi);
1779
1780                 if (progressBarDisplay) {
1781                         counter--;
1782                         if (counter <= 0) {
1783                                 counter = progressBarStep;
1784                                 _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
1785                         }
1786                 }
1787         }
1788 }
1789
1790 void ViewMapBuilder::FindOccludee(FEdge *fe, Grid *iGrid, real epsilon, Polygon3r **oaPolygon, unsigned timestamp,
1791                                   Vec3r& u, Vec3r& A, Vec3r& origin, Vec3r& edge, vector<WVertex*>& faceVertices)
1792 {
1793         WFace *face = NULL;
1794         if (fe->isSmooth()) {
1795                 FEdgeSmooth *fes = dynamic_cast<FEdgeSmooth*>(fe);
1796                 face = (WFace *)fes->face();
1797         }
1798         OccludersSet occluders;
1799         WFace *oface;
1800         bool skipFace;
1801
1802         WVertex::incoming_edge_iterator ie;
1803         OccludersSet::iterator p, pend;
1804
1805         *oaPolygon = NULL;
1806         if (((fe)->getNature() & Nature::SILHOUETTE) || ((fe)->getNature() & Nature::BORDER)) {
1807                 occluders.clear();
1808                 // we cast a ray from A in the same direction but looking behind
1809                 Vec3r v(-u[0], -u[1], -u[2]);
1810                 iGrid->castInfiniteRay(A, v, occluders, timestamp);
1811
1812                 bool noIntersection = true;
1813                 real mint = FLT_MAX;
1814                 // we met some occluders, let us fill the aShape field with the first intersected occluder
1815                 for (p = occluders.begin(), pend = occluders.end(); p != pend; p++) {
1816                         // check whether the edge and the polygon plane are coincident:
1817                         //-------------------------------------------------------------
1818                         //first let us compute the plane equation.
1819                         oface = (WFace *)(*p)->userdata;
1820                         Vec3r v1(((*p)->getVertices())[0]);
1821                         Vec3r normal((*p)->getNormal());
1822                         real d = -(v1 * normal);
1823                         real t, t_u, t_v;
1824
1825                         if (face) {
1826                                 skipFace = false;
1827
1828                                 if (face == oface)
1829                                         continue;
1830
1831                                 if (faceVertices.empty())
1832                                         continue;
1833
1834                                 for (vector<WVertex*>::iterator fv = faceVertices.begin(), fvend = faceVertices.end();
1835                                      fv != fvend;
1836                                      ++fv)
1837                                 {
1838                                         if ((*fv)->isBoundary())
1839                                                 continue;
1840                                         WVertex::incoming_edge_iterator iebegin = (*fv)->incoming_edges_begin();
1841                                         WVertex::incoming_edge_iterator ieend = (*fv)->incoming_edges_end();
1842                                         for (ie = iebegin; ie != ieend; ++ie) {
1843                                                 if ((*ie) == 0)
1844                                                         continue;
1845
1846                                                 WFace *sface = (*ie)->GetbFace();
1847                                                 if (sface == oface) {
1848                                                         skipFace = true;
1849                                                         break;
1850                                                 }
1851                                         }
1852                                         if (skipFace)
1853                                                 break;
1854                                 }
1855                                 if (skipFace)
1856                                         continue;
1857                         }
1858                         else {
1859                                 if (GeomUtils::COINCIDENT == GeomUtils::intersectRayPlane(origin, edge, normal, d, t, epsilon))
1860                                         continue;
1861                         }
1862                         if ((*p)->rayIntersect(A, v, t, t_u, t_v)) {
1863                                 if (fabs(v * normal) > 0.0001) {
1864                                         if (t > 0.0) { // && t < 1.0) {
1865                                                 if (t < mint) {
1866                                                         *oaPolygon = (*p);
1867                                                         mint = t;
1868                                                         noIntersection = false;
1869                                                         fe->setOccludeeIntersection(Vec3r(A + t * v));
1870                                                 }
1871                                         }
1872                                 }
1873                         }
1874                 }
1875
1876                 if (noIntersection)
1877                         *oaPolygon = NULL;
1878         }
1879 }
1880
1881 void ViewMapBuilder::FindOccludee(FEdge *fe, Grid *iGrid, real epsilon, Polygon3r **oaPolygon, unsigned timestamp)
1882 {
1883         OccludersSet occluders;
1884
1885         Vec3r A;
1886         Vec3r edge;
1887         Vec3r origin;
1888         A = Vec3r(((fe)->vertexA()->point3D() + (fe)->vertexB()->point3D()) / 2.0);
1889         edge = Vec3r((fe)->vertexB()->point3D() - (fe)->vertexA()->point3D());
1890         origin = Vec3r((fe)->vertexA()->point3D());
1891         Vec3r u;
1892         if (_orthographicProjection) {
1893                 u = Vec3r(0.0, 0.0, _viewpoint.z() - A.z());
1894         }
1895         else {
1896                 u = Vec3r(_viewpoint - A);
1897         }
1898         u.normalize();
1899         if (A < iGrid->getOrigin())
1900                 cerr << "Warning: point is out of the grid for fedge " << fe->getId().getFirst() << "-" <<
1901                         fe->getId().getSecond() << endl;
1902
1903         vector<WVertex*> faceVertices;
1904
1905         WFace *face = NULL;
1906         if (fe->isSmooth()) {
1907                 FEdgeSmooth *fes = dynamic_cast<FEdgeSmooth*>(fe);
1908                 face = (WFace *)fes->face();
1909         }
1910         if (face)
1911                 face->RetrieveVertexList(faceVertices);
1912
1913         return FindOccludee(fe, iGrid, epsilon, oaPolygon, timestamp, u, A, origin, edge, faceVertices);
1914 }
1915
1916 int ViewMapBuilder::ComputeRayCastingVisibility(FEdge *fe, Grid *iGrid, real epsilon, set<ViewShape*>& oOccluders,
1917                                                 Polygon3r **oaPolygon, unsigned timestamp)
1918 {
1919         OccludersSet occluders;
1920         int qi = 0;
1921
1922         Vec3r center;
1923         Vec3r edge;
1924         Vec3r origin;
1925
1926         center = fe->center3d();
1927         edge = Vec3r(fe->vertexB()->point3D() - fe->vertexA()->point3D());
1928         origin = Vec3r(fe->vertexA()->point3D());
1929         // Is the edge outside the view frustum ?
1930         Vec3r gridOrigin(iGrid->getOrigin());
1931         Vec3r gridExtremity(iGrid->getOrigin() + iGrid->gridSize());
1932
1933         if ((center.x() < gridOrigin.x()) || (center.y() < gridOrigin.y()) || (center.z() < gridOrigin.z()) ||
1934             (center.x() > gridExtremity.x()) || (center.y() > gridExtremity.y()) || (center.z() > gridExtremity.z()))
1935         {
1936                 cerr << "Warning: point is out of the grid for fedge " << fe->getId() << endl;
1937                 //return 0;
1938         }
1939
1940 #if 0
1941         Vec3r A(fe->vertexA()->point2d());
1942         Vec3r B(fe->vertexB()->point2d());
1943         int viewport[4];
1944         SilhouetteGeomEngine::retrieveViewport(viewport);
1945         if ((A.x() < viewport[0]) || (A.x() > viewport[2]) || (A.y() < viewport[1]) || (A.y() > viewport[3]) ||
1946             (B.x() < viewport[0]) || (B.x() > viewport[2]) || (B.y() < viewport[1]) || (B.y() > viewport[3])) {
1947                 cerr << "Warning: point is out of the grid for fedge " << fe->getId() << endl;
1948                 //return 0;
1949         }
1950 #endif
1951
1952         Vec3r vp;
1953         if (_orthographicProjection) {
1954                 vp = Vec3r(center.x(), center.y(), _viewpoint.z());
1955         }
1956         else {
1957                 vp = Vec3r(_viewpoint);
1958         }
1959         Vec3r u(vp - center);
1960         real raylength = u.norm();
1961         u.normalize();
1962 #if 0
1963         if (_global.debug & G_DEBUG_FREESTYLE) {
1964                 cout << "grid origin " << iGrid->getOrigin().x() << "," << iGrid->getOrigin().y() << ","
1965                      << iGrid->getOrigin().z() << endl;
1966                 cout << "center " << center.x() << "," << center.y() << "," << center.z() << endl;
1967         }
1968 #endif
1969
1970         iGrid->castRay(center, vp, occluders, timestamp);
1971
1972         WFace *face = NULL;
1973         if (fe->isSmooth()) {
1974                 FEdgeSmooth *fes = dynamic_cast<FEdgeSmooth *>(fe);
1975                 face = (WFace *)fes->face();
1976         }
1977         vector<WVertex *> faceVertices;
1978         WVertex::incoming_edge_iterator ie;
1979
1980         WFace *oface;
1981         bool skipFace;
1982         OccludersSet::iterator p, pend;
1983         if (face)
1984                 face->RetrieveVertexList(faceVertices);
1985
1986         for (p = occluders.begin(), pend = occluders.end(); p != pend; p++) {
1987                 // If we're dealing with an exact silhouette, check whether we must take care of this occluder of not.
1988                 // (Indeed, we don't consider the occluders that share at least one vertex with the face containing this edge).
1989                 //-----------
1990                 oface = (WFace *)(*p)->userdata;
1991 #if LOGGING
1992                 if (_global.debug & G_DEBUG_FREESTYLE) {
1993                         cout << "\t\tEvaluating intersection for occluder " << ((*p)->getVertices())[0] <<
1994                                 ((*p)->getVertices())[1] << ((*p)->getVertices())[2] << endl << "\t\t\tand ray " << vp <<
1995                                 " * " << u << " (center " << center << ")" << endl;
1996                 }
1997 #endif
1998                 Vec3r v1(((*p)->getVertices())[0]);
1999                 Vec3r normal((*p)->getNormal());
2000                 real d = -(v1 * normal);
2001                 real t, t_u, t_v;
2002
2003 #if LOGGING
2004                 if (_global.debug & G_DEBUG_FREESTYLE) {
2005                         cout << "\t\tp:  " << ((*p)->getVertices())[0] << ((*p)->getVertices())[1] << ((*p)->getVertices())[2] <<
2006                                 ", norm: " << (*p)->getNormal() << endl;
2007                 }
2008 #endif
2009
2010                 if (face) {
2011 #if LOGGING
2012                         if (_global.debug & G_DEBUG_FREESTYLE) {
2013                                 cout << "\t\tDetermining face adjacency...";
2014                         }
2015 #endif
2016                         skipFace = false;
2017
2018                         if (face == oface) {
2019 #if LOGGING
2020                                 if (_global.debug & G_DEBUG_FREESTYLE) {
2021                                         cout << "  Rejecting occluder for face concurrency." << endl;
2022                                 }
2023 #endif
2024                                 continue;
2025                         }
2026
2027                         for (vector<WVertex*>::iterator fv = faceVertices.begin(), fvend = faceVertices.end();
2028                              fv != fvend;
2029                              ++fv)
2030                         {
2031                                 if ((*fv)->isBoundary())
2032                                         continue;
2033
2034                                 WVertex::incoming_edge_iterator iebegin = (*fv)->incoming_edges_begin();
2035                                 WVertex::incoming_edge_iterator ieend = (*fv)->incoming_edges_end();
2036                                 for (ie = iebegin; ie != ieend; ++ie) {
2037                                         if ((*ie) == 0)
2038                                                 continue;
2039
2040                                         WFace *sface = (*ie)->GetbFace();
2041                                         //WFace *sfacea = (*ie)->GetaFace();
2042                                         //if ((sface == oface) || (sfacea == oface)) {
2043                                         if (sface == oface) {
2044                                                 skipFace = true;
2045                                                 break;
2046                                         }
2047                                 }
2048                                 if (skipFace)
2049                                         break;
2050                         }
2051                         if (skipFace) {
2052 #if LOGGING
2053                                 if (_global.debug & G_DEBUG_FREESTYLE) {
2054                                         cout << "  Rejecting occluder for face adjacency." << endl;
2055                                 }
2056 #endif
2057                                 continue;
2058                         }
2059                 }
2060                 else {
2061                         // check whether the edge and the polygon plane are coincident:
2062                         //-------------------------------------------------------------
2063                         //first let us compute the plane equation.
2064
2065                         if (GeomUtils::COINCIDENT == GeomUtils::intersectRayPlane(origin, edge, normal, d, t, epsilon)) {
2066 #if LOGGING
2067                                 if (_global.debug & G_DEBUG_FREESTYLE) {
2068                                         cout << "\t\tRejecting occluder for target coincidence." << endl;
2069                                 }
2070 #endif
2071                                 continue;
2072                         }
2073                 }
2074
2075                 if ((*p)->rayIntersect(center, u, t, t_u, t_v)) {
2076 #if LOGGING
2077                         if (_global.debug & G_DEBUG_FREESTYLE) {
2078                                 cout << "\t\tRay " << vp << " * " << u << " intersects at time " << t << " (raylength is " <<
2079                                         raylength << ")" << endl;
2080                                 cout << "\t\t(u * normal) == " << (u * normal) << " for normal " << normal << endl;
2081                         }
2082 #endif
2083                         if (fabs(u * normal) > 0.0001) {
2084                                 if ((t>0.0) && (t<raylength)) {
2085 #if LOGGING
2086                                         if (_global.debug & G_DEBUG_FREESTYLE) {
2087                                                 cout << "\t\tIs occluder" << endl;
2088                                         }
2089 #endif
2090                                         WFace *f = (WFace *)((*p)->userdata);
2091                                         ViewShape *vshape = _ViewMap->viewShape(f->GetVertex(0)->shape()->GetId());
2092                                         oOccluders.insert(vshape);
2093                                         ++qi;
2094                                         if (!_EnableQI)
2095                                                 break;
2096                                 }
2097                         }
2098                 }
2099         }
2100
2101         // Find occludee
2102         FindOccludee(fe, iGrid, epsilon, oaPolygon, timestamp, u, center, origin, edge, faceVertices);
2103
2104         return qi;
2105 }
2106
2107 void ViewMapBuilder::ComputeIntersections(ViewMap *ioViewMap, intersection_algo iAlgo, real epsilon)
2108 {
2109         switch (iAlgo) {
2110                 case sweep_line:
2111                         ComputeSweepLineIntersections(ioViewMap, epsilon);
2112                         break;
2113                 default:
2114                         break;
2115         }
2116 #if 0
2117         if (_global.debug & G_DEBUG_FREESTYLE) {
2118                 ViewMap::viewvertices_container& vvertices = ioViewMap->ViewVertices();
2119                 for (ViewMap::viewvertices_container::iterator vv = vvertices.begin(), vvend = vvertices.end();
2120                      vv != vvend; ++vv)
2121                 {
2122                         if ((*vv)->getNature() == Nature::T_VERTEX) {
2123                                 TVertex *tvertex = (TVertex *)(*vv);
2124                                 cout << "TVertex " << tvertex->getId() << " has :" << endl;
2125                                 cout << "FrontEdgeA: " << tvertex->frontEdgeA().first << endl;
2126                                 cout << "FrontEdgeB: " << tvertex->frontEdgeB().first << endl;
2127                                 cout << "BackEdgeA: " << tvertex->backEdgeA().first << endl;
2128                                 cout << "BackEdgeB: " << tvertex->backEdgeB().first << endl << endl;
2129                         }
2130                 }
2131         }
2132 #endif
2133 }
2134
2135 struct less_SVertex2D : public binary_function<SVertex *, SVertex *, bool>
2136 {
2137         real epsilon;
2138
2139         less_SVertex2D(real eps) : binary_function<SVertex *, SVertex *, bool>()
2140         {
2141                 epsilon = eps;
2142         }
2143
2144         bool operator()(SVertex *x, SVertex *y)
2145         {
2146                 Vec3r A = x->point2D();
2147                 Vec3r B = y->point2D();
2148                 for (unsigned int i = 0; i < 3; i++) {
2149                         if ((fabs(A[i] - B[i])) < epsilon)
2150                                 continue;
2151                         if (A[i] < B[i])
2152                                 return true;
2153                         if (A[i] > B[i])
2154                                 return false;
2155                 }
2156                 return false;
2157         }
2158 };
2159
2160 typedef Segment<FEdge *, Vec3r> segment;
2161 typedef Intersection<segment> intersection;
2162
2163 struct less_Intersection : public binary_function<intersection *, intersection *, bool>
2164 {
2165         segment *edge;
2166
2167         less_Intersection(segment *iEdge) : binary_function<intersection *, intersection *, bool>()
2168         {
2169                 edge = iEdge;
2170         }
2171
2172         bool operator()(intersection *x, intersection *y)
2173         {
2174                 real tx = x->getParameter(edge);
2175                 real ty = y->getParameter(edge);
2176                 if (tx > ty)
2177                         return true;
2178                 return false;
2179         }
2180 };
2181
2182 struct silhouette_binary_rule : public binary_rule<segment, segment>
2183 {
2184         silhouette_binary_rule() : binary_rule<segment, segment>() {}
2185
2186         virtual bool operator()(segment& s1, segment& s2)
2187         {
2188                 FEdge *f1 = s1.edge();
2189                 FEdge *f2 = s2.edge();
2190
2191                 if ((!(((f1)->getNature() & Nature::SILHOUETTE) || ((f1)->getNature() & Nature::BORDER))) &&
2192                     (!(((f2)->getNature() & Nature::SILHOUETTE) || ((f2)->getNature() & Nature::BORDER))))
2193                 {
2194                         return false;
2195                 }
2196
2197                 return true;
2198         }
2199 };
2200
2201 void ViewMapBuilder::ComputeSweepLineIntersections(ViewMap *ioViewMap, real epsilon)
2202 {
2203         vector<SVertex *>& svertices = ioViewMap->SVertices();
2204         bool progressBarDisplay = false;
2205         unsigned sVerticesSize = svertices.size();
2206         unsigned fEdgesSize = ioViewMap->FEdges().size();
2207 #if 0
2208         if (_global.debug & G_DEBUG_FREESTYLE) {
2209                 ViewMap::fedges_container& fedges = ioViewMap->FEdges();
2210                 for (ViewMap::fedges_container::const_iterator f = fedges.begin(), end = fedges.end(); f != end; ++f) {
2211                         cout << (*f)->aMaterialIndex() << "-" << (*f)->bMaterialIndex() << endl;
2212                 }
2213         }
2214 #endif
2215         unsigned progressBarStep = 0;
2216
2217         if (_pProgressBar != NULL && fEdgesSize > gProgressBarMinSize) {
2218                 unsigned progressBarSteps = min(gProgressBarMaxSteps, sVerticesSize);
2219                 progressBarStep = sVerticesSize / progressBarSteps;
2220                 _pProgressBar->reset();
2221                 _pProgressBar->setLabelText("Computing Sweep Line Intersections");
2222                 _pProgressBar->setTotalSteps(progressBarSteps);
2223                 _pProgressBar->setProgress(0);
2224                 progressBarDisplay = true;
2225         }
2226
2227         unsigned counter = progressBarStep;
2228
2229         sort(svertices.begin(), svertices.end(), less_SVertex2D(epsilon));
2230
2231         SweepLine<FEdge *, Vec3r> SL;
2232
2233         vector<FEdge *>& ioEdges = ioViewMap->FEdges();
2234
2235         vector<segment*> segments;
2236
2237         vector<FEdge*>::iterator fe, fend;
2238
2239         for (fe = ioEdges.begin(), fend = ioEdges.end(); fe != fend; fe++) {
2240                 segment *s = new segment((*fe), (*fe)->vertexA()->point2D(), (*fe)->vertexB()->point2D());
2241                 (*fe)->userdata = s;
2242                 segments.push_back(s);
2243         }
2244
2245         vector<segment*> vsegments;
2246         for (vector<SVertex*>::iterator sv = svertices.begin(), svend = svertices.end(); sv != svend; sv++) {
2247                 if (_pRenderMonitor && _pRenderMonitor->testBreak())
2248                         break;
2249
2250                 const vector<FEdge*>& vedges = (*sv)->fedges();
2251
2252                 for (vector<FEdge *>::const_iterator sve = vedges.begin(), sveend = vedges.end(); sve != sveend; sve++) {
2253                         vsegments.push_back((segment *)((*sve)->userdata));
2254                 }
2255
2256                 Vec3r evt((*sv)->point2D());
2257                 silhouette_binary_rule sbr;
2258                 SL.process(evt, vsegments, sbr, epsilon);
2259
2260                 if (progressBarDisplay) {
2261                         counter--;
2262                         if (counter <= 0) {
2263                                 counter = progressBarStep;
2264                                 _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
2265                         }
2266                 }
2267                 vsegments.clear();
2268         }
2269
2270         if (_pRenderMonitor && _pRenderMonitor->testBreak()) {
2271                 // delete segments
2272                 if (!segments.empty()) {
2273                         vector<segment*>::iterator s, send;
2274                         for (s = segments.begin(), send = segments.end(); s != send; s++) {
2275                                 delete *s;
2276                         }
2277                 }
2278                 return;
2279         }
2280
2281         // reset userdata:
2282         for (fe = ioEdges.begin(), fend = ioEdges.end(); fe != fend; fe++)
2283                 (*fe)->userdata = NULL;
2284
2285         // list containing the new edges resulting from splitting operations. 
2286         vector<FEdge*> newEdges;
2287
2288         // retrieve the intersected edges:
2289         vector<segment*>& iedges = SL.intersectedEdges();
2290         // retrieve the intersections:
2291         vector<intersection*>& intersections = SL.intersections();
2292
2293         int id = 0;
2294         // create a view vertex for each intersection and linked this one with the intersection object
2295         vector<intersection*>::iterator i, iend;
2296         for (i = intersections.begin(), iend = intersections.end(); i != iend; i++) {
2297                 FEdge *fA = (*i)->EdgeA->edge();
2298                 FEdge *fB = (*i)->EdgeB->edge();
2299
2300                 Vec3r A1 = fA->vertexA()->point3D();
2301                 Vec3r A2 = fA->vertexB()->point3D();
2302                 Vec3r B1 = fB->vertexA()->point3D();
2303                 Vec3r B2 = fB->vertexB()->point3D();
2304
2305                 Vec3r a1 = fA->vertexA()->point2D();
2306                 Vec3r a2 = fA->vertexB()->point2D();
2307                 Vec3r b1 = fB->vertexA()->point2D();
2308                 Vec3r b2 = fB->vertexB()->point2D();
2309
2310                 real ta = (*i)->tA;
2311                 real tb = (*i)->tB;
2312
2313                 if ((ta < -epsilon) || (ta > 1 + epsilon))
2314                         cerr << "Warning: 2D intersection out of range for edge " << fA->vertexA()->getId() << " - " <<
2315                                 fA->vertexB()->getId() << endl;
2316
2317                 if ((tb < -epsilon) || (tb > 1 + epsilon))
2318                         cerr << "Warning: 2D intersection out of range for edge " << fB->vertexA()->getId() << " - " <<
2319                                 fB->vertexB()->getId() << endl;
2320
2321                 real Ta = SilhouetteGeomEngine::ImageToWorldParameter(fA, ta);
2322                 real Tb = SilhouetteGeomEngine::ImageToWorldParameter(fB, tb);
2323
2324                 if ((Ta < -epsilon) || (Ta > 1 + epsilon))
2325                         cerr << "Warning: 3D intersection out of range for edge " << fA->vertexA()->getId() << " - " <<
2326                                 fA->vertexB()->getId() << endl;
2327
2328                 if ((Tb < -epsilon) || (Tb > 1 + epsilon))
2329                         cerr << "Warning: 3D intersection out of range for edge " << fB->vertexA()->getId() << " - " <<
2330                                 fB->vertexB()->getId() << endl;
2331
2332 #if 0
2333                 if (_global.debug & G_DEBUG_FREESTYLE) {
2334                         if ((Ta < -epsilon) || (Ta > 1 + epsilon) || (Tb < -epsilon) || (Tb > 1 + epsilon)) {
2335                                 printf("ta %.12e\n", ta);
2336                                 printf("tb %.12e\n", tb);
2337                                 printf("a1 %e, %e -- a2 %e, %e\n", a1[0], a1[1], a2[0], a2[1]);
2338                                 printf("b1 %e, %e -- b2 %e, %e\n", b1[0], b1[1], b2[0], b2[1]);
2339                                 //printf("line([%e, %e], [%e, %e]);\n", a1[0], a2[0], a1[1], a2[1]);
2340                                 //printf("line([%e, %e], [%e, %e]);\n", b1[0], b2[0], b1[1], b2[1]);
2341                                 if ((Ta < -epsilon) || (Ta > 1 + epsilon))
2342                                         printf("Ta %.12e\n", Ta);
2343                                 if ((Tb < -epsilon) || (Tb > 1 + epsilon))
2344                                         printf("Tb %.12e\n", Tb);
2345                                 printf("A1 %e, %e, %e -- A2 %e, %e, %e\n", A1[0], A1[1], A1[2], A2[0], A2[1], A2[2]);
2346                                 printf("B1 %e, %e, %e -- B2 %e, %e, %e\n", B1[0], B1[1], B1[2], B2[0], B2[1], B2[2]);
2347                         }
2348                 }
2349 #endif
2350
2351                 TVertex *tvertex = ioViewMap->CreateTVertex(Vec3r(A1 + Ta * (A2 - A1)), Vec3r(a1 + ta * (a2 - a1)), fA,
2352                                                             Vec3r(B1 + Tb * (B2 - B1)), Vec3r(b1 + tb * (b2 - b1)), fB, id);
2353
2354                 (*i)->userdata = tvertex;
2355                 ++id;
2356         }
2357
2358         progressBarStep = 0;
2359
2360         if (progressBarDisplay) {
2361                 unsigned iEdgesSize = iedges.size();
2362                 unsigned progressBarSteps = min(gProgressBarMaxSteps, iEdgesSize);
2363                 progressBarStep = iEdgesSize / progressBarSteps;
2364                 _pProgressBar->reset();
2365                 _pProgressBar->setLabelText("Splitting intersected edges");
2366                 _pProgressBar->setTotalSteps(progressBarSteps);
2367                 _pProgressBar->setProgress(0);
2368         }
2369
2370         counter = progressBarStep;
2371
2372         vector<TVertex*> edgeVVertices;
2373         vector<ViewEdge*> newVEdges;
2374         vector<segment*>::iterator s, send;
2375         for (s = iedges.begin(), send = iedges.end(); s != send; s++) {
2376                 edgeVVertices.clear();
2377                 newEdges.clear();
2378                 newVEdges.clear();
2379
2380                 FEdge *fedge = (*s)->edge();
2381                 ViewEdge *vEdge = fedge->viewedge();
2382                 ViewShape *shape = vEdge->viewShape();
2383
2384                 vector<intersection*>& eIntersections = (*s)->intersections();
2385                 // we first need to sort these intersections from farther to closer to A
2386                 sort(eIntersections.begin(), eIntersections.end(), less_Intersection(*s));
2387                 for (i = eIntersections.begin(), iend = eIntersections.end(); i != iend; i++)
2388                         edgeVVertices.push_back((TVertex *)(*i)->userdata);
2389
2390                 shape->SplitEdge(fedge, edgeVVertices, ioViewMap->FEdges(), ioViewMap->ViewEdges());
2391
2392                 if (progressBarDisplay) {
2393                         counter--;
2394                         if (counter <= 0) {
2395                                 counter = progressBarStep;
2396                                 _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
2397                         }
2398                 }
2399         }
2400
2401         // reset userdata:
2402         for (fe = ioEdges.begin(), fend = ioEdges.end(); fe != fend; fe++)
2403                 (*fe)->userdata = NULL;
2404
2405         // delete segments
2406         if (!segments.empty()) {
2407                 for (s = segments.begin(), send = segments.end(); s != send; s++) {
2408                         delete *s;
2409                 }
2410         }
2411 }
2412
2413 } /* namespace Freestyle */