OpenVDB  6.2.1
points/PointScatter.h
Go to the documentation of this file.
1 //
3 // Copyright (c) DreamWorks Animation LLC
4 //
5 // All rights reserved. This software is distributed under the
6 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
7 //
8 // Redistributions of source code must retain the above copyright
9 // and license notice and the following restrictions and disclaimer.
10 //
11 // * Neither the name of DreamWorks Animation nor the names of
12 // its contributors may be used to endorse or promote products derived
13 // from this software without specific prior written permission.
14 //
15 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
16 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
17 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
18 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
19 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY INDIRECT, INCIDENTAL,
20 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
21 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
22 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
23 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
25 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26 // IN NO EVENT SHALL THE COPYRIGHT HOLDERS' AND CONTRIBUTORS' AGGREGATE
27 // LIABILITY FOR ALL CLAIMS REGARDLESS OF THEIR BASIS EXCEED US$250.00.
28 //
30 
40 
41 #ifndef OPENVDB_POINTS_POINT_SCATTER_HAS_BEEN_INCLUDED
42 #define OPENVDB_POINTS_POINT_SCATTER_HAS_BEEN_INCLUDED
43 
44 #include <type_traits>
45 #include <algorithm>
46 #include <thread>
47 #include <random>
48 
49 #include <openvdb/openvdb.h>
50 #include <openvdb/Types.h>
52 #include <openvdb/tools/Prune.h>
54 
55 #include "AttributeArray.h"
56 #include "PointCount.h"
57 #include "PointDataGrid.h"
58 
59 #include <tbb/parallel_sort.h>
60 #include <tbb/parallel_for.h>
61 
62 namespace openvdb {
64 namespace OPENVDB_VERSION_NAME {
65 namespace points {
66 
84 
85 
98 template<
99  typename GridT,
100  typename RandGenT = std::mt19937,
101  typename PositionArrayT = TypedAttributeArray<Vec3f, NullCodec>,
102  typename PointDataGridT = Grid<
103  typename points::TreeConverter<typename GridT::TreeType>::Type>,
104  typename InterrupterT = util::NullInterrupter>
105 inline typename PointDataGridT::Ptr
106 uniformPointScatter(const GridT& grid,
107  const Index64 count,
108  const unsigned int seed = 0,
109  const float spread = 1.0f,
110  InterrupterT* interrupter = nullptr);
111 
126 
127 template<
128  typename GridT,
129  typename RandGenT = std::mt19937,
130  typename PositionArrayT = TypedAttributeArray<Vec3f, NullCodec>,
131  typename PointDataGridT = Grid<
132  typename points::TreeConverter<typename GridT::TreeType>::Type>,
133  typename InterrupterT = util::NullInterrupter>
134 inline typename PointDataGridT::Ptr
135 denseUniformPointScatter(const GridT& grid,
136  const float pointsPerVoxel,
137  const unsigned int seed = 0,
138  const float spread = 1.0f,
139  InterrupterT* interrupter = nullptr);
140 
159 template<
160  typename GridT,
161  typename RandGenT = std::mt19937,
162  typename PositionArrayT = TypedAttributeArray<Vec3f, NullCodec>,
163  typename PointDataGridT = Grid<
164  typename points::TreeConverter<typename GridT::TreeType>::Type>,
165  typename InterrupterT = util::NullInterrupter>
166 inline typename PointDataGridT::Ptr
167 nonUniformPointScatter(const GridT& grid,
168  const float pointsPerVoxel,
169  const unsigned int seed = 0,
170  const float spread = 1.0f,
171  InterrupterT* interrupter = nullptr);
172 
173 
175 
176 
177 namespace point_scatter_internal
178 {
179 
183 template<typename PointDataGridT, typename GridT>
184 inline typename PointDataGridT::Ptr
185 initialisePointTopology(const GridT& grid)
186 {
187  typename PointDataGridT::Ptr points(new PointDataGridT);
188  points->setTransform(grid.transform().copy());
189  points->topologyUnion(grid);
190  if (points->tree().hasActiveTiles()) {
191  points->tree().voxelizeActiveTiles();
192  }
193 
194  return points;
195 }
196 
204 template<typename PositionType,
205  typename CodecT,
206  typename RandGenT,
207  typename LeafNodeT>
208 inline void
209 generatePositions(LeafNodeT& leaf,
210  const AttributeSet::Descriptor::Ptr& descriptor,
211  const Index64& count,
212  const float spread,
213  RandGenT& rand01)
214 {
215  using PositionTraits = VecTraits<PositionType>;
216  using ValueType = typename PositionTraits::ElementType;
217  using PositionWriteHandle = AttributeWriteHandle<PositionType, CodecT>;
218 
219  leaf.initializeAttributes(descriptor, static_cast<Index>(count));
220 
221  // directly expand to avoid needlessly setting uniform values in the
222  // write handle
223  auto& array = leaf.attributeArray(0);
224  array.expand(/*fill*/false);
225 
226  PositionWriteHandle pHandle(array, /*expand*/false);
227  PositionType P;
228  for (Index64 index = 0; index < count; ++index) {
229  P[0] = (spread * (rand01() - ValueType(0.5)));
230  P[1] = (spread * (rand01() - ValueType(0.5)));
231  P[2] = (spread * (rand01() - ValueType(0.5)));
232  pHandle.set(static_cast<Index>(index), P);
233  }
234 }
235 
236 } // namespace point_scatter_internal
237 
238 
240 
241 
242 template<
243  typename GridT,
244  typename RandGenT,
245  typename PositionArrayT,
246  typename PointDataGridT,
247  typename InterrupterT>
248 inline typename PointDataGridT::Ptr
249 uniformPointScatter(const GridT& grid,
250  const Index64 count,
251  const unsigned int seed,
252  const float spread,
253  InterrupterT* interrupter)
254 {
255  using PositionType = typename PositionArrayT::ValueType;
256  using PositionTraits = VecTraits<PositionType>;
257  using ValueType = typename PositionTraits::ElementType;
258  using CodecType = typename PositionArrayT::Codec;
259 
260  using RandomGenerator = math::Rand01<ValueType, RandGenT>;
261 
262  using TreeType = typename PointDataGridT::TreeType;
263  using LeafNodeType = typename TreeType::LeafNodeType;
264  using LeafManagerT = tree::LeafManager<TreeType>;
265 
266  struct Local
267  {
271  static void getPrefixSum(LeafManagerT& leafManager,
272  std::vector<Index64>& offsets)
273  {
274  Index64 offset = 0;
275  offsets.reserve(leafManager.leafCount() + 1);
276  offsets.push_back(0);
277  const auto leafRange = leafManager.leafRange();
278  for (auto leaf = leafRange.begin(); leaf; ++leaf) {
279  offset += leaf->onVoxelCount();
280  offsets.push_back(offset);
281  }
282  }
283  };
284 
285  static_assert(PositionTraits::IsVec && PositionTraits::Size == 3,
286  "Invalid Position Array type.");
287 
288  if (spread < 0.0f || spread > 1.0f) {
289  OPENVDB_THROW(ValueError, "Spread must be between 0 and 1.");
290  }
291 
292  if (interrupter) interrupter->start("Uniform scattering with fixed point count");
293 
294  typename PointDataGridT::Ptr points =
295  point_scatter_internal::initialisePointTopology<PointDataGridT>(grid);
296  TreeType& tree = points->tree();
297  if (!tree.cbeginLeaf()) return points;
298 
299  LeafManagerT leafManager(tree);
300  const Index64 voxelCount = leafManager.activeLeafVoxelCount();
301  assert(voxelCount != 0);
302 
303  const double pointsPerVolume = double(count) / double(voxelCount);
304  const Index32 pointsPerVoxel = static_cast<Index32>(math::RoundDown(pointsPerVolume));
305  const Index64 remainder = count - (pointsPerVoxel * voxelCount);
306 
307  if (remainder == 0) {
309  GridT, RandGenT, PositionArrayT, PointDataGridT, InterrupterT>(
310  grid, float(pointsPerVoxel), seed, spread, interrupter);
311  }
312 
313  std::vector<Index64> voxelOffsets, values;
314  std::thread worker(&Local::getPrefixSum, std::ref(leafManager), std::ref(voxelOffsets));
315 
316  {
317  math::RandInt<Index64, RandGenT> gen(seed, 0, voxelCount-1);
318  values.reserve(remainder);
319  for (Index64 i = 0; i < remainder; ++i) values.emplace_back(gen());
320  }
321 
322  worker.join();
323 
324  if (util::wasInterrupted<InterrupterT>(interrupter)) {
325  tree.clear();
326  return points;
327  }
328 
329  tbb::parallel_sort(values.begin(), values.end());
330  const bool fractionalOnly(pointsPerVoxel == 0);
331 
332  leafManager.foreach([&voxelOffsets, &values, fractionalOnly]
333  (LeafNodeType& leaf, const size_t idx)
334  {
335  const Index64 lowerOffset = voxelOffsets[idx]; // inclusive
336  const Index64 upperOffset = voxelOffsets[idx + 1]; // exclusive
337  assert(upperOffset > lowerOffset);
338 
339  const auto valuesEnd = values.end();
340  auto lower = std::lower_bound(values.begin(), valuesEnd, lowerOffset);
341 
342  auto* const data = leaf.buffer().data();
343  auto iter = leaf.beginValueOn();
344 
345  Index32 currentOffset(0);
346  bool addedPoints(!fractionalOnly);
347  while (lower != valuesEnd) {
348  const Index64 vId = *lower;
349  if (vId >= upperOffset) break;
350 
351  const Index32 nextOffset = Index32(vId - lowerOffset);
352  iter.increment(nextOffset - currentOffset);
353  currentOffset = nextOffset;
354  assert(iter);
355 
356  auto& value = data[iter.pos()];
357  value = value + 1; // no += operator support
358  addedPoints = true;
359  ++lower;
360  }
361 
362  // deactivate this leaf if no points were added. This will speed up
363  // the unthreaded rng
364  if (!addedPoints) leaf.setValuesOff();
365  });
366 
367  voxelOffsets.clear();
368  values.clear();
369 
370  if (fractionalOnly) {
371  tools::pruneInactive(tree);
372  leafManager.rebuild();
373  }
374 
375  const AttributeSet::Descriptor::Ptr descriptor =
376  AttributeSet::Descriptor::create(PositionArrayT::attributeType());
377  RandomGenerator rand01(seed);
378 
379  const auto leafRange = leafManager.leafRange();
380  auto leaf = leafRange.begin();
381  for (; leaf; ++leaf) {
382  if (util::wasInterrupted<InterrupterT>(interrupter)) break;
383  Index32 offset(0);
384  for (auto iter = leaf->beginValueAll(); iter; ++iter) {
385  if (iter.isValueOn()) {
386  const Index32 value = Index32(pointsPerVolume + Index32(*iter));
387  if (value == 0) leaf->setValueOff(iter.pos());
388  else offset += value;
389  }
390  // @note can't use iter.setValue(offset) on point grids
391  leaf->setOffsetOnly(iter.pos(), offset);
392  }
393 
394  // offset should always be non zero
395  assert(offset != 0);
396  point_scatter_internal::generatePositions<PositionType, CodecType>
397  (*leaf, descriptor, offset, spread, rand01);
398  }
399 
400  // if interrupted, remove remaining leaf nodes
401  if (leaf) {
402  for (; leaf; ++leaf) leaf->setValuesOff();
403  tools::pruneInactive(tree);
404  }
405 
406  if (interrupter) interrupter->end();
407  return points;
408 }
409 
410 
412 
413 
414 template<
415  typename GridT,
416  typename RandGenT,
417  typename PositionArrayT,
418  typename PointDataGridT,
419  typename InterrupterT>
420 inline typename PointDataGridT::Ptr
421 denseUniformPointScatter(const GridT& grid,
422  const float pointsPerVoxel,
423  const unsigned int seed,
424  const float spread,
425  InterrupterT* interrupter)
426 {
427  using PositionType = typename PositionArrayT::ValueType;
428  using PositionTraits = VecTraits<PositionType>;
429  using ValueType = typename PositionTraits::ElementType;
430  using CodecType = typename PositionArrayT::Codec;
431 
432  using RandomGenerator = math::Rand01<ValueType, RandGenT>;
433 
434  using TreeType = typename PointDataGridT::TreeType;
435 
436  static_assert(PositionTraits::IsVec && PositionTraits::Size == 3,
437  "Invalid Position Array type.");
438 
439  if (pointsPerVoxel < 0.0f) {
440  OPENVDB_THROW(ValueError, "Points per voxel must not be less than zero.");
441  }
442 
443  if (spread < 0.0f || spread > 1.0f) {
444  OPENVDB_THROW(ValueError, "Spread must be between 0 and 1.");
445  }
446 
447  if (interrupter) interrupter->start("Dense uniform scattering with fixed point count");
448 
449  typename PointDataGridT::Ptr points =
450  point_scatter_internal::initialisePointTopology<PointDataGridT>(grid);
451  TreeType& tree = points->tree();
452  auto leafIter = tree.beginLeaf();
453  if (!leafIter) return points;
454 
455  const Index32 pointsPerVoxelInt = math::Floor(pointsPerVoxel);
456  const double delta = pointsPerVoxel - float(pointsPerVoxelInt);
457  const bool fractional = !math::isApproxZero(delta, 1.0e-6);
458  const bool fractionalOnly = pointsPerVoxelInt == 0;
459 
460  const AttributeSet::Descriptor::Ptr descriptor =
461  AttributeSet::Descriptor::create(PositionArrayT::attributeType());
462  RandomGenerator rand01(seed);
463 
464  for (; leafIter; ++leafIter) {
465  if (util::wasInterrupted<InterrupterT>(interrupter)) break;
466  Index32 offset(0);
467  for (auto iter = leafIter->beginValueAll(); iter; ++iter) {
468  if (iter.isValueOn()) {
469  offset += pointsPerVoxelInt;
470  if (fractional && rand01() < delta) ++offset;
471  else if (fractionalOnly) leafIter->setValueOff(iter.pos());
472  }
473  // @note can't use iter.setValue(offset) on point grids
474  leafIter->setOffsetOnly(iter.pos(), offset);
475  }
476 
477  if (offset != 0) {
478  point_scatter_internal::generatePositions<PositionType, CodecType>
479  (*leafIter, descriptor, offset, spread, rand01);
480  }
481  }
482 
483  // if interrupted, remove remaining leaf nodes
484  const bool prune(leafIter || fractionalOnly);
485  for (; leafIter; ++leafIter) leafIter->setValuesOff();
486 
487  if (prune) tools::pruneInactive(tree);
488  if (interrupter) interrupter->end();
489  return points;
490 }
491 
492 
494 
495 
496 template<
497  typename GridT,
498  typename RandGenT,
499  typename PositionArrayT,
500  typename PointDataGridT,
501  typename InterrupterT>
502 inline typename PointDataGridT::Ptr
503 nonUniformPointScatter(const GridT& grid,
504  const float pointsPerVoxel,
505  const unsigned int seed,
506  const float spread,
507  InterrupterT* interrupter)
508 {
509  using PositionType = typename PositionArrayT::ValueType;
510  using PositionTraits = VecTraits<PositionType>;
511  using ValueType = typename PositionTraits::ElementType;
512  using CodecType = typename PositionArrayT::Codec;
513 
514  using RandomGenerator = math::Rand01<ValueType, RandGenT>;
515 
516  using TreeType = typename PointDataGridT::TreeType;
517 
518  static_assert(PositionTraits::IsVec && PositionTraits::Size == 3,
519  "Invalid Position Array type.");
520  static_assert(std::is_arithmetic<typename GridT::ValueType>::value,
521  "Scalar grid type required for weighted voxel scattering.");
522 
523  if (pointsPerVoxel < 0.0f) {
524  OPENVDB_THROW(ValueError, "Points per voxel must not be less than zero.");
525  }
526 
527  if (spread < 0.0f || spread > 1.0f) {
528  OPENVDB_THROW(ValueError, "Spread must be between 0 and 1.");
529  }
530 
531  if (interrupter) interrupter->start("Non-uniform scattering with local point density");
532 
533  typename PointDataGridT::Ptr points =
534  point_scatter_internal::initialisePointTopology<PointDataGridT>(grid);
535  TreeType& tree = points->tree();
536  auto leafIter = tree.beginLeaf();
537  if (!leafIter) return points;
538 
539  const AttributeSet::Descriptor::Ptr descriptor =
540  AttributeSet::Descriptor::create(PositionArrayT::attributeType());
541  RandomGenerator rand01(seed);
542  const auto accessor = grid.getConstAccessor();
543 
544  for (; leafIter; ++leafIter) {
545  if (util::wasInterrupted<InterrupterT>(interrupter)) break;
546  Index32 offset(0);
547  for (auto iter = leafIter->beginValueAll(); iter; ++iter) {
548  if (iter.isValueOn()) {
549  double fractional =
550  double(accessor.getValue(iter.getCoord())) * pointsPerVoxel;
551  fractional = std::max(0.0, fractional);
552  int count = int(fractional);
553  if (rand01() < (fractional - double(count))) ++count;
554  else if (count == 0) leafIter->setValueOff(iter.pos());
555  offset += count;
556  }
557  // @note can't use iter.setValue(offset) on point grids
558  leafIter->setOffsetOnly(iter.pos(), offset);
559  }
560 
561  if (offset != 0) {
562  point_scatter_internal::generatePositions<PositionType, CodecType>
563  (*leafIter, descriptor, offset, spread, rand01);
564  }
565  }
566 
567  // if interrupted, remove remaining leaf nodes
568  for (; leafIter; ++leafIter) leafIter->setValuesOff();
569 
570  tools::pruneInactive(points->tree());
571  if (interrupter) interrupter->end();
572  return points;
573 }
574 
575 
576 } // namespace points
577 } // namespace OPENVDB_VERSION_NAME
578 } // namespace openvdb
579 
580 
581 #endif // OPENVDB_POINTS_POINT_SCATTER_HAS_BEEN_INCLUDED
582 
583 // Copyright (c) DreamWorks Animation LLC
584 // All rights reserved. This software is distributed under the
585 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:133
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:109
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:249
Attribute Array storage templated on type and compression codec.
Methods for counting points in VDB Point grids.
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:361
LeafRange leafRange(size_t grainsize=1) const
Return a TBB-compatible LeafRange.
Definition: LeafManager.h:385
Write-able version of AttributeHandle.
Definition: AttributeArray.h:944
Definition: Types.h:231
Simple random integer generator.
Definition: Math.h:171
Defined various multi-threaded utility functions for trees.
openvdb::GridBase Grid
Definition: Utils.h:60
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:128
bool isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance...
Definition: Math.h:320
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:381
PointDataGridT::Ptr initialisePointTopology(const GridT &grid)
initialise the topology of a PointDataGrid and ensure everything is voxelized
Definition: points/PointScatter.h:185
Simple generator of random numbers over the range [0, 1)
Definition: Math.h:135
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:209
Definition: Exceptions.h:40
uint64_t Index64
Definition: Types.h:60
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:421
This class manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional auxil...
Definition: LeafManager.h:109
Attribute-owned data structure for points. Point attributes are stored in leaf nodes and ordered by v...
Definition: Exceptions.h:92
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:180
A LeafManager manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional au...
float RoundDown(float x)
Return x rounded down to the nearest integer.
Definition: Math.h:777
uint32_t Index32
Definition: Types.h:59
int Floor(float x)
Return the floor of x.
Definition: Math.h:822
void expand(bool fill=true)
If this array is uniform, replace it with an array of length size().
Definition: AttributeArray.h:2286
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:503