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