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