OpenVDB 11.0.0
Loading...
Searching...
No Matches
PointScatterImpl.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 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
12namespace openvdb {
14namespace OPENVDB_VERSION_NAME {
15namespace points {
16
17/// @cond OPENVDB_DOCS_INTERNAL
18
19namespace 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
25template<typename PointDataGridT, typename GridT>
26inline typename PointDataGridT::Ptr
27initialisePointTopology(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.
46template<typename PositionType,
47 typename CodecT,
48 typename RandGenT,
49 typename LeafNodeT>
50inline void
51generatePositions(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
85template<
86 typename GridT,
87 typename RandGenT,
88 typename PositionArrayT,
89 typename PointDataGridT,
90 typename InterrupterT>
91inline typename PointDataGridT::Ptr
92uniformPointScatter(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 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 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 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 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
257template<
258 typename GridT,
259 typename RandGenT,
260 typename PositionArrayT,
261 typename PointDataGridT,
262 typename InterrupterT>
263inline typename PointDataGridT::Ptr
264denseUniformPointScatter(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
339template<
340 typename GridT,
341 typename RandGenT,
342 typename PositionArrayT,
343 typename PointDataGridT,
344 typename InterrupterT>
345inline typename PointDataGridT::Ptr
346nonUniformPointScatter(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
Definition Exceptions.h:65
Simple generator of random numbers over the range [0, 1)
Definition Math.h:167
Simple random integer generator.
Definition Math.h:203
std::shared_ptr< Descriptor > Ptr
Definition AttributeSet.h:314
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition LeafManager.h:85
LeafRange leafRange(size_t grainsize=1) const
Return a TBB-compatible LeafRange.
Definition LeafManager.h:345
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
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
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
Index32 Index
Definition Types.h:54
uint32_t Index32
Definition Types.h:52
uint64_t Index64
Definition Types.h:53
Definition Exceptions.h:13
#define OPENVDB_THROW(exception, message)
Definition Exceptions.h:74
Definition Types.h:244
#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:212