OpenVDB  9.0.1
PointScatter.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 //
4 /// @author Ken Museth
5 ///
6 /// @file tools/PointScatter.h
7 ///
8 /// @brief We offer three different algorithms (each in its own class)
9 /// for scattering of points in active voxels:
10 ///
11 /// 1) UniformPointScatter. Has two modes: Either randomly distributes
12 /// a fixed number of points into the active voxels, or the user can
13 /// specify a fixed probability of having a points per unit of volume.
14 ///
15 /// 2) DenseUniformPointScatter. Randomly distributes points into active
16 /// voxels using a fixed number of points per voxel.
17 ///
18 /// 3) NonIniformPointScatter. Define the local probability of having
19 /// a point in a voxel as the product of a global density and the
20 /// value of the voxel itself.
21 
22 #ifndef OPENVDB_TOOLS_POINT_SCATTER_HAS_BEEN_INCLUDED
23 #define OPENVDB_TOOLS_POINT_SCATTER_HAS_BEEN_INCLUDED
24 
25 #include <openvdb/Types.h>
26 #include <openvdb/Grid.h>
27 #include <openvdb/math/Math.h>
29 #include <tbb/parallel_sort.h>
30 #include <tbb/parallel_for.h>
31 #include <iostream>
32 #include <memory>
33 #include <string>
34 
35 namespace openvdb {
37 namespace OPENVDB_VERSION_NAME {
38 namespace tools {
39 
40 /// Forward declaration of base class
41 template<typename PointAccessorType,
42  typename RandomGenerator,
43  typename InterruptType = util::NullInterrupter>
45 
46 /// @brief The two point scatters UniformPointScatter and
47 /// NonUniformPointScatter depend on the following two classes:
48 ///
49 /// The @c PointAccessorType template argument below refers to any class
50 /// with the following interface:
51 /// @code
52 /// class PointAccessor {
53 /// ...
54 /// public:
55 /// void add(const openvdb::Vec3R &pos);// appends point with world positions pos
56 /// };
57 /// @endcode
58 ///
59 ///
60 /// The @c InterruptType template argument below refers to any class
61 /// with the following interface:
62 /// @code
63 /// class Interrupter {
64 /// ...
65 /// public:
66 /// void start(const char* name = nullptr) // called when computations begin
67 /// void end() // called when computations end
68 /// bool wasInterrupted(int percent=-1) // return true to break computation
69 ///};
70 /// @endcode
71 ///
72 /// @note If no template argument is provided for this InterruptType
73 /// the util::NullInterrupter is used which implies that all
74 /// interrupter calls are no-ops (i.e. incurs no computational overhead).
75 
76 
77 /// @brief Uniformly scatters points in the active voxels.
78 /// The point count is either explicitly defined or implicitly
79 /// through the specification of a global density (=points-per-volume)
80 ///
81 /// @note This uniform scattering technique assumes that the number of
82 /// points is generally smaller than the number of active voxels
83 /// (including virtual active voxels in active tiles).
84 template<typename PointAccessorType,
85  typename RandomGenerator,
86  typename InterruptType = util::NullInterrupter>
87 class UniformPointScatter : public BasePointScatter<PointAccessorType,
88  RandomGenerator,
89  InterruptType>
90 {
91 public:
93 
94  UniformPointScatter(PointAccessorType& points,
96  RandomGenerator& randGen,
97  double spread = 1.0,
98  InterruptType* interrupt = nullptr)
99  : BaseT(points, randGen, spread, interrupt)
100  , mTargetPointCount(pointCount)
101  , mPointsPerVolume(0.0f)
102  {
103  }
104  UniformPointScatter(PointAccessorType& points,
105  float pointsPerVolume,
106  RandomGenerator& randGen,
107  double spread = 1.0,
108  InterruptType* interrupt = nullptr)
109  : BaseT(points, randGen, spread, interrupt)
110  , mTargetPointCount(0)
111  , mPointsPerVolume(pointsPerVolume)
112  {
113  }
114 
115  /// This is the main functor method implementing the actual scattering of points.
116  template<typename GridT>
117  bool operator()(const GridT& grid)
118  {
119  mVoxelCount = grid.activeVoxelCount();
120  if (mVoxelCount == 0) return false;
121 
122  const auto voxelVolume = grid.transform().voxelVolume();
123  if (mPointsPerVolume > 0) {
124  BaseT::start("Uniform scattering with fixed point density");
125  mTargetPointCount = Index64(mPointsPerVolume * voxelVolume * double(mVoxelCount));
126  } else if (mTargetPointCount > 0) {
127  BaseT::start("Uniform scattering with fixed point count");
128  mPointsPerVolume = float(mTargetPointCount) / float(voxelVolume * double(mVoxelCount));
129  } else {
130  return false;
131  }
132 
133  std::unique_ptr<Index64[]> idList{new Index64[mTargetPointCount]};
134  math::RandInt<Index64, RandomGenerator> rand(BaseT::mRand01.engine(), 0, mVoxelCount-1);
135  for (Index64 i=0; i<mTargetPointCount; ++i) idList[i] = rand();
136  tbb::parallel_sort(idList.get(), idList.get() + mTargetPointCount);
137 
138  CoordBBox bbox;
139  const Vec3R offset(0.5, 0.5, 0.5);
140  typename GridT::ValueOnCIter valueIter = grid.cbeginValueOn();
141 
142  for (Index64 i=0, n=valueIter.getVoxelCount() ; i != mTargetPointCount; ++i) {
143  if (BaseT::interrupt()) return false;
144  const Index64 voxelId = idList[i];
145  while ( n <= voxelId ) {
146  ++valueIter;
147  n += valueIter.getVoxelCount();
148  }
149  if (valueIter.isVoxelValue()) {// a majority is expected to be voxels
150  BaseT::addPoint(grid, valueIter.getCoord() - offset);
151  } else {// tiles contain multiple (virtual) voxels
152  valueIter.getBoundingBox(bbox);
153  BaseT::addPoint(grid, bbox.min() - offset, bbox.extents());
154  }
155  }//loop over all the active voxels and tiles
156  //}
157 
158  BaseT::end();
159  return true;
160  }
161 
162  // The following methods should only be called after the
163  // the operator() method was called
164  void print(const std::string &name, std::ostream& os = std::cout) const
165  {
166  os << "Uniformly scattered " << mPointCount << " points into " << mVoxelCount
167  << " active voxels in \"" << name << "\" corresponding to "
168  << mPointsPerVolume << " points per volume." << std::endl;
169  }
170 
171  float getPointsPerVolume() const { return mPointsPerVolume; }
172  Index64 getTargetPointCount() const { return mTargetPointCount; }
173 
174 private:
175 
176  using BaseT::mPointCount;
177  using BaseT::mVoxelCount;
178  Index64 mTargetPointCount;
179  float mPointsPerVolume;
180 
181 }; // class UniformPointScatter
182 
183 /// @brief Scatters a fixed (and integer) number of points in all
184 /// active voxels and tiles.
185 template<typename PointAccessorType,
186  typename RandomGenerator,
187  typename InterruptType = util::NullInterrupter>
188 class DenseUniformPointScatter : public BasePointScatter<PointAccessorType,
189  RandomGenerator,
190  InterruptType>
191 {
192 public:
194 
195  DenseUniformPointScatter(PointAccessorType& points,
196  float pointsPerVoxel,
197  RandomGenerator& randGen,
198  double spread = 1.0,
199  InterruptType* interrupt = nullptr)
200  : BaseT(points, randGen, spread, interrupt)
201  , mPointsPerVoxel(pointsPerVoxel)
202  {
203  }
204 
205  /// This is the main functor method implementing the actual scattering of points.
206  template<typename GridT>
207  bool operator()(const GridT& grid)
208  {
209  using ValueIter = typename GridT::ValueOnCIter;
210  if (mPointsPerVoxel < 1.0e-6) return false;
211  mVoxelCount = grid.activeVoxelCount();
212  if (mVoxelCount == 0) return false;
213  BaseT::start("Dense uniform scattering with fixed point count");
214  CoordBBox bbox;
215  const Vec3R offset(0.5, 0.5, 0.5);
216 
217  const int ppv = math::Floor(mPointsPerVoxel);
218  const double delta = mPointsPerVoxel - float(ppv);
219  const bool fractional = !math::isApproxZero(delta, 1.0e-6);
220 
221  for (ValueIter iter = grid.cbeginValueOn(); iter; ++iter) {
222  if (BaseT::interrupt()) return false;
223  if (iter.isVoxelValue()) {// a majority is expected to be voxels
224  const Vec3R dmin = iter.getCoord() - offset;
225  for (int n = 0; n != ppv; ++n) BaseT::addPoint(grid, dmin);
226  if (fractional && BaseT::getRand01() < delta) BaseT::addPoint(grid, dmin);
227  } else {// tiles contain multiple (virtual) voxels
228  iter.getBoundingBox(bbox);
229  const Coord size(bbox.extents());
230  const Vec3R dmin = bbox.min() - offset;
231  const double d = mPointsPerVoxel * float(iter.getVoxelCount());
232  const int m = math::Floor(d);
233  for (int n = 0; n != m; ++n) BaseT::addPoint(grid, dmin, size);
234  if (BaseT::getRand01() < d - m) BaseT::addPoint(grid, dmin, size);
235  }
236  }//loop over all the active voxels and tiles
237  //}
238  BaseT::end();
239  return true;
240  }
241 
242  // The following methods should only be called after the
243  // the operator() method was called
244  void print(const std::string &name, std::ostream& os = std::cout) const
245  {
246  os << "Dense uniformly scattered " << mPointCount << " points into " << mVoxelCount
247  << " active voxels in \"" << name << "\" corresponding to "
248  << mPointsPerVoxel << " points per voxel." << std::endl;
249  }
250 
251  float getPointsPerVoxel() const { return mPointsPerVoxel; }
252 
253 private:
254  using BaseT::mPointCount;
255  using BaseT::mVoxelCount;
256  float mPointsPerVoxel;
257 }; // class DenseUniformPointScatter
258 
259 /// @brief Non-uniform scatters of point in the active voxels.
260 /// The local point count is implicitly defined as a product of
261 /// of a global density (called pointsPerVolume) and the local voxel
262 /// (or tile) value.
263 ///
264 /// @note This scattering technique can be significantly slower
265 /// than a uniform scattering since its computational complexity
266 /// is proportional to the active voxel (and tile) count.
267 template<typename PointAccessorType,
268  typename RandomGenerator,
269  typename InterruptType = util::NullInterrupter>
270 class NonUniformPointScatter : public BasePointScatter<PointAccessorType,
271  RandomGenerator,
272  InterruptType>
273 {
274 public:
276 
277  NonUniformPointScatter(PointAccessorType& points,
278  float pointsPerVolume,
279  RandomGenerator& randGen,
280  double spread = 1.0,
281  InterruptType* interrupt = nullptr)
282  : BaseT(points, randGen, spread, interrupt)
283  , mPointsPerVolume(pointsPerVolume)//note this is merely a
284  //multiplier for the local point density
285  {
286  }
287 
288  /// This is the main functor method implementing the actual scattering of points.
289  template<typename GridT>
290  bool operator()(const GridT& grid)
291  {
292  if (mPointsPerVolume <= 0.0f) return false;
293  mVoxelCount = grid.activeVoxelCount();
294  if (mVoxelCount == 0) return false;
295  BaseT::start("Non-uniform scattering with local point density");
296  const Vec3d dim = grid.voxelSize();
297  const double volumePerVoxel = dim[0]*dim[1]*dim[2],
298  pointsPerVoxel = mPointsPerVolume * volumePerVoxel;
299  CoordBBox bbox;
300  const Vec3R offset(0.5, 0.5, 0.5);
301  for (typename GridT::ValueOnCIter iter = grid.cbeginValueOn(); iter; ++iter) {
302  if (BaseT::interrupt()) return false;
303  const double d = double(*iter) * pointsPerVoxel * double(iter.getVoxelCount());
304  const int n = int(d);
305  if (iter.isVoxelValue()) { // a majority is expected to be voxels
306  const Vec3R dmin =iter.getCoord() - offset;
307  for (int i = 0; i < n; ++i) BaseT::addPoint(grid, dmin);
308  if (BaseT::getRand01() < (d - n)) BaseT::addPoint(grid, dmin);
309  } else { // tiles contain multiple (virtual) voxels
310  iter.getBoundingBox(bbox);
311  const Coord size(bbox.extents());
312  const Vec3R dmin = bbox.min() - offset;
313  for (int i = 0; i < n; ++i) BaseT::addPoint(grid, dmin, size);
314  if (BaseT::getRand01() < (d - n)) BaseT::addPoint(grid, dmin, size);
315  }
316  }//loop over all the active voxels and tiles
317  BaseT::end();
318  return true;
319  }
320 
321  // The following methods should only be called after the
322  // the operator() method was called
323  void print(const std::string &name, std::ostream& os = std::cout) const
324  {
325  os << "Non-uniformly scattered " << mPointCount << " points into " << mVoxelCount
326  << " active voxels in \"" << name << "\"." << std::endl;
327  }
328 
329  float getPointPerVolume() const { return mPointsPerVolume; }
330 
331 private:
332  using BaseT::mPointCount;
333  using BaseT::mVoxelCount;
334  float mPointsPerVolume;
335 
336 }; // class NonUniformPointScatter
337 
338 /// Base class of all the point scattering classes defined above
339 template<typename PointAccessorType,
340  typename RandomGenerator,
341  typename InterruptType>
342 class BasePointScatter
343 {
344 public:
345 
346  Index64 getPointCount() const { return mPointCount; }
347  Index64 getVoxelCount() const { return mVoxelCount; }
348 
349 protected:
350 
351  PointAccessorType& mPoints;
352  InterruptType* mInterrupter;
356  const double mSpread;
358 
359  /// This is a base class so the constructor is protected
360  BasePointScatter(PointAccessorType& points,
361  RandomGenerator& randGen,
362  double spread,
363  InterruptType* interrupt = nullptr)
364  : mPoints(points)
365  , mInterrupter(interrupt)
366  , mPointCount(0)
367  , mVoxelCount(0)
368  , mInterruptCount(0)
369  , mSpread(math::Clamp01(spread))
370  , mRand01(randGen)
371  {
372  }
373 
374  inline void start(const char* name)
375  {
376  if (mInterrupter) mInterrupter->start(name);
377  }
378 
379  inline void end()
380  {
381  if (mInterrupter) mInterrupter->end();
382  }
383 
384  inline bool interrupt()
385  {
386  //only check interrupter for every 32'th call
387  return !(mInterruptCount++ & ((1<<5)-1)) && util::wasInterrupted(mInterrupter);
388  }
389 
390  /// @brief Return a random floating point number between zero and one
391  inline double getRand01() { return mRand01(); }
392 
393  /// @brief Return a random floating point number between 0.5 -+ mSpread/2
394  inline double getRand() { return 0.5 + mSpread * (mRand01() - 0.5); }
395 
396  template <typename GridT>
397  inline void addPoint(const GridT &grid, const Vec3R &dmin)
398  {
399  const Vec3R pos(dmin[0] + this->getRand(),
400  dmin[1] + this->getRand(),
401  dmin[2] + this->getRand());
402  mPoints.add(grid.indexToWorld(pos));
403  ++mPointCount;
404  }
405 
406  template <typename GridT>
407  inline void addPoint(const GridT &grid, const Vec3R &dmin, const Coord &size)
408  {
409  const Vec3R pos(dmin[0] + size[0]*this->getRand(),
410  dmin[1] + size[1]*this->getRand(),
411  dmin[2] + size[2]*this->getRand());
412  mPoints.add(grid.indexToWorld(pos));
413  ++mPointCount;
414  }
415 };// class BasePointScatter
416 
417 } // namespace tools
418 } // namespace OPENVDB_VERSION_NAME
419 } // namespace openvdb
420 
421 #endif // OPENVDB_TOOLS_POINT_SCATTER_HAS_BEEN_INCLUDED
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:49
Simple random integer generator.
Definition: Math.h:201
Index64 pointCount(const PointDataTreeT &tree, const FilterT &filter=NullFilter(), const bool inCoreOnly=false, const bool threaded=true)
Count the total number of points in a PointDataTree.
Definition: PointCount.h:88
void print(const std::string &name, std::ostream &os=std::cout) const
Definition: PointScatter.h:244
Type Clamp01(Type x)
Return x clamped to [0, 1].
Definition: Math.h:270
Index64 getTargetPointCount() const
Definition: PointScatter.h:172
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
float getPointsPerVoxel() const
Definition: PointScatter.h:251
bool interrupt()
Definition: PointScatter.h:384
Base class for interrupters.
Definition: NullInterrupter.h:25
Index64 mVoxelCount
Definition: PointScatter.h:354
void print(const std::string &name, std::ostream &os=std::cout) const
Definition: PointScatter.h:323
Index64 mInterruptCount
Definition: PointScatter.h:355
bool isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance...
Definition: Math.h:350
void print(const std::string &name, std::ostream &os=std::cout) const
Definition: PointScatter.h:164
double getRand01()
Return a random floating point number between zero and one.
Definition: PointScatter.h:391
int Floor(float x)
Return the floor of x.
Definition: Math.h:851
bool operator()(const GridT &grid)
This is the main functor method implementing the actual scattering of points.
Definition: PointScatter.h:290
Scatters a fixed (and integer) number of points in all active voxels and tiles.
Definition: PointScatter.h:188
BBox< Coord > CoordBBox
Definition: NanoVDB.h:1658
void end()
Definition: PointScatter.h:379
InterruptType * mInterrupter
Definition: PointScatter.h:352
Index64 getVoxelCount() const
Definition: PointScatter.h:347
The two point scatters UniformPointScatter and NonUniformPointScatter depend on the following two cla...
Definition: PointScatter.h:87
double getRand()
Return a random floating point number between 0.5 -+ mSpread/2.
Definition: PointScatter.h:394
bool operator()(const GridT &grid)
This is the main functor method implementing the actual scattering of points.
Definition: PointScatter.h:117
bool operator()(const GridT &grid)
This is the main functor method implementing the actual scattering of points.
Definition: PointScatter.h:207
Non-uniform scatters of point in the active voxels. The local point count is implicitly defined as a ...
Definition: PointScatter.h:270
UniformPointScatter(PointAccessorType &points, float pointsPerVolume, RandomGenerator &randGen, double spread=1.0, InterruptType *interrupt=nullptr)
Definition: PointScatter.h:104
uint64_t Index64
Definition: Types.h:53
float getPointsPerVolume() const
Definition: PointScatter.h:171
Definition: Exceptions.h:13
PointAccessorType & mPoints
Definition: PointScatter.h:351
BasePointScatter(PointAccessorType &points, RandomGenerator &randGen, double spread, InterruptType *interrupt=nullptr)
This is a base class so the constructor is protected.
Definition: PointScatter.h:360
Forward declaration of base class.
Definition: PointScatter.h:44
void addPoint(const GridT &grid, const Vec3R &dmin, const Coord &size)
Definition: PointScatter.h:407
float getPointPerVolume() const
Definition: PointScatter.h:329
void addPoint(const GridT &grid, const Vec3R &dmin)
Definition: PointScatter.h:397
void start(const char *name)
Definition: PointScatter.h:374
const double mSpread
Definition: PointScatter.h:356
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:116
Index64 getPointCount() const
Definition: PointScatter.h:346
NonUniformPointScatter(PointAccessorType &points, float pointsPerVolume, RandomGenerator &randGen, double spread=1.0, InterruptType *interrupt=nullptr)
Definition: PointScatter.h:277
math::Rand01< double, RandomGenerator > mRand01
Definition: PointScatter.h:357
Index64 mPointCount
Definition: PointScatter.h:353
DenseUniformPointScatter(PointAccessorType &points, float pointsPerVoxel, RandomGenerator &randGen, double spread=1.0, InterruptType *interrupt=nullptr)
Definition: PointScatter.h:195
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:202
UniformPointScatter(PointAccessorType &points, Index64 pointCount, RandomGenerator &randGen, double spread=1.0, InterruptType *interrupt=nullptr)
Definition: PointScatter.h:94