Freestyle: minor optimization for space from mesh importing to feature edge detection.
[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  * ***** END GPL LICENSE BLOCK *****
19  */
20
21 /** \file blender/freestyle/intern/view_map/FEdgeXDetector.cpp
22  *  \ingroup freestyle
23  *  \brief Detects/flags/builds extended features edges on the WXEdge structure
24  *  \author Stephane Grabli
25  *  \date 26/10/2003
26  */
27
28 #include <float.h>
29
30 #include "FEdgeXDetector.h"
31
32 #include "../geometry/GeomUtils.h"
33 #include "../geometry/normal_cycle.h"
34
35 #include "BKE_global.h"
36
37 namespace Freestyle {
38
39 void FEdgeXDetector::processShapes(WingedEdge& we)
40 {
41         bool progressBarDisplay = false;
42         Vec3r Min, Max;
43         vector<WShape*> wshapes = we.getWShapes();
44         WXShape *wxs;
45
46         if (_pProgressBar != NULL) {
47                 _pProgressBar->reset();
48                 _pProgressBar->setLabelText("Detecting feature lines");
49                 _pProgressBar->setTotalSteps(wshapes.size() * 3);
50                 _pProgressBar->setProgress(0);
51                 progressBarDisplay = true;
52         }
53
54         for (vector<WShape*>::const_iterator it = wshapes.begin(); it != wshapes.end(); it++) {
55                 if (_pRenderMonitor && _pRenderMonitor->testBreak())
56                         break;
57                 wxs = dynamic_cast<WXShape*>(*it);
58 #if 0
59                 wxs->bbox(Min, Max);
60                 _bbox_diagonal = (Max - Min).norm();
61 #endif
62                 if (_changes) {
63                         vector<WFace*>& wfaces = wxs->GetFaceList();
64                         for (vector<WFace*>::iterator wf = wfaces.begin(), wfend = wfaces.end(); wf != wfend; ++wf) {
65                                 WXFace *wxf = dynamic_cast<WXFace*>(*wf);
66                                 wxf->Clear();
67                         }
68                         _computeViewIndependent = true;
69                 }
70                 else if (!(wxs)->getComputeViewIndependentFlag()) {
71                         wxs->Reset();
72                         _computeViewIndependent = false;
73                 }
74                 else {
75                         _computeViewIndependent = true;
76                 }
77                 preProcessShape(wxs);
78                 if (progressBarDisplay)
79                         _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
80                 processBorderShape(wxs);
81                 if (_computeMaterialBoundaries)
82                         processMaterialBoundaryShape(wxs);
83                 processCreaseShape(wxs);
84                 if (_computeRidgesAndValleys)
85                         processRidgesAndValleysShape(wxs);
86                 if (_computeSuggestiveContours)
87                         processSuggestiveContourShape(wxs);
88                 processSilhouetteShape(wxs);
89                 processEdgeMarksShape(wxs);
90                 if (progressBarDisplay)
91                         _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
92
93                 // build smooth edges:
94                 buildSmoothEdges(wxs);
95
96                 // Post processing for suggestive contours
97                 if (_computeSuggestiveContours)
98                         postProcessSuggestiveContourShape(wxs);
99                 if (progressBarDisplay)
100                         _pProgressBar->setProgress(_pProgressBar->getProgress() + 1);
101
102                 wxs->setComputeViewIndependentFlag(false);
103                 _computeViewIndependent = false;
104                 _changes = false;
105
106                 // reset user data
107                 (*it)->ResetUserData();
108         }
109 }
110
111 // GENERAL STUFF
112 ////////////////
113 void FEdgeXDetector::preProcessShape(WXShape *iWShape)
114 {
115         _meanK1 = 0;
116         _meanKr = 0;
117         _minK1 = FLT_MAX;
118         _maxK1 = -FLT_MAX;
119         _minKr = FLT_MAX;
120         _maxKr = -FLT_MAX;
121         _nPoints = 0;
122 #if 0
123         _meanEdgeSize = iWShape->getMeanEdgeSize();
124 #else
125         _meanEdgeSize = iWShape->ComputeMeanEdgeSize();
126 #endif
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 independent stuff
197         if (_computeViewIndependent) {
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
325                 // for 1 vertex 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 (!_computeViewIndependent)
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 (!_computeViewIndependent)
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 (!_computeViewIndependent)
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 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
593         // of the smooth edge.
594         // If the derivative is smaller than a given threshold _kr_derivative_epsilon, discard the edge.
595
596         // Find the suggestive contour layer of the face (zero or one edge).
597         vector<WXFaceLayer*> sc_layers;
598         iFace->retrieveSmoothEdgesLayers(Nature::SUGGESTIVE_CONTOUR, sc_layers);
599         if (sc_layers.empty())
600                 return;
601
602         WXFaceLayer *sc_layer;
603         sc_layer = sc_layers[0];
604
605         // Compute the derivative value at each vertex of the face, and add it in a vector.
606         vector<real> kr_derivatives;
607
608         unsigned vertices_nb = iFace->numberOfVertices();
609         WXVertex *v, *opposite_vertex_a, *opposite_vertex_b;
610         WXFace *wxf;
611         WOEdge *opposite_edge;
612         Vec3r normal_vec, radial_normal_vec, er_vec, v_vec, inter, inter1, inter2, tmp_vec;
613         GeomUtils::intersection_test res;
614         real kr(0), kr1(0), kr2(0), t;
615
616         for (unsigned int i = 0; i < vertices_nb; ++i) {
617                 v = (WXVertex *)(iFace->GetVertex(i));
618
619                 // v is a singular vertex, skip it.
620                 if (v->isBoundary()) {
621                         kr_derivatives.push_back(0);
622                         continue;
623                 }
624
625                 v_vec = v->GetVertex();
626                 er_vec = v->curvatures()->er;
627
628                 // For each vertex, iterate on its adjacent faces.
629                 for (WVertex::face_iterator fit = v->faces_begin(), fitend = v->faces_end(); fit != fitend; ++fit) {
630                         wxf = dynamic_cast<WXFace*>(*fit);
631                         if (!wxf->getOppositeEdge(v, opposite_edge))
632                                 continue;
633
634                         opposite_vertex_a = (WXVertex *)opposite_edge->GetaVertex();
635                         opposite_vertex_b = (WXVertex *)opposite_edge->GetbVertex();
636                         normal_vec = wxf->GetVertexNormal(v); // FIXME: what about e1 ^ e2 ?
637                         radial_normal_vec = er_vec ^ normal_vec;
638
639                         // Test wether the radial plan intersects with the edge at the opposite of v.
640                         res = GeomUtils::intersectRayPlane(opposite_vertex_a->GetVertex(), opposite_edge->GetVec(),
641                                                            radial_normal_vec, -(v_vec * radial_normal_vec),
642                                                            t, 1.0e-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                                 }
654                                 else {
655                                         kr1 = kr;
656                                         inter1 = inter;
657                                 }
658                         }
659                 }
660
661                 // Now we have kr1 and kr2 along the radial direction, for one vertex of iFace.
662                 // We have to compute the derivative of kr for that vertex, equal to:
663                 // (kr2 - kr1) / dist(inter1, inter2).
664                 // Then we add it to the vector of derivatives.
665                 v->curvatures()->dKr = (kr2 - kr1) / (inter2 - inter1).norm();
666                 kr_derivatives.push_back(v->curvatures()->dKr);
667         }
668
669         // At that point, we have the derivatives for each vertex of iFace.
670         // All we have to do now is to use linear interpolation to compute the values at 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         {
677                 sc_layer->removeSmoothEdge();
678                 return;
679         }
680         sc_oedge = sc_edge->woeb();
681         t = sc_edge->tb();
682         if (t * kr_derivatives[iFace->GetIndex(sc_oedge->GetaVertex())] +
683             (1 - t) * kr_derivatives[iFace->GetIndex(sc_oedge->GetbVertex())] < _kr_derivative_epsilon)
684         {
685                 sc_layer->removeSmoothEdge();
686         }
687 }
688
689 // MATERIAL_BOUNDARY
690 ////////////////////
691 void FEdgeXDetector::processMaterialBoundaryShape(WXShape *iWShape)
692 {
693         if (!_computeViewIndependent)
694                 return;
695         // Make a pass on the edges to detect material boundaries
696         vector<WEdge*>::iterator we, weend;
697         vector<WEdge*> &wedges = iWShape->getEdgeList();
698         for (we = wedges.begin(), weend = wedges.end(); we != weend; ++we) {
699                 ProcessMaterialBoundaryEdge((WXEdge *)(*we));
700         }
701 }
702
703 void FEdgeXDetector::ProcessMaterialBoundaryEdge(WXEdge *iEdge)
704 {
705         // check whether the edge is a material boundary?
706         WFace *aFace = iEdge->GetaFace();
707         WFace *bFace = iEdge->GetbFace();
708         if (aFace && bFace && aFace->frs_materialIndex() != bFace->frs_materialIndex()) {
709                 iEdge->AddNature(Nature::MATERIAL_BOUNDARY);
710         }
711 }
712
713 // EDGE MARKS
714 /////////////
715 void FEdgeXDetector::processEdgeMarksShape(WXShape *iShape)
716 {
717         // Make a pass on the edges to detect material boundaries
718         vector<WEdge*>::iterator we, weend;
719         vector<WEdge*> &wedges = iShape->getEdgeList();
720         for (we = wedges.begin(), weend = wedges.end(); we != weend; ++we) {
721                 ProcessEdgeMarks((WXEdge *)(*we));
722         }
723 }
724
725 void FEdgeXDetector::ProcessEdgeMarks(WXEdge *iEdge)
726 {
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 {
736         bool hasSmoothEdges = false;
737
738         // Make a last pass to build smooth edges from the previous stored values:
739         //--------------------------------------------------------------------------
740         vector<WFace*>& wfaces = iShape->GetFaceList();
741         for (vector<WFace *>::iterator f = wfaces.begin(), fend = wfaces.end(); f != fend; ++f) {
742                 vector<WXFaceLayer *>& faceLayers = ((WXFace *)(*f))->getSmoothLayers();
743                 for (vector<WXFaceLayer *>::iterator wxfl = faceLayers.begin(), wxflend = faceLayers.end();
744                      wxfl != wxflend;
745                      ++wxfl)
746                 {
747                         if ((*wxfl)->BuildSmoothEdge())
748                                 hasSmoothEdges = true;
749                 }
750         }
751
752         if (hasSmoothEdges && !_computeRidgesAndValleys && !_computeSuggestiveContours) {
753                 vector<WVertex *>& wvertices = iShape->getVertexList();
754                 for (vector<WVertex*>::iterator wv = wvertices.begin(), wvend = wvertices.end(); wv != wvend; ++wv) {
755                         // Compute curvatures
756                         WXVertex *wxv = dynamic_cast<WXVertex *>(*wv);
757                         computeCurvatures(wxv);
758                 }
759                 _meanK1 /= (real)(_nPoints);
760                 _meanKr /= (real)(_nPoints);
761         }
762 }
763
764 } /* namespace Freestyle */