Fixed a bug in SilhouetteGeomEngine::ImageToWorldParameter() that caused
[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 <algorithm>
23 #include "ViewMapBuilder.h"
24
25 using namespace std;
26
27 ViewMap* ViewMapBuilder::BuildViewMap(WingedEdge& we, visibility_algo iAlgo, real epsilon) {
28   _ViewMap = new ViewMap;
29   _currentId = 1;
30   _currentFId = 0;
31   _currentSVertexId = 0;
32   
33   // Builds initial view edges
34   computeInitialViewEdges(we);
35   
36   // Detects cusps
37   computeCusps(_ViewMap); 
38   
39   // Compute intersections
40   ComputeIntersections(_ViewMap, sweep_line, epsilon);
41   
42   // Compute visibility
43   ComputeEdgesVisibility(_ViewMap, iAlgo, _Grid, epsilon);
44   
45   return _ViewMap;
46 }
47
48 void ViewMapBuilder::computeInitialViewEdges(WingedEdge& we)
49 {
50   vector<WShape*> wshapes = we.getWShapes();
51   SShape* psShape;
52
53   for (vector<WShape*>::const_iterator it = wshapes.begin();
54        it != wshapes.end();
55        it++) {
56     // create the embedding
57     psShape = new SShape;
58     psShape->setId((*it)->GetId());
59     psShape->setFrsMaterials((*it)->frs_materials()); // FIXME
60
61     // create the view shape
62     ViewShape * vshape = new ViewShape(psShape);
63     // add this view shape to the view map:
64     _ViewMap->AddViewShape(vshape);
65   
66     _pViewEdgeBuilder->setCurrentViewId(_currentId); // we want to number the view edges in a unique way for the while scene.
67     _pViewEdgeBuilder->setCurrentFId(_currentFId); // we want to number the feature edges in a unique way for the while scene.
68     _pViewEdgeBuilder->setCurrentSVertexId(_currentFId); // we want to number the SVertex in a unique way for the while scene.
69     _pViewEdgeBuilder->BuildViewEdges(dynamic_cast<WXShape*>(*it), vshape, 
70                                       _ViewMap->ViewEdges(), 
71                                       _ViewMap->ViewVertices(), 
72                                       _ViewMap->FEdges(), 
73                                       _ViewMap->SVertices());
74   
75     _currentId = _pViewEdgeBuilder->currentViewId()+1;
76     _currentFId = _pViewEdgeBuilder->currentFId()+1;
77     _currentSVertexId = _pViewEdgeBuilder->currentSVertexId()+1;
78
79     psShape->ComputeBBox();
80   }
81 }
82
83 void ViewMapBuilder::computeCusps(ViewMap *ioViewMap){
84   vector<ViewVertex*> newVVertices;
85   vector<ViewEdge*> newVEdges;
86   ViewMap::viewedges_container& vedges = ioViewMap->ViewEdges();
87   ViewMap::viewedges_container::iterator ve=vedges.begin(), veend=vedges.end();
88   for(;
89   ve!=veend;
90   ++ve){
91     if((!((*ve)->getNature() & Nature::SILHOUETTE)) || (!((*ve)->fedgeA()->isSmooth())))
92       continue;
93     FEdge *fe = (*ve)->fedgeA();
94     FEdge * fefirst = fe;
95     bool first = true;
96     bool positive = true;
97     do{
98       FEdgeSmooth * fes = dynamic_cast<FEdgeSmooth*>(fe);
99       Vec3r A((fes)->vertexA()->point3d());
100       Vec3r B((fes)->vertexB()->point3d());
101       Vec3r AB(B-A);
102       AB.normalize();
103       Vec3r m((A+B)/2.0);
104       Vec3r crossP(AB^(fes)->normal()); 
105       crossP.normalize();
106       Vec3r viewvector(m-_viewpoint);
107       viewvector.normalize();
108       if(first){
109         if(((crossP)*(viewvector)) > 0)
110           positive = true;
111         else
112           positive = false;
113         first = false;
114       }
115       // If we're in a positive part, we need 
116       // a stronger negative value to change
117       NonTVertex *cusp = 0;
118       if(positive){
119         if(((crossP)*(viewvector)) < -0.1){
120           // state changes
121           positive = false;
122           // creates and insert cusp
123           cusp = dynamic_cast<NonTVertex*>(ioViewMap->InsertViewVertex(fes->vertexA(), newVEdges));
124           if(cusp!=0)
125             cusp->setNature(cusp->getNature()|Nature::CUSP);
126         }
127           
128       }else{
129         // If we're in a negative part, we need 
130         // a stronger negative value to change
131         if(((crossP)*(viewvector)) > 0.1){
132           positive = true;
133           cusp = dynamic_cast<NonTVertex*>(ioViewMap->InsertViewVertex(fes->vertexA(), newVEdges));
134           if(cusp!=0)
135             cusp->setNature(cusp->getNature()|Nature::CUSP);
136         }
137       }
138       fe = fe->nextEdge();
139     }while((fe!=0) && (fe!=fefirst));
140   }
141   for(ve=newVEdges.begin(), veend=newVEdges.end();
142   ve!=veend;
143   ++ve){
144     (*ve)->viewShape()->AddEdge(*ve);
145     vedges.push_back(*ve);
146   }
147 }
148 void ViewMapBuilder::ComputeEdgesVisibility(ViewMap *ioViewMap, visibility_algo iAlgo,  Grid *iGrid, real epsilon)
149 {
150   if((iAlgo == ray_casting ||
151       iAlgo == ray_casting_fast ||
152       iAlgo == ray_casting_very_fast) && (NULL == iGrid))
153   {
154     cerr << "Error: can't cast ray, no grid defined" << endl;
155     return;
156   }
157
158   switch(iAlgo)
159   {
160   case ray_casting:
161     ComputeRayCastingVisibility(ioViewMap, iGrid, epsilon);
162     break;
163   case ray_casting_fast:
164     ComputeFastRayCastingVisibility(ioViewMap, iGrid, epsilon);
165     break;
166   case ray_casting_very_fast:
167     ComputeVeryFastRayCastingVisibility(ioViewMap, iGrid, epsilon);
168     break;
169   default:
170     break;
171   }
172 }
173
174 static const unsigned gProgressBarMaxSteps = 10;
175 static const unsigned gProgressBarMinSize = 2000;
176
177 void ViewMapBuilder::ComputeRayCastingVisibility(ViewMap *ioViewMap, Grid* iGrid, real epsilon)
178 {
179   vector<ViewEdge*>& vedges = ioViewMap->ViewEdges();
180   bool progressBarDisplay = false;
181   unsigned progressBarStep = 0;
182   unsigned vEdgesSize = vedges.size();
183   unsigned fEdgesSize = ioViewMap->FEdges().size();
184
185   if(_pProgressBar != NULL && fEdgesSize > gProgressBarMinSize) {
186     unsigned progressBarSteps = min(gProgressBarMaxSteps, vEdgesSize);
187     progressBarStep = vEdgesSize / progressBarSteps;
188     _pProgressBar->reset();
189     _pProgressBar->setLabelText("Computing Ray casting Visibility");
190     _pProgressBar->setTotalSteps(progressBarSteps);
191     _pProgressBar->setProgress(0);
192     progressBarDisplay = true;
193   }
194   
195   unsigned counter = progressBarStep;
196   FEdge * fe, *festart;
197   int nSamples = 0;
198   vector<Polygon3r*> aFaces;
199   Polygon3r *aFace = 0;
200   unsigned tmpQI = 0;
201   unsigned qiClasses[256];
202   unsigned maxIndex, maxCard;
203   unsigned qiMajority;
204   static unsigned timestamp = 1;
205   for(vector<ViewEdge*>::iterator ve=vedges.begin(), veend=vedges.end();
206   ve!=veend;
207   ve++)
208   {
209     festart = (*ve)->fedgeA();
210     fe = (*ve)->fedgeA();
211     qiMajority = 1;
212     do {
213        qiMajority++;
214        fe = fe->nextEdge();
215     } while (fe && fe != festart);
216     qiMajority >>= 1;
217
218     tmpQI = 0;
219     maxIndex = 0;
220     maxCard = 0;
221     nSamples = 0;
222     fe = (*ve)->fedgeA();
223     memset(qiClasses, 0, 256 * sizeof(*qiClasses));
224     set<ViewShape*> occluders;
225     do
226     {
227       if((maxCard < qiMajority)) {
228         tmpQI = ComputeRayCastingVisibility(fe, iGrid, epsilon, occluders, &aFace, timestamp++);
229         
230         if(tmpQI >= 256)
231           cerr << "Warning: too many occluding levels" << endl;
232
233         if (++qiClasses[tmpQI] > maxCard) {
234           maxCard = qiClasses[tmpQI];
235           maxIndex = tmpQI;
236         }
237       }
238       else
239         FindOccludee(fe, iGrid, epsilon, &aFace, timestamp++);
240
241       if(aFace) { 
242         fe->setaFace(*aFace);
243         aFaces.push_back(aFace);
244         fe->setOccludeeEmpty(false);
245       }
246       else
247         fe->setOccludeeEmpty(true);
248
249       ++nSamples;
250       fe = fe->nextEdge();
251     }
252     while((maxCard < qiMajority) && (0!=fe) && (fe!=festart));
253
254     // ViewEdge
255     // qi --
256     (*ve)->setQI(maxIndex);
257     // occluders --
258     for(set<ViewShape*>::iterator o=occluders.begin(), oend=occluders.end();
259     o!=oend;
260     ++o)
261       (*ve)->AddOccluder((*o));
262     // occludee --
263     if(!aFaces.empty())
264     {
265       if(aFaces.size() <= (float)nSamples/2.f)
266       {
267         (*ve)->setaShape(0);
268       }
269       else
270       {
271         vector<Polygon3r*>::iterator p = aFaces.begin();
272         WFace * wface = (WFace*)((*p)->userdata);
273         ViewShape *vshape = ioViewMap->viewShape(wface->GetVertex(0)->shape()->GetId());
274         ++p;
275         (*ve)->setaShape(vshape);
276       }
277     }
278     
279     if(progressBarDisplay) {  
280       counter--;
281       if (counter <= 0) {
282         counter = progressBarStep;
283         _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
284       }
285     }   
286     aFaces.clear();
287   }    
288 }
289
290 void ViewMapBuilder::ComputeFastRayCastingVisibility(ViewMap *ioViewMap, Grid* iGrid, real epsilon)
291 {
292   vector<ViewEdge*>& vedges = ioViewMap->ViewEdges();
293   bool progressBarDisplay = false;
294   unsigned progressBarStep = 0;
295   unsigned vEdgesSize = vedges.size();
296   unsigned fEdgesSize = ioViewMap->FEdges().size();
297
298   if(_pProgressBar != NULL && fEdgesSize > gProgressBarMinSize) {
299     unsigned progressBarSteps = min(gProgressBarMaxSteps, vEdgesSize);
300     progressBarStep = vEdgesSize / progressBarSteps;
301     _pProgressBar->reset();
302     _pProgressBar->setLabelText("Computing Ray casting Visibility");
303     _pProgressBar->setTotalSteps(progressBarSteps);
304     _pProgressBar->setProgress(0);
305     progressBarDisplay = true;
306   }
307
308   unsigned counter = progressBarStep;  
309   FEdge * fe, *festart;
310   unsigned nSamples = 0;
311   vector<Polygon3r*> aFaces;
312   Polygon3r *aFace = 0;
313   unsigned tmpQI = 0;
314   unsigned qiClasses[256];
315   unsigned maxIndex, maxCard;
316   unsigned qiMajority;
317   static unsigned timestamp = 1;
318   bool even_test;
319   for(vector<ViewEdge*>::iterator ve=vedges.begin(), veend=vedges.end();
320   ve!=veend;
321   ve++)
322   {
323     festart = (*ve)->fedgeA();
324     fe = (*ve)->fedgeA();
325     qiMajority = 1;
326     do {
327        qiMajority++;
328        fe = fe->nextEdge();
329     } while (fe && fe != festart);
330     if (qiMajority >= 4)
331       qiMajority >>= 2;
332     else
333       qiMajority = 1;
334
335     set<ViewShape*> occluders;
336
337     even_test = true;
338     maxIndex = 0;
339     maxCard = 0;
340     nSamples = 0;
341     memset(qiClasses, 0, 256 * sizeof(*qiClasses));
342     fe = (*ve)->fedgeA();
343     do
344     {
345       if (even_test)
346       {
347         if((maxCard < qiMajority)) {
348           tmpQI = ComputeRayCastingVisibility(fe, iGrid, epsilon, occluders, &aFace, timestamp++);
349           
350           if(tmpQI >= 256)
351             cerr << "Warning: too many occluding levels" << endl;
352
353           if (++qiClasses[tmpQI] > maxCard) {
354             maxCard = qiClasses[tmpQI];
355             maxIndex = tmpQI;
356           }
357         }
358         else
359           FindOccludee(fe, iGrid, epsilon, &aFace, timestamp++);
360
361         if(aFace)
362         { 
363           fe->setaFace(*aFace);
364           aFaces.push_back(aFace);
365         } 
366         ++nSamples;
367         even_test = false;
368       }
369       else
370         even_test = true;
371       fe = fe->nextEdge();
372     } while ((maxCard < qiMajority) && (0!=fe) && (fe!=festart));
373
374     (*ve)->setQI(maxIndex);
375
376     if(!aFaces.empty())
377     {
378       if(aFaces.size() < nSamples / 2)
379       {
380         (*ve)->setaShape(0);
381       }
382       else
383       {
384         vector<Polygon3r*>::iterator p = aFaces.begin();
385         WFace * wface = (WFace*)((*p)->userdata);
386         ViewShape *vshape = ioViewMap->viewShape(wface->GetVertex(0)->shape()->GetId());
387         ++p;
388         //        for(;
389         //        p!=pend;
390         //        ++p)
391         //        { 
392         //          WFace *f = (WFace*)((*p)->userdata);
393         //          ViewShape *vs = ioViewMap->viewShape(f->GetVertex(0)->shape()->GetId());
394         //          if(vs != vshape)
395         //          { 
396         //            sameShape = false;
397         //            break;
398         //          } 
399         //        } 
400         //        if(sameShape)
401         (*ve)->setaShape(vshape);
402       }
403     }
404     
405     //(*ve)->setaFace(aFace);
406     
407     if(progressBarDisplay) {  
408       counter--;
409       if (counter <= 0) {
410         counter = progressBarStep;
411         _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
412       }
413     }
414     aFaces.clear();
415   }
416 }
417
418 void ViewMapBuilder::ComputeVeryFastRayCastingVisibility(ViewMap *ioViewMap, Grid* iGrid, real epsilon)
419 {
420   vector<ViewEdge*>& vedges = ioViewMap->ViewEdges();
421   bool progressBarDisplay = false;
422   unsigned progressBarStep = 0;
423   unsigned vEdgesSize = vedges.size();
424   unsigned fEdgesSize = ioViewMap->FEdges().size();
425
426   if(_pProgressBar != NULL && fEdgesSize > gProgressBarMinSize) {
427     unsigned progressBarSteps = min(gProgressBarMaxSteps, vEdgesSize);
428     progressBarStep = vEdgesSize / progressBarSteps;
429     _pProgressBar->reset();
430     _pProgressBar->setLabelText("Computing Ray casting Visibility");
431     _pProgressBar->setTotalSteps(progressBarSteps);
432     _pProgressBar->setProgress(0);
433     progressBarDisplay = true;
434   }
435
436   unsigned counter = progressBarStep;
437   FEdge* fe;
438   unsigned qi = 0;
439   Polygon3r *aFace = 0;
440   static unsigned timestamp = 1;
441   for(vector<ViewEdge*>::iterator ve=vedges.begin(), veend=vedges.end();
442   ve!=veend;
443   ve++)
444   {
445     set<ViewShape*> occluders;
446
447     fe = (*ve)->fedgeA();
448     qi = ComputeRayCastingVisibility(fe, iGrid, epsilon, occluders, &aFace, timestamp++);
449     if(aFace)
450     { 
451       fe->setaFace(*aFace);
452       WFace * wface = (WFace*)(aFace->userdata);
453       ViewShape *vshape = ioViewMap->viewShape(wface->GetVertex(0)->shape()->GetId());
454       (*ve)->setaShape(vshape);
455     } 
456     else
457     {
458       (*ve)->setaShape(0);
459     }
460
461     (*ve)->setQI(qi);
462     
463     if(progressBarDisplay) {  
464       counter--;
465       if (counter <= 0) {
466         counter = progressBarStep;
467         _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
468       }
469     }
470   }  
471 }
472
473
474 void ViewMapBuilder::FindOccludee(FEdge *fe, Grid* iGrid, real epsilon, Polygon3r** oaPolygon, unsigned timestamp, 
475                                   Vec3r& u, Vec3r& A, Vec3r& origin, Vec3r& edge, vector<WVertex*>& faceVertices)
476 {
477   WFace *face = 0;
478   if(fe->isSmooth()){
479     FEdgeSmooth * fes = dynamic_cast<FEdgeSmooth*>(fe);
480     face = (WFace*)fes->face();
481   }
482   OccludersSet occluders;
483   WFace * oface;
484   bool skipFace;
485   
486   WVertex::incoming_edge_iterator ie;
487   OccludersSet::iterator p, pend;
488   
489   *oaPolygon = 0;
490   if(((fe)->getNature() & Nature::SILHOUETTE) || ((fe)->getNature() & Nature::BORDER))
491   {
492     occluders.clear();
493     // we cast a ray from A in the same direction but looking behind
494     Vec3r v(-u[0],-u[1],-u[2]);
495     iGrid->castInfiniteRay(A, v, occluders, timestamp);
496     
497     bool noIntersection = true;
498     real mint=FLT_MAX;
499     // we met some occluders, let us fill the aShape field 
500     // with the first intersected occluder
501     for(p=occluders.begin(),pend=occluders.end();
502     p!=pend;
503     p++)
504     { 
505       // check whether the edge and the polygon plane are coincident:
506       //-------------------------------------------------------------
507       //first let us compute the plane equation.
508       oface = (WFace*)(*p)->userdata;
509       Vec3r v1(((*p)->getVertices())[0]);
510       Vec3r normal((*p)->getNormal());
511       real d = -(v1 * normal);
512       real t,t_u,t_v;
513       
514       if(0 != face)
515       { 
516         skipFace = false;
517         
518         if(face == oface)
519           continue;
520         
521         if(faceVertices.empty())
522           continue;
523         
524         for(vector<WVertex*>::iterator fv=faceVertices.begin(), fvend=faceVertices.end();
525         fv!=fvend;
526         ++fv)
527         { 
528           if((*fv)->isBoundary())
529             continue;
530           WVertex::incoming_edge_iterator iebegin=(*fv)->incoming_edges_begin();
531           WVertex::incoming_edge_iterator ieend=(*fv)->incoming_edges_end();
532           for(ie=iebegin;ie!=ieend; ++ie)
533           {  
534             if((*ie) == 0)
535               continue;
536             
537             WFace * sface = (*ie)->GetbFace();
538             if(sface == oface)
539             { 
540               skipFace = true;
541               break;
542             } 
543           }  
544           if(skipFace)
545             break;
546         } 
547         if(skipFace)
548           continue;
549       }  
550       else
551       {
552         if(GeomUtils::COINCIDENT == GeomUtils::intersectRayPlane(origin, edge, normal, d, t, epsilon))
553           continue;
554       }
555       if((*p)->rayIntersect(A, v, t,t_u,t_v))
556       {
557         if (fabs(v * normal) > 0.0001)
558           if ((t>0.0)) // && (t<1.0))
559           {
560             if (t<mint)
561             {
562               *oaPolygon = (*p);
563               mint = t;
564               noIntersection = false;
565               fe->setOccludeeIntersection(Vec3r(A+t*v));
566             }
567           }
568       }  
569     }
570     
571     if(noIntersection)
572       *oaPolygon = 0;
573   }
574 }
575
576 void ViewMapBuilder::FindOccludee(FEdge *fe, Grid* iGrid, real epsilon, Polygon3r** oaPolygon, unsigned timestamp)
577 {
578   OccludersSet occluders;
579
580   Vec3r A;
581   Vec3r edge;
582   Vec3r origin;
583   A = Vec3r(((fe)->vertexA()->point3D() + (fe)->vertexB()->point3D())/2.0);
584   edge = Vec3r((fe)->vertexB()->point3D()-(fe)->vertexA()->point3D());
585   origin = Vec3r((fe)->vertexA()->point3D());
586   Vec3r u(_viewpoint-A);
587   u.normalize();
588   if(A < iGrid->getOrigin())
589     cerr << "Warning: point is out of the grid for fedge " << fe->getId().getFirst() << "-" << fe->getId().getSecond() << endl;
590
591   vector<WVertex*> faceVertices;
592      
593   WFace *face = 0;
594   if(fe->isSmooth()){
595     FEdgeSmooth * fes = dynamic_cast<FEdgeSmooth*>(fe);
596     face = (WFace*)fes->face();
597   }
598   if(0 != face)
599     face->RetrieveVertexList(faceVertices);
600   
601   return FindOccludee(fe,iGrid, epsilon, oaPolygon, timestamp, 
602                       u, A, origin, edge, faceVertices);
603 }
604
605 int ViewMapBuilder::ComputeRayCastingVisibility(FEdge *fe, Grid* iGrid, real epsilon, set<ViewShape*>& oOccluders,
606                                                 Polygon3r** oaPolygon, unsigned timestamp)
607 {
608   OccludersSet occluders;
609   int qi = 0;
610
611   Vec3r center;
612   Vec3r edge;
613   Vec3r origin;
614
615   center = fe->center3d();
616   edge = Vec3r(fe->vertexB()->point3D() - fe->vertexA()->point3D());
617   origin = Vec3r(fe->vertexA()->point3D());
618   //
619   //   // Is the edge outside the view frustum ?
620   Vec3r gridOrigin(iGrid->getOrigin());
621   Vec3r gridExtremity(iGrid->getOrigin()+iGrid->gridSize());
622   
623   if( (center.x() < gridOrigin.x()) || (center.y() < gridOrigin.y()) || (center.z() < gridOrigin.z())
624     ||(center.x() > gridExtremity.x()) || (center.y() > gridExtremity.y()) || (center.z() > gridExtremity.z())){
625      cerr << "Warning: point is out of the grid for fedge " << fe->getId() << endl;
626     //return 0;
627   }
628
629  
630   //  Vec3r A(fe->vertexA()->point2d());
631   //  Vec3r B(fe->vertexB()->point2d());
632   //  int viewport[4];
633   //  SilhouetteGeomEngine::retrieveViewport(viewport);
634   //  if( (A.x() < viewport[0]) || (A.x() > viewport[2]) || (A.y() < viewport[1]) || (A.y() > viewport[3])
635   //    ||(B.x() < viewport[0]) || (B.x() > viewport[2]) || (B.y() < viewport[1]) || (B.y() > viewport[3])){
636   //    cerr << "Warning: point is out of the grid for fedge " << fe->getId() << endl;
637   //    //return 0;
638   //  }
639
640   Vec3r u(_viewpoint - center);
641   real raylength = u.norm();
642   u.normalize();
643   //cout << "grid origin " << iGrid->getOrigin().x() << "," << iGrid->getOrigin().y() << "," << iGrid->getOrigin().z() << endl;
644   //cout << "center " << center.x() << "," << center.y() << "," << center.z() << endl;
645   
646   iGrid->castRay(center, Vec3r(_viewpoint), occluders, timestamp);
647
648   WFace *face = 0;
649   if(fe->isSmooth()){
650     FEdgeSmooth * fes = dynamic_cast<FEdgeSmooth*>(fe);
651     face = (WFace*)fes->face();
652   }
653   vector<WVertex*> faceVertices;
654   WVertex::incoming_edge_iterator ie;
655             
656   WFace * oface;
657   bool skipFace;
658   OccludersSet::iterator p, pend;
659   if(face)
660     face->RetrieveVertexList(faceVertices);
661
662   for(p=occluders.begin(),pend=occluders.end();
663       p!=pend;
664       p++)
665     { 
666       // If we're dealing with an exact silhouette, check whether 
667       // we must take care of this occluder of not.
668       // (Indeed, we don't consider the occluders that 
669       // share at least one vertex with the face containing 
670       // this edge).
671       //-----------
672       oface = (WFace*)(*p)->userdata;
673       Vec3r v1(((*p)->getVertices())[0]);
674       Vec3r normal((*p)->getNormal());
675       real d = -(v1 * normal);
676       real t, t_u, t_v;
677             
678       if(0 != face)
679         {
680           skipFace = false;
681               
682           if(face == oface)
683             continue;
684               
685               
686           for(vector<WVertex*>::iterator fv=faceVertices.begin(), fvend=faceVertices.end();
687               fv!=fvend;
688               ++fv)
689             {
690                 if((*fv)->isBoundary())
691                   continue;
692
693               WVertex::incoming_edge_iterator iebegin=(*fv)->incoming_edges_begin();
694               WVertex::incoming_edge_iterator ieend=(*fv)->incoming_edges_end();
695               for(ie=iebegin;ie!=ieend; ++ie)
696                 { 
697                   if((*ie) == 0)
698                     continue;
699                   
700                   WFace * sface = (*ie)->GetbFace();
701                   //WFace * sfacea = (*ie)->GetaFace();
702                   //if((sface == oface) || (sfacea == oface))
703                   if(sface == oface)
704                     {
705                       skipFace = true;
706                       break;
707                     }
708                 }
709               if(skipFace)
710                 break;
711             }
712           if(skipFace)
713             continue;
714         }
715       else
716         {
717           // check whether the edge and the polygon plane are coincident:
718           //-------------------------------------------------------------
719           //first let us compute the plane equation.
720            
721           if(GeomUtils::COINCIDENT == GeomUtils::intersectRayPlane(origin, edge, normal, d, t, epsilon))
722             continue;
723         }
724
725       if((*p)->rayIntersect(center, u, t, t_u, t_v))
726         {
727           if (fabs(u * normal) > 0.0001)
728             if ((t>0.0) && (t<raylength))
729               {
730                 WFace *f = (WFace*)((*p)->userdata);
731                 ViewShape *vshape = _ViewMap->viewShape(f->GetVertex(0)->shape()->GetId());
732                 oOccluders.insert(vshape);
733                 ++qi;
734     if(!_EnableQI)
735       break;
736               }
737         }
738     }
739
740   // Find occludee
741   FindOccludee(fe,iGrid, epsilon, oaPolygon, timestamp, 
742                u, center, edge, origin, faceVertices);
743
744   return qi;
745 }
746
747 void ViewMapBuilder::ComputeIntersections(ViewMap *ioViewMap, intersection_algo iAlgo, real epsilon)
748 {
749   switch(iAlgo)
750   {
751   case sweep_line:
752     ComputeSweepLineIntersections(ioViewMap, epsilon);
753     break;
754   default:
755     break;
756   }
757   ViewMap::viewvertices_container& vvertices = ioViewMap->ViewVertices();
758   for(ViewMap::viewvertices_container::iterator vv=vvertices.begin(), vvend=vvertices.end();
759   vv!=vvend;
760   ++vv)
761   {
762     if((*vv)->getNature() == Nature::T_VERTEX)
763     {
764       TVertex *tvertex = (TVertex*)(*vv);
765       cout << "TVertex " << tvertex->getId() << " has :" << endl;
766       cout << "FrontEdgeA: " << tvertex->frontEdgeA().first << endl;
767       cout << "FrontEdgeB: " << tvertex->frontEdgeB().first << endl;
768       cout << "BackEdgeA: " << tvertex->backEdgeA().first << endl;
769       cout << "BackEdgeB: " << tvertex->backEdgeB().first << endl << endl;
770     }
771   }
772 }
773
774 struct less_SVertex2D : public binary_function<SVertex*, SVertex*, bool> 
775 {
776   real epsilon;
777   less_SVertex2D(real eps)
778     : binary_function<SVertex*,SVertex*,bool>()
779   {
780     epsilon = eps;
781   }
782         bool operator()(SVertex* x, SVertex* y) 
783   {
784     Vec3r A = x->point2D();
785     Vec3r B = y->point2D();
786     for(unsigned int i=0; i<3; i++)
787     {
788       if((fabs(A[i] - B[i])) < epsilon)
789         continue;
790       if(A[i] < B[i])
791         return true;
792       if(A[i] > B[i])
793         return false;
794     }
795     
796     return false;
797   }
798 };
799
800 typedef Segment<FEdge*,Vec3r > segment;
801 typedef Intersection<segment> intersection;
802
803 struct less_Intersection : public binary_function<intersection*, intersection*, bool> 
804 {
805   segment *edge;
806   less_Intersection(segment *iEdge)
807     : binary_function<intersection*,intersection*,bool>()
808   {
809     edge = iEdge;
810   }
811         bool operator()(intersection* x, intersection* y) 
812   {
813     real tx = x->getParameter(edge);
814     real ty = y->getParameter(edge);
815     if(tx > ty)
816       return true;
817     return false;
818   }
819 };
820
821 struct silhouette_binary_rule : public binary_rule<segment,segment>
822 {
823   silhouette_binary_rule() : binary_rule<segment,segment>() {}
824   virtual bool operator() (segment& s1, segment& s2)
825   {
826     FEdge * f1 = s1.edge();
827     FEdge * f2 = s2.edge();
828
829     if((!(((f1)->getNature() & Nature::SILHOUETTE) || ((f1)->getNature() & Nature::BORDER))) && (!(((f2)->getNature() & Nature::SILHOUETTE) || ((f2)->getNature() & Nature::BORDER))))
830       return false;
831       
832     return true;
833   }
834 };
835
836 void ViewMapBuilder::ComputeSweepLineIntersections(ViewMap *ioViewMap, real epsilon)
837 {
838   vector<SVertex*>& svertices = ioViewMap->SVertices();
839   bool progressBarDisplay = false;
840   unsigned sVerticesSize = svertices.size();
841   unsigned fEdgesSize = ioViewMap->FEdges().size();
842   //  ViewMap::fedges_container& fedges = ioViewMap->FEdges();
843   //  for(ViewMap::fedges_container::const_iterator f=fedges.begin(), end=fedges.end();
844   //  f!=end;
845   //  ++f){
846   //    cout << (*f)->aMaterialIndex() << "-" << (*f)->bMaterialIndex() << endl;
847   //  }
848     
849   unsigned progressBarStep = 0;
850   
851   if(_pProgressBar != NULL && fEdgesSize > gProgressBarMinSize) {
852     unsigned progressBarSteps = min(gProgressBarMaxSteps, sVerticesSize);
853     progressBarStep = sVerticesSize / progressBarSteps;
854     _pProgressBar->reset();
855     _pProgressBar->setLabelText("Computing Sweep Line Intersections");
856     _pProgressBar->setTotalSteps(progressBarSteps);
857     _pProgressBar->setProgress(0);
858     progressBarDisplay = true;
859   }
860   
861   unsigned counter = progressBarStep;
862
863   sort(svertices.begin(), svertices.end(), less_SVertex2D(epsilon));
864
865   SweepLine<FEdge*,Vec3r> SL;
866
867   vector<FEdge*>& ioEdges = ioViewMap->FEdges();
868
869   vector<segment* > segments;
870
871   vector<FEdge*>::iterator fe,fend;
872
873   for(fe=ioEdges.begin(), fend=ioEdges.end();
874   fe!=fend;
875   fe++)
876   {
877     segment * s = new segment((*fe), (*fe)->vertexA()->point2D(), (*fe)->vertexB()->point2D());
878     (*fe)->userdata = s;    
879     segments.push_back(s);
880   }
881
882   vector<segment*> vsegments;
883   for(vector<SVertex*>::iterator sv=svertices.begin(),svend=svertices.end();
884   sv!=svend;
885   sv++)
886   {
887     const vector<FEdge*>& vedges = (*sv)->fedges();
888     
889     for(vector<FEdge*>::const_iterator sve=vedges.begin(), sveend=vedges.end();
890     sve!=sveend;
891     sve++)
892     {
893       vsegments.push_back((segment*)((*sve)->userdata));
894     }
895
896     Vec3r evt((*sv)->point2D());
897     silhouette_binary_rule sbr;
898     SL.process(evt, vsegments, sbr);
899
900     if(progressBarDisplay) {  
901       counter--;
902       if (counter <= 0) {
903         counter = progressBarStep;
904         _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
905       }
906     }
907     vsegments.clear();
908   }
909
910   // reset userdata:
911   for(fe=ioEdges.begin(), fend=ioEdges.end();
912   fe!=fend;
913   fe++)
914     (*fe)->userdata = NULL;
915  
916   // list containing the new edges resulting from splitting operations. 
917   vector<FEdge*> newEdges;
918   
919   // retrieve the intersected edges:
920   vector<segment* >& iedges = SL.intersectedEdges();
921   // retrieve the intersections:
922   vector<intersection*>& intersections = SL.intersections();
923   
924   int id=0;
925   // create a view vertex for each intersection and linked this one 
926   // with the intersection object
927   vector<intersection*>::iterator i, iend;
928   for(i=intersections.begin(),iend=intersections.end();
929   i!=iend;
930   i++)
931   {
932     FEdge *fA = (*i)->EdgeA->edge();
933     FEdge *fB = (*i)->EdgeB->edge();
934     
935     Vec3r A1 = fA->vertexA()->point3D();
936     Vec3r A2 = fA->vertexB()->point3D();
937     Vec3r B1 = fB->vertexA()->point3D();
938     Vec3r B2 = fB->vertexB()->point3D();
939
940     Vec3r a1 = fA->vertexA()->point2D();
941     Vec3r a2 = fA->vertexB()->point2D();
942     Vec3r b1 = fB->vertexA()->point2D();
943     Vec3r b2 = fB->vertexB()->point2D();
944
945     real ta = (*i)->tA;
946     real tb = (*i)->tB;
947
948     if((ta < -epsilon) || (ta > 1+epsilon))
949         cerr << "Warning: 2D intersection out of range for edge " << fA->vertexA()->getId() << " - " << fA->vertexB()->getId() << endl;
950     
951     if((tb < -epsilon) || (tb > 1+epsilon))
952         cerr << "Warning: 2D intersection out of range for edge " << fB->vertexA()->getId() << " - " << fB->vertexB()->getId() << endl;
953     
954     real Ta = SilhouetteGeomEngine::ImageToWorldParameter(fA, ta);
955     real Tb = SilhouetteGeomEngine::ImageToWorldParameter(fB, tb);
956
957     if((Ta < -epsilon) || (Ta > 1+epsilon))
958         cerr << "Warning: 3D intersection out of range for edge " << fA->vertexA()->getId() << " - " << fA->vertexB()->getId() << endl;
959     
960     if((Tb < -epsilon) || (Tb > 1+epsilon))
961         cerr << "Warning: 3D intersection out of range for edge " << fB->vertexA()->getId() << " - " << fB->vertexB()->getId() << endl;
962
963     TVertex * tvertex = ioViewMap->CreateTVertex(Vec3r(A1 + Ta*(A2-A1)), Vec3r(a1 + ta*(a2-a1)), fA, 
964                                                  Vec3r(B1 + Tb*(B2-B1)), Vec3r(b1 + tb*(b2-b1)), fB, id);
965      
966     (*i)->userdata = tvertex;
967     ++id;
968   }
969
970   progressBarStep = 0;
971
972   if(progressBarDisplay) {
973     unsigned iEdgesSize = iedges.size();
974     unsigned progressBarSteps = min(gProgressBarMaxSteps, iEdgesSize);
975     progressBarStep = iEdgesSize / progressBarSteps;
976     _pProgressBar->reset();
977     _pProgressBar->setLabelText("Splitting intersected edges");
978     _pProgressBar->setTotalSteps(progressBarSteps);
979     _pProgressBar->setProgress(0);
980   }
981   
982   counter = progressBarStep;
983
984   vector<TVertex*> edgeVVertices;
985   vector<ViewEdge*> newVEdges;
986   vector<segment* >::iterator s, send;
987   for(s=iedges.begin(),send=iedges.end();
988   s!=send;
989   s++)
990   {
991     edgeVVertices.clear();
992     newEdges.clear();
993     newVEdges.clear();
994     
995     FEdge* fedge = (*s)->edge();
996     ViewEdge *vEdge = fedge->viewedge();
997     ViewShape *shape = vEdge->viewShape();
998     
999     vector<intersection*>& eIntersections = (*s)->intersections();
1000     // we first need to sort these intersections from farther to closer to A
1001     sort(eIntersections.begin(), eIntersections.end(), less_Intersection(*s));
1002     for(i=eIntersections.begin(),iend=eIntersections.end();
1003     i!=iend;
1004     i++)
1005       edgeVVertices.push_back((TVertex*)(*i)->userdata);
1006
1007     shape->SplitEdge(fedge, edgeVVertices, ioViewMap->FEdges(), ioViewMap->ViewEdges()); 
1008
1009     if(progressBarDisplay) {  
1010       counter--;
1011       if (counter <= 0) {
1012         counter = progressBarStep;
1013         _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
1014       }
1015     }
1016   }
1017
1018   // reset userdata:
1019   for(fe=ioEdges.begin(), fend=ioEdges.end();
1020   fe!=fend;
1021   fe++)
1022     (*fe)->userdata = NULL;
1023
1024   // delete segments
1025   //  if(!segments.empty()){
1026   //    for(s=segments.begin(),send=segments.end();
1027   //    s!=send;
1028   //    s++){
1029   //      delete *s;
1030   //    }
1031   //    segments.clear();
1032   //  }
1033 }