MSDM
 All Classes Namespaces Files Functions Variables Typedefs Pages
GL - Distance3D.h
1 /*------------------------------------------------------------------------
2 Copyright LIRIS, M2DisCo 2009, Licence QPL
3 Author: Guillaume Lavou¨¦
4 ------------------------------------------------------------------------*/
5 
6 #ifndef DISTANCE_3D_H
7 #define DISTANCE_3D_H
8 
9 
10 #include <CGAL/circulator.h>
11 #include <CGAL/HalfedgeDS_decorator.h>
12 #include <CGAL/HalfedgeDS_face_base.h>
13 #include<set>
14 #include<map>
15 #include<list>
16 #include<vector>
17 
18 #include<set>
19 #include<map>
20 #include<list>
21 #include<stack>
22 #include <algorithm>
23 
24 
25 
26 #include <time.h>
27 //#include"matrix.h"
28 
29 #define myPI 3.1415926535
30 
31 
32 
33 template <class _Poly>
34 class Distance3D
35 {
36  typedef _Poly Polyhedron;
37 
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;
42 
43  typedef typename Polyhedron::Vertex_handle Vertex_handle;
44  typedef typename Polyhedron::Halfedge_handle Halfedge_handle;
45  typedef typename Polyhedron::Facet_handle Facet_handle;
46 
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;
51 
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;
59 
60 
61 private:
62 
63 
64 
65 
66  bool m_IsLoaded,is_calculated;
67  FILE * times;
68 
69 public:
70 
71  Polyhedron *m_PolyOriginal;
72  Polyhedron m_PolyDegrad;
73 
74  bool IsMapDistanceCalculated;
75  Distance3D()
76  {
77  IsMapDistanceCalculated=false;
78 
79  m_IsLoaded=false;
80 
81 }
82 
83  Distance3D(Polyhedron *_PolyOriginal)
84  {
85  PolyOriginal(_PolyOriginal);
86  m_IsLoaded=false;
87 
88  }
89 
90 
91  const bool& IsLoaded() const { return m_IsLoaded; }
92  bool& IsLoaded() { return m_IsLoaded; }
93  void IsLoaded(const bool& p) { m_IsLoaded = p;}
94 
95  const Polyhedron& PolyOriginal() const { return *m_PolyOriginal; }
96  Polyhedron& PolyOriginal() { return *m_PolyOriginal; }
97  void PolyOriginal(Polyhedron* p) { m_PolyOriginal = p; }
98 
99  const Polyhedron& PolyDegrad() const { return m_PolyDegrad; }
100  Polyhedron& PolyDegrad() { return m_PolyDegrad; }
101  void PolyDegrad(const Polyhedron& p) {
102 
103  m_PolyDegrad = p;
104  for(Facet_iterator pFacet = m_PolyDegrad.facets_begin();pFacet!= m_PolyDegrad.facets_end();pFacet++)
105  {
106  Halfedge_around_facet_circulator pHalfedge = pFacet->facet_begin();
107  int size=CGAL::circulator_size(pHalfedge);
108 
109  if(size==4)
110  {
111  Halfedge_iterator V1=pHalfedge;
112  Halfedge_iterator V2=pHalfedge->next()->next();
113  m_PolyDegrad.split_facet(V1,V2);
114  }
115 
116  }
117 
118 
119  }
120 
121  double Processroughness_curve_Dual(double radius, double maxdim,bool IsGauss = true)
122  {
123 
124  m_PolyOriginal->MinNrmMeanCurvature(100000000);
125  m_PolyOriginal->MaxNrmMeanCurvature(0);
126 
127  m_PolyOriginal->MinNrmSalCurvature(100000000);
128  m_PolyOriginal->MaxNrmSalCurvature(0);
129 
130  m_PolyDegrad.MinNrmMeanCurvature(100000000);
131  m_PolyDegrad.MaxNrmMeanCurvature(0);
132 
133  m_PolyDegrad.MinNrmSalCurvature(100000000);
134  m_PolyDegrad.MaxNrmSalCurvature(0);
135 
136  double somme_roughness=0;
137  int NbVert=0;
138 
139  Vertex_iterator pVertexDeg=m_PolyDegrad.vertices_begin();
140  for(Vertex_iterator pVertex = m_PolyOriginal->vertices_begin();
141  pVertex != m_PolyOriginal->vertices_end();
142  pVertex++)
143  {
144  std::vector<double> TabDistance;
145  std::vector<Point> TabPoint;
146 
147  std::vector<double> TabDistanceDeg;
148  std::vector<Point> TabPointDeg;
149 
150  double moyenne,moyenneDeg;
151 
152  double var1=Processroughness_per_vertex_curve(m_PolyOriginal,(&(*pVertex)),radius,TabDistance,TabPoint,moyenne,maxdim,IsGauss);
153 
154  double var2=Processroughness_per_vertex_curve(&m_PolyDegrad,(&(*pVertexDeg)),radius,TabDistanceDeg,TabPointDeg,moyenneDeg,maxdim,IsGauss);
155 
156  float cov=static_cast<float>(ProcessCovariance((&(*pVertexDeg)),moyenne,moyenneDeg,TabDistance,TabPoint,TabDistanceDeg,TabPointDeg,maxdim,IsGauss));
157 
158  //double cov=1;
159  pVertexDeg->CourbureCoVariance=pVertex->CourbureCoVariance=cov;
160 
161 
162 
163  pVertexDeg++;
164 
165 
166  }
167  return 0;
168  }
169 
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)
171  {
172  double Cov1,Cov2;
173  double somme1=0;
174  double somme2=0;
175  //int Nb1;//Nb2;
176  double x1,x2;
177 
178  double buff,maxBuff;
179  Vector VBuff;
180 
181  //int IndMax;
182  int taille1=TabDistance.size();
183  int taille2=TabDistanceDeg.size();
184  if(taille1==0 || taille2==0)
185  int h=3;
186 
187  int ecart=abs(taille1-taille2);
188 
189  double sommewi=0;
190  for(int i=0;i<(int)TabDistance.size();i++)
191  {
192  int IndMin=i-2*ecart;
193  if(IndMin<0)
194  IndMin=0;
195 
196  int IndMax=i+2*ecart;
197  if(IndMax>(int)TabPointDeg.size()-1)
198  IndMax=TabPointDeg.size()-1;
199 
200 
201  x1=TabDistance[i];
202  x2=TabDistanceDeg[IndMin];
203  VBuff=TabPointDeg[IndMin]-TabPoint[i];
204  maxBuff=sqrt(VBuff*VBuff);
205 
206 
207 
208  for(int j=IndMin+1;j<=IndMax;j++)
209  {
210  VBuff=TabPointDeg[j]-TabPoint[i];
211  buff=sqrt(VBuff*VBuff);
212  if(buff<maxBuff)
213  {
214 
215  maxBuff=buff;
216  x2=TabDistanceDeg[j];
217 
218  }
219 
220  }
221 
222  if(IsGauss==false)
223  {
224  somme1+=(x1-moyenne)*(x2-moyenneDeg);
225  sommewi+=1;
226  }
227  else
228  {
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);
232 
233  sommewi+=wi;
234  somme1+=wi*(x1-moyenne)*(x2-moyenneDeg);
235  }
236  }
237  Cov1=somme1/sommewi;
238 
239  sommewi=0;
240  for(int i=0;i<(int)TabDistanceDeg.size();i++)
241  {
242 
243  int IndMin=i-2*ecart;
244  if(IndMin<0)
245  IndMin=0;
246 
247  int IndMax=i+2*ecart;
248  if(IndMax>(int)TabPoint.size()-1)
249  IndMax=TabPoint.size()-1;
250 
251  x1=TabDistanceDeg[i];
252  x2=TabDistance[IndMin];
253  VBuff=TabPoint[IndMin]-TabPointDeg[i];
254  maxBuff=sqrt(VBuff*VBuff);
255 
256 
257 
258  for(int j=IndMin+1;j<=IndMax;j++)
259  {
260  VBuff=TabPoint[j]-TabPointDeg[i];
261  buff=sqrt(VBuff*VBuff);
262  if(buff<maxBuff)
263  {
264  maxBuff=buff;
265  x2=TabDistance[j];
266 
267  }
268 
269  }
270  if(IsGauss==false)
271  {
272  somme2+=(x1-moyenneDeg)*(x2-moyenne);
273  sommewi+=1;
274  }
275  else
276  {
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);
280 
281  sommewi+=wi;
282  somme2+=wi*(x1-moyenneDeg)*(x2-moyenne);
283  }
284  }
285  Cov2=somme2/sommewi;
286 
287  if(Cov2>100000 || Cov2>100000)
288  int klj=3;
289  return (Cov2+Cov1)/2;
290 
291 
292  }
293 
294 
295 
296 
297  bool sphere_clip_vector(Point &O, double r,const Point &P, Vector &V)
298  {
299 
300  Vector W = P - O ;
301  double a = (V*V);
302  double b = 2.0 * V * W ;
303  double c = (W*W) - r*r ;
304  double delta = b*b - 4*a*c ;
305 
306 
307 
308  if( a==0)
309  return true ;
310 
311  if(delta < 0) {
312  // Should not happen, but happens sometimes (numerical precision)
313 
314  return true ;
315  }
316  double t = (- b + std::sqrt(delta)) / (2.0 * a) ;
317  if(t < 0.0) {
318 
319  // Should not happen, but happens sometimes (numerical precision)
320  return true ;
321  }
322  if(t >= 1.0) {
323  // Inside the sphere
324  return false ;
325  }
326 
327  if(t==0)
328  {
329 
330  t=0.01;
331  }
332 
333  V=V*t;
334 
335  return true ;
336  }
337 
338 
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)
340  {
341  std::set<Vertex*> vertices ;
342  Point O = pVertex->point() ;
343  std::stack<Vertex*> S ;
344  S.push(pVertex) ;
345  vertices.insert(pVertex) ;
346 
347 
348  int NbSommetInSphere=0;
349  double SommeDistance=0;
350 
351 
352  while(!S.empty())
353  {
354  Vertex* v = S.top() ;
355  S.pop() ;
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)
360  {
361  Point p1 = h->vertex()->point();
362  Point p2 = h->opposite()->vertex()->point();
363  Vector V = (p2-p1);
364  if(v==pVertex || V * (P - O) > 0.0)
365  {
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);
369 
370  NbSommetInSphere++;
372 
373  double DistancePondere;
374  Point PPondere;
375  if(len_old!=0)
376  {
377  DistancePondere=(1-len_edge/len_old)*h->vertex()->Kmax()+len_edge/len_old*h->opposite()->vertex()->Kmax();
378  PPondere=p1+V;
379  }
380  else
381  {
382  DistancePondere=h->opposite()->vertex()->Kmax();
383  PPondere=p2;
384  }
385 
386  TabDistance.push_back(DistancePondere);
387  TabPoint.push_back(PPondere);
388 
389  SommeDistance+=DistancePondere;
390 
391  if(!isect)
392  {
393 
394  Vertex_iterator w=h->opposite()->vertex();
395  if(vertices.find(&(*w)) == vertices.end())
396  {
397  vertices.insert(&(*w)) ;
398  S.push(&(*w)) ;
399  }
400  }
401 
402  }
403 
404  }
405 
406  }
407 
408  double moyenne=0;
409  double variance=0;
410 
411  if(IsGauss==false)
412  {
413  moyenne=SommeDistance/(double)NbSommetInSphere;
414 
415 
416 
417 
418  variance=0;
419  if(NbSommetInSphere!=0)
420  {
421  for(int i=0;i<NbSommetInSphere;i++)
422  variance+=pow(TabDistance[i]-moyenne,2);
423 
424  variance=variance/(double)NbSommetInSphere;
425  variance=sqrt(variance);
426  }
427  }
428  else
429  {//variance = 0.008
430  SommeDistance=0;
431  double SommeWi=0;
432  for(int i=0;i<NbSommetInSphere;i++)
433  {
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);
437  SommeWi+=wi;
438  SommeDistance+=TabDistance[i]*wi;
439  }
440  moyenne=SommeDistance/(double)SommeWi;
441 
442 
443  for(int i=0;i<NbSommetInSphere;i++)
444  {
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);
448 
449  variance+=wi*pow(TabDistance[i]-moyenne,2);
450  }
451 
452  variance=variance/SommeWi;
453  variance=sqrt(variance);
454 
455  }
456  //variance=sqrt(variance);
457 
458  //variance=log(variance+10);
459 
460 
461 
462 
464  pVertex->Kh(variance);
465  pVertex->Kg(moyenne);
466 
467  pVertex->CourbureMoyenne=moyenne;
468  pVertex->CourbureVariance=variance;
469 
470  if(variance<PolyUsed->MinNrmMeanCurvature())
471  PolyUsed->MinNrmMeanCurvature(variance);
472 
473  if(variance>PolyUsed->MaxNrmMeanCurvature())
474  PolyUsed->MaxNrmMeanCurvature(variance);
475 
476  if(moyenne<PolyUsed->MinNrmSalCurvature())
477  PolyUsed->MinNrmSalCurvature(moyenne);
478 
479  if(moyenne>PolyUsed->MaxNrmSalCurvature())
480  PolyUsed->MaxNrmSalCurvature(moyenne);
481 
482  if(moyenne<PolyUsed->MinNrmGaussCurvature())
483  PolyUsed->MinNrmGaussCurvature(moyenne);
484 
485  if(moyenne>PolyUsed->MaxNrmGaussCurvature())
486  PolyUsed->MaxNrmGaussCurvature(moyenne);
487 
488 
489 
490 
491  moyenneRet=moyenne;
492 
493 
494 
495  return variance;
496 
497  }
498 
499 
500 
501 
502 
503 
504 
505  void ComputeDistanceEcartNormalRoughnessPonderate(double Param,double &L1,double &L2,double &L3,double &L4)
506  {
507  double SommeDist=0;
508  double SommeDistCarre=0;
509  double SommeDist3=0;
510  double SommeDist4=0;
511  double Max=0;
512 
513 
514  int NbVertex=0;
515 
516  double PolyDegradRoughMin=m_PolyDegrad.MinNrmMeanCurvature();
517  double PolyDegradRoughMax=m_PolyDegrad.MaxNrmMeanCurvature();
518 
519  double PolyOriginalRoughMin=m_PolyOriginal->MinNrmMeanCurvature();
520  double PolyOriginalRoughMax=m_PolyOriginal->MaxNrmMeanCurvature();
521 
522 
523  m_PolyOriginal->MinNrmGaussCurvature(10000);
524  m_PolyOriginal->MaxNrmGaussCurvature(0);
525 
526  int CMPT=0;
527  Vertex_iterator pVertexDeg = m_PolyDegrad.vertices_begin();
528  for(Vertex_iterator pVertex = m_PolyOriginal->vertices_begin();
529  pVertex != m_PolyOriginal->vertices_end();
530  pVertex++)
531  {
532 
533 
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;
540 
541  double angleRad;
542 
544  double fact1=(fabs(MoyX-MoyY))/(std::max(MoyX,MoyY)+1);//(2*MoyX*MoyY)/(MoyX*MoyX+MoyY*MoyY);
545  double fact2=(fabs(SigX-SigY))/(std::max(SigX,SigY)+1);//2*SigX*SigY/(SigX*SigX+SigY*SigY);
546  double fact3=(fabs(SigX*SigY-SigXY))/(SigX*SigY+1);//fabs(SigXY/(SigX*SigY));
547 
548 
549  angleRad=pow((pow(fact1,3)+pow(fact2,3)+Param*pow(fact3,3))/(2.+Param),1./3.);
550 
551 
552  if(MoyX<1)
553  int h=2;
554 
555  CMPT++;
556 
557 
558 
559 
560 
561  SommeDist+=angleRad;
562  SommeDistCarre+=angleRad*angleRad;
563  SommeDist4+=angleRad*angleRad*angleRad*angleRad;
564  SommeDist3+=angleRad*angleRad*angleRad;
565  if(Max<angleRad)
566  Max=angleRad;
567 
568 
569  double Indice=angleRad;
570 
571  pVertex->Kg(Indice);
572  if(Indice<m_PolyOriginal->MinNrmGaussCurvature())
573  m_PolyOriginal->MinNrmGaussCurvature(Indice);
574 
575  if(Indice>m_PolyOriginal->MaxNrmGaussCurvature())
576  m_PolyOriginal->MaxNrmGaussCurvature(Indice);
577 
578  NbVertex++;
579  pVertexDeg++;
580 
581 
582 
583 
584  }
585 
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);
590 
591 
592  }
593 
594 
595 
596 };
597 
598 #endif
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
K::Point_3 Point
Definition: Header.h:23
Polyhedron class for the MSDM computation.
Definition: enriched_polyhedron.h:189