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