OpenVDB  12.0.0
PointScatterImpl.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: Apache-2.0
3 
4 /// @author Nick Avramoussis
5 ///
6 /// @file PointScatterImpl.h
7 ///
8 
9 #ifndef OPENVDB_POINTS_POINT_SCATTER_IMPL_HAS_BEEN_INCLUDED
10 #define OPENVDB_POINTS_POINT_SCATTER_IMPL_HAS_BEEN_INCLUDED
11 
12 namespace openvdb {
14 namespace OPENVDB_VERSION_NAME {
15 namespace points {
16 
17 /// @cond OPENVDB_DOCS_INTERNAL
18 
19 namespace point_scatter_internal
20 {
21 
22 /// @brief initialise the topology of a PointDataGrid and ensure
23 /// everything is voxelized
24 /// @param grid The source grid from which to base the topology generation
25 template<typename PointDataGridT, typename GridT>
26 inline typename PointDataGridT::Ptr
27 initialisePointTopology(const GridT& grid)
28 {
29  typename PointDataGridT::Ptr points(new PointDataGridT);
30  points->setTransform(grid.transform().copy());
31  points->topologyUnion(grid);
32  if (points->tree().hasActiveTiles()) {
33  points->tree().voxelizeActiveTiles();
34  }
35 
36  return points;
37 }
38 
39 /// @brief Generate random point positions for a leaf node
40 /// @param leaf The leaf node to initialize
41 /// @param descriptor The descriptor containing the position type
42 /// @param count The number of points to generate
43 /// @param spread The spread of points from the voxel center
44 /// @param rand01 The random number generator, expected to produce floating point
45 /// values between 0 and 1.
46 template<typename PositionType,
47  typename CodecT,
48  typename RandGenT,
49  typename LeafNodeT>
50 inline void
51 generatePositions(LeafNodeT& leaf,
52  const AttributeSet::Descriptor::Ptr& descriptor,
53  const Index64& count,
54  const float spread,
55  RandGenT& rand01)
56 {
57  using PositionTraits = VecTraits<PositionType>;
58  using ValueType = typename PositionTraits::ElementType;
59  using PositionWriteHandle = AttributeWriteHandle<PositionType, CodecT>;
60 
61  leaf.initializeAttributes(descriptor, static_cast<Index>(count));
62 
63  // directly expand to avoid needlessly setting uniform values in the
64  // write handle
65  auto& array = leaf.attributeArray(0);
66  array.expand(/*fill*/false);
67 
68  PositionWriteHandle pHandle(array, /*expand*/false);
69  PositionType P;
70  for (Index64 index = 0; index < count; ++index) {
71  P[0] = (spread * (rand01() - ValueType(0.5)));
72  P[1] = (spread * (rand01() - ValueType(0.5)));
73  P[2] = (spread * (rand01() - ValueType(0.5)));
74  pHandle.set(static_cast<Index>(index), P);
75  }
76 }
77 
78 } // namespace point_scatter_internal
79 
80 /// @endcond
81 
82 ////////////////////////////////////////
83 
84 
85 template<
86  typename GridT,
87  typename RandGenT,
88  typename PositionArrayT,
89  typename PointDataGridT,
90  typename InterrupterT>
91 inline typename PointDataGridT::Ptr
92 uniformPointScatter(const GridT& grid,
93  const Index64 count,
94  const unsigned int seed,
95  const float spread,
96  InterrupterT* interrupter)
97 {
98  using PositionType = typename PositionArrayT::ValueType;
99  using PositionTraits = VecTraits<PositionType>;
100  using ValueType = typename PositionTraits::ElementType;
101  using CodecType = typename PositionArrayT::Codec;
102 
103  using RandomGenerator = math::Rand01<ValueType, RandGenT>;
104 
105  using TreeType = typename PointDataGridT::TreeType;
106  using LeafNodeType = typename TreeType::LeafNodeType;
107  using LeafManagerT = tree::LeafManager<TreeType>;
108 
109  struct Local
110  {
111  /// @brief Get the prefixed voxel counts for each leaf node with an
112  /// additional value to represent the end voxel count.
113  /// See also LeafManager::getPrefixSum()
114  static void getPrefixSum(LeafManagerT& leafManager,
115  std::vector<Index64>& offsets)
116  {
117  Index64 offset = 0;
118  offsets.reserve(leafManager.leafCount() + 1);
119  offsets.push_back(0);
120  const auto leafRange = leafManager.leafRange();
121  for (auto leaf = leafRange.begin(); leaf; ++leaf) {
122  offset += leaf->onVoxelCount();
123  offsets.push_back(offset);
124  }
125  }
126  };
127 
128  static_assert(PositionTraits::IsVec && PositionTraits::Size == 3,
129  "Invalid Position Array type.");
130 
131  if (spread < 0.0f || spread > 1.0f) {
132  OPENVDB_THROW(ValueError, "Spread must be between 0 and 1.");
133  }
134 
135  if (interrupter) interrupter->start("Uniform scattering with fixed point count");
136 
137  typename PointDataGridT::Ptr points =
138  point_scatter_internal::initialisePointTopology<PointDataGridT>(grid);
139  TreeType& tree = points->tree();
140  if (!tree.cbeginLeaf()) return points;
141 
142  LeafManagerT leafManager(tree);
143  const Index64 voxelCount = leafManager.activeLeafVoxelCount();
144  OPENVDB_ASSERT(voxelCount != 0);
145 
146  const double pointsPerVolume = double(count) / double(voxelCount);
147  const Index32 pointsPerVoxel = static_cast<Index32>(math::RoundDown(pointsPerVolume));
148  const Index64 remainder = count - (pointsPerVoxel * voxelCount);
149 
150  if (remainder == 0) {
152  GridT, RandGenT, PositionArrayT, PointDataGridT, InterrupterT>(
153  grid, float(pointsPerVoxel), seed, spread, interrupter);
154  }
155 
156  std::vector<Index64> voxelOffsets, values;
157  std::thread worker(&Local::getPrefixSum, std::ref(leafManager), std::ref(voxelOffsets));
158 
159  {
160  math::RandInt<Index64, RandGenT> gen(seed, 0, voxelCount-1);
161  values.reserve(remainder);
162  for (Index64 i = 0; i < remainder; ++i) values.emplace_back(gen());
163  }
164 
165  worker.join();
166 
167  if (util::wasInterrupted<InterrupterT>(interrupter)) {
168  tree.clear();
169  return points;
170  }
171 
172  tbb::parallel_sort(values.begin(), values.end());
173  const bool fractionalOnly(pointsPerVoxel == 0);
174 
175  leafManager.foreach([&voxelOffsets, &values, fractionalOnly]
176  (LeafNodeType& leaf, const size_t idx)
177  {
178  const Index64 lowerOffset = voxelOffsets[idx]; // inclusive
179  const Index64 upperOffset = voxelOffsets[idx + 1]; // exclusive
180  OPENVDB_ASSERT(upperOffset > lowerOffset);
181 
182  const auto valuesEnd = values.end();
183  auto lower = std::lower_bound(values.begin(), valuesEnd, lowerOffset);
184 
185  auto* const data = leaf.buffer().data();
186  auto iter = leaf.beginValueOn();
187 
188  Index32 currentOffset(0);
189  bool addedPoints(!fractionalOnly);
190  while (lower != valuesEnd) {
191  const Index64 vId = *lower;
192  if (vId >= upperOffset) break;
193 
194  const Index32 nextOffset = Index32(vId - lowerOffset);
195  iter.increment(nextOffset - currentOffset);
196  currentOffset = nextOffset;
197  OPENVDB_ASSERT(iter);
198 
199  auto& value = data[iter.pos()];
200  value = value + 1; // no += operator support
201  addedPoints = true;
202  ++lower;
203  }
204 
205  // deactivate this leaf if no points were added. This will speed up
206  // the unthreaded rng
207  if (!addedPoints) leaf.setValuesOff();
208  });
209 
210  voxelOffsets.clear();
211  values.clear();
212 
213  if (fractionalOnly) {
214  tools::pruneInactive(tree);
215  leafManager.rebuild();
216  }
217 
218  const AttributeSet::Descriptor::Ptr descriptor =
219  AttributeSet::Descriptor::create(PositionArrayT::attributeType());
220  RandomGenerator rand01(seed);
221 
222  const auto leafRange = leafManager.leafRange();
223  auto leaf = leafRange.begin();
224  for (; leaf; ++leaf) {
225  if (util::wasInterrupted<InterrupterT>(interrupter)) break;
226  Index32 offset(0);
227  for (auto iter = leaf->beginValueAll(); iter; ++iter) {
228  if (iter.isValueOn()) {
229  const Index32 value = Index32(pointsPerVolume + Index32(*iter));
230  if (value == 0) leaf->setValueOff(iter.pos());
231  else offset += value;
232  }
233  // @note can't use iter.setValue(offset) on point grids
234  leaf->setOffsetOnly(iter.pos(), offset);
235  }
236 
237  // offset should always be non zero
238  OPENVDB_ASSERT(offset != 0);
239  point_scatter_internal::generatePositions<PositionType, CodecType>
240  (*leaf, descriptor, offset, spread, rand01);
241  }
242 
243  // if interrupted, remove remaining leaf nodes
244  if (leaf) {
245  for (; leaf; ++leaf) leaf->setValuesOff();
246  tools::pruneInactive(tree);
247  }
248 
249  if (interrupter) interrupter->end();
250  return points;
251 }
252 
253 
254 ////////////////////////////////////////
255 
256 
257 template<
258  typename GridT,
259  typename RandGenT,
260  typename PositionArrayT,
261  typename PointDataGridT,
262  typename InterrupterT>
263 inline typename PointDataGridT::Ptr
264 denseUniformPointScatter(const GridT& grid,
265  const float pointsPerVoxel,
266  const unsigned int seed,
267  const float spread,
268  InterrupterT* interrupter)
269 {
270  using PositionType = typename PositionArrayT::ValueType;
271  using PositionTraits = VecTraits<PositionType>;
272  using ValueType = typename PositionTraits::ElementType;
273  using CodecType = typename PositionArrayT::Codec;
274 
275  using RandomGenerator = math::Rand01<ValueType, RandGenT>;
276 
277  using TreeType = typename PointDataGridT::TreeType;
278 
279  static_assert(PositionTraits::IsVec && PositionTraits::Size == 3,
280  "Invalid Position Array type.");
281 
282  if (pointsPerVoxel < 0.0f) {
283  OPENVDB_THROW(ValueError, "Points per voxel must not be less than zero.");
284  }
285 
286  if (spread < 0.0f || spread > 1.0f) {
287  OPENVDB_THROW(ValueError, "Spread must be between 0 and 1.");
288  }
289 
290  if (interrupter) interrupter->start("Dense uniform scattering with fixed point count");
291 
292  typename PointDataGridT::Ptr points =
293  point_scatter_internal::initialisePointTopology<PointDataGridT>(grid);
294  TreeType& tree = points->tree();
295  auto leafIter = tree.beginLeaf();
296  if (!leafIter) return points;
297 
298  const Index32 pointsPerVoxelInt = math::Floor(pointsPerVoxel);
299  const double delta = pointsPerVoxel - float(pointsPerVoxelInt);
300  const bool fractional = !math::isApproxZero(delta, 1.0e-6);
301  const bool fractionalOnly = pointsPerVoxelInt == 0;
302 
303  const AttributeSet::Descriptor::Ptr descriptor =
304  AttributeSet::Descriptor::create(PositionArrayT::attributeType());
305  RandomGenerator rand01(seed);
306 
307  for (; leafIter; ++leafIter) {
308  if (util::wasInterrupted<InterrupterT>(interrupter)) break;
309  Index32 offset(0);
310  for (auto iter = leafIter->beginValueAll(); iter; ++iter) {
311  if (iter.isValueOn()) {
312  offset += pointsPerVoxelInt;
313  if (fractional && rand01() < delta) ++offset;
314  else if (fractionalOnly) leafIter->setValueOff(iter.pos());
315  }
316  // @note can't use iter.setValue(offset) on point grids
317  leafIter->setOffsetOnly(iter.pos(), offset);
318  }
319 
320  if (offset != 0) {
321  point_scatter_internal::generatePositions<PositionType, CodecType>
322  (*leafIter, descriptor, offset, spread, rand01);
323  }
324  }
325 
326  // if interrupted, remove remaining leaf nodes
327  const bool prune(leafIter || fractionalOnly);
328  for (; leafIter; ++leafIter) leafIter->setValuesOff();
329 
330  if (prune) tools::pruneInactive(tree);
331  if (interrupter) interrupter->end();
332  return points;
333 }
334 
335 
336 ////////////////////////////////////////
337 
338 
339 template<
340  typename GridT,
341  typename RandGenT,
342  typename PositionArrayT,
343  typename PointDataGridT,
344  typename InterrupterT>
345 inline typename PointDataGridT::Ptr
346 nonUniformPointScatter(const GridT& grid,
347  const float pointsPerVoxel,
348  const unsigned int seed,
349  const float spread,
350  InterrupterT* interrupter)
351 {
352  using PositionType = typename PositionArrayT::ValueType;
353  using PositionTraits = VecTraits<PositionType>;
354  using ValueType = typename PositionTraits::ElementType;
355  using CodecType = typename PositionArrayT::Codec;
356 
357  using RandomGenerator = math::Rand01<ValueType, RandGenT>;
358 
359  using TreeType = typename PointDataGridT::TreeType;
360 
361  static_assert(PositionTraits::IsVec && PositionTraits::Size == 3,
362  "Invalid Position Array type.");
363  static_assert(std::is_arithmetic<typename GridT::ValueType>::value,
364  "Scalar grid type required for weighted voxel scattering.");
365 
366  if (pointsPerVoxel < 0.0f) {
367  OPENVDB_THROW(ValueError, "Points per voxel must not be less than zero.");
368  }
369 
370  if (spread < 0.0f || spread > 1.0f) {
371  OPENVDB_THROW(ValueError, "Spread must be between 0 and 1.");
372  }
373 
374  if (interrupter) interrupter->start("Non-uniform scattering with local point density");
375 
376  typename PointDataGridT::Ptr points =
377  point_scatter_internal::initialisePointTopology<PointDataGridT>(grid);
378  TreeType& tree = points->tree();
379  auto leafIter = tree.beginLeaf();
380  if (!leafIter) return points;
381 
382  const AttributeSet::Descriptor::Ptr descriptor =
383  AttributeSet::Descriptor::create(PositionArrayT::attributeType());
384  RandomGenerator rand01(seed);
385  const auto accessor = grid.getConstAccessor();
386 
387  for (; leafIter; ++leafIter) {
388  if (util::wasInterrupted<InterrupterT>(interrupter)) break;
389  Index32 offset(0);
390  for (auto iter = leafIter->beginValueAll(); iter; ++iter) {
391  if (iter.isValueOn()) {
392  double fractional =
393  double(accessor.getValue(iter.getCoord())) * pointsPerVoxel;
394  fractional = std::max(0.0, fractional);
395  int count = int(fractional);
396  if (rand01() < (fractional - double(count))) ++count;
397  else if (count == 0) leafIter->setValueOff(iter.pos());
398  offset += count;
399  }
400  // @note can't use iter.setValue(offset) on point grids
401  leafIter->setOffsetOnly(iter.pos(), offset);
402  }
403 
404  if (offset != 0) {
405  point_scatter_internal::generatePositions<PositionType, CodecType>
406  (*leafIter, descriptor, offset, spread, rand01);
407  }
408  }
409 
410  // if interrupted, remove remaining leaf nodes
411  for (; leafIter; ++leafIter) leafIter->setValuesOff();
412 
413  tools::pruneInactive(points->tree());
414  if (interrupter) interrupter->end();
415  return points;
416 }
417 
418 
419 } // namespace points
420 } // namespace OPENVDB_VERSION_NAME
421 } // namespace openvdb
422 
423 
424 #endif // OPENVDB_POINTS_POINT_SCATTER_IMPL_HAS_BEEN_INCLUDED
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: PointScatterImpl.h:92
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
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: PointScatterImpl.h:346
float RoundDown(float x)
Return x rounded down to the nearest integer.
Definition: Math.h:803
uint64_t Index64
Definition: Types.h:53
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:110
Definition: Exceptions.h:65
bool isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance...
Definition: Math.h:349
#define OPENVDB_ASSERT(X)
Definition: Assert.h:41
LeafRange leafRange(size_t grainsize=1) const
Return a TBB-compatible LeafRange.
Definition: LeafManager.h:346
int Floor(float x)
Return the floor of x.
Definition: Math.h:848
Simple random integer generator.
Definition: Math.h:202
Definition: Exceptions.h:13
Simple generator of random numbers over the range [0, 1)
Definition: Math.h:166
uint32_t Index32
Definition: Types.h:52
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:335
Codec
Define compression codecs.
Definition: NanoVDB.h:6224
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: PointScatterImpl.h:264
Definition: Types.h:243
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:355
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:121
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:218