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