10 #include <CGAL/circulator.h>
11 #include <CGAL/HalfedgeDS_decorator.h>
12 #include <CGAL/HalfedgeDS_face_base.h>
29 #define myPI 3.1415926535
33 template <
class _Poly>
38 typedef typename Polyhedron::Halfedge_data_structure HDS;
39 typedef typename Polyhedron::Vertex Vertex;
40 typedef typename Polyhedron::Halfedge Halfedge;
41 typedef typename Polyhedron::Facet Facet;
43 typedef typename Polyhedron::Vertex_handle Vertex_handle;
44 typedef typename Polyhedron::Halfedge_handle Halfedge_handle;
45 typedef typename Polyhedron::Facet_handle Facet_handle;
47 typedef typename Polyhedron::Vertex_iterator Vertex_iterator;
48 typedef typename Polyhedron::Halfedge_iterator Halfedge_iterator;
49 typedef typename Polyhedron::Edge_iterator Edge_iterator;
50 typedef typename Polyhedron::Facet_iterator Facet_iterator;
52 typedef typename Polyhedron::Halfedge_around_facet_circulator
53 Halfedge_around_facet_circulator;
54 typedef typename Polyhedron::Halfedge_around_vertex_circulator
55 Halfedge_around_vertex_circulator;
56 typedef typename Polyhedron::Point_3
Point;
57 typedef typename Polyhedron::FT FT;
58 typedef typename Polyhedron::Vector Vector;
66 bool m_IsLoaded,is_calculated;
71 Polyhedron *m_PolyOriginal;
72 Polyhedron m_PolyDegrad;
74 bool IsMapDistanceCalculated;
77 IsMapDistanceCalculated=
false;
85 PolyOriginal(_PolyOriginal);
91 const bool& IsLoaded()
const {
return m_IsLoaded; }
92 bool& IsLoaded() {
return m_IsLoaded; }
93 void IsLoaded(
const bool& p) { m_IsLoaded = p;}
95 const Polyhedron& PolyOriginal()
const {
return *m_PolyOriginal; }
96 Polyhedron& PolyOriginal() {
return *m_PolyOriginal; }
97 void PolyOriginal(Polyhedron* p) { m_PolyOriginal = p; }
99 const Polyhedron& PolyDegrad()
const {
return m_PolyDegrad; }
100 Polyhedron& PolyDegrad() {
return m_PolyDegrad; }
101 void PolyDegrad(
const Polyhedron& p) {
104 for(Facet_iterator pFacet = m_PolyDegrad.facets_begin();pFacet!= m_PolyDegrad.facets_end();pFacet++)
106 Halfedge_around_facet_circulator pHalfedge = pFacet->facet_begin();
107 int size=CGAL::circulator_size(pHalfedge);
111 Halfedge_iterator V1=pHalfedge;
112 Halfedge_iterator V2=pHalfedge->next()->next();
113 m_PolyDegrad.split_facet(V1,V2);
121 double Processroughness_curve_Dual(
double radius,
double maxdim,
bool IsGauss =
true)
124 m_PolyOriginal->MinNrmMeanCurvature(100000000);
125 m_PolyOriginal->MaxNrmMeanCurvature(0);
127 m_PolyOriginal->MinNrmSalCurvature(100000000);
128 m_PolyOriginal->MaxNrmSalCurvature(0);
130 m_PolyDegrad.MinNrmMeanCurvature(100000000);
131 m_PolyDegrad.MaxNrmMeanCurvature(0);
133 m_PolyDegrad.MinNrmSalCurvature(100000000);
134 m_PolyDegrad.MaxNrmSalCurvature(0);
136 double somme_roughness=0;
139 Vertex_iterator pVertexDeg=m_PolyDegrad.vertices_begin();
140 for(Vertex_iterator pVertex = m_PolyOriginal->vertices_begin();
141 pVertex != m_PolyOriginal->vertices_end();
144 std::vector<double> TabDistance;
145 std::vector<Point> TabPoint;
147 std::vector<double> TabDistanceDeg;
148 std::vector<Point> TabPointDeg;
150 double moyenne,moyenneDeg;
156 float cov=
static_cast<float>(ProcessCovariance((&(*pVertexDeg)),moyenne,moyenneDeg,TabDistance,TabPoint,TabDistanceDeg,TabPointDeg,maxdim,IsGauss));
159 pVertexDeg->CourbureCoVariance=pVertex->CourbureCoVariance=cov;
170 double ProcessCovariance(Vertex* pVertex,
double moyenne,
double moyenneDeg,std::vector<double> TabDistance,std::vector<Point> TabPoint,std::vector<double> TabDistanceDeg,std::vector<Point> TabPointDeg,
double dim,
bool IsGauss)
182 int taille1=TabDistance.size();
183 int taille2=TabDistanceDeg.size();
184 if(taille1==0 || taille2==0)
187 int ecart=abs(taille1-taille2);
190 for(
int i=0;i<(int)TabDistance.size();i++)
192 int IndMin=i-2*ecart;
196 int IndMax=i+2*ecart;
197 if(IndMax>(
int)TabPointDeg.size()-1)
198 IndMax=TabPointDeg.size()-1;
202 x2=TabDistanceDeg[IndMin];
203 VBuff=TabPointDeg[IndMin]-TabPoint[i];
204 maxBuff=sqrt(VBuff*VBuff);
208 for(
int j=IndMin+1;j<=IndMax;j++)
210 VBuff=TabPointDeg[j]-TabPoint[i];
211 buff=sqrt(VBuff*VBuff);
216 x2=TabDistanceDeg[j];
224 somme1+=(x1-moyenne)*(x2-moyenneDeg);
229 Vector DistancePt=TabPoint[i]-pVertex->point();
230 double distPt=sqrt(DistancePt*DistancePt);
231 double wi=1/0.008/dim/sqrt(2*3.141592)*exp(-(distPt*distPt)/2/0.008/0.008/dim/dim);
234 somme1+=wi*(x1-moyenne)*(x2-moyenneDeg);
240 for(
int i=0;i<(int)TabDistanceDeg.size();i++)
243 int IndMin=i-2*ecart;
247 int IndMax=i+2*ecart;
248 if(IndMax>(
int)TabPoint.size()-1)
249 IndMax=TabPoint.size()-1;
251 x1=TabDistanceDeg[i];
252 x2=TabDistance[IndMin];
253 VBuff=TabPoint[IndMin]-TabPointDeg[i];
254 maxBuff=sqrt(VBuff*VBuff);
258 for(
int j=IndMin+1;j<=IndMax;j++)
260 VBuff=TabPoint[j]-TabPointDeg[i];
261 buff=sqrt(VBuff*VBuff);
272 somme2+=(x1-moyenneDeg)*(x2-moyenne);
277 Vector DistancePt=TabPointDeg[i]-pVertex->point();
278 double distPt=sqrt(DistancePt*DistancePt);
279 double wi=1/0.008/dim/sqrt(2*3.141592)*exp(-(distPt*distPt)/2/0.008/0.008/dim/dim);
282 somme2+=wi*(x1-moyenneDeg)*(x2-moyenne);
287 if(Cov2>100000 || Cov2>100000)
289 return (Cov2+Cov1)/2;
297 bool sphere_clip_vector(Point &O,
double r,
const Point &P, Vector &V)
302 double b = 2.0 * V * W ;
303 double c = (W*W) - r*r ;
304 double delta = b*b - 4*a*c ;
316 double t = (- b + std::sqrt(delta)) / (2.0 * a) ;
339 double Processroughness_per_vertex_curve(Polyhedron * PolyUsed,Vertex* pVertex,
double radius,std::vector<double> &TabDistance,std::vector<Point> &TabPoint,
double & moyenneRet,
double dim,
bool IsGauss =
false)
341 std::set<Vertex*> vertices ;
342 Point O = pVertex->point() ;
343 std::stack<Vertex*> S ;
345 vertices.insert(pVertex) ;
348 int NbSommetInSphere=0;
349 double SommeDistance=0;
354 Vertex* v = S.top() ;
356 Point P = v->point() ;
357 Halfedge_around_vertex_circulator h = v->vertex_begin();
358 Halfedge_around_vertex_circulator pHalfedgeStart = h;
359 CGAL_For_all(h,pHalfedgeStart)
361 Point p1 = h->vertex()->point();
362 Point p2 = h->opposite()->vertex()->point();
364 if(v==pVertex || V * (P - O) > 0.0)
366 double len_old = std::sqrt(V*V);
367 bool isect = sphere_clip_vector(O, radius, P, V) ;
368 double len_edge = std::sqrt(V*V);
373 double DistancePondere;
377 DistancePondere=(1-len_edge/len_old)*h->vertex()->Kmax()+len_edge/len_old*h->opposite()->vertex()->Kmax();
382 DistancePondere=h->opposite()->vertex()->Kmax();
386 TabDistance.push_back(DistancePondere);
387 TabPoint.push_back(PPondere);
389 SommeDistance+=DistancePondere;
394 Vertex_iterator w=h->opposite()->vertex();
395 if(vertices.find(&(*w)) == vertices.end())
397 vertices.insert(&(*w)) ;
413 moyenne=SommeDistance/(double)NbSommetInSphere;
419 if(NbSommetInSphere!=0)
421 for(
int i=0;i<NbSommetInSphere;i++)
422 variance+=pow(TabDistance[i]-moyenne,2);
424 variance=variance/(double)NbSommetInSphere;
425 variance=sqrt(variance);
432 for(
int i=0;i<NbSommetInSphere;i++)
434 Vector DistancePt=TabPoint[i]-pVertex->point();
435 double distPt=sqrt(DistancePt*DistancePt);
436 double wi=1/0.008/dim/sqrt(2*3.141592)*exp(-(distPt*distPt)/2/0.008/0.008/dim/dim);
438 SommeDistance+=TabDistance[i]*wi;
440 moyenne=SommeDistance/(double)SommeWi;
443 for(
int i=0;i<NbSommetInSphere;i++)
445 Vector DistancePt=TabPoint[i]-pVertex->point();
446 double distPt=sqrt(DistancePt*DistancePt);
447 double wi=1/0.008/dim/sqrt(2*3.141592)*exp(-(distPt*distPt)/2/0.008/0.008/dim/dim);
449 variance+=wi*pow(TabDistance[i]-moyenne,2);
452 variance=variance/SommeWi;
453 variance=sqrt(variance);
464 pVertex->Kh(variance);
465 pVertex->Kg(moyenne);
467 pVertex->CourbureMoyenne=moyenne;
468 pVertex->CourbureVariance=variance;
470 if(variance<PolyUsed->MinNrmMeanCurvature())
471 PolyUsed->MinNrmMeanCurvature(variance);
473 if(variance>PolyUsed->MaxNrmMeanCurvature())
474 PolyUsed->MaxNrmMeanCurvature(variance);
476 if(moyenne<PolyUsed->MinNrmSalCurvature())
477 PolyUsed->MinNrmSalCurvature(moyenne);
479 if(moyenne>PolyUsed->MaxNrmSalCurvature())
480 PolyUsed->MaxNrmSalCurvature(moyenne);
482 if(moyenne<PolyUsed->MinNrmGaussCurvature())
483 PolyUsed->MinNrmGaussCurvature(moyenne);
485 if(moyenne>PolyUsed->MaxNrmGaussCurvature())
486 PolyUsed->MaxNrmGaussCurvature(moyenne);
505 void ComputeDistanceEcartNormalRoughnessPonderate(
double Param,
double &L1,
double &L2,
double &L3,
double &L4)
508 double SommeDistCarre=0;
516 double PolyDegradRoughMin=m_PolyDegrad.MinNrmMeanCurvature();
517 double PolyDegradRoughMax=m_PolyDegrad.MaxNrmMeanCurvature();
519 double PolyOriginalRoughMin=m_PolyOriginal->MinNrmMeanCurvature();
520 double PolyOriginalRoughMax=m_PolyOriginal->MaxNrmMeanCurvature();
523 m_PolyOriginal->MinNrmGaussCurvature(10000);
524 m_PolyOriginal->MaxNrmGaussCurvature(0);
527 Vertex_iterator pVertexDeg = m_PolyDegrad.vertices_begin();
528 for(Vertex_iterator pVertex = m_PolyOriginal->vertices_begin();
529 pVertex != m_PolyOriginal->vertices_end();
535 double MoyX=pVertex->CourbureMoyenne;
536 double MoyY=pVertexDeg->CourbureMoyenne;
537 double SigX=pVertex->CourbureVariance;
538 double SigY=pVertexDeg->CourbureVariance;
539 double SigXY=pVertex->CourbureCoVariance;
544 double fact1=(fabs(MoyX-MoyY))/(std::max(MoyX,MoyY)+1);
545 double fact2=(fabs(SigX-SigY))/(std::max(SigX,SigY)+1);
546 double fact3=(fabs(SigX*SigY-SigXY))/(SigX*SigY+1);
549 angleRad=pow((pow(fact1,3)+pow(fact2,3)+Param*pow(fact3,3))/(2.+Param),1./3.);
562 SommeDistCarre+=angleRad*angleRad;
563 SommeDist4+=angleRad*angleRad*angleRad*angleRad;
564 SommeDist3+=angleRad*angleRad*angleRad;
569 double Indice=angleRad;
572 if(Indice<m_PolyOriginal->MinNrmGaussCurvature())
573 m_PolyOriginal->MinNrmGaussCurvature(Indice);
575 if(Indice>m_PolyOriginal->MaxNrmGaussCurvature())
576 m_PolyOriginal->MaxNrmGaussCurvature(Indice);
586 L2=sqrt(SommeDistCarre/(
double)NbVertex);
587 L1=SommeDist/(double)NbVertex;
588 L4=pow(SommeDist4/(
double)NbVertex,0.25);
589 L3=pow(SommeDist3/(
double)NbVertex,0.33333);
double Processroughness_per_vertex_curve(Polyhedron *PolyUsed, Vertex *pVertex, double radius, std::vector< double > &TabDistance, std::vector< Point > &TabPoint, double &moyenneRet, double dim, bool IsGauss=false)
Definition: GL - Distance3D.h:339
Estimation of the MSDM between 2 polyhedron with curvature information.
Definition: Distance3D.h:16
Polyhedron class for the MSDM computation.
Definition: enriched_polyhedron.h:189