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