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