Fix for a static variable in BlenderStrokeRenderer::RenderStrokeRep() left after
[blender.git] / source / blender / freestyle / intern / view_map / FEdgeXDetector.cpp
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * The Original Code is Copyright (C) 2010 Blender Foundation.
19  * All rights reserved.
20  *
21  * The Original Code is: all of this file.
22  *
23  * Contributor(s): none yet.
24  *
25  * ***** END GPL LICENSE BLOCK *****
26  */
27
28 /** \file blender/freestyle/intern/view_map/FEdgeXDetector.cpp
29  *  \ingroup freestyle
30  *  \brief Detects/flags/builds extended features edges on the WXEdge structure
31  *  \author Stephane Grabli
32  *  \date 26/10/2003
33  */
34
35 #include <float.h>
36 #include <math.h>
37
38 #include "FEdgeXDetector.h"
39
40 #include "../geometry/GeomUtils.h"
41 #include "../geometry/normal_cycle.h"
42
43 #include "BKE_global.h"
44
45 void FEdgeXDetector::processShapes(WingedEdge& we)
46 {
47         bool progressBarDisplay = false;
48         Vec3r Min, Max;
49         vector<WShape*> wshapes = we.getWShapes();
50         WXShape *wxs;
51
52         if (_pProgressBar != NULL) {
53                 _pProgressBar->reset();
54                 _pProgressBar->setLabelText("Detecting feature lines");
55                 _pProgressBar->setTotalSteps(wshapes.size() * 3);
56                 _pProgressBar->setProgress(0);
57                 progressBarDisplay = true;
58         }
59
60         for (vector<WShape*>::const_iterator it = wshapes.begin(); it != wshapes.end(); it++) {
61                 if (_pRenderMonitor && _pRenderMonitor->testBreak())
62                         break;
63                 wxs = dynamic_cast<WXShape*>(*it);
64                 wxs->bbox(Min, Max);
65                 _bbox_diagonal = (Max - Min).norm();
66                 if (_changes) {
67                         vector<WFace*>& wfaces = wxs->GetFaceList();
68                         for (vector<WFace*>::iterator wf = wfaces.begin(), wfend = wfaces.end(); wf != wfend; ++wf) {
69                                 WXFace *wxf = dynamic_cast<WXFace*>(*wf);
70                                 wxf->Clear();
71                         }
72                         _computeViewIndependant = true;
73                 }
74                 else if (!(wxs)->getComputeViewIndependantFlag()) {
75                         wxs->Reset();
76                         _computeViewIndependant = false;
77                 }
78                 else {
79                         _computeViewIndependant = true;
80                 }
81                 preProcessShape(wxs);
82                 if (progressBarDisplay)
83                         _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
84                 processBorderShape(wxs);
85                 if (_computeMaterialBoundaries)
86                         processMaterialBoundaryShape(wxs);
87                 processCreaseShape(wxs);
88                 if (_computeRidgesAndValleys)
89                         processRidgesAndValleysShape(wxs);
90                 if (_computeSuggestiveContours)
91                         processSuggestiveContourShape(wxs);
92                 processSilhouetteShape(wxs);
93                 processEdgeMarksShape(wxs);
94                 if (progressBarDisplay)
95                         _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
96
97                 // build smooth edges:
98                 buildSmoothEdges(wxs);
99
100                 // Post processing for suggestive contours
101                 if (_computeSuggestiveContours)
102                         postProcessSuggestiveContourShape(wxs);
103                 if (progressBarDisplay)
104                         _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
105
106                 wxs->setComputeViewIndependantFlag(false);
107                 _computeViewIndependant = false;
108                 _changes = false;
109
110                 // reset user data
111                 (*it)->ResetUserData();
112         }
113 }
114
115 // GENERAL STUFF
116 ////////////////
117 void FEdgeXDetector::preProcessShape(WXShape *iWShape)
118 {
119         _meanK1 = 0;
120         _meanKr = 0;
121         _minK1 = FLT_MAX;
122         _maxK1 = -FLT_MAX;
123         _minKr = FLT_MAX;
124         _maxKr = -FLT_MAX;
125         _nPoints = 0;
126         _meanEdgeSize = iWShape->getMeanEdgeSize();
127
128         vector<WFace*>& wfaces = iWShape->GetFaceList();
129         vector<WFace*>::iterator f, fend;
130         // view dependant stuff
131         for (f = wfaces.begin(), fend = wfaces.end(); f != fend; ++f) {
132                 preProcessFace((WXFace*)(*f));
133         }
134
135         if (_computeRidgesAndValleys || _computeSuggestiveContours) {
136                 vector<WVertex*>& wvertices = iWShape->getVertexList();
137                 for (vector<WVertex*>::iterator wv = wvertices.begin(), wvend = wvertices.end();  wv != wvend; ++wv) {
138                         // Compute curvatures
139                         WXVertex *wxv = dynamic_cast<WXVertex*>(*wv);
140                         computeCurvatures(wxv);
141                 }
142                 _meanK1 /= (real)(_nPoints);
143                 _meanKr /= (real)(_nPoints);
144         }
145 }
146
147 void FEdgeXDetector::preProcessFace(WXFace *iFace)
148 {
149         Vec3r firstPoint = iFace->GetVertex(0)->GetVertex();
150         Vec3r N = iFace->GetNormal();
151
152         // Compute the dot product between V (=_Viewpoint - firstPoint) and N:
153         Vec3r V;
154         if (_orthographicProjection) {
155                 V = Vec3r(0.0, 0.0, _Viewpoint.z() - firstPoint.z());
156         }
157         else {
158                 V = Vec3r(_Viewpoint - firstPoint);
159         }
160         N.normalize();
161         V.normalize();
162         iFace->setDotP(N * V);
163
164         // compute the distance between the face center and the viewpoint:
165         if (_orthographicProjection) {
166                 iFace->setZ(iFace->center().z() - _Viewpoint.z());
167         }
168         else {
169                 Vec3r dist_vec(iFace->center() - _Viewpoint);
170                 iFace->setZ(dist_vec.norm());
171         }
172 }
173
174 void FEdgeXDetector::computeCurvatures(WXVertex *vertex)
175 {
176         // TODO: for some reason, the 'vertex' may have no associated edges
177         // (i.e., WVertex::_EdgeList is empty), which causes a crash due to
178         // a subsequent call of WVertex::_EdgeList.front().
179         if (vertex->GetEdges().empty()) {
180                 if (G.debug & G_DEBUG_FREESTYLE) {
181                         printf("Warning: WVertex %d has no associated edges.\n", vertex->GetId());
182                 }
183                 return;
184         }
185
186         // CURVATURE LAYER
187         // store all the curvature datas for each vertex
188
189         //soc unused - real K1, K2
190         real cos2theta, sin2theta;
191         Vec3r e1, n, v;
192         // one vertex curvature info :
193         CurvatureInfo *C;
194         float radius = _sphereRadius * _meanEdgeSize; 
195
196         // view independant stuff
197         if (_computeViewIndependant) {
198                 C = new CurvatureInfo();
199                 vertex->setCurvatures(C);
200                 OGF::NormalCycle ncycle;
201                 ncycle.begin();
202                 if (radius > 0) {
203                         OGF::compute_curvature_tensor(vertex, radius, ncycle);
204                 }
205                 else {
206                         OGF::compute_curvature_tensor_one_ring(vertex, ncycle);
207                 }
208                 ncycle.end();
209                 C->K1 = ncycle.kmin();
210                 C->K2 = ncycle.kmax();
211                 C->e1 = ncycle.Kmax(); //ncycle.kmin() * ncycle.Kmax();
212                 C->e2 = ncycle.Kmin(); //ncycle.kmax() * ncycle.Kmin();
213
214                 real absK1 = fabs(C->K1);
215                 _meanK1 += absK1;
216                 if (absK1 > _maxK1)
217                         _maxK1 = absK1;
218                 if (absK1 < _minK1)
219                         _minK1 = absK1;
220         }
221         // view dependant
222         C = vertex->curvatures();
223         if (C == 0)
224                 return;
225
226         // compute radial curvature :
227         n = C->e1 ^ C->e2;
228         if (_orthographicProjection) {
229                 v = Vec3r(0.0, 0.0, _Viewpoint.z() - vertex->GetVertex().z());
230         }
231         else {
232                 v = Vec3r(_Viewpoint - vertex->GetVertex());
233         }
234         C->er = v - (v * n) * n;
235         C->er.normalize();
236         e1 = C->e1;
237         e1.normalize();
238         cos2theta = C->er * e1;
239         cos2theta *= cos2theta;
240         sin2theta = 1 - cos2theta;
241         C->Kr = C->K1 * cos2theta + C->K2 * sin2theta;
242         real absKr = fabs(C->Kr);
243         _meanKr += absKr;
244         if (absKr > _maxKr)
245                 _maxKr = absKr;
246         if (absKr < _minKr)
247                 _minKr = absKr;
248
249         ++_nPoints;
250 }
251
252 // SILHOUETTE
253 /////////////
254 void FEdgeXDetector::processSilhouetteShape(WXShape *iWShape)
255 {
256         // Make a first pass on every polygons in order to compute all their silhouette relative values:
257         vector<WFace*>& wfaces = iWShape->GetFaceList();
258         vector<WFace*>::iterator f, fend;
259         for (f = wfaces.begin(), fend = wfaces.end(); f != fend; ++f) {
260                 ProcessSilhouetteFace((WXFace*)(*f));
261         }
262
263         // Make a pass on the edges to detect the silhouette edges that are not smooth
264         vector<WEdge*>::iterator we, weend;
265         vector<WEdge*> &wedges = iWShape->getEdgeList();
266         for (we = wedges.begin(), weend = wedges.end(); we != weend; ++we) {
267                 ProcessSilhouetteEdge((WXEdge*)(*we));
268         }
269 }
270
271 void FEdgeXDetector::ProcessSilhouetteFace(WXFace *iFace)
272 {
273         // SILHOUETTE LAYER
274         Vec3r normal;
275         // Compute the dot products between View direction and N at each vertex of the face:
276         Vec3r point;
277         int closestPointId = 0;
278         real dist, minDist = FLT_MAX;
279         int numVertices = iFace->numberOfVertices();
280         WXFaceLayer *faceLayer = new WXFaceLayer(iFace, Nature::SILHOUETTE, true);
281         for (int i = 0; i < numVertices; i++) {
282                 point = iFace->GetVertex(i)->GetVertex();
283                 normal = iFace->GetVertexNormal(i);
284                 normal.normalize();
285                 Vec3r V;
286                 if (_orthographicProjection) {
287                         V = Vec3r(0.0, 0.0, _Viewpoint.z() - point.z());
288                 }
289                 else {
290                         V = Vec3r(_Viewpoint - point);
291                 }
292                 V.normalize();
293                 real d = normal * V;
294                 faceLayer->PushDotP(d);
295                 // Find the point the closest to the viewpoint
296                 if (_orthographicProjection) {
297                         dist = point.z() - _Viewpoint.z();
298                 }
299                 else {
300                         Vec3r dist_vec(point - _Viewpoint);
301                         dist = dist_vec.norm();
302                 }
303                 if (dist < minDist) {
304                         minDist = dist;
305                         closestPointId = i;
306                 }
307         }
308         // Set the closest point id:
309         faceLayer->setClosestPointIndex(closestPointId);
310         // Add this layer to the face:
311         iFace->AddSmoothLayer(faceLayer);
312 }
313
314 void FEdgeXDetector::ProcessSilhouetteEdge(WXEdge *iEdge)
315 {
316         if (iEdge->nature() & Nature::BORDER)
317                 return;
318         // SILHOUETTE ?
319         //-------------
320         WXFace *fA = (WXFace *)iEdge->GetaOEdge()->GetaFace();
321         WXFace *fB = (WXFace *)iEdge->GetaOEdge()->GetbFace();
322
323         if ((fA->front()) ^ (fB->front())) { // fA->visible XOR fB->visible (true if one is 0 and the other is 1)
324                 // The only edges we want to set as silhouette edges in this way are the ones with 2 different normals for 1 vertex
325                 // for these two faces
326                 //--------------------
327                 // In reality we only test the normals for 1 of the 2 vertices.
328                 if (fA->GetVertexNormal(iEdge->GetaVertex()) == fB->GetVertexNormal(iEdge->GetaVertex()))
329                         return;
330                 iEdge->AddNature(Nature::SILHOUETTE);
331                 if (fB->front())
332                         iEdge->setOrder(1);
333                 else
334                         iEdge->setOrder(-1);
335         }
336 }
337
338 // BORDER
339 /////////
340 void FEdgeXDetector::processBorderShape(WXShape *iWShape)
341 {
342         if (!_computeViewIndependant)
343                 return;
344         // Make a pass on the edges to detect the BORDER
345         vector<WEdge*>::iterator we, weend;
346         vector<WEdge*> &wedges = iWShape->getEdgeList();
347         for (we = wedges.begin(), weend = wedges.end(); we != weend; ++we) {
348                 ProcessBorderEdge((WXEdge *)(*we));
349         }
350 }
351
352 void FEdgeXDetector::ProcessBorderEdge(WXEdge *iEdge)
353 {
354         // first check whether it is a border edge: BORDER ?
355         //---------
356         if (iEdge->GetaFace() == 0) {
357                 // it is a border edge
358                 iEdge->AddNature(Nature::BORDER);
359         }
360 }
361
362
363 // CREASE
364 /////////
365 void FEdgeXDetector::processCreaseShape(WXShape *iWShape)
366 {
367         if (!_computeViewIndependant)
368                 return;
369
370         // Make a pass on the edges to detect the CREASE 
371         vector<WEdge*>::iterator we, weend;
372         vector<WEdge*> &wedges = iWShape->getEdgeList();
373         for (we = wedges.begin(), weend = wedges.end(); we != weend; ++we) {
374                 ProcessCreaseEdge((WXEdge *)(*we));
375         }
376 }
377
378 void FEdgeXDetector::ProcessCreaseEdge(WXEdge *iEdge)
379 {
380         // CREASE ?
381         //---------
382         if (iEdge->nature() & Nature::BORDER)
383                 return;
384         WXFace *fA = (WXFace *)iEdge->GetaOEdge()->GetaFace();
385         WXFace *fB = (WXFace *)iEdge->GetaOEdge()->GetbFace();
386
387         WVertex *aVertex = iEdge->GetaVertex();
388         if ((fA->GetVertexNormal(aVertex) * fB->GetVertexNormal(aVertex)) <= _creaseAngle)
389                 iEdge->AddNature(Nature::CREASE);
390 }
391
392 // RIDGES AND VALLEYS
393 /////////////////////
394 void FEdgeXDetector::processRidgesAndValleysShape(WXShape *iWShape)
395 {
396         // Don't forget to add the built layer to the face at the end of the ProcessFace:
397         //iFace->AddSmoothLayer(faceLayer);
398
399         if (!_computeViewIndependant)
400                 return;
401
402         // Here the curvatures must already have been computed
403         vector<WFace*>& wfaces = iWShape->GetFaceList();
404         vector<WFace*>::iterator f, fend;
405         for (f = wfaces.begin(), fend = wfaces.end(); f != fend; ++f) {
406                 ProcessRidgeFace((WXFace*)(*f));
407         }
408 }
409
410
411 // RIDGES
412 /////////
413 void FEdgeXDetector::ProcessRidgeFace(WXFace *iFace)
414 {
415         WXFaceLayer *flayer = new WXFaceLayer(iFace, Nature::RIDGE|Nature::VALLEY, false);
416         iFace->AddSmoothLayer(flayer);
417
418         unsigned int numVertices = iFace->numberOfVertices();
419         for (unsigned int i = 0; i < numVertices; ++i) {
420                 WVertex *wv = iFace->GetVertex(i);
421                 WXVertex *wxv = dynamic_cast<WXVertex*>(wv);
422                 flayer->PushDotP(wxv->curvatures()->K1);
423         }
424
425         real threshold = 0;
426         //real threshold = _maxK1 - (_maxK1 - _meanK1) / 20.0;
427
428         if (flayer->nPosDotP() != numVertices) {
429                 if ((fabs(flayer->dotP(0)) < threshold) && (fabs(flayer->dotP(1)) < threshold) &&
430                     (fabs(flayer->dotP(2)) < threshold))
431                 {
432                         flayer->ReplaceDotP(0, 0);
433                         flayer->ReplaceDotP(1, 0);
434                         flayer->ReplaceDotP(2, 0);
435                 }
436         }
437 }
438
439 #if 0
440 void FEdgeXDetector::ProcessRidgeFace(WXFace *iFace)
441 {
442         // RIDGE LAYER
443         // Compute the RidgeFunction, that is the derivative of the ppal curvature along e1 at each vertex of the face
444         WVertex *v;
445         Vec3r v1v2;
446         real t;
447         vector<WXFaceLayer*> SmoothLayers;
448         WXFaceLayer *faceLayer;
449         Face_Curvature_Info *layer_info;
450         real K1_a(0), K1_b(0);
451         Vec3r Inter_a, Inter_b;
452
453         // find the ridge layer of the face
454         iFace->retrieveSmoothLayers(Nature::RIDGE, SmoothLayers);
455         if ( SmoothLayers.size()!=1 )
456                 return;
457         faceLayer = SmoothLayers[0];
458         // retrieve the curvature info of this layer
459         layer_info = (Face_Curvature_Info *)faceLayer->userdata;
460
461         int numVertices = iFace->numberOfVertices();
462         for (int i = 0; i < numVertices; i++) {
463                 v = iFace->GetVertex(i);
464                 // vec_curvature_info[i] contains the curvature info of this vertex
465                 Vec3r e2 = layer_info->vec_curvature_info[i]->K2*layer_info->vec_curvature_info[i]->e2;
466                 Vec3r e1 = layer_info->vec_curvature_info[i]->K1*layer_info->vec_curvature_info[i]->e1;
467                 e2.normalize();
468
469                 WVertex::face_iterator fit = v->faces_begin();
470                 WVertex::face_iterator fitend = v->faces_end();
471                 for (; fit != fitend; ++fit) {
472                         WXFace *wxf = dynamic_cast<WXFace*>(*fit);
473                         WOEdge *oppositeEdge;
474                         if (!(wxf->getOppositeEdge(v, oppositeEdge)))
475                                 continue;
476                         v1v2 = oppositeEdge->GetbVertex()->GetVertex() - oppositeEdge->GetaVertex()->GetVertex();
477                         GeomUtils::intersection_test res;
478                         res = GeomUtils::intersectRayPlane(oppositeEdge->GetaVertex()->GetVertex(), v1v2, e2, -(v->GetVertex()*e2),
479                                                            t, 1.0e-06);
480                         if ((res == GeomUtils::DO_INTERSECT) && (t >= 0.0) && (t <= 1.0)) {
481                                 vector<WXFaceLayer*> second_ridge_layer;
482                                 wxf->retrieveSmoothLayers(Nature::RIDGE, second_ridge_layer);
483                                 if (second_ridge_layer.size() != 1)
484                                         continue;
485                                 Face_Curvature_Info *second_layer_info = (Face_Curvature_Info*)second_ridge_layer[0]->userdata;
486
487                                 unsigned index1 = wxf->GetIndex(oppositeEdge->GetaVertex());
488                                 unsigned index2 = wxf->GetIndex(oppositeEdge->GetbVertex());
489                                 real K1_1 = second_layer_info->vec_curvature_info[index1]->K1;
490                                 real K1_2 = second_layer_info->vec_curvature_info[index2]->K1;
491                                 real K1 = (1.0 - t) * K1_1 + t * K1_2;
492                                 Vec3r inter((1.0 - t) * oppositeEdge->GetaVertex()->GetVertex() +
493                                             t * oppositeEdge->GetbVertex()->GetVertex());
494                                 Vec3r vtmp(inter - v->GetVertex());
495                                 // is it K1_a or K1_b ?
496                                 if (vtmp * e1 > 0) {
497                                         K1_b = K1;
498                                         Inter_b = inter;
499                                 }
500                                 else {
501                                         K1_a = K1;
502                                         Inter_a = inter;
503                                 }
504                         }
505                 }
506                 // Once we have K1 along the the ppal direction compute the derivative : K1b - K1a put it in DotP
507                 //real d = fabs(K1_b) - fabs(K1_a);
508                 real d = 0;
509                 real threshold = _meanK1 + (_maxK1 - _meanK1) / 7.0;
510                 //real threshold = _meanK1;
511                 //if ((fabs(K1_b) > threshold) || ((fabs(K1_a) > threshold)))
512                 d = (K1_b) - (K1_a) / (Inter_b - Inter_a).norm();
513                 faceLayer->PushDotP(d);
514                 //faceLayer->PushDotP(layer_info->vec_curvature_info[i]->K1);
515         }
516
517         // Make the values relevant by checking whether all principal directions have the "same" direction:
518         Vec3r e0((layer_info->vec_curvature_info[0]->K1 * layer_info->vec_curvature_info[0]->e1));
519         e0.normalize();
520         Vec3r e1((layer_info->vec_curvature_info[1]->K1 * layer_info->vec_curvature_info[1]->e1));
521         e1.normalize();
522         Vec3r e2((layer_info->vec_curvature_info[2]->K1 * layer_info->vec_curvature_info[2]->e1));
523         e2.normalize();
524         if (e0 * e1 < 0)
525                 // invert dotP[1]
526                 faceLayer->ReplaceDotP(1, -faceLayer->dotP(1));
527         if (e0 * e2 < 0)
528                 // invert dotP[2]
529                 faceLayer->ReplaceDotP(2, -faceLayer->dotP(2));
530
531 #if 0 // remove the weakest values;
532         real minDiff = (_maxK1 - _minK1) / 10.0;
533         real minDiff = _meanK1;
534         if ((faceLayer->dotP(0) < minDiff) && (faceLayer->dotP(1) < minDiff) && (faceLayer->dotP(2) < minDiff)) {
535                 faceLayer->ReplaceDotP(0, 0);
536                 faceLayer->ReplaceDotP(1, 0);
537                 faceLayer->ReplaceDotP(2, 0);
538         }
539 #endif
540 }
541 #endif
542
543 // SUGGESTIVE CONTOURS
544 //////////////////////
545
546 void FEdgeXDetector::processSuggestiveContourShape(WXShape *iWShape)
547 {
548         // Here the curvatures must already have been computed
549         vector<WFace*>& wfaces = iWShape->GetFaceList();
550         vector<WFace*>::iterator f, fend;
551         for (f = wfaces.begin(), fend = wfaces.end(); f != fend; ++f) {
552                 ProcessSuggestiveContourFace((WXFace*)(*f));
553         }
554 }
555
556 void FEdgeXDetector::ProcessSuggestiveContourFace(WXFace *iFace)
557 {
558         WXFaceLayer *faceLayer = new WXFaceLayer(iFace, Nature::SUGGESTIVE_CONTOUR, true);
559         iFace->AddSmoothLayer(faceLayer);
560
561         unsigned int numVertices = iFace->numberOfVertices();
562         for (unsigned int i = 0; i < numVertices; ++i) {
563                 WVertex *wv = iFace->GetVertex(i);
564                 WXVertex *wxv = dynamic_cast<WXVertex*>(wv);
565                 faceLayer->PushDotP(wxv->curvatures()->Kr);
566         }
567
568 #if 0 // FIXME: find a more clever way to compute the threshold
569         real threshold = _meanKr;
570         if (faceLayer->nPosDotP()!=numVertices) {
571                 if ((fabs(faceLayer->dotP(0)) < threshold) && (fabs(faceLayer->dotP(1)) < threshold) &&
572                     (fabs(faceLayer->dotP(2)) < threshold)) {
573                         faceLayer->ReplaceDotP(0, 0);
574                         faceLayer->ReplaceDotP(1, 0);
575                         faceLayer->ReplaceDotP(2, 0);
576                 }
577         }
578 #endif
579 }
580
581 void FEdgeXDetector::postProcessSuggestiveContourShape(WXShape *iShape)
582 {
583         vector<WFace*>& wfaces = iShape->GetFaceList();
584         vector<WFace*>::iterator f, fend;
585         for (f = wfaces.begin(), fend = wfaces.end(); f != fend; ++f) {
586                 postProcessSuggestiveContourFace((WXFace*)(*f));
587         }
588 }
589
590 void FEdgeXDetector::postProcessSuggestiveContourFace(WXFace *iFace)
591 {
592         // Compute the derivative of the radial curvature in the radial direction, at the two extremities of the smooth edge.
593         // If the derivative is smaller than a given threshold _kr_derivative_epsilon, 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 int 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(); fit != fitend; ++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, 1.0e-06);
642
643                         // If there is an intersection, compute the value of the derivative ath that point.
644                         if ((res == GeomUtils::DO_INTERSECT) && (t >= 0) && (t <= 1)) {
645                                 kr = t * opposite_vertex_a->curvatures()->Kr + (1 - t) * opposite_vertex_b->curvatures()->Kr;
646                                 inter = opposite_vertex_a->GetVertex() + t * opposite_edge->GetVec();
647                                 tmp_vec = inter - v->GetVertex();
648                                 // Is it kr1 or kr2?
649                                 if (tmp_vec * er_vec > 0) {
650                                         kr2 = kr;
651                                         inter2 = inter;
652                                 }
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 the extremities of the smooth edge.
670         WXSmoothEdge *sc_edge = sc_layer->getSmoothEdge();
671         WOEdge *sc_oedge = sc_edge->woea();
672         t = sc_edge->ta();
673         if (t * kr_derivatives[iFace->GetIndex(sc_oedge->GetaVertex())] +
674             (1 - t) * kr_derivatives[iFace->GetIndex(sc_oedge->GetbVertex())] < _kr_derivative_epsilon)
675         {
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         {
684                 sc_layer->removeSmoothEdge();
685         }
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(); we != weend; ++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 void FEdgeXDetector::processEdgeMarksShape(WXShape *iShape)
715 {
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(); we != weend; ++we) {
720                 ProcessEdgeMarks((WXEdge*)(*we));
721         }
722 }
723
724 void FEdgeXDetector::ProcessEdgeMarks(WXEdge *iEdge)
725 {
726         if (iEdge->GetMark()) {
727                 iEdge->AddNature(Nature::EDGE_MARK);
728         }
729 }
730
731 // Build Smooth edges
732 /////////////////////
733 void FEdgeXDetector::buildSmoothEdges(WXShape *iShape)
734 {
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(); f != fend; ++f) {
741                 vector<WXFaceLayer*>& faceLayers = ((WXFace*)(*f))->getSmoothLayers();
742                 for (vector<WXFaceLayer*>::iterator wxfl = faceLayers.begin(), wxflend = faceLayers.end();
743                      wxfl != wxflend;
744                      ++wxfl)
745                 {
746                         if ((*wxfl)->BuildSmoothEdge())
747                                 hasSmoothEdges = true;
748                 }
749         }
750
751         if (hasSmoothEdges && !_computeRidgesAndValleys && !_computeSuggestiveContours) {
752                 vector<WVertex*>& wvertices = iShape->getVertexList();
753                 for (vector<WVertex*>::iterator wv = wvertices.begin(), wvend = wvertices.end(); wv != wvend; ++wv) {
754                         // Compute curvatures
755                         WXVertex *wxv = dynamic_cast<WXVertex*>(*wv);
756                         computeCurvatures(wxv);
757                 }
758                 _meanK1 /= (real)(_nPoints);
759                 _meanKr /= (real)(_nPoints);
760         }
761 }