Merged changes in the trunk up to revision 44266 (including BMesh).
[blender.git] / source / blender / freestyle / intern / view_map / FEdgeXDetector.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 "FEdgeXDetector.h"
23 #include "float.h"
24 #include "../geometry/GeomUtils.h"
25 #include <math.h>
26 #include "../geometry/normal_cycle.h"
27
28 void FEdgeXDetector::processShapes(WingedEdge& we) {
29   bool progressBarDisplay = false;
30   Vec3r Min, Max;
31   vector<WShape*> wshapes = we.getWShapes();
32   WXShape * wxs;
33
34   if(_pProgressBar != NULL) {
35     _pProgressBar->reset();
36     _pProgressBar->setLabelText("Detecting feature lines");
37     _pProgressBar->setTotalSteps(wshapes.size() * 3);
38     _pProgressBar->setProgress(0);
39     progressBarDisplay = true;
40   }
41
42   for(vector<WShape*>::const_iterator it = wshapes.begin();
43   it != wshapes.end();
44   it++){
45     wxs = dynamic_cast<WXShape*>(*it);
46     wxs->bbox(Min, Max);
47     _bbox_diagonal = (Max-Min).norm();
48     if(_changes){
49       vector<WFace*>& wfaces = wxs->GetFaceList();
50       for(vector<WFace*>::iterator wf=wfaces.begin(), wfend=wfaces.end();
51       wf!=wfend;
52       ++wf){
53         WXFace* wxf = dynamic_cast<WXFace*>(*wf);
54         wxf->Clear();
55       }
56       _computeViewIndependant = true;
57     } else if (!(wxs)->getComputeViewIndependantFlag()) {
58       wxs->Reset();
59       _computeViewIndependant = false;
60     } else {
61       _computeViewIndependant = true;
62     }
63     preProcessShape(wxs);
64     if (progressBarDisplay)
65       _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
66     processBorderShape(wxs);
67     if(_computeMaterialBoundaries)
68       processMaterialBoundaryShape(wxs);
69     processCreaseShape(wxs);
70     if(_computeRidgesAndValleys)
71       processRidgesAndValleysShape(wxs);
72     if(_computeSuggestiveContours)
73       processSuggestiveContourShape(wxs);
74     processSilhouetteShape(wxs);
75     processEdgeMarksShape(wxs);
76     if (progressBarDisplay)
77       _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
78
79     // build smooth edges:
80     buildSmoothEdges(wxs);
81
82     // Post processing for suggestive contours
83     if(_computeSuggestiveContours)
84       postProcessSuggestiveContourShape(wxs);
85     if (progressBarDisplay)
86       _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
87
88     wxs->setComputeViewIndependantFlag(false);
89     _computeViewIndependant = false;
90     _changes = false;
91
92     // reset user data
93     (*it)->ResetUserData();
94   }
95 }
96
97 // GENERAL STUFF
98 ////////////////
99 void FEdgeXDetector::preProcessShape(WXShape* iWShape) {
100   _meanK1 = 0;
101   _meanKr = 0;
102   _minK1 = FLT_MAX;
103   _maxK1 = -FLT_MAX;
104   _minKr = FLT_MAX;
105   _maxKr = -FLT_MAX;
106   _nPoints = 0;
107   _meanEdgeSize = iWShape->getMeanEdgeSize();
108   
109   vector<WFace*>& wfaces = iWShape->GetFaceList();
110   vector<WFace*>::iterator f,fend;
111   // view dependant stuff
112   for(f=wfaces.begin(), fend=wfaces.end();
113   f!=fend;
114   ++f){
115     preProcessFace((WXFace*)(*f));
116   }
117
118   if(_computeRidgesAndValleys || _computeSuggestiveContours ) {
119     vector<WVertex*>& wvertices = iWShape->getVertexList();
120     for(vector<WVertex*>::iterator wv=wvertices.begin(), wvend=wvertices.end(); 
121     wv!=wvend;
122     ++wv){
123       // Compute curvatures
124       WXVertex * wxv = dynamic_cast<WXVertex*>(*wv);
125       computeCurvatures(wxv);
126     }
127     _meanK1 /= (real)(_nPoints);
128     _meanKr /= (real)(_nPoints);
129   }
130 }
131
132 void FEdgeXDetector::preProcessFace(WXFace *iFace){
133   Vec3r firstPoint = iFace->GetVertex(0)->GetVertex();
134   Vec3r N = iFace->GetNormal();
135
136   // Compute the dot product between V (=_Viewpoint - firstPoint) and N:
137   Vec3r V;
138   if (_orthographicProjection) {
139     V = Vec3r(0.0, 0.0, _Viewpoint.z() - firstPoint.z());
140   } else {
141     V = Vec3r(_Viewpoint - firstPoint);
142   }
143   N.normalize();
144   V.normalize();
145   iFace->setDotP(N * V);
146   
147   // compute the distance between the face center and the viewpoint:
148   if (_orthographicProjection) {
149     iFace->setZ(iFace->center().z() - _Viewpoint.z());
150   } else {
151     Vec3r dist_vec(iFace->center() - _Viewpoint);
152     iFace->setZ(dist_vec.norm());
153   }
154 }
155
156 void FEdgeXDetector::computeCurvatures(WXVertex *vertex){
157    // CURVATURE LAYER
158   // store all the curvature datas for each vertex
159
160   //soc unused - real K1, K2
161   real cos2theta, sin2theta;
162   Vec3r e1, n, v;
163   // one vertex curvature info :
164   CurvatureInfo *C;
165   float radius = _sphereRadius*_meanEdgeSize; 
166   
167   // view independant stuff
168   if(_computeViewIndependant){
169     C = new CurvatureInfo();   
170     vertex->setCurvatures(C);
171     OGF::NormalCycle ncycle ;
172     ncycle.begin() ;
173     if(radius > 0) {
174       OGF::compute_curvature_tensor(vertex, radius, ncycle) ;
175     } else {
176       OGF::compute_curvature_tensor_one_ring(vertex, ncycle) ;
177     }
178     ncycle.end() ;
179     C->K1 = ncycle.kmin();
180     C->K2 = ncycle.kmax();
181     C->e1 = ncycle.Kmax(); //ncycle.kmin() * ncycle.Kmax();
182     C->e2 = ncycle.Kmin(); //ncycle.kmax() * ncycle.Kmin() ;
183
184     real absK1 = fabs(C->K1);
185     _meanK1 += absK1;
186     if(absK1 > _maxK1)
187       _maxK1 = absK1;
188     if(absK1 < _minK1)
189       _minK1 = absK1;
190   }
191   // view dependant
192   C = vertex->curvatures();
193   if(C == 0)
194     return;
195     
196   // compute radial curvature :
197   n = C->e1 ^ C->e2;
198   if (_orthographicProjection) {
199     v = Vec3r(0.0, 0.0, _Viewpoint.z() - vertex->GetVertex().z());
200   } else {
201     v = Vec3r(_Viewpoint - vertex->GetVertex());
202   }
203   C->er = v - (v * n) * n;
204   C->er.normalize();
205   e1 = C->e1;
206   e1.normalize();
207   cos2theta = C->er * e1;
208   cos2theta *= cos2theta;
209   sin2theta = 1 - cos2theta;
210   C->Kr = C->K1 * cos2theta + C->K2 * sin2theta;
211   real absKr = fabs(C->Kr);
212   _meanKr += absKr;
213   if(absKr > _maxKr)
214     _maxKr = absKr;
215   if(absKr < _minKr)
216     _minKr = absKr;
217
218   ++_nPoints;
219 }
220
221 // SILHOUETTE
222 /////////////
223 void FEdgeXDetector::processSilhouetteShape(WXShape* iWShape) {
224   // Make a first pass on every polygons in order 
225   // to compute all their silhouette relative values:
226   //------------------------------------------------
227   vector<WFace*>& wfaces = iWShape->GetFaceList();
228   vector<WFace*>::iterator f,fend;
229   for(f=wfaces.begin(), fend=wfaces.end();
230   f!=fend;
231   ++f)
232   {
233     ProcessSilhouetteFace((WXFace*)(*f));
234   }
235
236   // Make a pass on the edges to detect 
237   // the silhouette edges that are not smooth
238   // --------------------
239   vector<WEdge*>::iterator we, weend;
240   vector<WEdge*> &wedges = iWShape->getEdgeList();
241   for(we=wedges.begin(), weend=wedges.end();
242   we!=weend;
243   ++we)
244   {
245     ProcessSilhouetteEdge((WXEdge*)(*we));
246   }
247 }
248
249 void FEdgeXDetector::ProcessSilhouetteFace(WXFace *iFace)
250 {
251   
252   // SILHOUETTE LAYER
253   Vec3r normal;
254   // Compute the dot products between View direction and N at each vertex
255   // of the face:
256   Vec3r point;
257   int closestPointId = 0;
258   real dist, minDist = FLT_MAX;
259   int numVertices = iFace->numberOfVertices();
260   WXFaceLayer * faceLayer = new WXFaceLayer(iFace, Nature::SILHOUETTE, true);
261   for(int i=0; i<numVertices; i++){
262     point = iFace->GetVertex(i)->GetVertex();
263     normal = iFace->GetVertexNormal(i);
264     normal.normalize();
265     Vec3r V;
266     if (_orthographicProjection) {
267       V = Vec3r(0.0, 0.0, _Viewpoint.z() - point.z());
268     } else {
269       V = Vec3r(_Viewpoint - point);
270     }
271     V.normalize();
272     real d = normal * V;
273     faceLayer->PushDotP(d);
274     // Find the point the closest to the viewpoint
275     if (_orthographicProjection) {
276       dist = point.z() - _Viewpoint.z();
277     } else {
278       Vec3r dist_vec(point - _Viewpoint);
279       dist = dist_vec.norm();
280     }
281     if(dist < minDist) {
282       minDist = dist;
283       closestPointId = i;
284     }
285   }
286   // Set the closest point id:
287   faceLayer->setClosestPointIndex(closestPointId);
288   // Add this layer to the face:
289   iFace->AddSmoothLayer(faceLayer);
290 }
291
292 void FEdgeXDetector::ProcessSilhouetteEdge(WXEdge *iEdge)
293 {
294   if(iEdge->nature() & Nature::BORDER)
295     return;
296   // SILHOUETTE ?
297   //-------------
298   WXFace * fA = (WXFace *)iEdge->GetaOEdge()->GetaFace();
299   WXFace * fB = (WXFace *)iEdge->GetaOEdge()->GetbFace();
300
301   if((fA->front())^(fB->front())){ // fA->visible XOR fB->visible (true if one is 0 and the other is 1)
302     // The only edges we want to set as silhouette edges in this 
303     // way are the ones with 2 different normals for 1 vertex 
304     // for these two faces
305     //--------------------
306     // In reality we only test the normals for 1 of the 2 vertices.
307     if(fA->GetVertexNormal(iEdge->GetaVertex()) == fB->GetVertexNormal(iEdge->GetaVertex()))
308       return;
309     iEdge->AddNature(Nature::SILHOUETTE);
310     if(fB->front())
311       iEdge->setOrder(1);
312     else
313       iEdge->setOrder(-1);
314   }
315 }
316
317
318 // BORDER
319 /////////
320 void FEdgeXDetector::processBorderShape(WXShape* iWShape) {
321
322   if(!_computeViewIndependant)
323     return;
324   // Make a pass on the edges to detect 
325   // the BORDER
326   // --------------------
327   vector<WEdge*>::iterator we, weend;
328   vector<WEdge*> &wedges = iWShape->getEdgeList();
329   for(we=wedges.begin(), weend=wedges.end();
330   we!=weend;
331   ++we){
332     ProcessBorderEdge((WXEdge*)(*we));
333   }
334 }
335
336 void FEdgeXDetector::ProcessBorderEdge(WXEdge *iEdge)
337 {
338   // first check whether it is a border edge:
339   // BORDER ?
340   //---------
341   if(iEdge->GetaFace() == 0){
342     // it is a border edge
343     iEdge->AddNature(Nature::BORDER);
344   }
345 }
346
347
348 // CREASE
349 /////////
350 void FEdgeXDetector::processCreaseShape(WXShape* iWShape) {
351   if(!_computeViewIndependant)
352     return;
353
354   // Make a pass on the edges to detect 
355   // the CREASE 
356   // --------------------
357   vector<WEdge*>::iterator we, weend;
358   vector<WEdge*> &wedges = iWShape->getEdgeList();
359   for(we=wedges.begin(), weend=wedges.end();
360   we!=weend;
361   ++we){
362     ProcessCreaseEdge((WXEdge*)(*we));
363   }
364 }
365
366 void FEdgeXDetector::ProcessCreaseEdge(WXEdge *iEdge)
367 {
368   // CREASE ?
369   //---------
370   if(iEdge->nature() & Nature::BORDER)
371     return;
372   WXFace * fA = (WXFace *)iEdge->GetaOEdge()->GetaFace();
373   WXFace * fB = (WXFace *)iEdge->GetaOEdge()->GetbFace();
374
375   WVertex * aVertex = iEdge->GetaVertex();
376   if((fA->GetVertexNormal(aVertex) * fB->GetVertexNormal(aVertex)) <= _creaseAngle)
377     iEdge->AddNature(Nature::CREASE);
378 }
379
380 // RIDGES AND VALLEYS
381 /////////////////////
382
383 void FEdgeXDetector::processRidgesAndValleysShape(WXShape* iWShape) {
384   // Don't forget to add the built layer to the face at the end
385   // of the ProcessFace:
386   //iFace->AddSmoothLayer(faceLayer);
387
388   if((!_computeViewIndependant))
389     return;
390   
391   // Here the curvatures must already have been computed
392   vector<WFace*>& wfaces = iWShape->GetFaceList();
393   vector<WFace*>::iterator f, fend;
394   for(f=wfaces.begin(), fend=wfaces.end();
395         f!=fend;
396         ++f)
397   { 
398         ProcessRidgeFace((WXFace*)(*f));
399   }
400 }
401
402
403 // RIDGES
404 /////////
405
406 void FEdgeXDetector::ProcessRidgeFace(WXFace *iFace)
407 {
408   WXFaceLayer * flayer = new WXFaceLayer(iFace, Nature::RIDGE|Nature::VALLEY, false);
409   iFace->AddSmoothLayer(flayer);
410   
411   unsigned int numVertices = iFace->numberOfVertices();
412   for(unsigned int i=0; i<numVertices; ++i){
413     WVertex *wv = iFace->GetVertex(i);
414     WXVertex * wxv = dynamic_cast<WXVertex*>(wv);
415     flayer->PushDotP(wxv->curvatures()->K1);
416   }
417
418   real threshold = 0;
419   //real threshold = _maxK1 - (_maxK1-_meanK1)/20.0;
420   
421   if(flayer->nPosDotP()!=numVertices){
422     if((fabs(flayer->dotP(0)) < threshold) && (fabs(flayer->dotP(1)) < threshold) && (fabs(flayer->dotP(2)) < threshold)){
423       flayer->ReplaceDotP(0, 0);
424       flayer->ReplaceDotP(1, 0);
425       flayer->ReplaceDotP(2, 0);
426     }
427   }
428 }
429
430 // void FEdgeXDetector::ProcessRidgeFace(WXFace *iFace)
431 // {
432  
433 //  //  RIDGE LAYER
434 //  //    compute the RidgeFunction, that is the derivative of the ppal curvature
435 //  //    along e1 at each vertex of the face
436  
437 //  WVertex *v;
438 //  Vec3r v1v2;
439 //  real t;
440 //  vector<WXFaceLayer*> SmoothLayers;
441 //  WXFaceLayer *faceLayer;
442 //  Face_Curvature_Info *layer_info;
443 //  real K1_a(0), K1_b(0);
444 //  Vec3r Inter_a, Inter_b;
445
446 //  //  find the ridge layer of the face
447 //  iFace->retrieveSmoothLayers(Nature::RIDGE, SmoothLayers);
448 //  if(SmoothLayers.size()!=1)
449 //    return;
450 //  faceLayer = SmoothLayers[0];
451 //  //  retrieve the curvature info of this layer
452 //  layer_info = (Face_Curvature_Info *)faceLayer->userdata;
453  
454 //  int numVertices = iFace->numberOfVertices();
455 //  for(int i=0; i<numVertices; i++){
456 //    v = iFace->GetVertex(i);
457 //    //    vec_curvature_info[i] contains the curvature info of this vertex
458 //    Vec3r e2 = layer_info->vec_curvature_info[i]->K2*layer_info->vec_curvature_info[i]->e2;
459 //    Vec3r e1 = layer_info->vec_curvature_info[i]->K1*layer_info->vec_curvature_info[i]->e1;
460 //    e2.normalize();
461    
462 //    WVertex::face_iterator fit = v->faces_begin();
463 //    WVertex::face_iterator fitend = v->faces_end();
464 //    for(; fit!=fitend; ++fit){  
465 //      WXFace * wxf = dynamic_cast<WXFace*>(*fit);
466 //      WOEdge * oppositeEdge;
467 //      if(!(wxf->getOppositeEdge(v, oppositeEdge)))
468 //        continue;
469 //      v1v2 = oppositeEdge->GetbVertex()->GetVertex() - oppositeEdge->GetaVertex()->GetVertex();
470 //      GeomUtils::intersection_test res;
471 //      res = GeomUtils::intersectRayPlane(oppositeEdge->GetaVertex()->GetVertex(), v1v2,
472 //        e2, -(v->GetVertex()*e2),
473 //        t,1.e-06);
474 //      if((res == GeomUtils::DO_INTERSECT) && (t>=0.0) && (t<=1.0)){
475 //        vector<WXFaceLayer*> second_ridge_layer;
476 //        wxf->retrieveSmoothLayers(Nature::RIDGE, second_ridge_layer);
477 //        if(second_ridge_layer.size()!=1)
478 //          continue;
479 //        Face_Curvature_Info *second_layer_info = (Face_Curvature_Info*)second_ridge_layer[0]->userdata;
480        
481 //        unsigned index1 = wxf->GetIndex(oppositeEdge->GetaVertex());
482 //        unsigned index2 = wxf->GetIndex(oppositeEdge->GetbVertex());
483 //        real K1_1 = second_layer_info->vec_curvature_info[index1]->K1;
484 //        real K1_2 = second_layer_info->vec_curvature_info[index2]->K1;
485 //        real K1 = (1.0-t)*K1_1 + t*K1_2;
486 //        Vec3r inter((1.0-t)*oppositeEdge->GetaVertex()->GetVertex() + t*oppositeEdge->GetbVertex()->GetVertex());
487 //        Vec3r vtmp(inter - v->GetVertex());
488 //        // is it K1_a or K1_b ?
489 //        if(vtmp*e1 > 0){
490 //          K1_b = K1;
491 //          Inter_b = inter;
492 //        }else{
493 //          K1_a = K1;
494 //          Inter_a = inter;
495 //        }
496 //      }  
497 //    }  
498 //    //    Once we have K1 along the the ppal direction
499 //    //      compute the derivative : K1b - K1a
500 //    //      put it in DotP
501 //    //real d = fabs(K1_b)-fabs(K1_a);
502 //    real d = 0;
503 //    real threshold = _meanK1 + (_maxK1-_meanK1)/7.0;
504 //    //real threshold = _meanK1;
505 //    //if((fabs(K1_b) > threshold) || ((fabs(K1_a) > threshold)))
506 //      d = (K1_b)-(K1_a)/(Inter_b-Inter_a).norm();
507 //    faceLayer->PushDotP(d);
508 //    //faceLayer->PushDotP(layer_info->vec_curvature_info[i]->K1);
509 //  } 
510
511 //  // Make the values relevant by checking whether all principal 
512 //  // directions have the "same" direction:
513 //  Vec3r e0((layer_info->vec_curvature_info[0]->K1*layer_info->vec_curvature_info[0]->e1));
514 //  e0.normalize();
515 //  Vec3r e1((layer_info->vec_curvature_info[1]->K1*layer_info->vec_curvature_info[1]->e1));
516 //  e1.normalize();
517 //  Vec3r e2((layer_info->vec_curvature_info[2]->K1*layer_info->vec_curvature_info[2]->e1));
518 //  e2.normalize();
519 //  if (e0 * e1 < 0)
520 //    // invert dotP[1]
521 //    faceLayer->ReplaceDotP(1, -faceLayer->dotP(1));
522 //  if (e0 * e2 < 0)
523 //    // invert dotP[2]
524 //    faceLayer->ReplaceDotP(2, -faceLayer->dotP(2));
525
526 //  // remove the weakest values;
527 //  //real minDiff = (_maxK1 - _minK1)/10.0;
528 //  //  real minDiff = _meanK1;
529 //  //  if((faceLayer->dotP(0) < minDiff) && (faceLayer->dotP(1) < minDiff) && (faceLayer->dotP(2) < minDiff)){
530 //  //    faceLayer->ReplaceDotP(0, 0);
531 //  //    faceLayer->ReplaceDotP(1, 0);
532 //  //    faceLayer->ReplaceDotP(2, 0);
533 //  //  } 
534 // }
535
536 // SUGGESTIVE CONTOURS
537 //////////////////////
538
539 void FEdgeXDetector::processSuggestiveContourShape(WXShape* iWShape) {
540   
541   // Here the curvatures must already have been computed
542   vector<WFace*>& wfaces = iWShape->GetFaceList();
543   vector<WFace*>::iterator f, fend;
544   for(f=wfaces.begin(), fend=wfaces.end();
545         f!=fend;
546         ++f)
547   {
548           ProcessSuggestiveContourFace((WXFace*)(*f));
549   }
550 }
551
552 void FEdgeXDetector::ProcessSuggestiveContourFace(WXFace *iFace)
553 {
554   WXFaceLayer *faceLayer = new WXFaceLayer(iFace, Nature::SUGGESTIVE_CONTOUR, true);
555   iFace->AddSmoothLayer(faceLayer);
556
557   unsigned int numVertices = iFace->numberOfVertices();
558   for(unsigned int i=0; i<numVertices; ++i){
559     WVertex *wv = iFace->GetVertex(i);
560     WXVertex * wxv = dynamic_cast<WXVertex*>(wv);
561     faceLayer->PushDotP(wxv->curvatures()->Kr);
562   }
563   
564   // FIXME: find a more clever way to compute the threshold
565 //   real threshold = _meanKr;
566 //   if(faceLayer->nPosDotP()!=numVertices){
567 //     if((fabs(faceLayer->dotP(0)) < threshold) && (fabs(faceLayer->dotP(1)) < threshold) && (fabs(faceLayer->dotP(2)) < threshold)){
568 //       faceLayer->ReplaceDotP(0, 0);
569 //       faceLayer->ReplaceDotP(1, 0);
570 //       faceLayer->ReplaceDotP(2, 0);
571 //     }
572 //   }
573 }
574
575 void FEdgeXDetector::postProcessSuggestiveContourShape(WXShape* iShape) {
576   vector<WFace*>& wfaces = iShape->GetFaceList();
577   vector<WFace*>::iterator f, fend;
578   for(f=wfaces.begin(), fend=wfaces.end();
579       f!=fend;
580       ++f)
581     {
582       postProcessSuggestiveContourFace((WXFace*)(*f));
583     }
584 }
585
586 void FEdgeXDetector::postProcessSuggestiveContourFace(WXFace *iFace) {
587
588   // Compute the derivative of the radial curvature in the radial direction,
589   // at the two extremities of the smooth edge.
590   // If the derivative is smaller than a given threshold _kr_derivative_epsilon,
591   // discard the edge.
592  
593   // Find the suggestive contour layer of the face (zero or one edge).
594   vector<WXFaceLayer*> sc_layers;
595   iFace->retrieveSmoothEdgesLayers(Nature::SUGGESTIVE_CONTOUR, sc_layers);
596   if(sc_layers.empty())
597     return;
598
599   WXFaceLayer *sc_layer;
600   sc_layer = sc_layers[0];
601
602   // Compute the derivative value at each vertex of the face, and add it in a vector.
603   vector<real> kr_derivatives;
604
605   unsigned vertices_nb = iFace->numberOfVertices();
606   WXVertex *v, *opposite_vertex_a, *opposite_vertex_b;
607   WXFace *wxf;
608   WOEdge *opposite_edge;
609   Vec3r normal_vec, radial_normal_vec, er_vec, v_vec, inter, inter1, inter2, tmp_vec;
610   GeomUtils::intersection_test res;
611   real kr(0), kr1(0), kr2(0), t;
612
613   for (unsigned i = 0; i < vertices_nb; ++i) {
614     v = (WXVertex*)(iFace->GetVertex(i));
615
616     // v is a singular vertex, skip it.
617     if (v->isBoundary()) {
618       kr_derivatives.push_back(0);
619       continue;
620     }
621
622     v_vec = v->GetVertex();
623     er_vec = v->curvatures()->er;
624     
625     // For each vertex, iterate on its adjacent faces.
626     for (WVertex::face_iterator fit = v->faces_begin(), fitend = v->faces_end();
627          fit != fitend;
628          ++fit) {  
629       wxf = dynamic_cast<WXFace*>(*fit);
630       if(!(wxf->getOppositeEdge(v, opposite_edge)))
631         continue;
632       
633       opposite_vertex_a = (WXVertex*)opposite_edge->GetaVertex();
634       opposite_vertex_b = (WXVertex*)opposite_edge->GetbVertex();
635       normal_vec = wxf->GetVertexNormal(v); // FIXME: what about e1 ^ e2 ?
636       radial_normal_vec = er_vec ^ normal_vec;
637
638       // Test wether the radial plan intersects with the edge at the opposite of v.
639       res = GeomUtils::intersectRayPlane(opposite_vertex_a->GetVertex(), opposite_edge->GetVec(),
640                                          radial_normal_vec, -(v_vec * radial_normal_vec),
641                                          t,
642                                          1.e-06);
643       
644       // If there is an intersection, compute the value of the derivative ath that point.
645       if ((res == GeomUtils::DO_INTERSECT) && (t >= 0) && (t <= 1)) {
646         kr = t * opposite_vertex_a->curvatures()->Kr + (1 - t) * opposite_vertex_b->curvatures()->Kr;
647         inter = opposite_vertex_a->GetVertex() + t * opposite_edge->GetVec();
648         tmp_vec = inter - v->GetVertex();
649         // Is it kr1 or kr2?
650         if (tmp_vec * er_vec > 0) {
651           kr2 = kr;
652           inter2 = inter;
653         } else {
654           kr1 = kr;
655           inter1 = inter;
656         }
657       }
658     }
659
660     // Now we have kr1 and kr2 along the radial direction, for one vertex of iFace.
661     // We have to compute the derivative of kr for that vertex, equal to:
662     // (kr2 - kr1) / dist(inter1, inter2).
663     // Then we add it to the vector of derivatives.
664     v->curvatures()->dKr = (kr2 - kr1) / (inter2 - inter1).norm();
665     kr_derivatives.push_back(v->curvatures()->dKr);
666   }
667
668   // At that point, we have the derivatives for each vertex of iFace.
669   // All we have to do now is to use linear interpolation to compute the values at
670   // the extremities of the smooth edge.
671   WXSmoothEdge *sc_edge = sc_layer->getSmoothEdge();
672   WOEdge *sc_oedge = sc_edge->woea();
673   t = sc_edge->ta();
674   if (t * kr_derivatives[iFace->GetIndex(sc_oedge->GetaVertex())] +
675       (1 - t) * kr_derivatives[iFace->GetIndex(sc_oedge->GetbVertex())] < _kr_derivative_epsilon) {
676     sc_layer->removeSmoothEdge();
677     return;
678   }
679   sc_oedge = sc_edge->woeb();
680   t = sc_edge->tb();
681   if (t * kr_derivatives[iFace->GetIndex(sc_oedge->GetaVertex())] +
682       (1 - t) * kr_derivatives[iFace->GetIndex(sc_oedge->GetbVertex())] < _kr_derivative_epsilon)
683     sc_layer->removeSmoothEdge();
684 }
685
686 // MATERIAL_BOUNDARY
687 ////////////////////
688 void FEdgeXDetector::processMaterialBoundaryShape(WXShape* iWShape) {
689
690   if(!_computeViewIndependant)
691     return;
692   // Make a pass on the edges to detect material boundaries
693   vector<WEdge*>::iterator we, weend;
694   vector<WEdge*> &wedges = iWShape->getEdgeList();
695   for(we=wedges.begin(), weend=wedges.end();
696   we!=weend;
697   ++we){
698     ProcessMaterialBoundaryEdge((WXEdge*)(*we));
699   }
700 }
701
702 void FEdgeXDetector::ProcessMaterialBoundaryEdge(WXEdge *iEdge)
703 {
704   // check whether the edge is a material boundary?
705   WFace *aFace = iEdge->GetaFace();
706   WFace *bFace = iEdge->GetbFace();
707   if(aFace && bFace && aFace->frs_materialIndex() != bFace->frs_materialIndex()){
708     iEdge->AddNature(Nature::MATERIAL_BOUNDARY);
709   }
710 }
711
712 // EDGE MARKS
713 /////////////
714   
715 void FEdgeXDetector::processEdgeMarksShape(WXShape* iShape) {
716   // Make a pass on the edges to detect material boundaries
717   vector<WEdge*>::iterator we, weend;
718   vector<WEdge*> &wedges = iShape->getEdgeList();
719   for(we=wedges.begin(), weend=wedges.end();
720   we!=weend;
721   ++we){
722     ProcessEdgeMarks((WXEdge*)(*we));
723   }
724 }
725
726 void FEdgeXDetector::ProcessEdgeMarks(WXEdge *iEdge) {
727   if (iEdge->GetMark()) {
728     iEdge->AddNature(Nature::EDGE_MARK);
729   }
730 }
731
732 // Build Smooth edges
733 /////////////////////
734 void FEdgeXDetector::buildSmoothEdges(WXShape* iShape){
735   bool hasSmoothEdges = false;
736
737   // Make a last pass to build smooth edges from the previous stored values:
738   //--------------------------------------------------------------------------
739   vector<WFace*>& wfaces = iShape->GetFaceList();
740   for(vector<WFace*>::iterator f=wfaces.begin(), fend=wfaces.end();
741   f!=fend;
742   ++f)
743   {
744     vector<WXFaceLayer*>& faceLayers = ((WXFace*)(*f))->getSmoothLayers();
745     for(vector<WXFaceLayer*>::iterator wxfl = faceLayers.begin(), wxflend=faceLayers.end();
746     wxfl!=wxflend;
747     ++wxfl){
748       if ((*wxfl)->BuildSmoothEdge())
749             hasSmoothEdges = true;
750     }
751   }
752
753   if (hasSmoothEdges && !_computeRidgesAndValleys && !_computeSuggestiveContours) {
754     vector<WVertex*>& wvertices = iShape->getVertexList();
755     for(vector<WVertex*>::iterator wv=wvertices.begin(), wvend=wvertices.end(); 
756     wv!=wvend;
757     ++wv){
758       // Compute curvatures
759       WXVertex * wxv = dynamic_cast<WXVertex*>(*wv);
760       computeCurvatures(wxv);
761     }
762     _meanK1 /= (real)(_nPoints);
763     _meanKr /= (real)(_nPoints);
764   }
765 }