$darkmode
VCG Library
point_sampling.h
1 /****************************************************************************
2 * VCGLib o o *
3 * Visual and Computer Graphics Library o o *
4 * _ O _ *
5 * Copyright(C) 2004-2016 \/)\/ *
6 * Visual Computing Lab /\/| *
7 * ISTI - Italian National Research Council | *
8 * \ *
9 * All rights reserved. *
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 * This program is distributed in the hope that it will be useful, *
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
19 * GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
20 * for more details. *
21 * *
22 ****************************************************************************/
23 /****************************************************************************
24 
25 The sampling Class has a set of static functions, that you can call to sample the surface of a mesh.
26 Each function is templated on the mesh and on a Sampler object s.
27 Each function calls many time the sample object with the sampling point as parameter.
28 
29 Sampler Classes and Sampling algorithms are independent.
30 Sampler classes exploits the sample that are generated with various algorithms.
31 For example, you can compute Hausdorff distance (that is a sampler) using various
32 sampling strategies (montecarlo, stratified etc).
33 
34 ****************************************************************************/
35 #ifndef __VCGLIB_POINT_SAMPLING
36 #define __VCGLIB_POINT_SAMPLING
37 
38 #include <random>
39 
40 #include <vcg/math/random_generator.h>
41 #include <vcg/complex/algorithms/closest.h>
42 #include <vcg/space/index/spatial_hashing.h>
43 #include <vcg/complex/algorithms/hole.h>
44 #include <vcg/complex/algorithms/stat.h>
45 #include <vcg/complex/algorithms/create/platonic.h>
46 #include <vcg/complex/algorithms/update/normal.h>
47 #include <vcg/complex/algorithms/update/bounding.h>
48 #include <vcg/space/segment2.h>
49 #include <vcg/space/index/grid_static_ptr.h>
50 
51 namespace vcg
52 {
53 namespace tri
54 {
57 
70 template <class MeshType>
72 {
73 public:
74  typedef typename MeshType::ScalarType ScalarType;
75  typedef typename MeshType::CoordType CoordType;
76  typedef typename MeshType::VertexType VertexType;
77  typedef typename MeshType::EdgeType EdgeType;
78  typedef typename MeshType::FaceType FaceType;
79 
80  void reset()
81  {
82  sampleVec->clear();
83  }
84 
86  {
87  sampleVec = new std::vector<CoordType>();
88  vectorOwner=true;
89  }
90 
91  TrivialSampler(std::vector<CoordType> &Vec)
92  {
93  sampleVec = &Vec;
94  vectorOwner=false;
95  reset();
96  }
97 
99  {
100  if(vectorOwner) delete sampleVec;
101  }
102 
103 private:
104  std::vector<CoordType> *sampleVec;
105  bool vectorOwner;
106 public:
107 
108  std::vector<CoordType> &SampleVec()
109  {
110  return *sampleVec;
111  }
112 
113  void AddVert(const VertexType &p)
114  {
115  sampleVec->push_back(p.cP());
116  }
117  void AddEdge(const EdgeType& e, ScalarType u ) // u==0 -> v(0) u==1 -> v(1);
118  {
119  sampleVec->push_back(e.cV(0)->cP()*(1.0-u)+e.cV(1)->cP()*u);
120  }
121 
122  void AddFace(const FaceType &f, const CoordType &p)
123  {
124  sampleVec->push_back(f.cP(0)*p[0] + f.cP(1)*p[1] +f.cP(2)*p[2] );
125  }
126 
127  void AddTextureSample(const FaceType &, const CoordType &, const Point2i &, float )
128  {
129  // Retrieve the color of the sample from the face f using the barycentric coord p
130  // and write that color in a texture image at position <tp[0], texHeight-tp[1]>
131  // if edgeDist is > 0 then the corrisponding point is affecting face color even if outside the face area (in texture space)
132  }
133 }; // end class TrivialSampler
134 
135 template <class MeshType>
137 {
138 public:
139  typedef typename MeshType::ScalarType ScalarType;
140  typedef typename MeshType::CoordType CoordType;
141  typedef typename MeshType::VertexType VertexType;
142  typedef typename MeshType::EdgeType EdgeType;
143  typedef typename MeshType::FaceType FaceType;
144 
147 
148  void reset()
149  {
150  sampleVec.clear();
151  }
152 
153 public:
154  std::vector<VertexType *> sampleVec;
155 
156  void AddVert(VertexType &p)
157  {
158  sampleVec.push_back(&p);
159  }
160 
161  void AddEdge(const EdgeType& e, ScalarType u ) // u==0 -> v(0) u==1 -> v(1);
162  {
163  if( u < 0.5 )
164  sampleVec.push_back(e.cV(0));
165  else
166  sampleVec.push_back(e.cV(1));
167  }
168 
169  // This sampler should be used only for getting vertex pointers. Meaningless in other case.
170  void AddFace(const FaceType &, const CoordType &) { assert(0); }
171  void AddTextureSample(const FaceType &, const CoordType &, const Point2i &, float ) { assert(0); }
172 }; // end class TrivialSampler
173 
174 
175 template <class MeshType>
177 {
178 public:
179  typedef typename MeshType::VertexType VertexType;
180  typedef typename MeshType::FaceType FaceType;
181  typedef typename MeshType::CoordType CoordType;
182 
183  MeshSampler(MeshType &_m):m(_m){
184  perFaceNormal = false;
185  }
186  MeshType &m;
187 
188  bool perFaceNormal; // default false; if true the sample normal is the face normal, otherwise it is interpolated
189 
190  void reset()
191  {
192  m.Clear();
193  }
194 
195  void AddVert(const VertexType &p)
196  {
198  m.vert.back().ImportData(p);
199  }
200 
201  void AddFace(const FaceType &f, CoordType p)
202  {
204  m.vert.back().P() = f.cP(0)*p[0] + f.cP(1)*p[1] +f.cP(2)*p[2];
205  if(perFaceNormal) m.vert.back().N() = f.cN();
206  else m.vert.back().N() = f.cV(0)->N()*p[0] + f.cV(1)->N()*p[1] + f.cV(2)->N()*p[2];
207  if(tri::HasPerVertexQuality(m) )
208  m.vert.back().Q() = f.cV(0)->Q()*p[0] + f.cV(1)->Q()*p[1] + f.cV(2)->Q()*p[2];
209  }
210 }; // end class MeshSampler
211 
212 
213 
214 /* This sampler is used to perform compute the Hausdorff measuring.
215  * It keep internally the spatial indexing structure used to find the closest point
216  * and the partial integration results needed to compute the average and rms error values.
217  * Averaged values assume that the samples are equi-distributed (e.g. a good unbiased montecarlo sampling of the surface).
218  */
219 template <class MeshType>
221 {
222  typedef typename MeshType::FaceType FaceType;
223  typedef typename MeshType::VertexType VertexType;
224  typedef typename MeshType::CoordType CoordType;
225  typedef typename MeshType::ScalarType ScalarType;
226  typedef GridStaticPtr<FaceType, ScalarType > MetroMeshFaceGrid;
227  typedef GridStaticPtr<VertexType, ScalarType > MetroMeshVertexGrid;
228 
229 public:
230 
231  HausdorffSampler(MeshType* _m, MeshType* _sampleMesh=0, MeshType* _closestMesh=0 ) :markerFunctor(_m)
232  {
233  m=_m;
234  init(_sampleMesh,_closestMesh);
235  }
236 
237  MeshType *m;
238  MeshType *samplePtMesh;
239  MeshType *closestPtMesh;
240 
241  MetroMeshVertexGrid unifGridVert;
242  MetroMeshFaceGrid unifGridFace;
243 
244  // Parameters
245  double min_dist;
246  double max_dist;
247  double mean_dist;
248  double RMS_dist;
249  double volume;
250  double area_S1;
251  Histogramf hist;
252  // globals parameters driving the samples.
253  int n_total_samples;
254  int n_samples;
255  bool useVertexSampling;
256  ScalarType dist_upper_bound; // samples that have a distance beyond this threshold distance are not considered.
257  typedef typename tri::FaceTmark<MeshType> MarkerFace;
258  MarkerFace markerFunctor;
259 
260 
261  float getMeanDist() const { return mean_dist / n_total_samples; }
262  float getMinDist() const { return min_dist ; }
263  float getMaxDist() const { return max_dist ; }
264  float getRMSDist() const { return sqrt(RMS_dist / n_total_samples); }
265 
266  void init(MeshType* _sampleMesh=0, MeshType* _closestMesh=0 )
267  {
268  samplePtMesh =_sampleMesh;
269  closestPtMesh = _closestMesh;
270  if(m)
271  {
273  if(m->fn==0) useVertexSampling = true;
274  else useVertexSampling = false;
275 
276  if(useVertexSampling) unifGridVert.Set(m->vert.begin(),m->vert.end());
277  else unifGridFace.Set(m->face.begin(),m->face.end());
278  markerFunctor.SetMesh(m);
279  hist.SetRange(0.0, m->bbox.Diag()/100.0, 100);
280  }
281  min_dist = std::numeric_limits<double>::max();
282  max_dist = 0;
283  mean_dist =0;
284  RMS_dist = 0;
285  n_total_samples = 0;
286  }
287 
288  void AddFace(const FaceType &f, CoordType interp)
289  {
290  CoordType startPt = f.cP(0)*interp[0] + f.cP(1)*interp[1] +f.cP(2)*interp[2]; // point to be sampled
291  CoordType startN = f.cV(0)->cN()*interp[0] + f.cV(1)->cN()*interp[1] +f.cV(2)->cN()*interp[2]; // Normal of the interpolated point
292  AddSample(startPt,startN); // point to be sampled);
293  }
294 
295  void AddVert(VertexType &p)
296  {
297  p.Q()=AddSample(p.cP(),p.cN());
298  }
299 
300 
301  float AddSample(const CoordType &startPt,const CoordType &startN)
302  {
303  // the results
304  CoordType closestPt;
305  ScalarType dist = dist_upper_bound;
306 
307  // compute distance between startPt and the mesh S2
308  FaceType *nearestF=0;
309  VertexType *nearestV=0;
310  vcg::face::PointDistanceBaseFunctor<ScalarType> PDistFunct;
311  dist=dist_upper_bound;
312  if(useVertexSampling)
313  nearestV = tri::GetClosestVertex<MeshType,MetroMeshVertexGrid>(*m,unifGridVert,startPt,dist_upper_bound,dist);
314  else
315  nearestF = unifGridFace.GetClosest(PDistFunct,markerFunctor,startPt,dist_upper_bound,dist,closestPt);
316 
317  // update distance measures
318  if(dist == dist_upper_bound)
319  return dist;
320 
321  if(dist > max_dist) max_dist = dist; // L_inf
322  if(dist < min_dist) min_dist = dist; // L_inf
323 
324  mean_dist += dist; // L_1
325  RMS_dist += dist*dist; // L_2
326  n_total_samples++;
327 
328  hist.Add((float)fabs(dist));
329  if(samplePtMesh)
330  {
332  samplePtMesh->vert.back().P() = startPt;
333  samplePtMesh->vert.back().Q() = dist;
334  samplePtMesh->vert.back().N() = startN;
335  }
336  if(closestPtMesh)
337  {
339  closestPtMesh->vert.back().P() = closestPt;
340  closestPtMesh->vert.back().Q() = dist;
341  closestPtMesh->vert.back().N() = startN;
342  }
343  return dist;
344  }
345 }; // end class HausdorffSampler
346 
347 
348 
349 /* This sampler is used to transfer the detail of a mesh onto another one.
350  * It keep internally the spatial indexing structure used to find the closest point
351  */
352 template <class MeshType>
354 {
355  typedef typename MeshType::FaceType FaceType;
356  typedef typename MeshType::VertexType VertexType;
357  typedef typename MeshType::CoordType CoordType;
358  typedef typename MeshType::ScalarType ScalarType;
359  typedef GridStaticPtr<FaceType, ScalarType > MetroMeshGrid;
360  typedef GridStaticPtr<VertexType, ScalarType > VertexMeshGrid;
361 
362 public:
363 
364  RedetailSampler():m(0) {}
365 
366  MeshType *m;
367  CallBackPos *cb;
368  int sampleNum; // the expected number of samples. Used only for the callback
369  int sampleCnt;
370  MetroMeshGrid unifGridFace;
371  VertexMeshGrid unifGridVert;
372  bool useVertexSampling;
373 
374  // Parameters
375  typedef tri::FaceTmark<MeshType> MarkerFace;
376  MarkerFace markerFunctor;
377 
378  bool coordFlag;
379  bool colorFlag;
380  bool normalFlag;
381  bool qualityFlag;
382  bool selectionFlag;
383  bool storeDistanceAsQualityFlag;
384  float dist_upper_bound;
385  void init(MeshType *_m, CallBackPos *_cb=0, int targetSz=0)
386  {
387  coordFlag=false;
388  colorFlag=false;
389  qualityFlag=false;
390  selectionFlag=false;
391  storeDistanceAsQualityFlag=false;
392  m=_m;
394  if(m->fn==0) useVertexSampling = true;
395  else useVertexSampling = false;
396 
397  if(useVertexSampling) unifGridVert.Set(m->vert.begin(),m->vert.end());
398  else unifGridFace.Set(m->face.begin(),m->face.end());
399  markerFunctor.SetMesh(m);
400  // sampleNum and sampleCnt are used only for the progress callback.
401  cb=_cb;
402  sampleNum = targetSz;
403  sampleCnt = 0;
404  }
405 
406  // this function is called for each vertex of the target mesh.
407  // and retrieve the closest point on the source mesh.
408  void AddVert(VertexType &p)
409  {
410  assert(m);
411  // the results
412  CoordType closestPt, normf, bestq, ip;
413  ScalarType dist = dist_upper_bound;
414  const CoordType &startPt= p.cP();
415  // compute distance between startPt and the mesh S2
416  if(useVertexSampling)
417  {
418  VertexType *nearestV=0;
419  nearestV = tri::GetClosestVertex<MeshType,VertexMeshGrid>(*m,unifGridVert,startPt,dist_upper_bound,dist); //(PDistFunct,markerFunctor,startPt,dist_upper_bound,dist,closestPt);
420  if(cb) cb(sampleCnt++*100/sampleNum,"Resampling Vertex attributes");
421  if(storeDistanceAsQualityFlag) p.Q() = dist;
422  if(dist == dist_upper_bound) return ;
423 
424  if(coordFlag) p.P()=nearestV->P();
425  if(colorFlag) p.C() = nearestV->C();
426  if(normalFlag) p.N() = nearestV->N();
427  if(qualityFlag) p.Q()= nearestV->Q();
428  if(selectionFlag) if(nearestV->IsS()) p.SetS();
429  }
430  else
431  {
432  FaceType *nearestF=0;
433  vcg::face::PointDistanceBaseFunctor<ScalarType> PDistFunct;
434  dist=dist_upper_bound;
435  if(cb) cb(sampleCnt++*100/sampleNum,"Resampling Vertex attributes");
436  nearestF = unifGridFace.GetClosest(PDistFunct,markerFunctor,startPt,dist_upper_bound,dist,closestPt);
437  if(dist == dist_upper_bound) return ;
438 
439  CoordType interp;
440  InterpolationParameters(*nearestF,(*nearestF).cN(),closestPt, interp);
441  interp[2]=1.0-interp[1]-interp[0];
442 
443  if(coordFlag) p.P()=closestPt;
444  if(colorFlag) p.C().lerp(nearestF->V(0)->C(),nearestF->V(1)->C(),nearestF->V(2)->C(),interp);
445  if(normalFlag) p.N() = nearestF->V(0)->N()*interp[0] + nearestF->V(1)->N()*interp[1] + nearestF->V(2)->N()*interp[2];
446  if(qualityFlag) p.Q()= nearestF->V(0)->Q()*interp[0] + nearestF->V(1)->Q()*interp[1] + nearestF->V(2)->Q()*interp[2];
447  if(selectionFlag) if(nearestF->IsS()) p.SetS();
448  }
449  }
450 }; // end class RedetailSampler
451 
452 
453 
454 
465 template <class MeshType, class VertexSampler = TrivialSampler< MeshType> >
467 {
468  typedef typename MeshType::CoordType CoordType;
469  typedef typename MeshType::BoxType BoxType;
470  typedef typename MeshType::ScalarType ScalarType;
471  typedef typename MeshType::VertexType VertexType;
472  typedef typename MeshType::VertexPointer VertexPointer;
473  typedef typename MeshType::VertexIterator VertexIterator;
474  typedef typename MeshType::EdgeType EdgeType;
475  typedef typename MeshType::EdgeIterator EdgeIterator;
476  typedef typename MeshType::FaceType FaceType;
477  typedef typename MeshType::FacePointer FacePointer;
478  typedef typename MeshType::FaceIterator FaceIterator;
479  typedef typename MeshType::FaceContainer FaceContainer;
480 
481  typedef typename vcg::SpatialHashTable<FaceType, ScalarType> MeshSHT;
482  typedef typename vcg::SpatialHashTable<FaceType, ScalarType>::CellIterator MeshSHTIterator;
483  typedef typename vcg::SpatialHashTable<VertexType, ScalarType> MontecarloSHT;
484  typedef typename vcg::SpatialHashTable<VertexType, ScalarType>::CellIterator MontecarloSHTIterator;
485  typedef typename vcg::SpatialHashTable<VertexType, ScalarType> SampleSHT;
486  typedef typename vcg::SpatialHashTable<VertexType, ScalarType>::CellIterator SampleSHTIterator;
487 
488  typedef typename MeshType::template PerVertexAttributeHandle<float> PerVertexFloatAttribute;
489 
490 public:
491 
492 static math::MarsenneTwisterRNG &SamplingRandomGenerator()
493 {
494  static math::MarsenneTwisterRNG rnd;
495  return rnd;
496 }
497 
498 // Returns an integer random number in the [0,i-1] interval using the improve Marsenne-Twister method.
499 // this functor is needed for passing it to the std functions.
500 static unsigned int RandomInt(unsigned int i)
501 {
502  return (SamplingRandomGenerator().generate(i));
503 }
504 
506 {
507 public:
508  typedef unsigned int result_type;
509  MarsenneTwisterURBG(result_type max){_max = max;}
510  static constexpr result_type min() {return 0;}
511  static constexpr result_type max() {return std::numeric_limits<result_type>::max();}
512  result_type operator()() {return SamplingRandomGenerator().generate(_max);}
513 private:
514  result_type _max;
515 };
516 
517 // Returns a random number in the [0,1) real interval using the improved Marsenne-Twister method.
518 static double RandomDouble01()
519 {
520  return SamplingRandomGenerator().generate01();
521 }
522 
523 #define FAK_LEN 1024
524 static double LnFac(int n) {
525  // Tabled log factorial function. gives natural logarithm of n!
526 
527  // define constants
528  static const double // coefficients in Stirling approximation
529  C0 = 0.918938533204672722, // ln(sqrt(2*pi))
530  C1 = 1./12.,
531  C3 = -1./360.;
532  // C5 = 1./1260., // use r^5 term if FAK_LEN < 50
533  // C7 = -1./1680.; // use r^7 term if FAK_LEN < 20
534  // static variables
535  static double fac_table[FAK_LEN]; // table of ln(n!):
536  static bool initialized = false; // remember if fac_table has been initialized
537 
538 
539  if (n < FAK_LEN) {
540  if (n <= 1) {
541  if (n < 0) assert(0);//("Parameter negative in LnFac function");
542  return 0;
543  }
544  if (!initialized) { // first time. Must initialize table
545  // make table of ln(n!)
546  double sum = fac_table[0] = 0.;
547  for (int i=1; i<FAK_LEN; i++) {
548  sum += log(double(i));
549  fac_table[i] = sum;
550  }
551  initialized = true;
552  }
553  return fac_table[n];
554  }
555  // not found in table. use Stirling approximation
556  double n1, r;
557  n1 = n; r = 1. / n1;
558  return (n1 + 0.5)*log(n1) - n1 + C0 + r*(C1 + r*r*C3);
559 }
560 
561 static int PoissonRatioUniforms(double L) {
562  /*
563 
564  This subfunction generates a integer with the poisson
565  distribution using the ratio-of-uniforms rejection method (PRUAt).
566  This approach is STABLE even for large L (e.g. it does not suffer from the overflow limit of the classical Knuth implementation)
567  Execution time does not depend on L, except that it matters whether
568  is within the range where ln(n!) is tabulated.
569 
570  Reference:
571 
572  E. Stadlober
573  "The ratio of uniforms approach for generating discrete random variates".
574  Journal of Computational and Applied Mathematics,
575  vol. 31, no. 1, 1990, pp. 181-189.
576 
577  Partially adapted/inspired from some subfunctions of the Agner Fog stocc library ( www.agner.org/random )
578  Same licensing scheme.
579 
580  */
581  // constants
582 
583  const double SHAT1 = 2.943035529371538573; // 8/e
584  const double SHAT2 = 0.8989161620588987408; // 3-sqrt(12/e)
585  double u; // uniform random
586  double lf; // ln(f(x))
587  double x; // real sample
588  int k; // integer sample
589 
590  double pois_a = L + 0.5; // hat center
591  int mode = (int)L; // mode
592  double pois_g = log(L);
593  double pois_f0 = mode * pois_g - LnFac(mode); // value at mode
594  double pois_h = sqrt(SHAT1 * (L+0.5)) + SHAT2; // hat width
595  double pois_bound = (int)(pois_a + 6.0 * pois_h); // safety-bound
596 
597  while(1) {
598  u = RandomDouble01();
599  if (u == 0) continue; // avoid division by 0
600  x = pois_a + pois_h * (RandomDouble01() - 0.5) / u;
601  if (x < 0 || x >= pois_bound) continue; // reject if outside valid range
602  k = (int)(x);
603  lf = k * pois_g - LnFac(k) - pois_f0;
604  if (lf >= u * (4.0 - u) - 3.0) break; // quick acceptance
605  if (u * (u - lf) > 1.0) continue; // quick rejection
606  if (2.0 * log(u) <= lf) break; // final acceptance
607  }
608  return k;
609 }
610 
611 
623 static int Poisson(double lambda)
624 {
625  if(lambda>50) return PoissonRatioUniforms(lambda);
626  double L = exp(-lambda);
627  int k =0;
628  double p = 1.0;
629  do
630  {
631  k = k+1;
632  p = p*RandomDouble01();
633  } while (p>L);
634 
635  return k -1;
636 }
637 
638 
639 static void AllVertex(MeshType & m, VertexSampler &ps)
640 {
641  AllVertex(m, ps, false);
642 }
643 
644 static void AllVertex(MeshType & m, VertexSampler &ps, bool onlySelected)
645 {
646  VertexIterator vi;
647  for(vi=m.vert.begin();vi!=m.vert.end();++vi)
648  if(!(*vi).IsD())
649  if ((!onlySelected) || ((*vi).IsS()))
650  {
651  ps.AddVert(*vi);
652  }
653 }
654 
662 
663 static void VertexWeighted(MeshType & m, VertexSampler &ps, int sampleNum)
664 {
665  ScalarType qSum = 0;
666  VertexIterator vi;
667  for(vi = m.vert.begin(); vi != m.vert.end(); ++vi)
668  if(!(*vi).IsD())
669  qSum += (*vi).Q();
670 
671  ScalarType samplePerUnit = sampleNum/qSum;
672  ScalarType floatSampleNum =0;
673  std::vector<VertexPointer> vertVec;
674  FillAndShuffleVertexPointerVector(m,vertVec);
675 
676  std::vector<bool> vertUsed(m.vn,false);
677 
678  int i=0; int cnt=0;
679  while(cnt < sampleNum)
680  {
681  if(vertUsed[i])
682  {
683  floatSampleNum += vertVec[i]->Q() * samplePerUnit;
684  int vertSampleNum = (int) floatSampleNum;
685  floatSampleNum -= (float) vertSampleNum;
686 
687  // for every sample p_i in T...
688  if(vertSampleNum > 1)
689  {
690  ps.AddVert(*vertVec[i]);
691  cnt++;
692  vertUsed[i]=true;
693  }
694  }
695  i = (i+1)%m.vn;
696  }
697 }
698 
701 static void VertexAreaUniform(MeshType & m, VertexSampler &ps, int sampleNum)
702 {
703  VertexIterator vi;
704  for(vi = m.vert.begin(); vi != m.vert.end(); ++vi)
705  if(!(*vi).IsD())
706  (*vi).Q() = 0;
707 
708  FaceIterator fi;
709  for(fi = m.face.begin(); fi != m.face.end(); ++fi)
710  if(!(*fi).IsD())
711  {
712  ScalarType areaThird = DoubleArea(*fi)/6.0;
713  (*fi).V(0)->Q()+=areaThird;
714  (*fi).V(1)->Q()+=areaThird;
715  (*fi).V(2)->Q()+=areaThird;
716  }
717 
718  VertexWeighted(m,ps,sampleNum);
719 }
720 
721 static void FillAndShuffleFacePointerVector(MeshType & m, std::vector<FacePointer> &faceVec)
722 {
723  for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi)
724  if(!(*fi).IsD()) faceVec.push_back(&*fi);
725 
726  assert((int)faceVec.size()==m.fn);
727 
728  //unsigned int (*p_myrandom)(unsigned int) = RandomInt;
729  //std::random_device rd;
730  //std::mt19937 g(rd());
731  MarsenneTwisterURBG g(faceVec.size());
732  std::shuffle(faceVec.begin(),faceVec.end(), g);
733 }
734 static void FillAndShuffleVertexPointerVector(MeshType & m, std::vector<VertexPointer> &vertVec)
735 {
736  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi)
737  if(!(*vi).IsD()) vertVec.push_back(&*vi);
738 
739  assert((int)vertVec.size()==m.vn);
740 
741  //unsigned int (*p_myrandom)(unsigned int) = RandomInt;
742  //std::random_device rd;
743  //std::mt19937 g(rd());
744  MarsenneTwisterURBG g(vertVec.size());
745  std::shuffle(vertVec.begin(),vertVec.end(), g);
746 }
747 
749 static void VertexUniform(MeshType & m, VertexSampler &ps, int sampleNum, bool onlySelected)
750 {
751  if (sampleNum >= m.vn) {
752  AllVertex(m, ps, onlySelected);
753  return;
754  }
755 
756  std::vector<VertexPointer> vertVec;
757  FillAndShuffleVertexPointerVector(m, vertVec);
758 
759  int added = 0;
760  for (int i = 0; ((i < m.vn) && (added < sampleNum)); ++i)
761  if (!(*vertVec[i]).IsD())
762  if ((!onlySelected) || (*vertVec[i]).IsS())
763  {
764  ps.AddVert(*vertVec[i]);
765  added++;
766  }
767 
768 }
769 
770 
771 static void VertexUniform(MeshType & m, VertexSampler &ps, int sampleNum)
772 {
773  VertexUniform(m, ps, sampleNum, false);
774 }
775 
776 
784 {
785  Floor = 0,
786  Round,
787  Ceil,
788 };
789 
797 static void EdgeMeshUniform(MeshType &m, VertexSampler &ps, float radius, EdgeSamplingStrategy strategy = Floor)
798 {
799  tri::RequireEEAdjacency(m);
800  tri::RequireCompactness(m);
801  tri::RequirePerEdgeFlags(m);
802  tri::RequirePerVertexFlags(m);
805  tri::MeshAssert<MeshType>::EEOneManifold(m);
806 
807  for (EdgeIterator ei = m.edge.begin(); ei != m.edge.end(); ++ei)
808  {
809  if (!ei->IsV())
810  {
811  edge::Pos<EdgeType> ep(&*ei,0);
812  edge::Pos<EdgeType> startep = ep;
813  do // first loop to search a boundary component.
814  {
815  ep.NextE();
816  if (ep.IsBorder())
817  break;
818  } while (startep != ep);
819  if (!ep.IsBorder())
820  {
821  // it's a loop
822  assert(ep == startep);
823 
824  // to keep the uniform resampling order-independent:
825  // 1) start from the 'lowest' point...
826  edge::Pos<EdgeType> altEp = ep;
827  altEp.NextE();
828  while (altEp != startep) {
829  if (altEp.V()->cP() < ep.V()->cP())
830  {
831  ep = altEp;
832  }
833  altEp.NextE();
834  }
835 
836  // 2) ... with consistent direction
837  const auto dir0 = ep.VFlip()->cP() - ep.V()->cP();
838  ep.FlipE();
839  const auto dir1 = ep.VFlip()->cP() - ep.V()->cP();
840  if (dir0 < dir1)
841  {
842  ep.FlipE();
843  }
844  }
845  else
846  {
847  // not a loop
848 
849  // to keep the uniform resampling order-independent
850  // start from the border with 'lowest' point
851  edge::Pos<EdgeType> altEp = ep;
852  do {
853  altEp.NextE();
854  } while (!altEp.IsBorder());
855 
856  if (altEp.V()->cP() < ep.V()->cP())
857  {
858  ep = altEp;
859  }
860  }
861 
862  ScalarType totalLen = 0;
863  ep.FlipV();
864  // second loop to compute the length of the chain.
865  do
866  {
867  ep.E()->SetV();
868  totalLen += Distance(ep.V()->cP(), ep.VFlip()->cP());
869  ep.NextE();
870  } while (!ep.E()->IsV() && !ep.IsBorder());
871  if (ep.IsBorder())
872  {
873  ep.E()->SetV();
874  totalLen += Distance(ep.V()->cP(), ep.VFlip()->cP());
875  }
876 
877  VertexPointer startVertex = ep.V();
878 
879  // Third loop actually performs the sampling.
880  int sampleNum = -1;
881  {
882  double div = double(totalLen) / radius;
883  switch (strategy) {
884  case Round:
885  sampleNum = int(round(div));
886  break;
887  case Ceil:
888  sampleNum = int(ceil(div));
889  break;
890  default: // Floor
891  sampleNum = int(floor(div));
892  break;
893  };
894  }
895  assert(sampleNum >= 0);
896 
897  ScalarType sampleDist = totalLen / sampleNum;
898 // printf("Found a chain of %f with %i samples every %f (%f)\n", totalLen, sampleNum, sampleDist, radius);
899 
900  ScalarType curLen = 0;
901  int sampleCnt = 1;
902  ps.AddEdge(*(ep.E()), ep.VInd() == 0 ? 0.0 : 1.0);
903 
904  do {
905  ep.NextE();
906  assert(ep.E()->IsV());
907  ScalarType edgeLen = Distance(ep.VFlip()->cP(), ep.V()->cP());
908  ScalarType d0 = curLen;
909  ScalarType d1 = d0 + edgeLen;
910 
911  while (d1 > sampleCnt * sampleDist && sampleCnt < sampleNum)
912  {
913  ScalarType off = (sampleCnt * sampleDist - d0) / edgeLen;
914 // printf("edgeLen %f off %f samplecnt %i\n", edgeLen, off, sampleCnt);
915  ps.AddEdge(*(ep.E()), ep.VInd() == 0 ? 1.0 - off : off);
916  sampleCnt++;
917  }
918  curLen += edgeLen;
919  } while(!ep.IsBorder() && ep.V() != startVertex);
920 
921  if(ep.V() != startVertex)
922  {
923  ps.AddEdge(*(ep.E()), ep.VInd() == 0 ? 0.0 : 1.0);
924  }
925  }
926  }
927 }
928 
929 
935 static void VertexBorderCorner(MeshType & m, VertexSampler &ps, ScalarType angleRad)
936 {
938  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi)
939  {
940  if(vi->IsS()) ps.AddVert(*vi);
941  }
942 }
943 
949 static void VertexBorder(MeshType & m, VertexSampler &ps)
950 {
951  VertexBorderCorner(m,ps,std::numeric_limits<ScalarType>::max());
952 }
953 
960 static void VertexCrease(MeshType & m, VertexSampler &ps)
961 {
962  typedef typename UpdateTopology<MeshType>::PEdge SimpleEdge;
963  std::vector< SimpleEdge > Edges;
964  typename std::vector< SimpleEdge >::iterator ei;
966 
967  typename MeshType::template PerVertexAttributeHandle <int> hv = tri::Allocator<MeshType>:: template GetPerVertexAttribute<int> (m);
968 
969  for(ei=Edges.begin(); ei!=Edges.end(); ++ei)
970  {
971  hv[ei->v[0]]++;
972  hv[ei->v[1]]++;
973  }
974 
975  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi)
976  {
977  if(hv[vi]>2)
978  ps.AddVert(*vi);
979  }
980 }
981 
982 
983 static void FaceUniform(MeshType & m, VertexSampler &ps, int sampleNum)
984 {
985  if(sampleNum>=m.fn) {
986  AllFace(m,ps);
987  return;
988  }
989 
990  std::vector<FacePointer> faceVec;
991  FillAndShuffleFacePointerVector(m,faceVec);
992 
993  for(int i =0; i< sampleNum; ++i)
994  ps.AddFace(*faceVec[i],Barycenter(*faceVec[i]));
995 }
996 
997 static void AllFace(MeshType & m, VertexSampler &ps)
998 {
999  FaceIterator fi;
1000  for(fi=m.face.begin();fi!=m.face.end();++fi)
1001  if(!(*fi).IsD())
1002  {
1003  ps.AddFace(*fi,Barycenter(*fi));
1004  }
1005 }
1006 
1007 
1008 static void AllEdge(MeshType & m, VertexSampler &ps)
1009 {
1010  // Edge sampling.
1011  typedef typename UpdateTopology<MeshType>::PEdge SimpleEdge;
1012  std::vector< SimpleEdge > Edges;
1013  typename std::vector< SimpleEdge >::iterator ei;
1014  UpdateTopology<MeshType>::FillUniqueEdgeVector(m,Edges);
1015 
1016  for(ei=Edges.begin(); ei!=Edges.end(); ++ei)
1017  ps.AddFace(*(*ei).f,ei->EdgeBarycentricToFaceBarycentric(0.5));
1018 }
1019 
1020 // Regular Uniform Edge sampling
1021 // Each edge is subdivided in a number of pieces proprtional to its length
1022 // Samples are chosen without touching the vertices.
1023 
1024 static void EdgeUniform(MeshType & m, VertexSampler &ps,int sampleNum, bool sampleFauxEdge=true)
1025 {
1026  typedef typename UpdateTopology<MeshType>::PEdge SimpleEdge;
1027 
1028  std::vector< SimpleEdge > Edges;
1029  UpdateTopology<MeshType>::FillUniqueEdgeVector(m,Edges,sampleFauxEdge);
1030  // First loop compute total edge length;
1031  float edgeSum=0;
1032  typename std::vector< SimpleEdge >::iterator ei;
1033  for(ei=Edges.begin(); ei!=Edges.end(); ++ei)
1034  edgeSum+=Distance((*ei).v[0]->P(),(*ei).v[1]->P());
1035 
1036  float sampleLen = edgeSum/sampleNum;
1037  float rest=0;
1038  for(ei=Edges.begin(); ei!=Edges.end(); ++ei)
1039  {
1040  float len = Distance((*ei).v[0]->P(),(*ei).v[1]->P());
1041  float samplePerEdge = floor((len+rest)/sampleLen);
1042  rest = (len+rest) - samplePerEdge * sampleLen;
1043  float step = 1.0/(samplePerEdge+1);
1044  for(int i=0;i<samplePerEdge;++i)
1045  {
1046  CoordType interp(0,0,0);
1047  interp[ (*ei).z ]=step*(i+1);
1048  interp[((*ei).z+1)%3]=1.0-step*(i+1);
1049  ps.AddFace(*(*ei).f,interp);
1050  }
1051  }
1052 }
1053 
1054 // Generate the barycentric coords of a random point over a single face,
1055 // with a uniform distribution over the triangle.
1056 // It uses the parallelogram folding trick.
1057 static CoordType RandomBarycentric()
1058 {
1059  return math::GenerateBarycentricUniform<ScalarType>(SamplingRandomGenerator());
1060 }
1061 
1062 // Given a triangle return a random point over it
1063 static CoordType RandomPointInTriangle(const FaceType &f)
1064 {
1065  CoordType u = RandomBarycentric();
1066  return f.cP(0)*u[0] + f.cP(1)*u[1] + f.cP(2)*u[2];
1067 }
1068 
1069 static void StratifiedMontecarlo(MeshType & m, VertexSampler &ps,int sampleNum)
1070 {
1071  ScalarType area = Stat<MeshType>::ComputeMeshArea(m);
1072  ScalarType samplePerAreaUnit = sampleNum/area;
1073  // Montecarlo sampling.
1074  double floatSampleNum = 0.0;
1075 
1076  FaceIterator fi;
1077  for(fi=m.face.begin(); fi != m.face.end(); fi++)
1078  if(!(*fi).IsD())
1079  {
1080  // compute # samples in the current face (taking into account of the remainders)
1081  floatSampleNum += 0.5*DoubleArea(*fi) * samplePerAreaUnit;
1082  int faceSampleNum = (int) floatSampleNum;
1083 
1084  // for every sample p_i in T...
1085  for(int i=0; i < faceSampleNum; i++)
1086  ps.AddFace(*fi,RandomBarycentric());
1087  floatSampleNum -= (double) faceSampleNum;
1088  }
1089 }
1090 
1105 static void MontecarloPoisson(MeshType & m, VertexSampler &ps,int sampleNum)
1106 {
1107  ScalarType area = Stat<MeshType>::ComputeMeshArea(m);
1108  ScalarType samplePerAreaUnit = sampleNum/area;
1109 
1110  FaceIterator fi;
1111  for(fi=m.face.begin(); fi != m.face.end(); fi++)
1112  if(!(*fi).IsD())
1113  {
1114  float areaT=DoubleArea(*fi) * 0.5f;
1115  int faceSampleNum = Poisson(areaT*samplePerAreaUnit);
1116 
1117  // for every sample p_i in T...
1118  for(int i=0; i < faceSampleNum; i++)
1119  ps.AddFace(*fi,RandomBarycentric());
1120 // SampleNum -= (double) faceSampleNum;
1121  }
1122 }
1123 
1124 
1131 static void EdgeMontecarlo(MeshType & m, VertexSampler &ps, int sampleNum, bool sampleAllEdges)
1132 {
1133  typedef typename UpdateTopology<MeshType>::PEdge SimpleEdge;
1134  std::vector< SimpleEdge > Edges;
1135  UpdateTopology<MeshType>::FillUniqueEdgeVector(m,Edges,sampleAllEdges);
1136 
1137  assert(!Edges.empty());
1138 
1139  typedef std::pair<ScalarType, SimpleEdge*> IntervalType;
1140  std::vector< IntervalType > intervals (Edges.size()+1);
1141  int i=0;
1142  intervals[i]=std::make_pair(0,(SimpleEdge*)(0));
1143  // First loop: build a sequence of consecutive segments proportional to the edge lenghts.
1144  typename std::vector< SimpleEdge >::iterator ei;
1145  for(ei=Edges.begin(); ei != Edges.end(); ei++)
1146  {
1147  intervals[i+1]=std::make_pair(intervals[i].first+Distance((*ei).v[0]->P(),(*ei).v[1]->P()), &*ei);
1148  ++i;
1149  }
1150 
1151  // Second Loop get a point on the line 0...Sum(edgeLen) to pick a point;
1152  ScalarType edgeSum = intervals.back().first;
1153  for(i=0;i<sampleNum;++i)
1154  {
1155  ScalarType val = edgeSum * RandomDouble01();
1156  // lower_bound returns the furthermost iterator i in [first, last) such that, for every iterator j in [first, i), *j < value.
1157  // E.g. An iterator pointing to the first element "not less than" val, or end() if every element is less than val.
1158  typename std::vector<IntervalType>::iterator it = lower_bound(intervals.begin(),intervals.end(),std::make_pair(val,(SimpleEdge*)(0)) );
1159  assert(it != intervals.end() && it != intervals.begin());
1160  assert( ( (*(it-1)).first < val ) && ((*(it)).first >= val) );
1161  SimpleEdge * ep=(*it).second;
1162  ps.AddFace( *(ep->f), ep->EdgeBarycentricToFaceBarycentric(RandomDouble01()) );
1163  }
1164 }
1165 
1172 static void Montecarlo(MeshType & m, VertexSampler &ps,int sampleNum)
1173 {
1174  typedef std::pair<ScalarType, FacePointer> IntervalType;
1175  std::vector< IntervalType > intervals (m.fn+1);
1176  FaceIterator fi;
1177  int i=0;
1178  intervals[i]=std::make_pair(0,FacePointer(0));
1179  // First loop: build a sequence of consecutive segments proportional to the triangle areas.
1180  for(fi=m.face.begin(); fi != m.face.end(); fi++)
1181  if(!(*fi).IsD())
1182  {
1183  intervals[i+1]=std::make_pair(intervals[i].first+0.5*DoubleArea(*fi), &*fi);
1184  ++i;
1185  }
1186  ScalarType meshArea = intervals.back().first;
1187  for(i=0;i<sampleNum;++i)
1188  {
1189  ScalarType val = meshArea * RandomDouble01();
1190  // lower_bound returns the furthermost iterator i in [first, last) such that, for every iterator j in [first, i), *j < value.
1191  // E.g. An iterator pointing to the first element "not less than" val, or end() if every element is less than val.
1192  typename std::vector<IntervalType>::iterator it = lower_bound(intervals.begin(),intervals.end(),std::make_pair(val,FacePointer(0)) );
1193  assert(it != intervals.end());
1194  assert(it != intervals.begin());
1195  assert( (*(it-1)).first <val );
1196  assert( (*(it)).first >= val);
1197  ps.AddFace( *(*it).second, RandomBarycentric() );
1198  }
1199 }
1200 
1201 static ScalarType WeightedArea(FaceType &f, PerVertexFloatAttribute &wH)
1202 {
1203  ScalarType averageQ = ( wH[f.V(0)] + wH[f.V(1)] + wH[f.V(2)] )/3.0;
1204  return averageQ*averageQ*DoubleArea(f)/2.0;
1205 }
1206 
1215 static void WeightedMontecarlo(MeshType & m, VertexSampler &ps,int sampleNum, float variance)
1216 {
1217  tri::RequirePerVertexQuality(m);
1218  tri::RequireCompactness(m);
1219  PerVertexFloatAttribute rH = tri::Allocator<MeshType>:: template GetPerVertexAttribute<float> (m,"radius");
1220  InitRadiusHandleFromQuality(m, rH, 1.0, variance, true);
1221 
1222  ScalarType weightedArea = 0;
1223  for(FaceIterator fi = m.face.begin(); fi != m.face.end(); ++fi)
1224  weightedArea += WeightedArea(*fi,rH);
1225 
1226  ScalarType samplePerAreaUnit = sampleNum/weightedArea;
1227  // Montecarlo sampling.
1228  double floatSampleNum = 0.0;
1229  for(FaceIterator fi=m.face.begin(); fi != m.face.end(); fi++)
1230  {
1231  // compute # samples in the current face (taking into account of the remainders)
1232  floatSampleNum += WeightedArea(*fi,rH) * samplePerAreaUnit;
1233  int faceSampleNum = (int) floatSampleNum;
1234 
1235  // for every sample p_i in T...
1236  for(int i=0; i < faceSampleNum; i++)
1237  ps.AddFace(*fi,RandomBarycentric());
1238 
1239  floatSampleNum -= (double) faceSampleNum;
1240  }
1241 }
1242 
1243 
1244 // Subdivision sampling of a single face.
1245 // return number of added samples
1246 
1247 static int SingleFaceSubdivision(int sampleNum, const CoordType & v0, const CoordType & v1, const CoordType & v2, VertexSampler &ps, FacePointer fp, bool randSample)
1248 {
1249  // recursive face subdivision.
1250  if(sampleNum == 1)
1251  {
1252  // ground case.
1253  CoordType SamplePoint;
1254  if(randSample)
1255  {
1256  CoordType rb=RandomBarycentric();
1257  SamplePoint=v0*rb[0]+v1*rb[1]+v2*rb[2];
1258  }
1259  else SamplePoint=((v0+v1+v2)*(1.0f/3.0f));
1260 
1261  ps.AddFace(*fp,SamplePoint);
1262  return 1;
1263  }
1264 
1265  int s0 = sampleNum /2;
1266  int s1 = sampleNum-s0;
1267  assert(s0>0);
1268  assert(s1>0);
1269 
1270  ScalarType w0 = ScalarType(s1)/ScalarType(sampleNum);
1271  ScalarType w1 = 1.0-w0;
1272  // compute the longest edge.
1273  ScalarType maxd01 = SquaredDistance(v0,v1);
1274  ScalarType maxd12 = SquaredDistance(v1,v2);
1275  ScalarType maxd20 = SquaredDistance(v2,v0);
1276  int res;
1277  if(maxd01 > maxd12)
1278  if(maxd01 > maxd20) res = 0;
1279  else res = 2;
1280  else
1281  if(maxd12 > maxd20) res = 1;
1282  else res = 2;
1283 
1284  int faceSampleNum=0;
1285  // break the input triangle along the midpoint of the longest edge.
1286  CoordType pp;
1287  switch(res)
1288  {
1289  case 0 : pp = v0*w0 + v1*w1;
1290  faceSampleNum+=SingleFaceSubdivision(s0,v0,pp,v2,ps,fp,randSample);
1291  faceSampleNum+=SingleFaceSubdivision(s1,pp,v1,v2,ps,fp,randSample);
1292  break;
1293  case 1 : pp = v1*w0 + v2*w1;
1294  faceSampleNum+=SingleFaceSubdivision(s0,v0,v1,pp,ps,fp,randSample);
1295  faceSampleNum+=SingleFaceSubdivision(s1,v0,pp,v2,ps,fp,randSample);
1296  break;
1297  case 2 : pp = v0*w0 + v2*w1;
1298  faceSampleNum+=SingleFaceSubdivision(s0,v0,v1,pp,ps,fp,randSample);
1299  faceSampleNum+=SingleFaceSubdivision(s1,pp,v1,v2,ps,fp,randSample);
1300  break;
1301  }
1302  return faceSampleNum;
1303 }
1304 
1305 
1307 static void FaceSubdivision(MeshType & m, VertexSampler &ps,int sampleNum, bool randSample)
1308 {
1309 
1310  ScalarType area = Stat<MeshType>::ComputeMeshArea(m);
1311  ScalarType samplePerAreaUnit = sampleNum/area;
1312  std::vector<FacePointer> faceVec;
1313  FillAndShuffleFacePointerVector(m,faceVec);
1315  double floatSampleNum = 0.0;
1316  int faceSampleNum;
1317  // Subdivision sampling.
1318  typename std::vector<FacePointer>::iterator fi;
1319  for(fi=faceVec.begin(); fi!=faceVec.end(); fi++)
1320  {
1321  const CoordType b0(1.0, 0.0, 0.0);
1322  const CoordType b1(0.0, 1.0, 0.0);
1323  const CoordType b2(0.0, 0.0, 1.0);
1324  // compute # samples in the current face.
1325  floatSampleNum += 0.5*DoubleArea(**fi) * samplePerAreaUnit;
1326  faceSampleNum = (int) floatSampleNum;
1327  if(faceSampleNum>0)
1328  faceSampleNum = SingleFaceSubdivision(faceSampleNum,b0,b1,b2,ps,*fi,randSample);
1329  floatSampleNum -= (double) faceSampleNum;
1330  }
1331 }
1332 //---------
1333 // Subdivision sampling of a single face.
1334 // return number of added samples
1335 
1336 static int SingleFaceSubdivisionOld(int sampleNum, const CoordType & v0, const CoordType & v1, const CoordType & v2, VertexSampler &ps, FacePointer fp, bool randSample)
1337 {
1338  // recursive face subdivision.
1339  if(sampleNum == 1)
1340  {
1341  // ground case.
1342  CoordType SamplePoint;
1343  if(randSample)
1344  {
1345  CoordType rb=RandomBarycentric();
1346  SamplePoint=v0*rb[0]+v1*rb[1]+v2*rb[2];
1347  }
1348  else SamplePoint=((v0+v1+v2)*(1.0f/3.0f));
1349 
1350  CoordType SampleBary;
1351  InterpolationParameters(*fp,SamplePoint,SampleBary);
1352  ps.AddFace(*fp,SampleBary);
1353  return 1;
1354  }
1355 
1356  int s0 = sampleNum /2;
1357  int s1 = sampleNum-s0;
1358  assert(s0>0);
1359  assert(s1>0);
1360 
1361  ScalarType w0 = ScalarType(s1)/ScalarType(sampleNum);
1362  ScalarType w1 = 1.0-w0;
1363  // compute the longest edge.
1364  ScalarType maxd01 = SquaredDistance(v0,v1);
1365  ScalarType maxd12 = SquaredDistance(v1,v2);
1366  ScalarType maxd20 = SquaredDistance(v2,v0);
1367  int res;
1368  if(maxd01 > maxd12)
1369  if(maxd01 > maxd20) res = 0;
1370  else res = 2;
1371  else
1372  if(maxd12 > maxd20) res = 1;
1373  else res = 2;
1374 
1375  int faceSampleNum=0;
1376  // break the input triangle along the midpoint of the longest edge.
1377  CoordType pp;
1378  switch(res)
1379  {
1380  case 0 : pp = v0*w0 + v1*w1;
1381  faceSampleNum+=SingleFaceSubdivision(s0,v0,pp,v2,ps,fp,randSample);
1382  faceSampleNum+=SingleFaceSubdivision(s1,pp,v1,v2,ps,fp,randSample);
1383  break;
1384  case 1 : pp = v1*w0 + v2*w1;
1385  faceSampleNum+=SingleFaceSubdivision(s0,v0,v1,pp,ps,fp,randSample);
1386  faceSampleNum+=SingleFaceSubdivision(s1,v0,pp,v2,ps,fp,randSample);
1387  break;
1388  case 2 : pp = v0*w0 + v2*w1;
1389  faceSampleNum+=SingleFaceSubdivision(s0,v0,v1,pp,ps,fp,randSample);
1390  faceSampleNum+=SingleFaceSubdivision(s1,pp,v1,v2,ps,fp,randSample);
1391  break;
1392  }
1393  return faceSampleNum;
1394 }
1395 
1396 
1398 static void FaceSubdivisionOld(MeshType & m, VertexSampler &ps,int sampleNum, bool randSample)
1399 {
1400 
1401  ScalarType area = Stat<MeshType>::ComputeMeshArea(m);
1402  ScalarType samplePerAreaUnit = sampleNum/area;
1403  std::vector<FacePointer> faceVec;
1404  FillAndShuffleFacePointerVector(m,faceVec);
1406  double floatSampleNum = 0.0;
1407  int faceSampleNum;
1408  // Subdivision sampling.
1409  typename std::vector<FacePointer>::iterator fi;
1410  for(fi=faceVec.begin(); fi!=faceVec.end(); fi++)
1411  {
1412  // compute # samples in the current face.
1413  floatSampleNum += 0.5*DoubleArea(**fi) * samplePerAreaUnit;
1414  faceSampleNum = (int) floatSampleNum;
1415  if(faceSampleNum>0)
1416  faceSampleNum = SingleFaceSubdivision(faceSampleNum,(**fi).V(0)->cP(), (**fi).V(1)->cP(), (**fi).V(2)->cP(),ps,*fi,randSample);
1417  floatSampleNum -= (double) faceSampleNum;
1418  }
1419 }
1420 
1421 
1422 //---------
1423 
1424 // Similar Triangles sampling.
1425 // Skip vertex and edges
1426 // Sample per edges includes vertexes, so here we should expect n_samples_per_edge >=4
1427 
1428 static int SingleFaceSimilar(FacePointer fp, VertexSampler &ps, int n_samples_per_edge)
1429 {
1430  int n_samples=0;
1431  int i, j;
1432  float segmentNum=n_samples_per_edge -1 ;
1433  float segmentLen = 1.0/segmentNum;
1434  // face sampling.
1435  for(i=1; i < n_samples_per_edge-1; i++)
1436  for(j=1; j < n_samples_per_edge-1-i; j++)
1437  {
1438  //AddSample( v0 + (V1*(double)i + V2*(double)j) );
1439  CoordType sampleBary(i*segmentLen,j*segmentLen, 1.0 - (i*segmentLen+j*segmentLen) ) ;
1440  n_samples++;
1441  ps.AddFace(*fp,sampleBary);
1442  }
1443  return n_samples;
1444 }
1445 static int SingleFaceSimilarDual(FacePointer fp, VertexSampler &ps, int n_samples_per_edge, bool randomFlag)
1446 {
1447  int n_samples=0;
1448  float i, j;
1449  float segmentNum=n_samples_per_edge -1 ;
1450  float segmentLen = 1.0/segmentNum;
1451  // face sampling.
1452  for(i=0; i < n_samples_per_edge-1; i++)
1453  for(j=0; j < n_samples_per_edge-1-i; j++)
1454  {
1455  //AddSample( v0 + (V1*(double)i + V2*(double)j) );
1456  CoordType V0((i+0)*segmentLen,(j+0)*segmentLen, 1.0 - ((i+0)*segmentLen+(j+0)*segmentLen) ) ;
1457  CoordType V1((i+1)*segmentLen,(j+0)*segmentLen, 1.0 - ((i+1)*segmentLen+(j+0)*segmentLen) ) ;
1458  CoordType V2((i+0)*segmentLen,(j+1)*segmentLen, 1.0 - ((i+0)*segmentLen+(j+1)*segmentLen) ) ;
1459  n_samples++;
1460  if(randomFlag) {
1461  CoordType rb=RandomBarycentric();
1462  ps.AddFace(*fp, V0*rb[0]+V1*rb[1]+V2*rb[2]);
1463  } else ps.AddFace(*fp,(V0+V1+V2)/3.0);
1464 
1465  if( j < n_samples_per_edge-i-2 )
1466  {
1467  CoordType V3((i+1)*segmentLen,(j+1)*segmentLen, 1.0 - ((i+1)*segmentLen+(j+1)*segmentLen) ) ;
1468  n_samples++;
1469  if(randomFlag) {
1470  CoordType rb=RandomBarycentric();
1471  ps.AddFace(*fp, V3*rb[0]+V1*rb[1]+V2*rb[2]);
1472  } else ps.AddFace(*fp,(V3+V1+V2)/3.0);
1473  }
1474  }
1475  return n_samples;
1476 }
1477 
1478 // Similar sampling
1479 // Each triangle is subdivided into similar triangles following a generalization of the classical 1-to-4 splitting rule of triangles.
1480 // According to the level of subdivision <k> you get 1, 4 , 9, 16 , <k^2> triangles.
1481 // Depending on the kind of the sampling strategies we can have two different approach to choosing the sample points.
1482 // 1) you have already sampled both edges and vertices
1483 // 2) you are not going to take samples on edges and vertices.
1484 //
1485 // In the first case you have to consider only internal vertices of the subdivided triangles (to avoid multiple sampling of edges and vertices).
1486 // Therefore the number of internal points is ((k-3)*(k-2))/2. where k is the number of points on an edge (vertex included)
1487 // E.g. for k=4 you get 3 segments on each edges and the original triangle is subdivided
1488 // into 9 smaller triangles and you get (1*2)/2 == 1 only a single internal point.
1489 // So if you want N samples in a triangle you have to solve k^2 -5k +6 - 2N = 0
1490 // from which you get:
1491 //
1492 // 5 + sqrt( 1 + 8N )
1493 // k = -------------------
1494 // 2
1495 //
1496 // In the second case if you are not interested to skip the sampling on edges and vertices you have to consider as sample number the number of triangles.
1497 // So if you want N samples in a triangle, the number <k> of points on an edge (vertex included) should be simply:
1498 // k = 1 + sqrt(N)
1499 // examples:
1500 // N = 4 -> k = 3
1501 // N = 9 -> k = 4
1502 
1503 
1504 
1505 //template <class MeshType>
1506 //void Sampling<MeshType>::SimilarFaceSampling()
1507 static void FaceSimilar(MeshType & m, VertexSampler &ps,int sampleNum, bool dualFlag, bool randomFlag)
1508 {
1509  ScalarType area = Stat<MeshType>::ComputeMeshArea(m);
1510  ScalarType samplePerAreaUnit = sampleNum/area;
1511 
1512  // Similar Triangles sampling.
1513  int n_samples_per_edge;
1514  double n_samples_decimal = 0.0;
1515  FaceIterator fi;
1516 
1517  for(fi=m.face.begin(); fi != m.face.end(); fi++)
1518  {
1519  // compute # samples in the current face.
1520  n_samples_decimal += 0.5*DoubleArea(*fi) * samplePerAreaUnit;
1521  int n_samples = (int) n_samples_decimal;
1522  if(n_samples>0)
1523  {
1524  // face sampling.
1525  if(dualFlag)
1526  {
1527  n_samples_per_edge = (int)((sqrt(1.0+8.0*(double)n_samples) +5.0)/2.0); // original for non dual case
1528  n_samples = SingleFaceSimilar(&*fi,ps, n_samples_per_edge);
1529  } else {
1530  n_samples_per_edge = (int)(sqrt((double)n_samples) +1.0);
1531  n_samples = SingleFaceSimilarDual(&*fi,ps, n_samples_per_edge,randomFlag);
1532  }
1533  }
1534  n_samples_decimal -= (double) n_samples;
1535  }
1536 }
1537 
1538 
1539  // Rasterization fuction
1540  // Take a triangle
1541  // T deve essere una classe funzionale che ha l'operatore ()
1542  // con due parametri x,y di tipo S esempio:
1543  // class Foo { public void operator()(int x, int y ) { ??? } };
1544 
1545 // This function does rasterization with a safety buffer area, thus accounting some points actually outside triangle area
1546 // The safety area samples are generated according to face flag BORDER which should be true for texture space border edges
1547 // Use correctSafePointsBaryCoords = true to map safety texels to closest point barycentric coords (on edge).
1548  static void SingleFaceRaster(typename MeshType::FaceType &f, VertexSampler &ps,
1549  const Point2<typename MeshType::ScalarType> & v0,
1550  const Point2<typename MeshType::ScalarType> & v1,
1551  const Point2<typename MeshType::ScalarType> & v2,
1552  bool correctSafePointsBaryCoords=true)
1553  {
1554  typedef typename MeshType::ScalarType S;
1555  // Calcolo bounding box
1556  Box2i bbox;
1557  Box2<S> bboxf;
1558  bboxf.Add(v0);
1559  bboxf.Add(v1);
1560  bboxf.Add(v2);
1561 
1562  bbox.min[0] = floor(bboxf.min[0]);
1563  bbox.min[1] = floor(bboxf.min[1]);
1564  bbox.max[0] = ceil(bboxf.max[0]);
1565  bbox.max[1] = ceil(bboxf.max[1]);
1566 
1567  // Calcolo versori degli spigoli
1568  Point2<S> d10 = v1 - v0;
1569  Point2<S> d21 = v2 - v1;
1570  Point2<S> d02 = v0 - v2;
1571 
1572  // Preparazione prodotti scalari
1573  S b0 = (bbox.min[0]-v0[0])*d10[1] - (bbox.min[1]-v0[1])*d10[0];
1574  S b1 = (bbox.min[0]-v1[0])*d21[1] - (bbox.min[1]-v1[1])*d21[0];
1575  S b2 = (bbox.min[0]-v2[0])*d02[1] - (bbox.min[1]-v2[1])*d02[0];
1576  // Preparazione degli steps
1577  S db0 = d10[1];
1578  S db1 = d21[1];
1579  S db2 = d02[1];
1580  // Preparazione segni
1581  S dn0 = -d10[0];
1582  S dn1 = -d21[0];
1583  S dn2 = -d02[0];
1584 
1585  //Calculating orientation
1586  bool flipped = !(d02 * vcg::Point2<S>(-d10[1], d10[0]) >= 0);
1587 
1588  // Calculating border edges
1589  Segment2<S> borderEdges[3];
1590  S edgeLength[3];
1591  unsigned char edgeMask = 0;
1592 
1593  if (f.IsB(0)) {
1594  borderEdges[0] = Segment2<S>(v0, v1);
1595  edgeLength[0] = borderEdges[0].Length();
1596  edgeMask |= 1;
1597  }
1598  if (f.IsB(1)) {
1599  borderEdges[1] = Segment2<S>(v1, v2);
1600  edgeLength[1] = borderEdges[1].Length();
1601  edgeMask |= 2;
1602  }
1603  if (f.IsB(2)) {
1604  borderEdges[2] = Segment2<S>(v2, v0);
1605  edgeLength[2] = borderEdges[2].Length();
1606  edgeMask |= 4;
1607  }
1608 
1609  // Rasterizzazione
1610  double de = v0[0]*v1[1]-v0[0]*v2[1]-v1[0]*v0[1]+v1[0]*v2[1]-v2[0]*v1[1]+v2[0]*v0[1];
1611 
1612  for(int x=bbox.min[0]-1;x<=bbox.max[0]+1;++x)
1613  {
1614  bool in = false;
1615  S n[3] = { b0-db0-dn0, b1-db1-dn1, b2-db2-dn2};
1616  for(int y=bbox.min[1]-1;y<=bbox.max[1]+1;++y)
1617  {
1618  if( ((n[0]>=0 && n[1]>=0 && n[2]>=0) || (n[0]<=0 && n[1]<=0 && n[2]<=0)) && (de != 0))
1619  {
1620  typename MeshType::CoordType baryCoord;
1621  baryCoord[0] = double(-y*v1[0]+v2[0]*y+v1[1]*x-v2[0]*v1[1]+v1[0]*v2[1]-x*v2[1])/de;
1622  baryCoord[1] = -double( x*v0[1]-x*v2[1]-v0[0]*y+v0[0]*v2[1]-v2[0]*v0[1]+v2[0]*y)/de;
1623  baryCoord[2] = 1-baryCoord[0]-baryCoord[1];
1624 
1625  ps.AddTextureSample(f, baryCoord, Point2i(x,y), 0);
1626  in = true;
1627  } else {
1628  // Check whether a pixel outside (on a border edge side) triangle affects color inside it
1629  Point2<S> px(x, y);
1630  Point2<S> closePoint;
1631  int closeEdge = -1;
1632  S minDst = FLT_MAX;
1633 
1634  // find the closest point (on some edge) that lies on the 2x2 squared neighborhood of the considered point
1635  for (int i=0; i<3; ++i)
1636  {
1637  if (edgeMask & (1 << i))
1638  {
1639  Point2<S> close;
1640  S dst;
1641  if ( ((!flipped) && (n[i]<0)) ||
1642  ( flipped && (n[i]>0)) )
1643  {
1644  dst = ((close = ClosestPoint(borderEdges[i], px)) - px).Norm();
1645  if(dst < minDst &&
1646  close.X() > px.X()-1 && close.X() < px.X()+1 &&
1647  close.Y() > px.Y()-1 && close.Y() < px.Y()+1)
1648  {
1649  minDst = dst;
1650  closePoint = close;
1651  closeEdge = i;
1652  }
1653  }
1654  }
1655  }
1656 
1657  if (closeEdge >= 0)
1658  {
1659  typename MeshType::CoordType baryCoord;
1660  if (correctSafePointsBaryCoords)
1661  {
1662  // Add x,y sample with closePoint barycentric coords (on edge)
1663  baryCoord[closeEdge] = (closePoint - borderEdges[closeEdge].P1()).Norm()/edgeLength[closeEdge];
1664  baryCoord[(closeEdge+1)%3] = 1 - baryCoord[closeEdge];
1665  baryCoord[(closeEdge+2)%3] = 0;
1666  } else {
1667  // Add x,y sample with his own barycentric coords (off edge)
1668  baryCoord[0] = double(-y*v1[0]+v2[0]*y+v1[1]*x-v2[0]*v1[1]+v1[0]*v2[1]-x*v2[1])/de;
1669  baryCoord[1] = -double( x*v0[1]-x*v2[1]-v0[0]*y+v0[0]*v2[1]-v2[0]*v0[1]+v2[0]*y)/de;
1670  baryCoord[2] = 1-baryCoord[0]-baryCoord[1];
1671  }
1672  ps.AddTextureSample(f, baryCoord, Point2i(x,y), minDst);
1673  in = true;
1674  }
1675  }
1676  n[0] += dn0;
1677  n[1] += dn1;
1678  n[2] += dn2;
1679  }
1680  b0 += db0;
1681  b1 += db1;
1682  b2 += db2;
1683  }
1684 }
1685 
1686 // check the radius constrain
1687 static bool checkPoissonDisk(SampleSHT & sht, const Point3<ScalarType> & p, ScalarType radius)
1688 {
1689  // get the samples closest to the given one
1690  std::vector<VertexType*> closests;
1691  typedef EmptyTMark<MeshType> MarkerVert;
1692  static MarkerVert mv;
1693 
1694  Box3f bb(p-Point3f(radius,radius,radius),p+Point3f(radius,radius,radius));
1695  GridGetInBox(sht, mv, bb, closests);
1696 
1697  ScalarType r2 = radius*radius;
1698  for(int i=0; i<closests.size(); ++i)
1699  if(SquaredDistance(p,closests[i]->cP()) < r2)
1700  return false;
1701 
1702  return true;
1703 }
1704 
1706 {
1708  {
1709  adaptiveRadiusFlag = false;
1710  bestSampleChoiceFlag = true;
1711  bestSamplePoolSize = 10;
1712  radiusVariance =1;
1713  MAXLEVELS = 5;
1714  invertQuality = false;
1715  preGenFlag = false;
1716  preGenMesh = NULL;
1717  geodesicDistanceFlag = false;
1718  randomSeed = 0;
1719  }
1720 
1721  struct Stat
1722  {
1723  int montecarloTime;
1724  int gridTime;
1725  int pruneTime;
1726  int totalTime;
1727  Point3i gridSize;
1728  int gridCellNum;
1729  size_t sampleNum;
1730  int montecarloSampleNum;
1731  };
1732 
1733  bool geodesicDistanceFlag;
1734  bool bestSampleChoiceFlag; // In poisson disk pruning when we choose a sample in a cell, we choose the sample that remove the minimal number of other samples. This previlege the "on boundary" samples.
1735  int bestSamplePoolSize;
1736  bool adaptiveRadiusFlag;
1737  float radiusVariance;
1738  bool invertQuality;
1739  bool preGenFlag; // when generating a poisson distribution, you can initialize the set of computed points with
1740  // ALL the vertices of another mesh. Useful for building progressive//prioritize refinements.
1741  MeshType *preGenMesh; // There are two ways of passing the pregen vertexes to the pruning, 1) is with a mesh pointer
1742  // 2) with a per vertex attribute.
1743  int MAXLEVELS;
1744  int randomSeed;
1745 
1746  Stat pds;
1747 };
1748 
1749 
1750 // generate Poisson-disk sample using a set of pre-generated samples (with the Montecarlo algorithm)
1751 // It always return a point.
1752 static VertexPointer getSampleFromCell(Point3i &cell, MontecarloSHT & samplepool)
1753 {
1754  MontecarloSHTIterator cellBegin, cellEnd;
1755  samplepool.Grid(cell, cellBegin, cellEnd);
1756  return *cellBegin;
1757 }
1758 
1759 // Given a cell of the grid it search the point that remove the minimum number of other samples
1760 // it linearly scan all the points of a cell.
1761 
1762 static VertexPointer getBestPrecomputedMontecarloSample(Point3i &cell, MontecarloSHT & samplepool, ScalarType diskRadius, const PoissonDiskParam &pp)
1763 {
1764  MontecarloSHTIterator cellBegin,cellEnd;
1765  samplepool.Grid(cell, cellBegin, cellEnd);
1766  VertexPointer bestSample=0;
1767  int minRemoveCnt = std::numeric_limits<int>::max();
1768  std::vector<typename MontecarloSHT::HashIterator> inSphVec;
1769  int i=0;
1770  for(MontecarloSHTIterator ci=cellBegin; ci!=cellEnd && i<pp.bestSamplePoolSize; ++ci,i++)
1771  {
1772  VertexPointer sp = *ci;
1773  if(pp.adaptiveRadiusFlag) diskRadius = sp->Q();
1774  int curRemoveCnt = samplepool.CountInSphere(sp->cP(),diskRadius,inSphVec);
1775  if(curRemoveCnt < minRemoveCnt)
1776  {
1777  bestSample = sp;
1778  minRemoveCnt = curRemoveCnt;
1779  }
1780  }
1781  return bestSample;
1782 }
1783 
1786 static ScalarType ComputePoissonDiskRadius(MeshType &origMesh, int sampleNum)
1787 {
1788  ScalarType meshArea = Stat<MeshType>::ComputeMeshArea(origMesh);
1789  // Manage approximately the PointCloud Case, use the half a area of the bbox.
1790  // TODO: If you had the radius a much better approximation could be done.
1791  if(meshArea ==0)
1792  {
1793  meshArea = (origMesh.bbox.DimX()*origMesh.bbox.DimY() +
1794  origMesh.bbox.DimX()*origMesh.bbox.DimZ() +
1795  origMesh.bbox.DimY()*origMesh.bbox.DimZ());
1796  }
1797  ScalarType diskRadius = sqrt(meshArea / (0.7 * M_PI * sampleNum)); // 0.7 is a density factor
1798  return diskRadius;
1799 }
1800 
1801 static int ComputePoissonSampleNum(MeshType &origMesh, ScalarType diskRadius)
1802 {
1803  ScalarType meshArea = Stat<MeshType>::ComputeMeshArea(origMesh);
1804  int sampleNum = meshArea / (diskRadius*diskRadius *M_PI *0.7) ; // 0.7 is a density factor
1805  return sampleNum;
1806 }
1807 
1814 
1815 static void InitRadiusHandleFromQuality(MeshType &sampleMesh, PerVertexFloatAttribute &rH, ScalarType diskRadius, ScalarType radiusVariance, bool invert)
1816 {
1817  std::pair<float,float> minmax = tri::Stat<MeshType>::ComputePerVertexQualityMinMax( sampleMesh);
1818  float minRad = diskRadius ;
1819  float maxRad = diskRadius * radiusVariance;
1820  float deltaQ = minmax.second-minmax.first;
1821  float deltaRad = maxRad-minRad;
1822  for (VertexIterator vi = sampleMesh.vert.begin(); vi != sampleMesh.vert.end(); vi++)
1823  {
1824  rH[*vi] = minRad + deltaRad*((invert ? minmax.second - (*vi).Q() : (*vi).Q() - minmax.first )/deltaQ);
1825  }
1826 }
1827 
1828 // initialize spatial hash table for searching
1829 // radius is the radius of empty disk centered over the samples (e.g. twice of the empty space disk)
1830 // This radius implies that when we pick a sample in a cell all that cell probably will not be touched again.
1831 // Howvever we must ensure that we do not put too many vertices inside each hash cell
1832 
1833 static void InitSpatialHashTable(MeshType &montecarloMesh, MontecarloSHT &montecarloSHT, ScalarType diskRadius,
1834  struct PoissonDiskParam pp=PoissonDiskParam())
1835 {
1836  ScalarType cellsize = 2.0f* diskRadius / sqrt(3.0);
1837  float occupancyRatio=0;
1838  do
1839  {
1840  // inflating
1841  BoxType bb=montecarloMesh.bbox;
1842  assert(!bb.IsNull());
1843  bb.Offset(cellsize);
1844 
1845  int sizeX = std::max(1,int(bb.DimX() / cellsize));
1846  int sizeY = std::max(1,int(bb.DimY() / cellsize));
1847  int sizeZ = std::max(1,int(bb.DimZ() / cellsize));
1848  Point3i gridsize(sizeX, sizeY, sizeZ);
1849 
1850  montecarloSHT.InitEmpty(bb, gridsize);
1851 
1852  for (VertexIterator vi = montecarloMesh.vert.begin(); vi != montecarloMesh.vert.end(); vi++)
1853  if(!(*vi).IsD())
1854  {
1855  montecarloSHT.Add(&(*vi));
1856  }
1857 
1858  montecarloSHT.UpdateAllocatedCells();
1859  pp.pds.gridSize = gridsize;
1860  pp.pds.gridCellNum = (int)montecarloSHT.AllocatedCells.size();
1861  cellsize/=2.0f;
1862  occupancyRatio = float(montecarloMesh.vn) / float(montecarloSHT.AllocatedCells.size());
1863  // qDebug(" %i / %i = %6.3f", montecarloMesh.vn , montecarloSHT.AllocatedCells.size(),occupancyRatio);
1864  }
1865  while( occupancyRatio> 100);
1866 }
1867 
1868 static void PoissonDiskPruningByNumber(VertexSampler &ps, MeshType &m,
1869  size_t sampleNum, ScalarType &diskRadius,
1870  PoissonDiskParam &pp,
1871  float tolerance=0.04,
1872  int maxIter=20)
1873 
1874 {
1875  size_t sampleNumMin = int(float(sampleNum)*(1.0f-tolerance));
1876  size_t sampleNumMax = int(float(sampleNum)*(1.0f+tolerance));
1877  float RangeMinRad = m.bbox.Diag()/50.0;
1878  float RangeMaxRad = m.bbox.Diag()/50.0;
1879  size_t RangeMinRadNum;
1880  size_t RangeMaxRadNum;
1881  // Note RangeMinRad < RangeMaxRad
1882  // but RangeMinRadNum > sampleNum > RangeMaxRadNum
1883  do {
1884  ps.reset();
1885  RangeMinRad/=2.0f;
1886  PoissonDiskPruning(ps, m ,RangeMinRad,pp);
1887  RangeMinRadNum = pp.pds.sampleNum;
1888 // qDebug("PoissonDiskPruning Iteratin Min (%6.3f:%5i) instead of %i",RangeMinRad,RangeMinRadNum,sampleNum);
1889  } while(RangeMinRadNum < sampleNum); // if the number of sample is still smaller you have to make radius larger.
1890 
1891  do {
1892  ps.reset();
1893  RangeMaxRad*=2.0f;
1894  PoissonDiskPruning(ps, m ,RangeMaxRad,pp);
1895  RangeMaxRadNum = pp.pds.sampleNum;
1896 // qDebug("PoissonDiskPruning Iteratin Max (%6.3f:%5i) instead of %i",RangeMaxRad,RangeMaxRadNum,sampleNum);
1897  } while(RangeMaxRadNum > sampleNum);
1898 
1899 
1900  float curRadius=RangeMaxRad;
1901  int iterCnt=0;
1902  while(iterCnt<maxIter &&
1903  (pp.pds.sampleNum < sampleNumMin || pp.pds.sampleNum > sampleNumMax) )
1904  {
1905  iterCnt++;
1906  ps.reset();
1907  curRadius=(RangeMaxRad+RangeMinRad)/2.0f;
1908  PoissonDiskPruning(ps, m ,curRadius,pp);
1909 // qDebug("PoissonDiskPruning Iteratin (%6.3f:%5lu %6.3f:%5lu) Cur Radius %f -> %lu sample instead of %lu",RangeMinRad,RangeMinRadNum,RangeMaxRad,RangeMaxRadNum,curRadius,pp.pds.sampleNum,sampleNum);
1910  if(pp.pds.sampleNum > sampleNum){
1911  RangeMinRad = curRadius;
1912  RangeMinRadNum = pp.pds.sampleNum;
1913  }
1914  if(pp.pds.sampleNum < sampleNum){
1915  RangeMaxRad = curRadius;
1916  RangeMaxRadNum = pp.pds.sampleNum;
1917  }
1918  }
1919  diskRadius = curRadius;
1920 }
1921 
1922 
1932 static void PoissonDiskPruning(VertexSampler &ps, MeshType &montecarloMesh,
1933  ScalarType diskRadius, PoissonDiskParam &pp)
1934 {
1935  tri::RequireCompactness(montecarloMesh);
1936  if(pp.randomSeed) SamplingRandomGenerator().initialize(pp.randomSeed);
1937  if(pp.adaptiveRadiusFlag)
1938  tri::RequirePerVertexQuality(montecarloMesh);
1939  int t0 = clock();
1940  // spatial index of montecarlo samples - used to choose a new sample to insert
1941  MontecarloSHT montecarloSHT;
1942  InitSpatialHashTable(montecarloMesh,montecarloSHT,diskRadius,pp);
1943 
1944  // if we are doing variable density sampling we have to prepare the handle that keeps the the random samples expected radii.
1945  // At this point we just assume that there is the quality values as sampled from the base mesh
1946  PerVertexFloatAttribute rH = tri::Allocator<MeshType>:: template GetPerVertexAttribute<float> (montecarloMesh,"radius");
1947  if(pp.adaptiveRadiusFlag)
1948  InitRadiusHandleFromQuality(montecarloMesh, rH, diskRadius, pp.radiusVariance, pp.invertQuality);
1949 
1950  //unsigned int (*p_myrandom)(unsigned int) = RandomInt;
1951 // std::random_device rd;
1952 // std::mt19937 g(rd());
1953 // std::shuffle(montecarloSHT.AllocatedCells.begin(),montecarloSHT.AllocatedCells.end(), g);
1954  MarsenneTwisterURBG g(montecarloSHT.AllocatedCells.size());
1955  std::shuffle(montecarloSHT.AllocatedCells.begin(),montecarloSHT.AllocatedCells.end(), g);
1956  int t1 = clock();
1957  pp.pds.montecarloSampleNum = montecarloMesh.vn;
1958  pp.pds.sampleNum =0;
1959  int removedCnt=0;
1960  // Initial pass for pruning the Hashed grid with the an eventual pre initialized set of samples
1961  if(pp.preGenFlag)
1962  {
1963  if(pp.preGenMesh==0)
1964  {
1965  typename MeshType::template PerVertexAttributeHandle<bool> fixed;
1966  fixed = tri::Allocator<MeshType>:: template GetPerVertexAttribute<bool> (montecarloMesh,"fixed");
1967  for(VertexIterator vi=montecarloMesh.vert.begin();vi!=montecarloMesh.vert.end();++vi)
1968  if(fixed[*vi]) {
1969  pp.pds.sampleNum++;
1970  ps.AddVert(*vi);
1971  removedCnt += montecarloSHT.RemoveInSphere(vi->cP(),diskRadius);
1972  }
1973  }
1974  else
1975  {
1976  for(VertexIterator vi =pp.preGenMesh->vert.begin(); vi!=pp.preGenMesh->vert.end();++vi)
1977  {
1978  ps.AddVert(*vi);
1979  pp.pds.sampleNum++;
1980  removedCnt += montecarloSHT.RemoveInSphere(vi->cP(),diskRadius);
1981  }
1982  }
1983  montecarloSHT.UpdateAllocatedCells();
1984  }
1985  vertex::ApproximateGeodesicDistanceFunctor<VertexType> GDF;
1986  while(!montecarloSHT.AllocatedCells.empty())
1987  {
1988  removedCnt=0;
1989  for (size_t i = 0; i < montecarloSHT.AllocatedCells.size(); i++)
1990  {
1991  if( montecarloSHT.EmptyCell(montecarloSHT.AllocatedCells[i]) ) continue;
1992  ScalarType currentRadius =diskRadius;
1993  VertexPointer sp;
1994  if(pp.bestSampleChoiceFlag)
1995  sp = getBestPrecomputedMontecarloSample(montecarloSHT.AllocatedCells[i], montecarloSHT, diskRadius, pp);
1996  else
1997  sp = getSampleFromCell(montecarloSHT.AllocatedCells[i], montecarloSHT);
1998 
1999  if(pp.adaptiveRadiusFlag)
2000  currentRadius = rH[sp];
2001 
2002  ps.AddVert(*sp);
2003  pp.pds.sampleNum++;
2004  if(pp.geodesicDistanceFlag) removedCnt += montecarloSHT.RemoveInSphereNormal(sp->cP(),sp->cN(),GDF,currentRadius);
2005  else removedCnt += montecarloSHT.RemoveInSphere(sp->cP(),currentRadius);
2006  }
2007  montecarloSHT.UpdateAllocatedCells();
2008  }
2009  int t2 = clock();
2010  pp.pds.gridTime = t1-t0;
2011  pp.pds.pruneTime = t2-t1;
2012 }
2013 
2024 static void HierarchicalPoissonDisk(MeshType &origMesh, VertexSampler &ps, MeshType &montecarloMesh, ScalarType diskRadius, const struct PoissonDiskParam pp=PoissonDiskParam())
2025 {
2026 // int t0=clock();
2027  // spatial index of montecarlo samples - used to choose a new sample to insert
2028  MontecarloSHT montecarloSHTVec[5];
2029 
2030 
2031 
2032  // initialize spatial hash table for searching
2033  // radius is the radius of empty disk centered over the samples (e.g. twice of the empty space disk)
2034  // This radius implies that when we pick a sample in a cell all that cell will not be touched again.
2035  ScalarType cellsize = 2.0f* diskRadius / sqrt(3.0);
2036 
2037  // inflating
2038  BoxType bb=origMesh.bbox;
2039  bb.Offset(cellsize);
2040 
2041  int sizeX = std::max(1.0f,bb.DimX() / cellsize);
2042  int sizeY = std::max(1.0f,bb.DimY() / cellsize);
2043  int sizeZ = std::max(1.0f,bb.DimZ() / cellsize);
2044  Point3i gridsize(sizeX, sizeY, sizeZ);
2045 
2046  // spatial hash table of the generated samples - used to check the radius constrain
2047  SampleSHT checkSHT;
2048  checkSHT.InitEmpty(bb, gridsize);
2049 
2050 
2051  // sampling algorithm
2052  // ------------------
2053  //
2054  // - generate millions of samples using montecarlo algorithm
2055  // - extract a cell (C) from the active cell list (with probability proportional to cell's volume)
2056  // - generate a sample inside C by choosing one of the contained pre-generated samples
2057  // - if the sample violates the radius constrain discard it, and add the cell to the cells-to-subdivide list
2058  // - iterate until the active cell list is empty or a pre-defined number of subdivisions is reached
2059  //
2060 
2061  int level = 0;
2062 
2063  // initialize spatial hash to index pre-generated samples
2064  montecarloSHTVec[0].InitEmpty(bb, gridsize);
2065  // create active cell list
2066  for (VertexIterator vi = montecarloMesh.vert.begin(); vi != montecarloMesh.vert.end(); vi++)
2067  montecarloSHTVec[0].Add(&(*vi));
2068  montecarloSHTVec[0].UpdateAllocatedCells();
2069 
2070  // if we are doing variable density sampling we have to prepare the random samples quality with the correct expected radii.
2071  PerVertexFloatAttribute rH = tri::Allocator<MeshType>:: template GetPerVertexAttribute<float> (montecarloMesh,"radius");
2072  if(pp.adaptiveRadiusFlag)
2073  InitRadiusHandleFromQuality(montecarloMesh, rH, diskRadius, pp.radiusVariance, pp.invertQuality);
2074 
2075  do
2076  {
2077  MontecarloSHT &montecarloSHT = montecarloSHTVec[level];
2078 
2079  if(level>0)
2080  {// initialize spatial hash with the remaining points
2081  montecarloSHT.InitEmpty(bb, gridsize);
2082  // create active cell list
2083  for (typename MontecarloSHT::HashIterator hi = montecarloSHTVec[level-1].hash_table.begin(); hi != montecarloSHTVec[level-1].hash_table.end(); hi++)
2084  montecarloSHT.Add((*hi).second);
2085  montecarloSHT.UpdateAllocatedCells();
2086  }
2087  // shuffle active cells
2088  //unsigned int (*p_myrandom)(unsigned int) = RandomInt;
2089  std::random_device rd;
2090 // std::mt19937 g(rd());
2091 // std::shuffle(montecarloSHT.AllocatedCells.begin(),montecarloSHT.AllocatedCells.end(), g);
2092  MarsenneTwisterURBG g(montecarloSHT.AllocatedCells.size());
2093  std::shuffle(montecarloSHT.AllocatedCells.begin(),montecarloSHT.AllocatedCells.end(), g);
2094 
2095  // generate a sample inside C by choosing one of the contained pre-generated samples
2097  int removedCnt=montecarloSHT.hash_table.size();
2098  int addedCnt=checkSHT.hash_table.size();
2099  for (int i = 0; i < montecarloSHT.AllocatedCells.size(); i++)
2100  {
2101  for(int j=0;j<4;j++)
2102  {
2103  if( montecarloSHT.EmptyCell(montecarloSHT.AllocatedCells[i]) ) continue;
2104 
2105  // generate a sample chosen from the pre-generated one
2106  typename MontecarloSHT::HashIterator hi = montecarloSHT.hash_table.find(montecarloSHT.AllocatedCells[i]);
2107 
2108  if(hi==montecarloSHT.hash_table.end()) {break;}
2109  VertexPointer sp = (*hi).second;
2110  // vr spans between 3.0*r and r / 4.0 according to vertex quality
2111  ScalarType sampleRadius = diskRadius;
2112  if(pp.adaptiveRadiusFlag) sampleRadius = rH[sp];
2113  if (checkPoissonDisk(checkSHT, sp->cP(), sampleRadius))
2114  {
2115  ps.AddVert(*sp);
2116  montecarloSHT.RemoveCell(sp);
2117  checkSHT.Add(sp);
2118  break;
2119  }
2120  else
2121  montecarloSHT.RemovePunctual(sp);
2122  }
2123  }
2124  addedCnt = checkSHT.hash_table.size()-addedCnt;
2125  removedCnt = removedCnt-montecarloSHT.hash_table.size();
2126 
2127  // proceed to the next level of subdivision
2128  // increase grid resolution
2129  gridsize *= 2;
2130 
2131  //
2132  level++;
2133  } while(level < 5);
2134 }
2135 
2136 //template <class MeshType>
2137 //void Sampling<MeshType>::SimilarFaceSampling()
2138 
2139 // This function also generates samples outside faces if those affects faces in texture space.
2140 // Use correctSafePointsBaryCoords = true to map safety texels to closest point barycentric coords (on edge)
2141 // otherwise obtained samples will map to barycentric coord actually outside face
2142 //
2143 // If you don't need to get those extra points clear faces Border Flags
2144 // vcg::tri::UpdateFlags<Mesh>::FaceClearB(m);
2145 //
2146 // Else make sure to update border flags from texture space FFadj
2147 // vcg::tri::UpdateTopology<Mesh>::FaceFaceFromTexCoord(m);
2148 // vcg::tri::UpdateFlags<Mesh>::FaceBorderFromFF(m);
2149 static void Texture(MeshType & m, VertexSampler &ps, int textureWidth, int textureHeight, bool correctSafePointsBaryCoords=true)
2150 {
2151 typedef Point2<ScalarType> Point2x;
2152  printf("Similar Triangles face sampling\n");
2153  for(FaceIterator fi=m.face.begin(); fi != m.face.end(); fi++)
2154  if (!fi->IsD())
2155  {
2156  Point2x ti[3];
2157  for(int i=0;i<3;++i)
2158  ti[i]=Point2x((*fi).WT(i).U() * textureWidth - 0.5, (*fi).WT(i).V() * textureHeight - 0.5);
2159  // - 0.5 constants are used to obtain correct texture mapping
2160 
2161  SingleFaceRaster(*fi, ps, ti[0],ti[1],ti[2], correctSafePointsBaryCoords);
2162  }
2163 }
2164 
2165 typedef GridStaticPtr<FaceType, ScalarType > TriMeshGrid;
2166 
2167 class RRParam
2168 {
2169 public:
2170 float offset;
2171 float minDiag;
2172 tri::FaceTmark<MeshType> markerFunctor;
2173 TriMeshGrid gM;
2174 };
2175 
2176 static void RegularRecursiveOffset(MeshType & m, std::vector<CoordType> &pvec, ScalarType offset, float minDiag)
2177 {
2178  Box3<ScalarType> bb=m.bbox;
2179  bb.Offset(offset*2.0);
2180 
2181  RRParam rrp;
2182 
2183  rrp.markerFunctor.SetMesh(&m);
2184 
2185  rrp.gM.Set(m.face.begin(),m.face.end(),bb);
2186 
2187 
2188  rrp.offset=offset;
2189  rrp.minDiag=minDiag;
2190  SubdivideAndSample(m, pvec, bb, rrp, bb.Diag());
2191 }
2192 
2193 static void SubdivideAndSample(MeshType & m, std::vector<CoordType> &pvec, const Box3<ScalarType> bb, RRParam &rrp, float curDiag)
2194 {
2195  CoordType startPt = bb.Center();
2196 
2197  ScalarType dist;
2198  // Compute mesh point nearest to bb center
2199  FaceType *nearestF=0;
2200  ScalarType dist_upper_bound = curDiag+rrp.offset;
2201  CoordType closestPt;
2202  vcg::face::PointDistanceBaseFunctor<ScalarType> PDistFunct;
2203  dist=dist_upper_bound;
2204  nearestF = rrp.gM.GetClosest(PDistFunct,rrp.markerFunctor,startPt,dist_upper_bound,dist,closestPt);
2205  curDiag /=2;
2206  if(dist < dist_upper_bound)
2207  {
2208  if(curDiag/3 < rrp.minDiag) //store points only for the last level of recursion (?)
2209  {
2210  if(rrp.offset==0)
2211  pvec.push_back(closestPt);
2212  else
2213  {
2214  if(dist>rrp.offset) // points below the offset threshold cannot be displaced at the right offset distance, we can only make points nearer.
2215  {
2216  CoordType delta = startPt-closestPt;
2217  pvec.push_back(closestPt+delta*(rrp.offset/dist));
2218  }
2219  }
2220  }
2221  if(curDiag < rrp.minDiag) return;
2222  CoordType hs = (bb.max-bb.min)/2;
2223  for(int i=0;i<2;i++)
2224  for(int j=0;j<2;j++)
2225  for(int k=0;k<2;k++)
2226  SubdivideAndSample(m, pvec,
2227  BoxType(CoordType( bb.min[0]+i*hs[0], bb.min[1]+j*hs[1], bb.min[2]+k*hs[2]),
2228  CoordType(startPt[0]+i*hs[0], startPt[1]+j*hs[1], startPt[2]+k*hs[2]) ),
2229  rrp,curDiag
2230  );
2231 
2232  }
2233 }
2234 }; // end sampling class
2235 
2236 template <class MeshType>
2237 typename MeshType::ScalarType ComputePoissonDiskRadius(MeshType &origMesh, int sampleNum)
2238 {
2239  typedef typename MeshType::ScalarType ScalarType;
2240  ScalarType meshArea = Stat<MeshType>::ComputeMeshArea(origMesh);
2241  // Manage approximately the PointCloud Case, use the half a area of the bbox.
2242  // TODO: If you had the radius a much better approximation could be done.
2243  if(meshArea ==0)
2244  {
2245  meshArea = (origMesh.bbox.DimX()*origMesh.bbox.DimY() +
2246  origMesh.bbox.DimX()*origMesh.bbox.DimZ() +
2247  origMesh.bbox.DimY()*origMesh.bbox.DimZ());
2248  }
2249  ScalarType diskRadius = sqrt(meshArea / (0.7 * M_PI * sampleNum)); // 0.7 is a density factor
2250  return diskRadius;
2251 }
2252 
2253 
2254 
2255 template <class MeshType>
2256 void MontecarloSampling(MeshType &m, // the mesh that has to be sampled
2257  MeshType &mm, // the mesh that will contain the samples
2258  int sampleNum) // the desired number sample, if zero you must set the radius to the wanted value
2259 {
2260  typedef tri::MeshSampler<MeshType> BaseSampler;
2261  MeshSampler<MeshType> mcSampler(&mm);
2263 }
2264 
2265 
2266 template <class MeshType>
2267 void MontecarloSampling(MeshType &m, // the mesh that has to be sampled
2268  std::vector<Point3f> &montercarloSamples, // the vector that will contain the set of points
2269  int sampleNum) // the desired number sample, if zero you must set the radius to the wanted value
2270 {
2271  typedef tri::TrivialSampler<MeshType> BaseSampler;
2272  BaseSampler mcSampler(montercarloSamples);
2274 }
2275 
2276 // Yet another simpler wrapper for the generation of a poisson disk distribution over a mesh.
2277 //
2278 template <class MeshType>
2279 void PoissonSampling(MeshType &m, // the mesh that has to be sampled
2280  std::vector<typename MeshType::CoordType> &poissonSamples, // the vector that will contain the set of points
2281  int sampleNum, // the desired number sample, if zero you must set the radius to the wanted value
2282  typename MeshType::ScalarType &radius, // the Poisson Disk Radius (used if sampleNum==0, setted if sampleNum!=0)
2283  typename MeshType::ScalarType radiusVariance=1,
2284  typename MeshType::ScalarType PruningByNumberTolerance=0.04f,
2285  unsigned int randSeed=0)
2286 
2287 {
2288  typedef tri::TrivialSampler<MeshType> BaseSampler;
2289  typedef tri::MeshSampler<MeshType> MontecarloSampler;
2290 
2291  typename tri::SurfaceSampling<MeshType, BaseSampler>::PoissonDiskParam pp;
2292  int t0=clock();
2293 
2294 // if(sampleNum>0) radius = tri::SurfaceSampling<MeshType,BaseSampler>::ComputePoissonDiskRadius(m,sampleNum);
2295  if(radius>0 && sampleNum==0) sampleNum = tri::SurfaceSampling<MeshType,BaseSampler>::ComputePoissonSampleNum(m,radius);
2296 
2297  pp.pds.sampleNum = sampleNum;
2298  pp.randomSeed = randSeed;
2299  poissonSamples.clear();
2300 // std::vector<Point3f> MontecarloSamples;
2301  MeshType MontecarloMesh;
2302 
2303  // First step build the sampling
2304  MontecarloSampler mcSampler(MontecarloMesh);
2305  BaseSampler pdSampler(poissonSamples);
2306 
2307  if(randSeed) tri::SurfaceSampling<MeshType,MontecarloSampler>::SamplingRandomGenerator().initialize(randSeed);
2308  tri::SurfaceSampling<MeshType,MontecarloSampler>::Montecarlo(m, mcSampler, std::max(10000,sampleNum*40));
2309  tri::UpdateBounding<MeshType>::Box(MontecarloMesh);
2310 // tri::Build(MontecarloMesh, MontecarloSamples);
2311  int t1=clock();
2312  pp.pds.montecarloTime = t1-t0;
2313 
2314  if(radiusVariance !=1)
2315  {
2316  pp.adaptiveRadiusFlag=true;
2317  pp.radiusVariance=radiusVariance;
2318  }
2319  if(sampleNum==0) tri::SurfaceSampling<MeshType,BaseSampler>::PoissonDiskPruning(pdSampler, MontecarloMesh, radius,pp);
2320  else tri::SurfaceSampling<MeshType,BaseSampler>::PoissonDiskPruningByNumber(pdSampler, MontecarloMesh, sampleNum, radius,pp,PruningByNumberTolerance);
2321  int t2=clock();
2322  pp.pds.totalTime = t2-t0;
2323 }
2324 
2328 //
2329 template <class MeshType>
2330 void PoissonPruning(MeshType &m, // the mesh that has to be pruned
2331  std::vector<typename MeshType::VertexPointer> &poissonSamples, // the vector that will contain the chosen set of points
2332  float radius, unsigned int randSeed=0)
2333 {
2334  typedef tri::TrivialPointerSampler<MeshType> BaseSampler;
2336  pp.randomSeed = randSeed;
2337 
2339  BaseSampler pdSampler;
2341  poissonSamples = pdSampler.sampleVec;
2342 }
2343 
2344 
2350 template <class MeshType>
2351 void PoissonPruning(MeshType &m, // the mesh that has to be pruned
2352  std::vector<typename MeshType::CoordType> &poissonSamples, // the vector that will contain the chosen set of points
2353  float radius, unsigned int randSeed=0)
2354 {
2355  std::vector<typename MeshType::VertexPointer> poissonSamplesVP;
2356  PoissonPruning(m,poissonSamplesVP,radius,randSeed);
2357  poissonSamples.resize(poissonSamplesVP.size());
2358  for(size_t i=0;i<poissonSamplesVP.size();++i)
2359  poissonSamples[i]=poissonSamplesVP[i]->P();
2360 }
2361 
2362 
2363 
2369 template <class MeshType>
2370 void PoissonPruningExact(MeshType &m,
2371  std::vector<typename MeshType::VertexPointer> &poissonSamples,
2372  typename MeshType::ScalarType & radius,
2373  int sampleNum,
2374  float tolerance=0.04,
2375  int maxIter=20,
2376  unsigned int randSeed=0)
2377 {
2378  size_t sampleNumMin = int(float(sampleNum)*(1.0f-tolerance)); // the expected values range.
2379  size_t sampleNumMax = int(float(sampleNum)*(1.0f+tolerance)); // e.g. any sampling in [sampleNumMin, sampleNumMax] is OK
2380  float RangeMinRad = m.bbox.Diag()/10.0f;
2381  float RangeMaxRad = m.bbox.Diag()/10.0f;
2382  size_t RangeMinSampleNum;
2383  size_t RangeMaxSampleNum;
2384  std::vector<typename MeshType::VertexPointer> poissonSamplesTmp;
2385 
2386  do
2387  {
2388  RangeMinRad/=2.0f;
2389  PoissonPruning(m,poissonSamplesTmp,RangeMinRad,randSeed);
2390  RangeMinSampleNum = poissonSamplesTmp.size();
2391  } while(RangeMinSampleNum < sampleNumMin);
2392 
2393  do
2394  {
2395  RangeMaxRad*=2.0f;
2396  PoissonPruning(m,poissonSamplesTmp,RangeMaxRad,randSeed);
2397  RangeMaxSampleNum = poissonSamplesTmp.size();
2398  } while(RangeMaxSampleNum > sampleNumMax);
2399 
2400  float curRadius;
2401  int iterCnt=0;
2402  while(iterCnt<maxIter &&
2403  (poissonSamplesTmp.size() < sampleNumMin || poissonSamplesTmp.size() > sampleNumMax) )
2404  {
2405  curRadius=(RangeMaxRad+RangeMinRad)/2.0f;
2406  PoissonPruning(m,poissonSamplesTmp,curRadius,randSeed);
2407  //qDebug("(%6.3f:%5i %6.3f:%5i) Cur Radius %f -> %i sample instead of %i",RangeMinRad,RangeMinSampleNum,RangeMaxRad,RangeMaxSampleNum,curRadius,poissonSamplesTmp.size(),sampleNum);
2408  if(poissonSamplesTmp.size() > size_t(sampleNum))
2409  RangeMinRad = curRadius;
2410  if(poissonSamplesTmp.size() < size_t(sampleNum))
2411  RangeMaxRad = curRadius;
2412  }
2413 
2414  swap(poissonSamples,poissonSamplesTmp);
2415  radius = curRadius;
2416 }
2417 } // end namespace tri
2418 } // end namespace vcg
2419 
2420 #endif
2421 
Definition: box3.h:42
Point3< BoxScalarType > max
max coordinate point
Definition: box3.h:51
Point3< BoxScalarType > Center() const
Return the center of the box.
Definition: box3.h:250
void Offset(const BoxScalarType s)
Definition: box3.h:75
Point3< BoxScalarType > min
min coordinate point
Definition: box3.h:49
BoxScalarType Diag() const
Return the lenght of the diagonal of the box .
Definition: box3.h:240
Class to safely add and delete elements in a mesh.
Definition: allocate.h:97
static VertexIterator AddVertices(MeshType &m, size_t n, PointerUpdater< VertexPointer > &pu)
Add n vertices to the mesh. Function to add n vertices to the mesh. The elements are added always to ...
Definition: allocate.h:189
Definition: point_sampling.h:221
double volume
from the wikipedia defintion RMS DIST is sqrt(Sum(distances^2)/n), here we store Sum(distances^2)
Definition: point_sampling.h:249
MeshType * closestPtMesh
the mesh containing the sample points
Definition: point_sampling.h:239
MeshType * samplePtMesh
the mesh for which we search the closest points.
Definition: point_sampling.h:238
MetroMeshVertexGrid unifGridVert
the mesh containing the corresponding closest points that have been found
Definition: point_sampling.h:241
Definition: point_sampling.h:177
Definition: point_sampling.h:354
CallBackPos * cb
the source mesh for which we search the closest points (e.g. the mesh from which we take colors etc).
Definition: point_sampling.h:367
Definition: point_sampling.h:506
Definition: point_sampling.h:2168
Main Class of the Sampling framework.
Definition: point_sampling.h:467
static void Montecarlo(MeshType &m, VertexSampler &ps, int sampleNum)
Definition: point_sampling.h:1172
static void HierarchicalPoissonDisk(MeshType &origMesh, VertexSampler &ps, MeshType &montecarloMesh, ScalarType diskRadius, const struct PoissonDiskParam pp=PoissonDiskParam())
Definition: point_sampling.h:2024
static void VertexAreaUniform(MeshType &m, VertexSampler &ps, int sampleNum)
Definition: point_sampling.h:701
static void VertexWeighted(MeshType &m, VertexSampler &ps, int sampleNum)
Definition: point_sampling.h:663
static void VertexCrease(MeshType &m, VertexSampler &ps)
Definition: point_sampling.h:960
static void VertexUniform(MeshType &m, VertexSampler &ps, int sampleNum, bool onlySelected)
Sample the vertices in a uniform way. Each vertex has the same probabiltiy of being chosen.
Definition: point_sampling.h:749
static void FaceSubdivision(MeshType &m, VertexSampler &ps, int sampleNum, bool randSample)
Compute a sampling of the surface where the points are regularly scattered over the face surface usin...
Definition: point_sampling.h:1307
static void PoissonDiskPruning(VertexSampler &ps, MeshType &montecarloMesh, ScalarType diskRadius, PoissonDiskParam &pp)
Definition: point_sampling.h:1932
static ScalarType ComputePoissonDiskRadius(MeshType &origMesh, int sampleNum)
Estimate the radius r that you should give to get a certain number of samples in a Poissson Disk Dist...
Definition: point_sampling.h:1786
static void VertexBorder(MeshType &m, VertexSampler &ps)
Sample all the border vertices.
Definition: point_sampling.h:949
static void VertexBorderCorner(MeshType &m, VertexSampler &ps, ScalarType angleRad)
Sample all the border corner vertices.
Definition: point_sampling.h:935
static void FaceSubdivisionOld(MeshType &m, VertexSampler &ps, int sampleNum, bool randSample)
Compute a sampling of the surface where the points are regularly scattered over the face surface usin...
Definition: point_sampling.h:1398
static int Poisson(double lambda)
Definition: point_sampling.h:623
static void MontecarloPoisson(MeshType &m, VertexSampler &ps, int sampleNum)
Definition: point_sampling.h:1105
static void WeightedMontecarlo(MeshType &m, VertexSampler &ps, int sampleNum, float variance)
Definition: point_sampling.h:1215
static void InitRadiusHandleFromQuality(MeshType &sampleMesh, PerVertexFloatAttribute &rH, ScalarType diskRadius, ScalarType radiusVariance, bool invert)
Definition: point_sampling.h:1815
static void EdgeMeshUniform(MeshType &m, VertexSampler &ps, float radius, EdgeSamplingStrategy strategy=Floor)
Definition: point_sampling.h:797
static void EdgeMontecarlo(MeshType &m, VertexSampler &ps, int sampleNum, bool sampleAllEdges)
Definition: point_sampling.h:1131
EdgeSamplingStrategy
The EdgeSamplingStrategy enum determines the sampling strategy for edge meshes. Given a sampling radi...
Definition: point_sampling.h:784
Definition: point_sampling.h:137
A basic sampler class that show the required interface used by the SurfaceSampling class.
Definition: point_sampling.h:72
static void Box(ComputeMeshType &m)
Calculates the bounding box of the given mesh m.
Definition: bounding.h:45
Management, updating and computation of per-vertex and per-face flags (like border flags).
Definition: flag.h:44
static void PerFaceNormalized(ComputeMeshType &m)
Equivalent to PerFace() and NormalizePerFace()
Definition: normal.h:276
static size_t VertexCornerBorder(MeshType &m, ScalarType angleRad, bool preserveSelection=false)
Select the border vertices that form a corner along the border with an angle that is below a certain ...
Definition: selection.h:641
Auxiliary data structure for computing face face adjacency information.
Definition: topology.h:149
Generation of per-vertex and per-face topological information.
Definition: topology.h:43
void PoissonPruningExact(MeshType &m, std::vector< typename MeshType::VertexPointer > &poissonSamples, typename MeshType::ScalarType &radius, int sampleNum, float tolerance=0.04, int maxIter=20, unsigned int randSeed=0)
Very simple wrapping for the Exact Poisson Disk Pruning.
Definition: point_sampling.h:2370
void PoissonPruning(MeshType &m, std::vector< typename MeshType::VertexPointer > &poissonSamples, float radius, unsigned int randSeed=0)
Low level wrapper for Poisson Disk Pruning.
Definition: point_sampling.h:2330
Definition: namespaces.dox:6
Definition: point_sampling.h:1722
Definition: point_sampling.h:1706