OpenVDB  7.0.0
points/PointScatter.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 
13 
14 #ifndef OPENVDB_POINTS_POINT_SCATTER_HAS_BEEN_INCLUDED
15 #define OPENVDB_POINTS_POINT_SCATTER_HAS_BEEN_INCLUDED
16 
17 #include <type_traits>
18 #include <algorithm>
19 #include <thread>
20 #include <random>
21 
22 #include <openvdb/openvdb.h>
23 #include <openvdb/Types.h>
25 #include <openvdb/tools/Prune.h>
27 
28 #include "AttributeArray.h"
29 #include "PointCount.h"
30 #include "PointDataGrid.h"
31 
32 #include <tbb/parallel_sort.h>
33 #include <tbb/parallel_for.h>
34 
35 namespace openvdb {
37 namespace OPENVDB_VERSION_NAME {
38 namespace points {
39 
57 
58 
71 template<
72  typename GridT,
73  typename RandGenT = std::mt19937,
74  typename PositionArrayT = TypedAttributeArray<Vec3f, NullCodec>,
75  typename PointDataGridT = Grid<
76  typename points::TreeConverter<typename GridT::TreeType>::Type>,
77  typename InterrupterT = util::NullInterrupter>
78 inline typename PointDataGridT::Ptr
79 uniformPointScatter(const GridT& grid,
80  const Index64 count,
81  const unsigned int seed = 0,
82  const float spread = 1.0f,
83  InterrupterT* interrupter = nullptr);
84 
99 
100 template<
101  typename GridT,
102  typename RandGenT = std::mt19937,
103  typename PositionArrayT = TypedAttributeArray<Vec3f, NullCodec>,
104  typename PointDataGridT = Grid<
105  typename points::TreeConverter<typename GridT::TreeType>::Type>,
106  typename InterrupterT = util::NullInterrupter>
107 inline typename PointDataGridT::Ptr
108 denseUniformPointScatter(const GridT& grid,
109  const float pointsPerVoxel,
110  const unsigned int seed = 0,
111  const float spread = 1.0f,
112  InterrupterT* interrupter = nullptr);
113 
132 template<
133  typename GridT,
134  typename RandGenT = std::mt19937,
135  typename PositionArrayT = TypedAttributeArray<Vec3f, NullCodec>,
136  typename PointDataGridT = Grid<
137  typename points::TreeConverter<typename GridT::TreeType>::Type>,
138  typename InterrupterT = util::NullInterrupter>
139 inline typename PointDataGridT::Ptr
140 nonUniformPointScatter(const GridT& grid,
141  const float pointsPerVoxel,
142  const unsigned int seed = 0,
143  const float spread = 1.0f,
144  InterrupterT* interrupter = nullptr);
145 
146 
148 
149 
150 namespace point_scatter_internal
151 {
152 
156 template<typename PointDataGridT, typename GridT>
157 inline typename PointDataGridT::Ptr
158 initialisePointTopology(const GridT& grid)
159 {
160  typename PointDataGridT::Ptr points(new PointDataGridT);
161  points->setTransform(grid.transform().copy());
162  points->topologyUnion(grid);
163  if (points->tree().hasActiveTiles()) {
164  points->tree().voxelizeActiveTiles();
165  }
166 
167  return points;
168 }
169 
177 template<typename PositionType,
178  typename CodecT,
179  typename RandGenT,
180  typename LeafNodeT>
181 inline void
182 generatePositions(LeafNodeT& leaf,
183  const AttributeSet::Descriptor::Ptr& descriptor,
184  const Index64& count,
185  const float spread,
186  RandGenT& rand01)
187 {
188  using PositionTraits = VecTraits<PositionType>;
189  using ValueType = typename PositionTraits::ElementType;
190  using PositionWriteHandle = AttributeWriteHandle<PositionType, CodecT>;
191 
192  leaf.initializeAttributes(descriptor, static_cast<Index>(count));
193 
194  // directly expand to avoid needlessly setting uniform values in the
195  // write handle
196  auto& array = leaf.attributeArray(0);
197  array.expand(/*fill*/false);
198 
199  PositionWriteHandle pHandle(array, /*expand*/false);
200  PositionType P;
201  for (Index64 index = 0; index < count; ++index) {
202  P[0] = (spread * (rand01() - ValueType(0.5)));
203  P[1] = (spread * (rand01() - ValueType(0.5)));
204  P[2] = (spread * (rand01() - ValueType(0.5)));
205  pHandle.set(static_cast<Index>(index), P);
206  }
207 }
208 
209 } // namespace point_scatter_internal
210 
211 
213 
214 
215 template<
216  typename GridT,
217  typename RandGenT,
218  typename PositionArrayT,
219  typename PointDataGridT,
220  typename InterrupterT>
221 inline typename PointDataGridT::Ptr
222 uniformPointScatter(const GridT& grid,
223  const Index64 count,
224  const unsigned int seed,
225  const float spread,
226  InterrupterT* interrupter)
227 {
228  using PositionType = typename PositionArrayT::ValueType;
229  using PositionTraits = VecTraits<PositionType>;
230  using ValueType = typename PositionTraits::ElementType;
231  using CodecType = typename PositionArrayT::Codec;
232 
233  using RandomGenerator = math::Rand01<ValueType, RandGenT>;
234 
235  using TreeType = typename PointDataGridT::TreeType;
236  using LeafNodeType = typename TreeType::LeafNodeType;
237  using LeafManagerT = tree::LeafManager<TreeType>;
238 
239  struct Local
240  {
244  static void getPrefixSum(LeafManagerT& leafManager,
245  std::vector<Index64>& offsets)
246  {
247  Index64 offset = 0;
248  offsets.reserve(leafManager.leafCount() + 1);
249  offsets.push_back(0);
250  const auto leafRange = leafManager.leafRange();
251  for (auto leaf = leafRange.begin(); leaf; ++leaf) {
252  offset += leaf->onVoxelCount();
253  offsets.push_back(offset);
254  }
255  }
256  };
257 
258  static_assert(PositionTraits::IsVec && PositionTraits::Size == 3,
259  "Invalid Position Array type.");
260 
261  if (spread < 0.0f || spread > 1.0f) {
262  OPENVDB_THROW(ValueError, "Spread must be between 0 and 1.");
263  }
264 
265  if (interrupter) interrupter->start("Uniform scattering with fixed point count");
266 
267  typename PointDataGridT::Ptr points =
268  point_scatter_internal::initialisePointTopology<PointDataGridT>(grid);
269  TreeType& tree = points->tree();
270  if (!tree.cbeginLeaf()) return points;
271 
272  LeafManagerT leafManager(tree);
273  const Index64 voxelCount = leafManager.activeLeafVoxelCount();
274  assert(voxelCount != 0);
275 
276  const double pointsPerVolume = double(count) / double(voxelCount);
277  const Index32 pointsPerVoxel = static_cast<Index32>(math::RoundDown(pointsPerVolume));
278  const Index64 remainder = count - (pointsPerVoxel * voxelCount);
279 
280  if (remainder == 0) {
282  GridT, RandGenT, PositionArrayT, PointDataGridT, InterrupterT>(
283  grid, float(pointsPerVoxel), seed, spread, interrupter);
284  }
285 
286  std::vector<Index64> voxelOffsets, values;
287  std::thread worker(&Local::getPrefixSum, std::ref(leafManager), std::ref(voxelOffsets));
288 
289  {
290  math::RandInt<Index64, RandGenT> gen(seed, 0, voxelCount-1);
291  values.reserve(remainder);
292  for (Index64 i = 0; i < remainder; ++i) values.emplace_back(gen());
293  }
294 
295  worker.join();
296 
297  if (util::wasInterrupted<InterrupterT>(interrupter)) {
298  tree.clear();
299  return points;
300  }
301 
302  tbb::parallel_sort(values.begin(), values.end());
303  const bool fractionalOnly(pointsPerVoxel == 0);
304 
305  leafManager.foreach([&voxelOffsets, &values, fractionalOnly]
306  (LeafNodeType& leaf, const size_t idx)
307  {
308  const Index64 lowerOffset = voxelOffsets[idx]; // inclusive
309  const Index64 upperOffset = voxelOffsets[idx + 1]; // exclusive
310  assert(upperOffset > lowerOffset);
311 
312  const auto valuesEnd = values.end();
313  auto lower = std::lower_bound(values.begin(), valuesEnd, lowerOffset);
314 
315  auto* const data = leaf.buffer().data();
316  auto iter = leaf.beginValueOn();
317 
318  Index32 currentOffset(0);
319  bool addedPoints(!fractionalOnly);
320  while (lower != valuesEnd) {
321  const Index64 vId = *lower;
322  if (vId >= upperOffset) break;
323 
324  const Index32 nextOffset = Index32(vId - lowerOffset);
325  iter.increment(nextOffset - currentOffset);
326  currentOffset = nextOffset;
327  assert(iter);
328 
329  auto& value = data[iter.pos()];
330  value = value + 1; // no += operator support
331  addedPoints = true;
332  ++lower;
333  }
334 
335  // deactivate this leaf if no points were added. This will speed up
336  // the unthreaded rng
337  if (!addedPoints) leaf.setValuesOff();
338  });
339 
340  voxelOffsets.clear();
341  values.clear();
342 
343  if (fractionalOnly) {
344  tools::pruneInactive(tree);
345  leafManager.rebuild();
346  }
347 
348  const AttributeSet::Descriptor::Ptr descriptor =
349  AttributeSet::Descriptor::create(PositionArrayT::attributeType());
350  RandomGenerator rand01(seed);
351 
352  const auto leafRange = leafManager.leafRange();
353  auto leaf = leafRange.begin();
354  for (; leaf; ++leaf) {
355  if (util::wasInterrupted<InterrupterT>(interrupter)) break;
356  Index32 offset(0);
357  for (auto iter = leaf->beginValueAll(); iter; ++iter) {
358  if (iter.isValueOn()) {
359  const Index32 value = Index32(pointsPerVolume + Index32(*iter));
360  if (value == 0) leaf->setValueOff(iter.pos());
361  else offset += value;
362  }
363  // @note can't use iter.setValue(offset) on point grids
364  leaf->setOffsetOnly(iter.pos(), offset);
365  }
366 
367  // offset should always be non zero
368  assert(offset != 0);
369  point_scatter_internal::generatePositions<PositionType, CodecType>
370  (*leaf, descriptor, offset, spread, rand01);
371  }
372 
373  // if interrupted, remove remaining leaf nodes
374  if (leaf) {
375  for (; leaf; ++leaf) leaf->setValuesOff();
376  tools::pruneInactive(tree);
377  }
378 
379  if (interrupter) interrupter->end();
380  return points;
381 }
382 
383 
385 
386 
387 template<
388  typename GridT,
389  typename RandGenT,
390  typename PositionArrayT,
391  typename PointDataGridT,
392  typename InterrupterT>
393 inline typename PointDataGridT::Ptr
394 denseUniformPointScatter(const GridT& grid,
395  const float pointsPerVoxel,
396  const unsigned int seed,
397  const float spread,
398  InterrupterT* interrupter)
399 {
400  using PositionType = typename PositionArrayT::ValueType;
401  using PositionTraits = VecTraits<PositionType>;
402  using ValueType = typename PositionTraits::ElementType;
403  using CodecType = typename PositionArrayT::Codec;
404 
405  using RandomGenerator = math::Rand01<ValueType, RandGenT>;
406 
407  using TreeType = typename PointDataGridT::TreeType;
408 
409  static_assert(PositionTraits::IsVec && PositionTraits::Size == 3,
410  "Invalid Position Array type.");
411 
412  if (pointsPerVoxel < 0.0f) {
413  OPENVDB_THROW(ValueError, "Points per voxel must not be less than zero.");
414  }
415 
416  if (spread < 0.0f || spread > 1.0f) {
417  OPENVDB_THROW(ValueError, "Spread must be between 0 and 1.");
418  }
419 
420  if (interrupter) interrupter->start("Dense uniform scattering with fixed point count");
421 
422  typename PointDataGridT::Ptr points =
423  point_scatter_internal::initialisePointTopology<PointDataGridT>(grid);
424  TreeType& tree = points->tree();
425  auto leafIter = tree.beginLeaf();
426  if (!leafIter) return points;
427 
428  const Index32 pointsPerVoxelInt = math::Floor(pointsPerVoxel);
429  const double delta = pointsPerVoxel - float(pointsPerVoxelInt);
430  const bool fractional = !math::isApproxZero(delta, 1.0e-6);
431  const bool fractionalOnly = pointsPerVoxelInt == 0;
432 
433  const AttributeSet::Descriptor::Ptr descriptor =
434  AttributeSet::Descriptor::create(PositionArrayT::attributeType());
435  RandomGenerator rand01(seed);
436 
437  for (; leafIter; ++leafIter) {
438  if (util::wasInterrupted<InterrupterT>(interrupter)) break;
439  Index32 offset(0);
440  for (auto iter = leafIter->beginValueAll(); iter; ++iter) {
441  if (iter.isValueOn()) {
442  offset += pointsPerVoxelInt;
443  if (fractional && rand01() < delta) ++offset;
444  else if (fractionalOnly) leafIter->setValueOff(iter.pos());
445  }
446  // @note can't use iter.setValue(offset) on point grids
447  leafIter->setOffsetOnly(iter.pos(), offset);
448  }
449 
450  if (offset != 0) {
451  point_scatter_internal::generatePositions<PositionType, CodecType>
452  (*leafIter, descriptor, offset, spread, rand01);
453  }
454  }
455 
456  // if interrupted, remove remaining leaf nodes
457  const bool prune(leafIter || fractionalOnly);
458  for (; leafIter; ++leafIter) leafIter->setValuesOff();
459 
460  if (prune) tools::pruneInactive(tree);
461  if (interrupter) interrupter->end();
462  return points;
463 }
464 
465 
467 
468 
469 template<
470  typename GridT,
471  typename RandGenT,
472  typename PositionArrayT,
473  typename PointDataGridT,
474  typename InterrupterT>
475 inline typename PointDataGridT::Ptr
476 nonUniformPointScatter(const GridT& grid,
477  const float pointsPerVoxel,
478  const unsigned int seed,
479  const float spread,
480  InterrupterT* interrupter)
481 {
482  using PositionType = typename PositionArrayT::ValueType;
483  using PositionTraits = VecTraits<PositionType>;
484  using ValueType = typename PositionTraits::ElementType;
485  using CodecType = typename PositionArrayT::Codec;
486 
487  using RandomGenerator = math::Rand01<ValueType, RandGenT>;
488 
489  using TreeType = typename PointDataGridT::TreeType;
490 
491  static_assert(PositionTraits::IsVec && PositionTraits::Size == 3,
492  "Invalid Position Array type.");
493  static_assert(std::is_arithmetic<typename GridT::ValueType>::value,
494  "Scalar grid type required for weighted voxel scattering.");
495 
496  if (pointsPerVoxel < 0.0f) {
497  OPENVDB_THROW(ValueError, "Points per voxel must not be less than zero.");
498  }
499 
500  if (spread < 0.0f || spread > 1.0f) {
501  OPENVDB_THROW(ValueError, "Spread must be between 0 and 1.");
502  }
503 
504  if (interrupter) interrupter->start("Non-uniform scattering with local point density");
505 
506  typename PointDataGridT::Ptr points =
507  point_scatter_internal::initialisePointTopology<PointDataGridT>(grid);
508  TreeType& tree = points->tree();
509  auto leafIter = tree.beginLeaf();
510  if (!leafIter) return points;
511 
512  const AttributeSet::Descriptor::Ptr descriptor =
513  AttributeSet::Descriptor::create(PositionArrayT::attributeType());
514  RandomGenerator rand01(seed);
515  const auto accessor = grid.getConstAccessor();
516 
517  for (; leafIter; ++leafIter) {
518  if (util::wasInterrupted<InterrupterT>(interrupter)) break;
519  Index32 offset(0);
520  for (auto iter = leafIter->beginValueAll(); iter; ++iter) {
521  if (iter.isValueOn()) {
522  double fractional =
523  double(accessor.getValue(iter.getCoord())) * pointsPerVoxel;
524  fractional = std::max(0.0, fractional);
525  int count = int(fractional);
526  if (rand01() < (fractional - double(count))) ++count;
527  else if (count == 0) leafIter->setValueOff(iter.pos());
528  offset += count;
529  }
530  // @note can't use iter.setValue(offset) on point grids
531  leafIter->setOffsetOnly(iter.pos(), offset);
532  }
533 
534  if (offset != 0) {
535  point_scatter_internal::generatePositions<PositionType, CodecType>
536  (*leafIter, descriptor, offset, spread, rand01);
537  }
538  }
539 
540  // if interrupted, remove remaining leaf nodes
541  for (; leafIter; ++leafIter) leafIter->setValuesOff();
542 
543  tools::pruneInactive(points->tree());
544  if (interrupter) interrupter->end();
545  return points;
546 }
547 
548 
549 } // namespace points
550 } // namespace OPENVDB_VERSION_NAME
551 } // namespace openvdb
552 
553 
554 #endif // OPENVDB_POINTS_POINT_SCATTER_HAS_BEEN_INCLUDED
PointDataGridT::Ptr denseUniformPointScatter(const GridT &grid, const float pointsPerVoxel, const unsigned int seed=0, const float spread=1.0f, InterrupterT *interrupter=nullptr)
Uniformly scatter a fixed number of points per active voxel. If the pointsPerVoxel value provided is ...
Definition: points/PointScatter.h:394
void pruneInactive(TreeT &tree, bool threaded=true, size_t grainSize=1)
Reduce the memory footprint of a tree by replacing with background tiles any nodes whose values are a...
Definition: Prune.h:354
Simple generator of random numbers over the range [0, 1)
Definition: Math.h:108
uint32_t Index32
Definition: Types.h:29
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:82
Attribute Array storage templated on type and compression codec.
Methods for counting points in VDB Point grids.
void generatePositions(LeafNodeT &leaf, const AttributeSet::Descriptor::Ptr &descriptor, const Index64 &count, const float spread, RandGenT &rand01)
Generate random point positions for a leaf node.
Definition: points/PointScatter.h:182
float RoundDown(float x)
Return x rounded down to the nearest integer.
Definition: Math.h:750
void expand(bool fill=true)
If this array is uniform, replace it with an array of length size().
Definition: AttributeArray.h:2278
Defined various multi-threaded utility functions for trees.
Definition: Exceptions.h:65
openvdb::GridBase Grid
Definition: Utils.h:33
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:102
PointDataGridT::Ptr nonUniformPointScatter(const GridT &grid, const float pointsPerVoxel, const unsigned int seed=0, const float spread=1.0f, InterrupterT *interrupter=nullptr)
Non uniformly scatter points per active voxel. The pointsPerVoxel value is used to weight each grids ...
Definition: points/PointScatter.h:476
Definition: Types.h:181
Definition: Exceptions.h:13
int Floor(float x)
Return the floor of x.
Definition: Math.h:795
Write-able version of AttributeHandle.
Definition: AttributeArray.h:920
This class manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional auxil...
Definition: LeafManager.h:82
LeafRange leafRange(size_t grainsize=1) const
Return a TBB-compatible LeafRange.
Definition: LeafManager.h:358
PointDataGridT::Ptr initialisePointTopology(const GridT &grid)
initialise the topology of a PointDataGrid and ensure everything is voxelized
Definition: points/PointScatter.h:158
PointDataGridT::Ptr uniformPointScatter(const GridT &grid, const Index64 count, const unsigned int seed=0, const float spread=1.0f, InterrupterT *interrupter=nullptr)
The free functions depend on the following class:
Definition: points/PointScatter.h:222
Attribute-owned data structure for points. Point attributes are stored in leaf nodes and ordered by v...
bool isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance...
Definition: Math.h:293
uint64_t Index64
Definition: Types.h:30
Simple random integer generator.
Definition: Math.h:144
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:154
void prune(TreeT &tree, typename TreeT::ValueType tolerance=zeroVal< typename TreeT::ValueType >(), bool threaded=true, size_t grainSize=1)
Reduce the memory footprint of a tree by replacing with tiles any nodes whose values are all the same...
Definition: Prune.h:334
A LeafManager manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional au...
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:106