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 BaseSampler
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 static void EdgeMeshUniform(MeshType &m, VertexSampler &ps, float radius, bool conservative = true)
766 {
767  tri::RequireEEAdjacency(m);
768  tri::RequireCompactness(m);
769  tri::RequirePerEdgeFlags(m);
770  tri::RequirePerVertexFlags(m);
773  tri::MeshAssert<MeshType>::EEOneManifold(m);
774 
775  for (EdgeIterator ei = m.edge.begin(); ei != m.edge.end(); ++ei)
776  {
777  if (!ei->IsV())
778  {
779  edge::Pos<EdgeType> ep(&*ei,0);
780  edge::Pos<EdgeType> startep = ep;
781  do // first loop to search a boundary component.
782  {
783  ep.NextE();
784  if (ep.IsBorder())
785  break;
786  } while (startep != ep);
787  if (!ep.IsBorder())
788  {
789  // it's a loop
790  assert(ep == startep);
791 
792  // to keep the uniform resampling order-independent:
793  // 1) start from the 'lowest' point...
794  edge::Pos<EdgeType> altEp = ep;
795  altEp.NextE();
796  while (altEp != startep) {
797  if (altEp.V()->cP() < ep.V()->cP())
798  {
799  ep = altEp;
800  }
801  altEp.NextE();
802  }
803 
804  // 2) ... with consistent direction
805  const auto dir0 = ep.VFlip()->cP() - ep.V()->cP();
806  ep.FlipE();
807  const auto dir1 = ep.VFlip()->cP() - ep.V()->cP();
808  if (dir0 < dir1)
809  {
810  ep.FlipE();
811  }
812  }
813  else
814  {
815  // not a loop
816 
817  // to keep the uniform resampling order-independent
818  // start from the border with 'lowest' point
819  edge::Pos<EdgeType> altEp = ep;
820  do {
821  altEp.NextE();
822  } while (!altEp.IsBorder());
823 
824  if (altEp.V()->cP() < ep.V()->cP())
825  {
826  ep = altEp;
827  }
828  }
829 
830  ScalarType totalLen = 0;
831  ep.FlipV();
832  // second loop to compute the length of the chain.
833  do
834  {
835  ep.E()->SetV();
836  totalLen += Distance(ep.V()->cP(), ep.VFlip()->cP());
837  ep.NextE();
838  } while (!ep.E()->IsV() && !ep.IsBorder());
839  if (ep.IsBorder())
840  {
841  ep.E()->SetV();
842  totalLen += Distance(ep.V()->cP(), ep.VFlip()->cP());
843  }
844 
845  VertexPointer startVertex = ep.V();
846 
847  // Third loop actually performs the sampling.
848  int sampleNum = conservative ? floor(totalLen / radius) : ceil(totalLen / radius);
849 
850  ScalarType sampleDist = totalLen / sampleNum;
851 // printf("Found a chain of %f with %i samples every %f (%f)\n", totalLen, sampleNum, sampleDist, radius);
852 
853  ScalarType curLen = 0;
854  int sampleCnt = 1;
855  ps.AddEdge(*(ep.E()), ep.VInd() == 0 ? 0.0 : 1.0);
856 
857  do {
858  ep.NextE();
859  assert(ep.E()->IsV());
860  ScalarType edgeLen = Distance(ep.VFlip()->cP(), ep.V()->cP());
861  ScalarType d0 = curLen;
862  ScalarType d1 = d0 + edgeLen;
863 
864  while (d1 > sampleCnt * sampleDist && sampleCnt < sampleNum)
865  {
866  ScalarType off = (sampleCnt * sampleDist - d0) / edgeLen;
867 // printf("edgeLen %f off %f samplecnt %i\n", edgeLen, off, sampleCnt);
868  ps.AddEdge(*(ep.E()), ep.VInd() == 0 ? 1.0 - off : off);
869  sampleCnt++;
870  }
871  curLen += edgeLen;
872  } while(!ep.IsBorder() && ep.V() != startVertex);
873 
874  if(ep.V() != startVertex)
875  {
876  ps.AddEdge(*(ep.E()), ep.VInd() == 0 ? 0.0 : 1.0);
877  }
878  }
879  }
880 }
881 
882 
888 static void VertexBorderCorner(MeshType & m, VertexSampler &ps, ScalarType angleRad)
889 {
891  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi)
892  {
893  if(vi->IsS()) ps.AddVert(*vi);
894  }
895 }
896 
902 static void VertexBorder(MeshType & m, VertexSampler &ps)
903 {
904  VertexBorderCorner(m,ps,std::numeric_limits<ScalarType>::max());
905 }
906 
913 static void VertexCrease(MeshType & m, VertexSampler &ps)
914 {
915  typedef typename UpdateTopology<MeshType>::PEdge SimpleEdge;
916  std::vector< SimpleEdge > Edges;
917  typename std::vector< SimpleEdge >::iterator ei;
919 
920  typename MeshType::template PerVertexAttributeHandle <int> hv = tri::Allocator<MeshType>:: template GetPerVertexAttribute<int> (m);
921 
922  for(ei=Edges.begin(); ei!=Edges.end(); ++ei)
923  {
924  hv[ei->v[0]]++;
925  hv[ei->v[1]]++;
926  }
927 
928  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi)
929  {
930  if(hv[vi]>2)
931  ps.AddVert(*vi);
932  }
933 }
934 
935 
936 static void FaceUniform(MeshType & m, VertexSampler &ps, int sampleNum)
937 {
938  if(sampleNum>=m.fn) {
939  AllFace(m,ps);
940  return;
941  }
942 
943  std::vector<FacePointer> faceVec;
944  FillAndShuffleFacePointerVector(m,faceVec);
945 
946  for(int i =0; i< sampleNum; ++i)
947  ps.AddFace(*faceVec[i],Barycenter(*faceVec[i]));
948 }
949 
950 static void AllFace(MeshType & m, VertexSampler &ps)
951 {
952  FaceIterator fi;
953  for(fi=m.face.begin();fi!=m.face.end();++fi)
954  if(!(*fi).IsD())
955  {
956  ps.AddFace(*fi,Barycenter(*fi));
957  }
958 }
959 
960 
961 static void AllEdge(MeshType & m, VertexSampler &ps)
962 {
963  // Edge sampling.
964  typedef typename UpdateTopology<MeshType>::PEdge SimpleEdge;
965  std::vector< SimpleEdge > Edges;
966  typename std::vector< SimpleEdge >::iterator ei;
967  UpdateTopology<MeshType>::FillUniqueEdgeVector(m,Edges);
968 
969  for(ei=Edges.begin(); ei!=Edges.end(); ++ei)
970  ps.AddFace(*(*ei).f,ei->EdgeBarycentricToFaceBarycentric(0.5));
971 }
972 
973 // Regular Uniform Edge sampling
974 // Each edge is subdivided in a number of pieces proprtional to its length
975 // Sample are choosen without touching the vertices.
976 
977 static void EdgeUniform(MeshType & m, VertexSampler &ps,int sampleNum, bool sampleFauxEdge=true)
978 {
979  typedef typename UpdateTopology<MeshType>::PEdge SimpleEdge;
980 
981  std::vector< SimpleEdge > Edges;
982  UpdateTopology<MeshType>::FillUniqueEdgeVector(m,Edges,sampleFauxEdge);
983  // First loop compute total edge length;
984  float edgeSum=0;
985  typename std::vector< SimpleEdge >::iterator ei;
986  for(ei=Edges.begin(); ei!=Edges.end(); ++ei)
987  edgeSum+=Distance((*ei).v[0]->P(),(*ei).v[1]->P());
988 
989  float sampleLen = edgeSum/sampleNum;
990  float rest=0;
991  for(ei=Edges.begin(); ei!=Edges.end(); ++ei)
992  {
993  float len = Distance((*ei).v[0]->P(),(*ei).v[1]->P());
994  float samplePerEdge = floor((len+rest)/sampleLen);
995  rest = (len+rest) - samplePerEdge * sampleLen;
996  float step = 1.0/(samplePerEdge+1);
997  for(int i=0;i<samplePerEdge;++i)
998  {
999  CoordType interp(0,0,0);
1000  interp[ (*ei).z ]=step*(i+1);
1001  interp[((*ei).z+1)%3]=1.0-step*(i+1);
1002  ps.AddFace(*(*ei).f,interp);
1003  }
1004  }
1005 }
1006 
1007 // Generate the barycentric coords of a random point over a single face,
1008 // with a uniform distribution over the triangle.
1009 // It uses the parallelogram folding trick.
1010 static CoordType RandomBarycentric()
1011 {
1012  return math::GenerateBarycentricUniform<ScalarType>(SamplingRandomGenerator());
1013 }
1014 
1015 // Given a triangle return a random point over it
1016 static CoordType RandomPointInTriangle(const FaceType &f)
1017 {
1018  CoordType u = RandomBarycentric();
1019  return f.cP(0)*u[0] + f.cP(1)*u[1] + f.cP(2)*u[2];
1020 }
1021 
1022 static void StratifiedMontecarlo(MeshType & m, VertexSampler &ps,int sampleNum)
1023 {
1024  ScalarType area = Stat<MeshType>::ComputeMeshArea(m);
1025  ScalarType samplePerAreaUnit = sampleNum/area;
1026  // Montecarlo sampling.
1027  double floatSampleNum = 0.0;
1028 
1029  FaceIterator fi;
1030  for(fi=m.face.begin(); fi != m.face.end(); fi++)
1031  if(!(*fi).IsD())
1032  {
1033  // compute # samples in the current face (taking into account of the remainders)
1034  floatSampleNum += 0.5*DoubleArea(*fi) * samplePerAreaUnit;
1035  int faceSampleNum = (int) floatSampleNum;
1036 
1037  // for every sample p_i in T...
1038  for(int i=0; i < faceSampleNum; i++)
1039  ps.AddFace(*fi,RandomBarycentric());
1040  floatSampleNum -= (double) faceSampleNum;
1041  }
1042 }
1043 
1058 static void MontecarloPoisson(MeshType & m, VertexSampler &ps,int sampleNum)
1059 {
1060  ScalarType area = Stat<MeshType>::ComputeMeshArea(m);
1061  ScalarType samplePerAreaUnit = sampleNum/area;
1062 
1063  FaceIterator fi;
1064  for(fi=m.face.begin(); fi != m.face.end(); fi++)
1065  if(!(*fi).IsD())
1066  {
1067  float areaT=DoubleArea(*fi) * 0.5f;
1068  int faceSampleNum = Poisson(areaT*samplePerAreaUnit);
1069 
1070  // for every sample p_i in T...
1071  for(int i=0; i < faceSampleNum; i++)
1072  ps.AddFace(*fi,RandomBarycentric());
1073 // SampleNum -= (double) faceSampleNum;
1074  }
1075 }
1076 
1077 
1084 static void EdgeMontecarlo(MeshType & m, VertexSampler &ps, int sampleNum, bool sampleAllEdges)
1085 {
1086  typedef typename UpdateTopology<MeshType>::PEdge SimpleEdge;
1087  std::vector< SimpleEdge > Edges;
1088  UpdateTopology<MeshType>::FillUniqueEdgeVector(m,Edges,sampleAllEdges);
1089 
1090  assert(!Edges.empty());
1091 
1092  typedef std::pair<ScalarType, SimpleEdge*> IntervalType;
1093  std::vector< IntervalType > intervals (Edges.size()+1);
1094  int i=0;
1095  intervals[i]=std::make_pair(0,(SimpleEdge*)(0));
1096  // First loop: build a sequence of consecutive segments proportional to the edge lenghts.
1097  typename std::vector< SimpleEdge >::iterator ei;
1098  for(ei=Edges.begin(); ei != Edges.end(); ei++)
1099  {
1100  intervals[i+1]=std::make_pair(intervals[i].first+Distance((*ei).v[0]->P(),(*ei).v[1]->P()), &*ei);
1101  ++i;
1102  }
1103 
1104  // Second Loop get a point on the line 0...Sum(edgeLen) to pick a point;
1105  ScalarType edgeSum = intervals.back().first;
1106  for(i=0;i<sampleNum;++i)
1107  {
1108  ScalarType val = edgeSum * RandomDouble01();
1109  // lower_bound returns the furthermost iterator i in [first, last) such that, for every iterator j in [first, i), *j < value.
1110  // E.g. An iterator pointing to the first element "not less than" val, or end() if every element is less than val.
1111  typename std::vector<IntervalType>::iterator it = lower_bound(intervals.begin(),intervals.end(),std::make_pair(val,(SimpleEdge*)(0)) );
1112  assert(it != intervals.end() && it != intervals.begin());
1113  assert( ( (*(it-1)).first < val ) && ((*(it)).first >= val) );
1114  SimpleEdge * ep=(*it).second;
1115  ps.AddFace( *(ep->f), ep->EdgeBarycentricToFaceBarycentric(RandomDouble01()) );
1116  }
1117 }
1118 
1125 static void Montecarlo(MeshType & m, VertexSampler &ps,int sampleNum)
1126 {
1127  typedef std::pair<ScalarType, FacePointer> IntervalType;
1128  std::vector< IntervalType > intervals (m.fn+1);
1129  FaceIterator fi;
1130  int i=0;
1131  intervals[i]=std::make_pair(0,FacePointer(0));
1132  // First loop: build a sequence of consecutive segments proportional to the triangle areas.
1133  for(fi=m.face.begin(); fi != m.face.end(); fi++)
1134  if(!(*fi).IsD())
1135  {
1136  intervals[i+1]=std::make_pair(intervals[i].first+0.5*DoubleArea(*fi), &*fi);
1137  ++i;
1138  }
1139  ScalarType meshArea = intervals.back().first;
1140  for(i=0;i<sampleNum;++i)
1141  {
1142  ScalarType val = meshArea * RandomDouble01();
1143  // lower_bound returns the furthermost iterator i in [first, last) such that, for every iterator j in [first, i), *j < value.
1144  // E.g. An iterator pointing to the first element "not less than" val, or end() if every element is less than val.
1145  typename std::vector<IntervalType>::iterator it = lower_bound(intervals.begin(),intervals.end(),std::make_pair(val,FacePointer(0)) );
1146  assert(it != intervals.end());
1147  assert(it != intervals.begin());
1148  assert( (*(it-1)).first <val );
1149  assert( (*(it)).first >= val);
1150  ps.AddFace( *(*it).second, RandomBarycentric() );
1151  }
1152 }
1153 
1154 static ScalarType WeightedArea(FaceType &f, PerVertexFloatAttribute &wH)
1155 {
1156  ScalarType averageQ = ( wH[f.V(0)] + wH[f.V(1)] + wH[f.V(2)] )/3.0;
1157  return averageQ*averageQ*DoubleArea(f)/2.0;
1158 }
1159 
1168 static void WeightedMontecarlo(MeshType & m, VertexSampler &ps,int sampleNum, float variance)
1169 {
1170  tri::RequirePerVertexQuality(m);
1171  tri::RequireCompactness(m);
1172  PerVertexFloatAttribute rH = tri::Allocator<MeshType>:: template GetPerVertexAttribute<float> (m,"radius");
1173  InitRadiusHandleFromQuality(m, rH, 1.0, variance, true);
1174 
1175  ScalarType weightedArea = 0;
1176  for(FaceIterator fi = m.face.begin(); fi != m.face.end(); ++fi)
1177  weightedArea += WeightedArea(*fi,rH);
1178 
1179  ScalarType samplePerAreaUnit = sampleNum/weightedArea;
1180  // Montecarlo sampling.
1181  double floatSampleNum = 0.0;
1182  for(FaceIterator fi=m.face.begin(); fi != m.face.end(); fi++)
1183  {
1184  // compute # samples in the current face (taking into account of the remainders)
1185  floatSampleNum += WeightedArea(*fi,rH) * samplePerAreaUnit;
1186  int faceSampleNum = (int) floatSampleNum;
1187 
1188  // for every sample p_i in T...
1189  for(int i=0; i < faceSampleNum; i++)
1190  ps.AddFace(*fi,RandomBarycentric());
1191 
1192  floatSampleNum -= (double) faceSampleNum;
1193  }
1194 }
1195 
1196 
1197 // Subdivision sampling of a single face.
1198 // return number of added samples
1199 
1200 static int SingleFaceSubdivision(int sampleNum, const CoordType & v0, const CoordType & v1, const CoordType & v2, VertexSampler &ps, FacePointer fp, bool randSample)
1201 {
1202  // recursive face subdivision.
1203  if(sampleNum == 1)
1204  {
1205  // ground case.
1206  CoordType SamplePoint;
1207  if(randSample)
1208  {
1209  CoordType rb=RandomBarycentric();
1210  SamplePoint=v0*rb[0]+v1*rb[1]+v2*rb[2];
1211  }
1212  else SamplePoint=((v0+v1+v2)*(1.0f/3.0f));
1213 
1214  ps.AddFace(*fp,SamplePoint);
1215  return 1;
1216  }
1217 
1218  int s0 = sampleNum /2;
1219  int s1 = sampleNum-s0;
1220  assert(s0>0);
1221  assert(s1>0);
1222 
1223  ScalarType w0 = ScalarType(s1)/ScalarType(sampleNum);
1224  ScalarType w1 = 1.0-w0;
1225  // compute the longest edge.
1226  ScalarType maxd01 = SquaredDistance(v0,v1);
1227  ScalarType maxd12 = SquaredDistance(v1,v2);
1228  ScalarType maxd20 = SquaredDistance(v2,v0);
1229  int res;
1230  if(maxd01 > maxd12)
1231  if(maxd01 > maxd20) res = 0;
1232  else res = 2;
1233  else
1234  if(maxd12 > maxd20) res = 1;
1235  else res = 2;
1236 
1237  int faceSampleNum=0;
1238  // break the input triangle along the midpoint of the longest edge.
1239  CoordType pp;
1240  switch(res)
1241  {
1242  case 0 : pp = v0*w0 + v1*w1;
1243  faceSampleNum+=SingleFaceSubdivision(s0,v0,pp,v2,ps,fp,randSample);
1244  faceSampleNum+=SingleFaceSubdivision(s1,pp,v1,v2,ps,fp,randSample);
1245  break;
1246  case 1 : pp = v1*w0 + v2*w1;
1247  faceSampleNum+=SingleFaceSubdivision(s0,v0,v1,pp,ps,fp,randSample);
1248  faceSampleNum+=SingleFaceSubdivision(s1,v0,pp,v2,ps,fp,randSample);
1249  break;
1250  case 2 : pp = v0*w0 + v2*w1;
1251  faceSampleNum+=SingleFaceSubdivision(s0,v0,v1,pp,ps,fp,randSample);
1252  faceSampleNum+=SingleFaceSubdivision(s1,pp,v1,v2,ps,fp,randSample);
1253  break;
1254  }
1255  return faceSampleNum;
1256 }
1257 
1258 
1260 static void FaceSubdivision(MeshType & m, VertexSampler &ps,int sampleNum, bool randSample)
1261 {
1262 
1263  ScalarType area = Stat<MeshType>::ComputeMeshArea(m);
1264  ScalarType samplePerAreaUnit = sampleNum/area;
1265  std::vector<FacePointer> faceVec;
1266  FillAndShuffleFacePointerVector(m,faceVec);
1268  double floatSampleNum = 0.0;
1269  int faceSampleNum;
1270  // Subdivision sampling.
1271  typename std::vector<FacePointer>::iterator fi;
1272  for(fi=faceVec.begin(); fi!=faceVec.end(); fi++)
1273  {
1274  const CoordType b0(1.0, 0.0, 0.0);
1275  const CoordType b1(0.0, 1.0, 0.0);
1276  const CoordType b2(0.0, 0.0, 1.0);
1277  // compute # samples in the current face.
1278  floatSampleNum += 0.5*DoubleArea(**fi) * samplePerAreaUnit;
1279  faceSampleNum = (int) floatSampleNum;
1280  if(faceSampleNum>0)
1281  faceSampleNum = SingleFaceSubdivision(faceSampleNum,b0,b1,b2,ps,*fi,randSample);
1282  floatSampleNum -= (double) faceSampleNum;
1283  }
1284 }
1285 //---------
1286 // Subdivision sampling of a single face.
1287 // return number of added samples
1288 
1289 static int SingleFaceSubdivisionOld(int sampleNum, const CoordType & v0, const CoordType & v1, const CoordType & v2, VertexSampler &ps, FacePointer fp, bool randSample)
1290 {
1291  // recursive face subdivision.
1292  if(sampleNum == 1)
1293  {
1294  // ground case.
1295  CoordType SamplePoint;
1296  if(randSample)
1297  {
1298  CoordType rb=RandomBarycentric();
1299  SamplePoint=v0*rb[0]+v1*rb[1]+v2*rb[2];
1300  }
1301  else SamplePoint=((v0+v1+v2)*(1.0f/3.0f));
1302 
1303  CoordType SampleBary;
1304  InterpolationParameters(*fp,SamplePoint,SampleBary);
1305  ps.AddFace(*fp,SampleBary);
1306  return 1;
1307  }
1308 
1309  int s0 = sampleNum /2;
1310  int s1 = sampleNum-s0;
1311  assert(s0>0);
1312  assert(s1>0);
1313 
1314  ScalarType w0 = ScalarType(s1)/ScalarType(sampleNum);
1315  ScalarType w1 = 1.0-w0;
1316  // compute the longest edge.
1317  ScalarType maxd01 = SquaredDistance(v0,v1);
1318  ScalarType maxd12 = SquaredDistance(v1,v2);
1319  ScalarType maxd20 = SquaredDistance(v2,v0);
1320  int res;
1321  if(maxd01 > maxd12)
1322  if(maxd01 > maxd20) res = 0;
1323  else res = 2;
1324  else
1325  if(maxd12 > maxd20) res = 1;
1326  else res = 2;
1327 
1328  int faceSampleNum=0;
1329  // break the input triangle along the midpoint of the longest edge.
1330  CoordType pp;
1331  switch(res)
1332  {
1333  case 0 : pp = v0*w0 + v1*w1;
1334  faceSampleNum+=SingleFaceSubdivision(s0,v0,pp,v2,ps,fp,randSample);
1335  faceSampleNum+=SingleFaceSubdivision(s1,pp,v1,v2,ps,fp,randSample);
1336  break;
1337  case 1 : pp = v1*w0 + v2*w1;
1338  faceSampleNum+=SingleFaceSubdivision(s0,v0,v1,pp,ps,fp,randSample);
1339  faceSampleNum+=SingleFaceSubdivision(s1,v0,pp,v2,ps,fp,randSample);
1340  break;
1341  case 2 : pp = v0*w0 + v2*w1;
1342  faceSampleNum+=SingleFaceSubdivision(s0,v0,v1,pp,ps,fp,randSample);
1343  faceSampleNum+=SingleFaceSubdivision(s1,pp,v1,v2,ps,fp,randSample);
1344  break;
1345  }
1346  return faceSampleNum;
1347 }
1348 
1349 
1351 static void FaceSubdivisionOld(MeshType & m, VertexSampler &ps,int sampleNum, bool randSample)
1352 {
1353 
1354  ScalarType area = Stat<MeshType>::ComputeMeshArea(m);
1355  ScalarType samplePerAreaUnit = sampleNum/area;
1356  std::vector<FacePointer> faceVec;
1357  FillAndShuffleFacePointerVector(m,faceVec);
1359  double floatSampleNum = 0.0;
1360  int faceSampleNum;
1361  // Subdivision sampling.
1362  typename std::vector<FacePointer>::iterator fi;
1363  for(fi=faceVec.begin(); fi!=faceVec.end(); fi++)
1364  {
1365  // compute # samples in the current face.
1366  floatSampleNum += 0.5*DoubleArea(**fi) * samplePerAreaUnit;
1367  faceSampleNum = (int) floatSampleNum;
1368  if(faceSampleNum>0)
1369  faceSampleNum = SingleFaceSubdivision(faceSampleNum,(**fi).V(0)->cP(), (**fi).V(1)->cP(), (**fi).V(2)->cP(),ps,*fi,randSample);
1370  floatSampleNum -= (double) faceSampleNum;
1371  }
1372 }
1373 
1374 
1375 //---------
1376 
1377 // Similar Triangles sampling.
1378 // Skip vertex and edges
1379 // Sample per edges includes vertexes, so here we should expect n_samples_per_edge >=4
1380 
1381 static int SingleFaceSimilar(FacePointer fp, VertexSampler &ps, int n_samples_per_edge)
1382 {
1383  int n_samples=0;
1384  int i, j;
1385  float segmentNum=n_samples_per_edge -1 ;
1386  float segmentLen = 1.0/segmentNum;
1387  // face sampling.
1388  for(i=1; i < n_samples_per_edge-1; i++)
1389  for(j=1; j < n_samples_per_edge-1-i; j++)
1390  {
1391  //AddSample( v0 + (V1*(double)i + V2*(double)j) );
1392  CoordType sampleBary(i*segmentLen,j*segmentLen, 1.0 - (i*segmentLen+j*segmentLen) ) ;
1393  n_samples++;
1394  ps.AddFace(*fp,sampleBary);
1395  }
1396  return n_samples;
1397 }
1398 static int SingleFaceSimilarDual(FacePointer fp, VertexSampler &ps, int n_samples_per_edge, bool randomFlag)
1399 {
1400  int n_samples=0;
1401  float i, j;
1402  float segmentNum=n_samples_per_edge -1 ;
1403  float segmentLen = 1.0/segmentNum;
1404  // face sampling.
1405  for(i=0; i < n_samples_per_edge-1; i++)
1406  for(j=0; j < n_samples_per_edge-1-i; j++)
1407  {
1408  //AddSample( v0 + (V1*(double)i + V2*(double)j) );
1409  CoordType V0((i+0)*segmentLen,(j+0)*segmentLen, 1.0 - ((i+0)*segmentLen+(j+0)*segmentLen) ) ;
1410  CoordType V1((i+1)*segmentLen,(j+0)*segmentLen, 1.0 - ((i+1)*segmentLen+(j+0)*segmentLen) ) ;
1411  CoordType V2((i+0)*segmentLen,(j+1)*segmentLen, 1.0 - ((i+0)*segmentLen+(j+1)*segmentLen) ) ;
1412  n_samples++;
1413  if(randomFlag) {
1414  CoordType rb=RandomBarycentric();
1415  ps.AddFace(*fp, V0*rb[0]+V1*rb[1]+V2*rb[2]);
1416  } else ps.AddFace(*fp,(V0+V1+V2)/3.0);
1417 
1418  if( j < n_samples_per_edge-i-2 )
1419  {
1420  CoordType V3((i+1)*segmentLen,(j+1)*segmentLen, 1.0 - ((i+1)*segmentLen+(j+1)*segmentLen) ) ;
1421  n_samples++;
1422  if(randomFlag) {
1423  CoordType rb=RandomBarycentric();
1424  ps.AddFace(*fp, V3*rb[0]+V1*rb[1]+V2*rb[2]);
1425  } else ps.AddFace(*fp,(V3+V1+V2)/3.0);
1426  }
1427  }
1428  return n_samples;
1429 }
1430 
1431 // Similar sampling
1432 // Each triangle is subdivided into similar triangles following a generalization of the classical 1-to-4 splitting rule of triangles.
1433 // According to the level of subdivision <k> you get 1, 4 , 9, 16 , <k^2> triangles.
1434 // Depending on the kind of the sampling strategies we can have two different approach to choosing the sample points.
1435 // 1) you have already sampled both edges and vertices
1436 // 2) you are not going to take samples on edges and vertices.
1437 //
1438 // In the first case you have to consider only internal vertices of the subdivided triangles (to avoid multiple sampling of edges and vertices).
1439 // Therefore the number of internal points is ((k-3)*(k-2))/2. where k is the number of points on an edge (vertex included)
1440 // E.g. for k=4 you get 3 segments on each edges and the original triangle is subdivided
1441 // into 9 smaller triangles and you get (1*2)/2 == 1 only a single internal point.
1442 // So if you want N samples in a triangle you have to solve k^2 -5k +6 - 2N = 0
1443 // from which you get:
1444 //
1445 // 5 + sqrt( 1 + 8N )
1446 // k = -------------------
1447 // 2
1448 //
1449 // 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.
1450 // So if you want N samples in a triangle, the number <k> of points on an edge (vertex included) should be simply:
1451 // k = 1 + sqrt(N)
1452 // examples:
1453 // N = 4 -> k = 3
1454 // N = 9 -> k = 4
1455 
1456 
1457 
1458 //template <class MeshType>
1459 //void Sampling<MeshType>::SimilarFaceSampling()
1460 static void FaceSimilar(MeshType & m, VertexSampler &ps,int sampleNum, bool dualFlag, bool randomFlag)
1461 {
1462  ScalarType area = Stat<MeshType>::ComputeMeshArea(m);
1463  ScalarType samplePerAreaUnit = sampleNum/area;
1464 
1465  // Similar Triangles sampling.
1466  int n_samples_per_edge;
1467  double n_samples_decimal = 0.0;
1468  FaceIterator fi;
1469 
1470  for(fi=m.face.begin(); fi != m.face.end(); fi++)
1471  {
1472  // compute # samples in the current face.
1473  n_samples_decimal += 0.5*DoubleArea(*fi) * samplePerAreaUnit;
1474  int n_samples = (int) n_samples_decimal;
1475  if(n_samples>0)
1476  {
1477  // face sampling.
1478  if(dualFlag)
1479  {
1480  n_samples_per_edge = (int)((sqrt(1.0+8.0*(double)n_samples) +5.0)/2.0); // original for non dual case
1481  n_samples = SingleFaceSimilar(&*fi,ps, n_samples_per_edge);
1482  } else {
1483  n_samples_per_edge = (int)(sqrt((double)n_samples) +1.0);
1484  n_samples = SingleFaceSimilarDual(&*fi,ps, n_samples_per_edge,randomFlag);
1485  }
1486  }
1487  n_samples_decimal -= (double) n_samples;
1488  }
1489 }
1490 
1491 
1492  // Rasterization fuction
1493  // Take a triangle
1494  // T deve essere una classe funzionale che ha l'operatore ()
1495  // con due parametri x,y di tipo S esempio:
1496  // class Foo { public void operator()(int x, int y ) { ??? } };
1497 
1498 // This function does rasterization with a safety buffer area, thus accounting some points actually outside triangle area
1499 // The safety area samples are generated according to face flag BORDER which should be true for texture space border edges
1500 // Use correctSafePointsBaryCoords = true to map safety texels to closest point barycentric coords (on edge).
1501  static void SingleFaceRaster(typename MeshType::FaceType &f, VertexSampler &ps,
1502  const Point2<typename MeshType::ScalarType> & v0,
1503  const Point2<typename MeshType::ScalarType> & v1,
1504  const Point2<typename MeshType::ScalarType> & v2,
1505  bool correctSafePointsBaryCoords=true)
1506  {
1507  typedef typename MeshType::ScalarType S;
1508  // Calcolo bounding box
1509  Box2i bbox;
1510  Box2<S> bboxf;
1511  bboxf.Add(v0);
1512  bboxf.Add(v1);
1513  bboxf.Add(v2);
1514 
1515  bbox.min[0] = floor(bboxf.min[0]);
1516  bbox.min[1] = floor(bboxf.min[1]);
1517  bbox.max[0] = ceil(bboxf.max[0]);
1518  bbox.max[1] = ceil(bboxf.max[1]);
1519 
1520  // Calcolo versori degli spigoli
1521  Point2<S> d10 = v1 - v0;
1522  Point2<S> d21 = v2 - v1;
1523  Point2<S> d02 = v0 - v2;
1524 
1525  // Preparazione prodotti scalari
1526  S b0 = (bbox.min[0]-v0[0])*d10[1] - (bbox.min[1]-v0[1])*d10[0];
1527  S b1 = (bbox.min[0]-v1[0])*d21[1] - (bbox.min[1]-v1[1])*d21[0];
1528  S b2 = (bbox.min[0]-v2[0])*d02[1] - (bbox.min[1]-v2[1])*d02[0];
1529  // Preparazione degli steps
1530  S db0 = d10[1];
1531  S db1 = d21[1];
1532  S db2 = d02[1];
1533  // Preparazione segni
1534  S dn0 = -d10[0];
1535  S dn1 = -d21[0];
1536  S dn2 = -d02[0];
1537 
1538  //Calculating orientation
1539  bool flipped = !(d02 * vcg::Point2<S>(-d10[1], d10[0]) >= 0);
1540 
1541  // Calculating border edges
1542  Segment2<S> borderEdges[3];
1543  S edgeLength[3];
1544  unsigned char edgeMask = 0;
1545 
1546  if (f.IsB(0)) {
1547  borderEdges[0] = Segment2<S>(v0, v1);
1548  edgeLength[0] = borderEdges[0].Length();
1549  edgeMask |= 1;
1550  }
1551  if (f.IsB(1)) {
1552  borderEdges[1] = Segment2<S>(v1, v2);
1553  edgeLength[1] = borderEdges[1].Length();
1554  edgeMask |= 2;
1555  }
1556  if (f.IsB(2)) {
1557  borderEdges[2] = Segment2<S>(v2, v0);
1558  edgeLength[2] = borderEdges[2].Length();
1559  edgeMask |= 4;
1560  }
1561 
1562  // Rasterizzazione
1563  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];
1564 
1565  for(int x=bbox.min[0]-1;x<=bbox.max[0]+1;++x)
1566  {
1567  bool in = false;
1568  S n[3] = { b0-db0-dn0, b1-db1-dn1, b2-db2-dn2};
1569  for(int y=bbox.min[1]-1;y<=bbox.max[1]+1;++y)
1570  {
1571  if( ((n[0]>=0 && n[1]>=0 && n[2]>=0) || (n[0]<=0 && n[1]<=0 && n[2]<=0)) && (de != 0))
1572  {
1573  typename MeshType::CoordType baryCoord;
1574  baryCoord[0] = double(-y*v1[0]+v2[0]*y+v1[1]*x-v2[0]*v1[1]+v1[0]*v2[1]-x*v2[1])/de;
1575  baryCoord[1] = -double( x*v0[1]-x*v2[1]-v0[0]*y+v0[0]*v2[1]-v2[0]*v0[1]+v2[0]*y)/de;
1576  baryCoord[2] = 1-baryCoord[0]-baryCoord[1];
1577 
1578  ps.AddTextureSample(f, baryCoord, Point2i(x,y), 0);
1579  in = true;
1580  } else {
1581  // Check whether a pixel outside (on a border edge side) triangle affects color inside it
1582  Point2<S> px(x, y);
1583  Point2<S> closePoint;
1584  int closeEdge = -1;
1585  S minDst = FLT_MAX;
1586 
1587  // find the closest point (on some edge) that lies on the 2x2 squared neighborhood of the considered point
1588  for (int i=0; i<3; ++i)
1589  {
1590  if (edgeMask & (1 << i))
1591  {
1592  Point2<S> close;
1593  S dst;
1594  if ( ((!flipped) && (n[i]<0)) ||
1595  ( flipped && (n[i]>0)) )
1596  {
1597  dst = ((close = ClosestPoint(borderEdges[i], px)) - px).Norm();
1598  if(dst < minDst &&
1599  close.X() > px.X()-1 && close.X() < px.X()+1 &&
1600  close.Y() > px.Y()-1 && close.Y() < px.Y()+1)
1601  {
1602  minDst = dst;
1603  closePoint = close;
1604  closeEdge = i;
1605  }
1606  }
1607  }
1608  }
1609 
1610  if (closeEdge >= 0)
1611  {
1612  typename MeshType::CoordType baryCoord;
1613  if (correctSafePointsBaryCoords)
1614  {
1615  // Add x,y sample with closePoint barycentric coords (on edge)
1616  baryCoord[closeEdge] = (closePoint - borderEdges[closeEdge].P1()).Norm()/edgeLength[closeEdge];
1617  baryCoord[(closeEdge+1)%3] = 1 - baryCoord[closeEdge];
1618  baryCoord[(closeEdge+2)%3] = 0;
1619  } else {
1620  // Add x,y sample with his own barycentric coords (off edge)
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), minDst);
1626  in = true;
1627  }
1628  }
1629  n[0] += dn0;
1630  n[1] += dn1;
1631  n[2] += dn2;
1632  }
1633  b0 += db0;
1634  b1 += db1;
1635  b2 += db2;
1636  }
1637 }
1638 
1639 // check the radius constrain
1640 static bool checkPoissonDisk(SampleSHT & sht, const Point3<ScalarType> & p, ScalarType radius)
1641 {
1642  // get the samples closest to the given one
1643  std::vector<VertexType*> closests;
1644  typedef EmptyTMark<MeshType> MarkerVert;
1645  static MarkerVert mv;
1646 
1647  Box3f bb(p-Point3f(radius,radius,radius),p+Point3f(radius,radius,radius));
1648  GridGetInBox(sht, mv, bb, closests);
1649 
1650  ScalarType r2 = radius*radius;
1651  for(int i=0; i<closests.size(); ++i)
1652  if(SquaredDistance(p,closests[i]->cP()) < r2)
1653  return false;
1654 
1655  return true;
1656 }
1657 
1658 struct PoissonDiskParam
1659 {
1660  PoissonDiskParam()
1661  {
1662  adaptiveRadiusFlag = false;
1663  bestSampleChoiceFlag = true;
1664  bestSamplePoolSize = 10;
1665  radiusVariance =1;
1666  MAXLEVELS = 5;
1667  invertQuality = false;
1668  preGenFlag = false;
1669  preGenMesh = NULL;
1670  geodesicDistanceFlag = false;
1671  randomSeed = 0;
1672  }
1673 
1674  struct Stat
1675  {
1676  int montecarloTime;
1677  int gridTime;
1678  int pruneTime;
1679  int totalTime;
1680  Point3i gridSize;
1681  int gridCellNum;
1682  size_t sampleNum;
1683  int montecarloSampleNum;
1684  };
1685 
1686  bool geodesicDistanceFlag;
1687  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.
1688  int bestSamplePoolSize;
1689  bool adaptiveRadiusFlag;
1690  float radiusVariance;
1691  bool invertQuality;
1692  bool preGenFlag; // when generating a poisson distribution, you can initialize the set of computed points with
1693  // ALL the vertices of another mesh. Useful for building progressive//prioritize refinements.
1694  MeshType *preGenMesh; // There are two ways of passing the pregen vertexes to the pruning, 1) is with a mesh pointer
1695  // 2) with a per vertex attribute.
1696  int MAXLEVELS;
1697  int randomSeed;
1698 
1699  Stat pds;
1700 };
1701 
1702 
1703 // generate Poisson-disk sample using a set of pre-generated samples (with the Montecarlo algorithm)
1704 // It always return a point.
1705 static VertexPointer getSampleFromCell(Point3i &cell, MontecarloSHT & samplepool)
1706 {
1707  MontecarloSHTIterator cellBegin, cellEnd;
1708  samplepool.Grid(cell, cellBegin, cellEnd);
1709  return *cellBegin;
1710 }
1711 
1712 // Given a cell of the grid it search the point that remove the minimum number of other samples
1713 // it linearly scan all the points of a cell.
1714 
1715 static VertexPointer getBestPrecomputedMontecarloSample(Point3i &cell, MontecarloSHT & samplepool, ScalarType diskRadius, const PoissonDiskParam &pp)
1716 {
1717  MontecarloSHTIterator cellBegin,cellEnd;
1718  samplepool.Grid(cell, cellBegin, cellEnd);
1719  VertexPointer bestSample=0;
1720  int minRemoveCnt = std::numeric_limits<int>::max();
1721  std::vector<typename MontecarloSHT::HashIterator> inSphVec;
1722  int i=0;
1723  for(MontecarloSHTIterator ci=cellBegin; ci!=cellEnd && i<pp.bestSamplePoolSize; ++ci,i++)
1724  {
1725  VertexPointer sp = *ci;
1726  if(pp.adaptiveRadiusFlag) diskRadius = sp->Q();
1727  int curRemoveCnt = samplepool.CountInSphere(sp->cP(),diskRadius,inSphVec);
1728  if(curRemoveCnt < minRemoveCnt)
1729  {
1730  bestSample = sp;
1731  minRemoveCnt = curRemoveCnt;
1732  }
1733  }
1734  return bestSample;
1735 }
1736 
1739 static ScalarType ComputePoissonDiskRadius(MeshType &origMesh, int sampleNum)
1740 {
1741  ScalarType meshArea = Stat<MeshType>::ComputeMeshArea(origMesh);
1742  // Manage approximately the PointCloud Case, use the half a area of the bbox.
1743  // TODO: If you had the radius a much better approximation could be done.
1744  if(meshArea ==0)
1745  {
1746  meshArea = (origMesh.bbox.DimX()*origMesh.bbox.DimY() +
1747  origMesh.bbox.DimX()*origMesh.bbox.DimZ() +
1748  origMesh.bbox.DimY()*origMesh.bbox.DimZ());
1749  }
1750  ScalarType diskRadius = sqrt(meshArea / (0.7 * M_PI * sampleNum)); // 0.7 is a density factor
1751  return diskRadius;
1752 }
1753 
1754 static int ComputePoissonSampleNum(MeshType &origMesh, ScalarType diskRadius)
1755 {
1756  ScalarType meshArea = Stat<MeshType>::ComputeMeshArea(origMesh);
1757  int sampleNum = meshArea / (diskRadius*diskRadius *M_PI *0.7) ; // 0.7 is a density factor
1758  return sampleNum;
1759 }
1760 
1767 
1768 static void InitRadiusHandleFromQuality(MeshType &sampleMesh, PerVertexFloatAttribute &rH, ScalarType diskRadius, ScalarType radiusVariance, bool invert)
1769 {
1770  std::pair<float,float> minmax = tri::Stat<MeshType>::ComputePerVertexQualityMinMax( sampleMesh);
1771  float minRad = diskRadius ;
1772  float maxRad = diskRadius * radiusVariance;
1773  float deltaQ = minmax.second-minmax.first;
1774  float deltaRad = maxRad-minRad;
1775  for (VertexIterator vi = sampleMesh.vert.begin(); vi != sampleMesh.vert.end(); vi++)
1776  {
1777  rH[*vi] = minRad + deltaRad*((invert ? minmax.second - (*vi).Q() : (*vi).Q() - minmax.first )/deltaQ);
1778  }
1779 }
1780 
1781 // initialize spatial hash table for searching
1782 // radius is the radius of empty disk centered over the samples (e.g. twice of the empty space disk)
1783 // This radius implies that when we pick a sample in a cell all that cell probably will not be touched again.
1784 // Howvever we must ensure that we do not put too many vertices inside each hash cell
1785 
1786 static void InitSpatialHashTable(MeshType &montecarloMesh, MontecarloSHT &montecarloSHT, ScalarType diskRadius,
1787  struct PoissonDiskParam pp=PoissonDiskParam())
1788 {
1789  ScalarType cellsize = 2.0f* diskRadius / sqrt(3.0);
1790  float occupancyRatio=0;
1791  do
1792  {
1793  // inflating
1794  BoxType bb=montecarloMesh.bbox;
1795  assert(!bb.IsNull());
1796  bb.Offset(cellsize);
1797 
1798  int sizeX = std::max(1,int(bb.DimX() / cellsize));
1799  int sizeY = std::max(1,int(bb.DimY() / cellsize));
1800  int sizeZ = std::max(1,int(bb.DimZ() / cellsize));
1801  Point3i gridsize(sizeX, sizeY, sizeZ);
1802 
1803  montecarloSHT.InitEmpty(bb, gridsize);
1804 
1805  for (VertexIterator vi = montecarloMesh.vert.begin(); vi != montecarloMesh.vert.end(); vi++)
1806  if(!(*vi).IsD())
1807  {
1808  montecarloSHT.Add(&(*vi));
1809  }
1810 
1811  montecarloSHT.UpdateAllocatedCells();
1812  pp.pds.gridSize = gridsize;
1813  pp.pds.gridCellNum = (int)montecarloSHT.AllocatedCells.size();
1814  cellsize/=2.0f;
1815  occupancyRatio = float(montecarloMesh.vn) / float(montecarloSHT.AllocatedCells.size());
1816  // qDebug(" %i / %i = %6.3f", montecarloMesh.vn , montecarloSHT.AllocatedCells.size(),occupancyRatio);
1817  }
1818  while( occupancyRatio> 100);
1819 }
1820 
1821 static void PoissonDiskPruningByNumber(VertexSampler &ps, MeshType &m,
1822  size_t sampleNum, ScalarType &diskRadius,
1823  PoissonDiskParam &pp,
1824  float tolerance=0.04,
1825  int maxIter=20)
1826 
1827 {
1828  size_t sampleNumMin = int(float(sampleNum)*(1.0f-tolerance));
1829  size_t sampleNumMax = int(float(sampleNum)*(1.0f+tolerance));
1830  float RangeMinRad = m.bbox.Diag()/50.0;
1831  float RangeMaxRad = m.bbox.Diag()/50.0;
1832  size_t RangeMinRadNum;
1833  size_t RangeMaxRadNum;
1834  // Note RangeMinRad < RangeMaxRad
1835  // but RangeMinRadNum > sampleNum > RangeMaxRadNum
1836  do {
1837  ps.reset();
1838  RangeMinRad/=2.0f;
1839  PoissonDiskPruning(ps, m ,RangeMinRad,pp);
1840  RangeMinRadNum = pp.pds.sampleNum;
1841 // qDebug("PoissonDiskPruning Iteratin Min (%6.3f:%5i) instead of %i",RangeMinRad,RangeMinRadNum,sampleNum);
1842  } while(RangeMinRadNum < sampleNum); // if the number of sample is still smaller you have to make radius larger.
1843 
1844  do {
1845  ps.reset();
1846  RangeMaxRad*=2.0f;
1847  PoissonDiskPruning(ps, m ,RangeMaxRad,pp);
1848  RangeMaxRadNum = pp.pds.sampleNum;
1849 // qDebug("PoissonDiskPruning Iteratin Max (%6.3f:%5i) instead of %i",RangeMaxRad,RangeMaxRadNum,sampleNum);
1850  } while(RangeMaxRadNum > sampleNum);
1851 
1852 
1853  float curRadius=RangeMaxRad;
1854  int iterCnt=0;
1855  while(iterCnt<maxIter &&
1856  (pp.pds.sampleNum < sampleNumMin || pp.pds.sampleNum > sampleNumMax) )
1857  {
1858  iterCnt++;
1859  ps.reset();
1860  curRadius=(RangeMaxRad+RangeMinRad)/2.0f;
1861  PoissonDiskPruning(ps, m ,curRadius,pp);
1862 // 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);
1863  if(pp.pds.sampleNum > sampleNum){
1864  RangeMinRad = curRadius;
1865  RangeMinRadNum = pp.pds.sampleNum;
1866  }
1867  if(pp.pds.sampleNum < sampleNum){
1868  RangeMaxRad = curRadius;
1869  RangeMaxRadNum = pp.pds.sampleNum;
1870  }
1871  }
1872  diskRadius = curRadius;
1873 }
1874 
1875 
1885 static void PoissonDiskPruning(VertexSampler &ps, MeshType &montecarloMesh,
1886  ScalarType diskRadius, PoissonDiskParam &pp)
1887 {
1888  tri::RequireCompactness(montecarloMesh);
1889  if(pp.randomSeed) SamplingRandomGenerator().initialize(pp.randomSeed);
1890  if(pp.adaptiveRadiusFlag)
1891  tri::RequirePerVertexQuality(montecarloMesh);
1892  int t0 = clock();
1893  // spatial index of montecarlo samples - used to choose a new sample to insert
1894  MontecarloSHT montecarloSHT;
1895  InitSpatialHashTable(montecarloMesh,montecarloSHT,diskRadius,pp);
1896 
1897  // if we are doing variable density sampling we have to prepare the handle that keeps the the random samples expected radii.
1898  // At this point we just assume that there is the quality values as sampled from the base mesh
1899  PerVertexFloatAttribute rH = tri::Allocator<MeshType>:: template GetPerVertexAttribute<float> (montecarloMesh,"radius");
1900  if(pp.adaptiveRadiusFlag)
1901  InitRadiusHandleFromQuality(montecarloMesh, rH, diskRadius, pp.radiusVariance, pp.invertQuality);
1902 
1903  unsigned int (*p_myrandom)(unsigned int) = RandomInt;
1904  std::random_shuffle(montecarloSHT.AllocatedCells.begin(),montecarloSHT.AllocatedCells.end(), p_myrandom);
1905  int t1 = clock();
1906  pp.pds.montecarloSampleNum = montecarloMesh.vn;
1907  pp.pds.sampleNum =0;
1908  int removedCnt=0;
1909  // Initial pass for pruning the Hashed grid with the an eventual pre initialized set of samples
1910  if(pp.preGenFlag)
1911  {
1912  if(pp.preGenMesh==0)
1913  {
1914  typename MeshType::template PerVertexAttributeHandle<bool> fixed;
1915  fixed = tri::Allocator<MeshType>:: template GetPerVertexAttribute<bool> (montecarloMesh,"fixed");
1916  for(VertexIterator vi=montecarloMesh.vert.begin();vi!=montecarloMesh.vert.end();++vi)
1917  if(fixed[*vi]) {
1918  pp.pds.sampleNum++;
1919  ps.AddVert(*vi);
1920  removedCnt += montecarloSHT.RemoveInSphere(vi->cP(),diskRadius);
1921  }
1922  }
1923  else
1924  {
1925  for(VertexIterator vi =pp.preGenMesh->vert.begin(); vi!=pp.preGenMesh->vert.end();++vi)
1926  {
1927  ps.AddVert(*vi);
1928  pp.pds.sampleNum++;
1929  removedCnt += montecarloSHT.RemoveInSphere(vi->cP(),diskRadius);
1930  }
1931  }
1932  montecarloSHT.UpdateAllocatedCells();
1933  }
1934  vertex::ApproximateGeodesicDistanceFunctor<VertexType> GDF;
1935  while(!montecarloSHT.AllocatedCells.empty())
1936  {
1937  removedCnt=0;
1938  for (size_t i = 0; i < montecarloSHT.AllocatedCells.size(); i++)
1939  {
1940  if( montecarloSHT.EmptyCell(montecarloSHT.AllocatedCells[i]) ) continue;
1941  ScalarType currentRadius =diskRadius;
1942  VertexPointer sp;
1943  if(pp.bestSampleChoiceFlag)
1944  sp = getBestPrecomputedMontecarloSample(montecarloSHT.AllocatedCells[i], montecarloSHT, diskRadius, pp);
1945  else
1946  sp = getSampleFromCell(montecarloSHT.AllocatedCells[i], montecarloSHT);
1947 
1948  if(pp.adaptiveRadiusFlag)
1949  currentRadius = rH[sp];
1950 
1951  ps.AddVert(*sp);
1952  pp.pds.sampleNum++;
1953  if(pp.geodesicDistanceFlag) removedCnt += montecarloSHT.RemoveInSphereNormal(sp->cP(),sp->cN(),GDF,currentRadius);
1954  else removedCnt += montecarloSHT.RemoveInSphere(sp->cP(),currentRadius);
1955  }
1956  montecarloSHT.UpdateAllocatedCells();
1957  }
1958  int t2 = clock();
1959  pp.pds.gridTime = t1-t0;
1960  pp.pds.pruneTime = t2-t1;
1961 }
1962 
1973 static void HierarchicalPoissonDisk(MeshType &origMesh, VertexSampler &ps, MeshType &montecarloMesh, ScalarType diskRadius, const struct PoissonDiskParam pp=PoissonDiskParam())
1974 {
1975 // int t0=clock();
1976  // spatial index of montecarlo samples - used to choose a new sample to insert
1977  MontecarloSHT montecarloSHTVec[5];
1978 
1979 
1980 
1981  // initialize spatial hash table for searching
1982  // radius is the radius of empty disk centered over the samples (e.g. twice of the empty space disk)
1983  // This radius implies that when we pick a sample in a cell all that cell will not be touched again.
1984  ScalarType cellsize = 2.0f* diskRadius / sqrt(3.0);
1985 
1986  // inflating
1987  BoxType bb=origMesh.bbox;
1988  bb.Offset(cellsize);
1989 
1990  int sizeX = std::max(1.0f,bb.DimX() / cellsize);
1991  int sizeY = std::max(1.0f,bb.DimY() / cellsize);
1992  int sizeZ = std::max(1.0f,bb.DimZ() / cellsize);
1993  Point3i gridsize(sizeX, sizeY, sizeZ);
1994 
1995  // spatial hash table of the generated samples - used to check the radius constrain
1996  SampleSHT checkSHT;
1997  checkSHT.InitEmpty(bb, gridsize);
1998 
1999 
2000  // sampling algorithm
2001  // ------------------
2002  //
2003  // - generate millions of samples using montecarlo algorithm
2004  // - extract a cell (C) from the active cell list (with probability proportional to cell's volume)
2005  // - generate a sample inside C by choosing one of the contained pre-generated samples
2006  // - if the sample violates the radius constrain discard it, and add the cell to the cells-to-subdivide list
2007  // - iterate until the active cell list is empty or a pre-defined number of subdivisions is reached
2008  //
2009 
2010  int level = 0;
2011 
2012  // initialize spatial hash to index pre-generated samples
2013  montecarloSHTVec[0].InitEmpty(bb, gridsize);
2014  // create active cell list
2015  for (VertexIterator vi = montecarloMesh.vert.begin(); vi != montecarloMesh.vert.end(); vi++)
2016  montecarloSHTVec[0].Add(&(*vi));
2017  montecarloSHTVec[0].UpdateAllocatedCells();
2018 
2019  // if we are doing variable density sampling we have to prepare the random samples quality with the correct expected radii.
2020  PerVertexFloatAttribute rH = tri::Allocator<MeshType>:: template GetPerVertexAttribute<float> (montecarloMesh,"radius");
2021  if(pp.adaptiveRadiusFlag)
2022  InitRadiusHandleFromQuality(montecarloMesh, rH, diskRadius, pp.radiusVariance, pp.invertQuality);
2023 
2024  do
2025  {
2026  MontecarloSHT &montecarloSHT = montecarloSHTVec[level];
2027 
2028  if(level>0)
2029  {// initialize spatial hash with the remaining points
2030  montecarloSHT.InitEmpty(bb, gridsize);
2031  // create active cell list
2032  for (typename MontecarloSHT::HashIterator hi = montecarloSHTVec[level-1].hash_table.begin(); hi != montecarloSHTVec[level-1].hash_table.end(); hi++)
2033  montecarloSHT.Add((*hi).second);
2034  montecarloSHT.UpdateAllocatedCells();
2035  }
2036  // shuffle active cells
2037  unsigned int (*p_myrandom)(unsigned int) = RandomInt;
2038  std::random_shuffle(montecarloSHT.AllocatedCells.begin(),montecarloSHT.AllocatedCells.end(), p_myrandom);
2039 
2040  // generate a sample inside C by choosing one of the contained pre-generated samples
2042  int removedCnt=montecarloSHT.hash_table.size();
2043  int addedCnt=checkSHT.hash_table.size();
2044  for (int i = 0; i < montecarloSHT.AllocatedCells.size(); i++)
2045  {
2046  for(int j=0;j<4;j++)
2047  {
2048  if( montecarloSHT.EmptyCell(montecarloSHT.AllocatedCells[i]) ) continue;
2049 
2050  // generate a sample chosen from the pre-generated one
2051  typename MontecarloSHT::HashIterator hi = montecarloSHT.hash_table.find(montecarloSHT.AllocatedCells[i]);
2052 
2053  if(hi==montecarloSHT.hash_table.end()) {break;}
2054  VertexPointer sp = (*hi).second;
2055  // vr spans between 3.0*r and r / 4.0 according to vertex quality
2056  ScalarType sampleRadius = diskRadius;
2057  if(pp.adaptiveRadiusFlag) sampleRadius = rH[sp];
2058  if (checkPoissonDisk(checkSHT, sp->cP(), sampleRadius))
2059  {
2060  ps.AddVert(*sp);
2061  montecarloSHT.RemoveCell(sp);
2062  checkSHT.Add(sp);
2063  break;
2064  }
2065  else
2066  montecarloSHT.RemovePunctual(sp);
2067  }
2068  }
2069  addedCnt = checkSHT.hash_table.size()-addedCnt;
2070  removedCnt = removedCnt-montecarloSHT.hash_table.size();
2071 
2072  // proceed to the next level of subdivision
2073  // increase grid resolution
2074  gridsize *= 2;
2075 
2076  //
2077  level++;
2078  } while(level < 5);
2079 }
2080 
2081 //template <class MeshType>
2082 //void Sampling<MeshType>::SimilarFaceSampling()
2083 
2084 // This function also generates samples outside faces if those affects faces in texture space.
2085 // Use correctSafePointsBaryCoords = true to map safety texels to closest point barycentric coords (on edge)
2086 // otherwise obtained samples will map to barycentric coord actually outside face
2087 //
2088 // If you don't need to get those extra points clear faces Border Flags
2089 // vcg::tri::UpdateFlags<Mesh>::FaceClearB(m);
2090 //
2091 // Else make sure to update border flags from texture space FFadj
2092 // vcg::tri::UpdateTopology<Mesh>::FaceFaceFromTexCoord(m);
2093 // vcg::tri::UpdateFlags<Mesh>::FaceBorderFromFF(m);
2094 static void Texture(MeshType & m, VertexSampler &ps, int textureWidth, int textureHeight, bool correctSafePointsBaryCoords=true)
2095 {
2096 typedef Point2<ScalarType> Point2x;
2097  printf("Similar Triangles face sampling\n");
2098  for(FaceIterator fi=m.face.begin(); fi != m.face.end(); fi++)
2099  if (!fi->IsD())
2100  {
2101  Point2x ti[3];
2102  for(int i=0;i<3;++i)
2103  ti[i]=Point2x((*fi).WT(i).U() * textureWidth - 0.5, (*fi).WT(i).V() * textureHeight - 0.5);
2104  // - 0.5 constants are used to obtain correct texture mapping
2105 
2106  SingleFaceRaster(*fi, ps, ti[0],ti[1],ti[2], correctSafePointsBaryCoords);
2107  }
2108 }
2109 
2110 typedef GridStaticPtr<FaceType, ScalarType > TriMeshGrid;
2111 
2112 class RRParam
2113 {
2114 public:
2115 float offset;
2116 float minDiag;
2117 tri::FaceTmark<MeshType> markerFunctor;
2118 TriMeshGrid gM;
2119 };
2120 
2121 static void RegularRecursiveOffset(MeshType & m, std::vector<CoordType> &pvec, ScalarType offset, float minDiag)
2122 {
2123  Box3<ScalarType> bb=m.bbox;
2124  bb.Offset(offset*2.0);
2125 
2126  RRParam rrp;
2127 
2128  rrp.markerFunctor.SetMesh(&m);
2129 
2130  rrp.gM.Set(m.face.begin(),m.face.end(),bb);
2131 
2132 
2133  rrp.offset=offset;
2134  rrp.minDiag=minDiag;
2135  SubdivideAndSample(m, pvec, bb, rrp, bb.Diag());
2136 }
2137 
2138 static void SubdivideAndSample(MeshType & m, std::vector<CoordType> &pvec, const Box3<ScalarType> bb, RRParam &rrp, float curDiag)
2139 {
2140  CoordType startPt = bb.Center();
2141 
2142  ScalarType dist;
2143  // Compute mesh point nearest to bb center
2144  FaceType *nearestF=0;
2145  ScalarType dist_upper_bound = curDiag+rrp.offset;
2146  CoordType closestPt;
2147  vcg::face::PointDistanceBaseFunctor<ScalarType> PDistFunct;
2148  dist=dist_upper_bound;
2149  nearestF = rrp.gM.GetClosest(PDistFunct,rrp.markerFunctor,startPt,dist_upper_bound,dist,closestPt);
2150  curDiag /=2;
2151  if(dist < dist_upper_bound)
2152  {
2153  if(curDiag/3 < rrp.minDiag) //store points only for the last level of recursion (?)
2154  {
2155  if(rrp.offset==0)
2156  pvec.push_back(closestPt);
2157  else
2158  {
2159  if(dist>rrp.offset) // points below the offset threshold cannot be displaced at the right offset distance, we can only make points nearer.
2160  {
2161  CoordType delta = startPt-closestPt;
2162  pvec.push_back(closestPt+delta*(rrp.offset/dist));
2163  }
2164  }
2165  }
2166  if(curDiag < rrp.minDiag) return;
2167  CoordType hs = (bb.max-bb.min)/2;
2168  for(int i=0;i<2;i++)
2169  for(int j=0;j<2;j++)
2170  for(int k=0;k<2;k++)
2171  SubdivideAndSample(m, pvec,
2172  BoxType(CoordType( bb.min[0]+i*hs[0], bb.min[1]+j*hs[1], bb.min[2]+k*hs[2]),
2173  CoordType(startPt[0]+i*hs[0], startPt[1]+j*hs[1], startPt[2]+k*hs[2]) ),
2174  rrp,curDiag
2175  );
2176 
2177  }
2178 }
2179 }; // end sampling class
2180 
2181 template <class MeshType>
2182 typename MeshType::ScalarType ComputePoissonDiskRadius(MeshType &origMesh, int sampleNum)
2183 {
2184  typedef typename MeshType::ScalarType ScalarType;
2185  ScalarType meshArea = Stat<MeshType>::ComputeMeshArea(origMesh);
2186  // Manage approximately the PointCloud Case, use the half a area of the bbox.
2187  // TODO: If you had the radius a much better approximation could be done.
2188  if(meshArea ==0)
2189  {
2190  meshArea = (origMesh.bbox.DimX()*origMesh.bbox.DimY() +
2191  origMesh.bbox.DimX()*origMesh.bbox.DimZ() +
2192  origMesh.bbox.DimY()*origMesh.bbox.DimZ());
2193  }
2194  ScalarType diskRadius = sqrt(meshArea / (0.7 * M_PI * sampleNum)); // 0.7 is a density factor
2195  return diskRadius;
2196 }
2197 
2198 
2199 
2200 template <class MeshType>
2201 void MontecarloSampling(MeshType &m, // the mesh that has to be sampled
2202  MeshType &mm, // the mesh that will contain the samples
2203  int sampleNum) // the desired number sample, if zero you must set the radius to the wanted value
2204 {
2205  typedef tri::MeshSampler<MeshType> BaseSampler;
2206  MeshSampler<MeshType> mcSampler(&mm);
2208 }
2209 
2210 
2211 template <class MeshType>
2212 void MontecarloSampling(MeshType &m, // the mesh that has to be sampled
2213  std::vector<Point3f> &montercarloSamples, // the vector that will contain the set of points
2214  int sampleNum) // the desired number sample, if zero you must set the radius to the wanted value
2215 {
2216  typedef tri::TrivialSampler<MeshType> BaseSampler;
2217  BaseSampler mcSampler(montercarloSamples);
2219 }
2220 
2221 // Yet another simpler wrapper for the generation of a poisson disk distribution over a mesh.
2222 //
2223 template <class MeshType>
2224 void PoissonSampling(MeshType &m, // the mesh that has to be sampled
2225  std::vector<typename MeshType::CoordType> &poissonSamples, // the vector that will contain the set of points
2226  int sampleNum, // the desired number sample, if zero you must set the radius to the wanted value
2227  typename MeshType::ScalarType &radius, // the Poisson Disk Radius (used if sampleNum==0, setted if sampleNum!=0)
2228  typename MeshType::ScalarType radiusVariance=1,
2229  typename MeshType::ScalarType PruningByNumberTolerance=0.04f,
2230  unsigned int randSeed=0)
2231 
2232 {
2233  typedef tri::TrivialSampler<MeshType> BaseSampler;
2234  typedef tri::MeshSampler<MeshType> MontecarloSampler;
2235 
2236  typename tri::SurfaceSampling<MeshType, BaseSampler>::PoissonDiskParam pp;
2237  int t0=clock();
2238 
2239 // if(sampleNum>0) radius = tri::SurfaceSampling<MeshType,BaseSampler>::ComputePoissonDiskRadius(m,sampleNum);
2240  if(radius>0 && sampleNum==0) sampleNum = tri::SurfaceSampling<MeshType,BaseSampler>::ComputePoissonSampleNum(m,radius);
2241 
2242  pp.pds.sampleNum = sampleNum;
2243  pp.randomSeed = randSeed;
2244  poissonSamples.clear();
2245 // std::vector<Point3f> MontecarloSamples;
2246  MeshType MontecarloMesh;
2247 
2248  // First step build the sampling
2249  MontecarloSampler mcSampler(MontecarloMesh);
2250  BaseSampler pdSampler(poissonSamples);
2251 
2252  if(randSeed) tri::SurfaceSampling<MeshType,MontecarloSampler>::SamplingRandomGenerator().initialize(randSeed);
2253  tri::SurfaceSampling<MeshType,MontecarloSampler>::Montecarlo(m, mcSampler, std::max(10000,sampleNum*40));
2254  tri::UpdateBounding<MeshType>::Box(MontecarloMesh);
2255 // tri::Build(MontecarloMesh, MontecarloSamples);
2256  int t1=clock();
2257  pp.pds.montecarloTime = t1-t0;
2258 
2259  if(radiusVariance !=1)
2260  {
2261  pp.adaptiveRadiusFlag=true;
2262  pp.radiusVariance=radiusVariance;
2263  }
2264  if(sampleNum==0) tri::SurfaceSampling<MeshType,BaseSampler>::PoissonDiskPruning(pdSampler, MontecarloMesh, radius,pp);
2265  else tri::SurfaceSampling<MeshType,BaseSampler>::PoissonDiskPruningByNumber(pdSampler, MontecarloMesh, sampleNum, radius,pp,PruningByNumberTolerance);
2266  int t2=clock();
2267  pp.pds.totalTime = t2-t0;
2268 }
2269 
2273 //
2274 template <class MeshType>
2275 void PoissonPruning(MeshType &m, // the mesh that has to be pruned
2276  std::vector<typename MeshType::VertexPointer> &poissonSamples, // the vector that will contain the chosen set of points
2277  float radius, unsigned int randSeed=0)
2278 {
2279  typedef tri::TrivialPointerSampler<MeshType> BaseSampler;
2281  pp.randomSeed = randSeed;
2282 
2284  BaseSampler pdSampler;
2286  std::swap(pdSampler.sampleVec,poissonSamples);
2287 }
2288 
2289 
2295 template <class MeshType>
2296 void PoissonPruning(MeshType &m, // the mesh that has to be pruned
2297  std::vector<typename MeshType::CoordType> &poissonSamples, // the vector that will contain the chosen set of points
2298  float radius, unsigned int randSeed=0)
2299 {
2300  std::vector<typename MeshType::VertexPointer> poissonSamplesVP;
2301  PoissonPruning(m,poissonSamplesVP,radius,randSeed);
2302  poissonSamples.resize(poissonSamplesVP.size());
2303  for(size_t i=0;i<poissonSamplesVP.size();++i)
2304  poissonSamples[i]=poissonSamplesVP[i]->P();
2305 }
2306 
2307 
2308 
2314 template <class MeshType>
2315 void PoissonPruningExact(MeshType &m,
2316  std::vector<typename MeshType::VertexPointer> &poissonSamples,
2317  typename MeshType::ScalarType & radius,
2318  int sampleNum,
2319  float tolerance=0.04,
2320  int maxIter=20,
2321  unsigned int randSeed=0)
2322 {
2323  size_t sampleNumMin = int(float(sampleNum)*(1.0f-tolerance)); // the expected values range.
2324  size_t sampleNumMax = int(float(sampleNum)*(1.0f+tolerance)); // e.g. any sampling in [sampleNumMin, sampleNumMax] is OK
2325  float RangeMinRad = m.bbox.Diag()/10.0f;
2326  float RangeMaxRad = m.bbox.Diag()/10.0f;
2327  size_t RangeMinSampleNum;
2328  size_t RangeMaxSampleNum;
2329  std::vector<typename MeshType::VertexPointer> poissonSamplesTmp;
2330 
2331  do
2332  {
2333  RangeMinRad/=2.0f;
2334  PoissonPruning(m,poissonSamplesTmp,RangeMinRad,randSeed);
2335  RangeMinSampleNum = poissonSamplesTmp.size();
2336  } while(RangeMinSampleNum < sampleNumMin);
2337 
2338  do
2339  {
2340  RangeMaxRad*=2.0f;
2341  PoissonPruning(m,poissonSamplesTmp,RangeMaxRad,randSeed);
2342  RangeMaxSampleNum = poissonSamplesTmp.size();
2343  } while(RangeMaxSampleNum > sampleNumMax);
2344 
2345  float curRadius;
2346  int iterCnt=0;
2347  while(iterCnt<maxIter &&
2348  (poissonSamplesTmp.size() < sampleNumMin || poissonSamplesTmp.size() > sampleNumMax) )
2349  {
2350  curRadius=(RangeMaxRad+RangeMinRad)/2.0f;
2351  PoissonPruning(m,poissonSamplesTmp,curRadius,randSeed);
2352  //qDebug("(%6.3f:%5i %6.3f:%5i) Cur Radius %f -> %i sample instead of %i",RangeMinRad,RangeMinSampleNum,RangeMaxRad,RangeMaxSampleNum,curRadius,poissonSamplesTmp.size(),sampleNum);
2353  if(poissonSamplesTmp.size() > size_t(sampleNum))
2354  RangeMinRad = curRadius;
2355  if(poissonSamplesTmp.size() < size_t(sampleNum))
2356  RangeMaxRad = curRadius;
2357  }
2358 
2359  swap(poissonSamples,poissonSamplesTmp);
2360  radius = curRadius;
2361 }
2362 } // end namespace tri
2363 } // end namespace vcg
2364 
2365 #endif
2366