OpenVDB 11.0.0
Loading...
Searching...
No Matches
PointStatisticsImpl.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 PointStatisticsImpl.h
7///
8
9#ifndef OPENVDB_POINTS_STATISTICS_IMPL_HAS_BEEN_INCLUDED
10#define OPENVDB_POINTS_STATISTICS_IMPL_HAS_BEEN_INCLUDED
11
12namespace openvdb {
14namespace OPENVDB_VERSION_NAME {
15namespace points {
16
17/// @cond OPENVDB_DOCS_INTERNAL
18
19namespace statistics_internal
20{
21
22/// @brief Scalar extent op to evaluate the min/max values of a
23/// single integral or floating point attribute type
24template <typename ValueT>
25struct ScalarMinMax
26{
27 using ExtentT = std::pair<ValueT, ValueT>;
28 ScalarMinMax(const ValueT& init) : mMinMax(init, init) {}
29 ScalarMinMax(const ExtentT& init) : mMinMax(init) {}
30 inline void operator()(const ValueT& b)
31 {
32 mMinMax.first = std::min(mMinMax.first, b);
33 mMinMax.second = std::max(mMinMax.second, b);
34 }
35 inline void operator()(const ExtentT& b)
36 {
37 mMinMax.first = std::min(mMinMax.first, b.first);
38 mMinMax.second = std::max(mMinMax.second, b.second);
39 }
40 inline const ExtentT& get() const { return mMinMax; }
41 ExtentT mMinMax;
42};
43
44/// @brief Vector squared magnitude op to evaluate the min/max of a
45/// vector attribute and return the result as a scalar of the
46/// appropriate precision
47template <typename VecT, bool MagResult = true>
48struct MagnitudeExtent
49 : public ScalarMinMax<typename ValueTraits<VecT>::ElementType>
50{
51 using ElementT = typename ValueTraits<VecT>::ElementType;
52 using ExtentT = typename ScalarMinMax<ElementT>::ExtentT;
53 using BaseT = ScalarMinMax<ElementT>;
54 MagnitudeExtent(const VecT& init) : BaseT(init.lengthSqr()) {}
55 MagnitudeExtent(const ExtentT& init) : BaseT(init) {}
56 inline void operator()(const VecT& b) { this->BaseT::operator()(b.lengthSqr()); }
57 inline void operator()(const ExtentT& b) { this->BaseT::operator()(b); }
58};
59
60/// @brief Vector squared magnitude op to evaluate the min/max of a
61/// vector attribute and return the result as the original vector
62template <typename VecT>
63struct MagnitudeExtent<VecT, false>
64{
65 using ElementT = typename ValueTraits<VecT>::ElementType;
66 using ExtentT = std::pair<VecT, VecT>;
67 MagnitudeExtent(const VecT& init)
68 : mLengths(), mMinMax(init, init) {
69 mLengths.first = init.lengthSqr();
70 mLengths.second = mLengths.first;
71 }
72 MagnitudeExtent(const ExtentT& init)
73 : mLengths(), mMinMax(init) {
74 mLengths.first = init.first.lengthSqr();
75 mLengths.second = init.second.lengthSqr();
76 }
77 inline const ExtentT& get() const { return mMinMax; }
78 inline void operator()(const VecT& b)
79 {
80 const ElementT l = b.lengthSqr();
81 if (l < mLengths.first) {
82 mLengths.first = l;
83 mMinMax.first = b;
84 }
85 else if (l > mLengths.second) {
86 mLengths.second = l;
87 mMinMax.second = b;
88 }
89 }
90 inline void operator()(const ExtentT& b)
91 {
92 ElementT l = b.first.lengthSqr();
93 if (l < mLengths.first) {
94 mLengths.first = l;
95 mMinMax.first = b.first;
96 }
97 l = b.second.lengthSqr();
98 if (l > mLengths.second) {
99 mLengths.second = l;
100 mMinMax.second = b.second;
101 }
102 }
103
104 std::pair<ElementT, ElementT> mLengths;
105 ExtentT mMinMax;
106};
107
108/// @brief Vector component-wise op to evaluate the min/max of
109/// vector components and return the result as a vector of
110/// equal size and precision
111template <typename VecT>
112struct ComponentExtent
113{
114 using ExtentT = std::pair<VecT, VecT>;
115 ComponentExtent(const VecT& init) : mMinMax(init, init) {}
116 ComponentExtent(const ExtentT& init) : mMinMax(init) {}
117 inline const ExtentT& get() const { return mMinMax; }
118 inline void operator()(const VecT& b)
119 {
120 mMinMax.first = math::minComponent(mMinMax.first, b);
121 mMinMax.second = math::maxComponent(mMinMax.second, b);
122 }
123 inline void operator()(const ExtentT& b)
124 {
125 mMinMax.first = math::minComponent(mMinMax.first, b.first);
126 mMinMax.second = math::maxComponent(mMinMax.second, b.second);
127 }
128
129 ExtentT mMinMax;
130};
131
132template <typename ValueT,
133 typename CodecT,
134 typename FilterT,
135 typename ExtentOp,
136 typename PointDataTreeT>
137bool evalExtents(const PointDataTreeT& points,
138 const std::string& attribute,
139 typename ExtentOp::ExtentT& ext,
140 const FilterT& filter,
141 typename PointDataTreeT::template ValueConverter
142 <typename ExtentOp::ExtentT::first_type>::Type* const minTree = nullptr,
143 typename PointDataTreeT::template ValueConverter
144 <typename ExtentOp::ExtentT::second_type>::Type* const maxTree = nullptr)
145{
146 static_assert(std::is_base_of<TreeBase, PointDataTreeT>::value,
147 "PointDataTreeT in instantiation of evalExtents is not an openvdb Tree type");
148
149 struct ResultType {
150 typename ExtentOp::ExtentT ext;
151 bool data = false;
152 };
153
154 tree::LeafManager<const PointDataTreeT> manager(points);
155 if (manager.leafCount() == 0) return false;
156 const size_t idx = manager.leaf(0).attributeSet().find(attribute);
157 if (idx == AttributeSet::INVALID_POS) return false;
158
159 // track results per leaf for min/max trees
160 std::vector<std::unique_ptr<typename ExtentOp::ExtentT>> values;
161 if (minTree || maxTree) values.resize(manager.leafCount());
162
163 const ResultType result = tbb::parallel_reduce(manager.leafRange(),
164 ResultType(),
165 [idx, &filter, &values]
166 (const auto& range, ResultType in) -> ResultType
167 {
168 for (auto leaf = range.begin(); leaf; ++leaf) {
169 AttributeHandle<ValueT, CodecT> handle(leaf->constAttributeArray(idx));
170 if (handle.size() == 0) continue;
171 if (filter.state() == index::ALL) {
172 const size_t size = handle.isUniform() ? 1 : handle.size();
173 ExtentOp op(handle.get(0));
174 for (size_t i = 1; i < size; ++i) {
175 assert(i < size_t(std::numeric_limits<Index>::max()));
176 op(handle.get(Index(i)));
177 }
178 if (!values.empty()) {
179 values[leaf.pos()].reset(new typename ExtentOp::ExtentT(op.get()));
180 }
181 if (in.data) op(in.ext);
182 in.data = true;
183 in.ext = op.get();
184 }
185 else {
186 auto iter = leaf->beginIndexOn(filter);
187 if (!iter) continue;
188 ExtentOp op(handle.get(*iter));
189 ++iter;
190 for (; iter; ++iter) op(handle.get(*iter));
191 if (!values.empty()) {
192 values[leaf.pos()].reset(new typename ExtentOp::ExtentT(op.get()));
193 }
194 if (in.data) op(in.ext);
195 in.data = true;
196 in.ext = op.get();
197 }
198 }
199
200 return in;
201 },
202 [](const ResultType& a, const ResultType& b) -> ResultType {
203 if (!b.data) return a;
204 if (!a.data) return b;
205 ExtentOp op(a.ext); op(b.ext);
206 ResultType t;
207 t.ext = op.get();
208 t.data = true;
209 return t;
210 });
211
212 // set minmax trees only if a new value was set - if the value
213 // hasn't changed, leave it as inactive background (this is
214 // only possible if a point leaf exists with no points or if a
215 // filter is provided but is not hit for a given leaf)
216 if (minTree || maxTree) {
217 manager.foreach([minTree, maxTree, &values]
218 (const auto& leaf, const size_t idx) {
219 const auto& v = values[idx];
220 if (v == nullptr) return;
221 const Coord& origin = leaf.origin();
222 if (minTree) minTree->addTile(1, origin, v->first, true);
223 if (maxTree) maxTree->addTile(1, origin, v->second, true);
224 }, false);
225 }
226
227 if (result.data) ext = result.ext;
228 return result.data;
229}
230
231template <typename ValueT,
232 typename CodecT,
233 typename FilterT,
234 typename PointDataTreeT,
235 typename std::enable_if<ValueTraits<ValueT>::IsVec, int>::type = 0>
236bool evalExtents(const PointDataTreeT& points,
237 const std::string& attribute,
238 ValueT& min,
239 ValueT& max,
240 const FilterT& filter,
241 typename PointDataTreeT::template ValueConverter<ValueT>::Type* minTree,
242 typename PointDataTreeT::template ValueConverter<ValueT>::Type* maxTree)
243{
244 typename ComponentExtent<ValueT>::ExtentT ext;
245 const bool s = evalExtents<ValueT, CodecT, FilterT,
246 ComponentExtent<ValueT>, PointDataTreeT>
247 (points, attribute, ext, filter, minTree, maxTree);
248 if (s) min = ext.first, max = ext.second;
249 return s;
250}
251
252template <typename ValueT,
253 typename CodecT,
254 typename FilterT,
255 typename PointDataTreeT,
256 typename std::enable_if<!ValueTraits<ValueT>::IsVec, int>::type = 0>
257bool evalExtents(const PointDataTreeT& points,
258 const std::string& attribute,
259 ValueT& min,
260 ValueT& max,
261 const FilterT& filter,
262 typename PointDataTreeT::template ValueConverter<ValueT>::Type* minTree,
263 typename PointDataTreeT::template ValueConverter<ValueT>::Type* maxTree)
264{
265 typename ScalarMinMax<ValueT>::ExtentT ext;
266 const bool s = evalExtents<ValueT, CodecT, FilterT,
267 ScalarMinMax<ValueT>, PointDataTreeT>
268 (points, attribute, ext, filter, minTree, maxTree);
269 if (s) min = ext.first, max = ext.second;
270 return s;
271}
272
273} // namespace statistics_internal
274
275/// @endcond
276
277template <typename ValueT,
278 typename CodecT,
279 typename FilterT,
280 typename PointDataTreeT>
281bool evalMinMax(const PointDataTreeT& points,
282 const std::string& attribute,
283 ValueT& min,
284 ValueT& max,
285 const FilterT& filter,
286 typename PointDataTreeT::template ValueConverter<ValueT>::Type* minTree,
287 typename PointDataTreeT::template ValueConverter<ValueT>::Type* maxTree)
288{
289 return statistics_internal::evalExtents<ValueT, CodecT, FilterT, PointDataTreeT>
290 (points, attribute, min, max, filter, minTree, maxTree);
291}
292
293template <typename ValueT,
294 typename CodecT,
295 typename FilterT,
296 typename PointDataTreeT,
297 typename ResultTreeT>
298bool evalAverage(const PointDataTreeT& points,
299 const std::string& attribute,
301 const FilterT& filter,
302 typename PointDataTreeT::template ValueConverter<ResultTreeT>::Type* averageTree)
303{
304 using ResultT = typename ConvertElementType<ValueT, double>::Type;
305
306 struct Sample
307 {
308 Sample(const ResultT& _avg, size_t _size) : avg(_avg), size(_size) {}
309
310 void add(const ResultT& val)
311 {
312 ++size;
313 const ResultT delta = val - avg;
314 avg = avg + (delta / static_cast<double>(size));
315 }
316
317 void add(const Sample& other)
318 {
319 assert(other.size > 0);
320 const double denom = 1.0 / static_cast<double>(size + other.size);
321 const ResultT delta = other.avg - avg;
322 avg = avg + (denom * delta * static_cast<double>(other.size));
323 size += other.size;
324 }
325
326 ResultT avg; size_t size;
327 };
328
329 static_assert(std::is_base_of<TreeBase, PointDataTreeT>::value,
330 "PointDataTreeT in instantiation of evalAverage is not an openvdb Tree type");
331 static_assert(std::is_constructible<ResultT, ValueT>::value,
332 "Target value in points::evalAverage is not constructible from the source value type.");
333
335 if (manager.leafCount() == 0) return false;
336 const size_t idx = manager.leaf(0).attributeSet().find(attribute);
337 if (idx == AttributeSet::INVALID_POS) return false;
338
339 std::vector<std::unique_ptr<Sample>> values;
340 values.resize(manager.leafCount());
341 tbb::parallel_for(manager.leafRange(),
342 [idx, &filter, &values] (const auto& range) {
343 for (auto leaf = range.begin(); leaf; ++leaf) {
344 AttributeHandle<ValueT, CodecT> handle(leaf->constAttributeArray(idx));
345 size_t size = handle.size();
346 if (size == 0) continue;
347 if (filter.state() == index::ALL) {
348 std::unique_ptr<Sample> S(new Sample(ResultT(handle.get(0)), 1));
349 if (handle.isUniform()) {
350 S->avg = S->avg / static_cast<double>(size);
351 S->size = size;
352 }
353 else {
354 for (size_t i = 1; i < size; ++i) {
355 assert(i < size_t(std::numeric_limits<Index>::max()));
356 S->add(ResultT(handle.get(Index(i))));
357 }
358 }
359 values[leaf.pos()] = std::move(S);
360 }
361 else {
362 auto iter = leaf->beginIndexOn(filter);
363 if (!iter) continue;
364 std::unique_ptr<Sample> S(new Sample(ResultT(handle.get(*iter)), 1));
365 ++iter;
366 for (; iter; ++iter, ++size) {
367 S->add(ResultT(handle.get(*iter)));
368 }
369 values[leaf.pos()] = std::move(S);
370 }
371 }
372 });
373
374 auto iter = values.cbegin();
375 while (iter != values.cend() && !(*iter)) ++iter;
376 if (iter == values.cend()) return false;
377 assert(*iter);
378
379 // serial deterministic reduction of floating point samples
380 Sample S = **iter;
381 ++iter;
382 for (; iter != values.cend(); ++iter) {
383 if (*iter) S.add(**iter);
384 }
385 average = S.avg;
386
387 // set average tree only if a new value was set - if the value
388 // hasn't changed, leave it as inactive background (this is
389 // only possible if a point leaf exists with no points or if a
390 // filter is provided but is not hit for a given leaf)
391 if (averageTree) {
392 manager.foreach([averageTree, &values]
393 (const auto& leaf, const size_t idx) {
394 const auto& S = values[idx];
395 if (S == nullptr) return;
396 const Coord& origin = leaf.origin();
397 averageTree->addTile(1, origin, S->avg, true);
398 }, false);
399 }
400
401 return true;
402}
403
404template <typename ValueT,
405 typename CodecT,
406 typename FilterT,
407 typename PointDataTreeT,
408 typename ResultTreeT>
409bool accumulate(const PointDataTreeT& points,
410 const std::string& attribute,
411 typename PromoteType<ValueT>::Highest& total,
412 const FilterT& filter,
413 typename PointDataTreeT::template ValueConverter<ResultTreeT>::Type* totalTree)
414{
415 using ResultT = typename PromoteType<ValueT>::Highest;
416 using ElementT = typename ValueTraits<ResultT>::ElementType;
417
418 static_assert(std::is_base_of<TreeBase, PointDataTreeT>::value,
419 "PointDataTreeT in instantiation of accumulate is not an openvdb Tree type");
420 static_assert(std::is_constructible<ResultT, ValueT>::value,
421 "Target value in points::accumulate is not constructible from the source value type.");
422
424 if (manager.leafCount() == 0) return false;
425 const size_t idx = manager.leaf(0).attributeSet().find(attribute);
426 if (idx == AttributeSet::INVALID_POS) return false;
427
428 std::vector<std::unique_ptr<ResultT>> values;
429 values.resize(manager.leafCount());
430 tbb::parallel_for(manager.leafRange(),
431 [idx, &filter, &values](const auto& range) {
432 for (auto leaf = range.begin(); leaf; ++leaf) {
433 AttributeHandle<ValueT, CodecT> handle(leaf->constAttributeArray(idx));
434 if (handle.size() == 0) continue;
435 if (filter.state() == index::ALL) {
436 const size_t size = handle.isUniform() ? 1 : handle.size();
437 auto total = ResultT(handle.get(0));
438 for (size_t i = 1; i < size; ++i) {
439 assert(i < size_t(std::numeric_limits<Index>::max()));
440 total += ResultT(handle.get(Index(i)));
441 }
442 values[leaf.pos()].reset(new ResultT(total));
443 }
444 else {
445 auto iter = leaf->beginIndexOn(filter);
446 if (!iter) continue;
447 auto total = ResultT(handle.get(*iter));
448 ++iter;
449 for (; iter; ++iter) total += ResultT(handle.get(*iter));
450 values[leaf.pos()].reset(new ResultT(total));
451 }
452 }
453 });
454
455 auto iter = values.cbegin();
456 while (iter != values.cend() && !(*iter)) ++iter;
457 if (iter == values.cend()) return false;
458 assert(*iter);
459 total = **iter; ++iter;
460
461 if (std::is_integral<ElementT>::value) {
462 using RangeT = tbb::blocked_range<const std::unique_ptr<ResultT>*>;
463 // reasonable grain size for accumulation of single to matrix types
464 total = tbb::parallel_reduce(RangeT(&(*iter), (&values.back())+1, 32), total,
465 [](const RangeT& range, ResultT p) -> ResultT {
466 for (const auto& r : range) if (r) p += *r;
467 return p;
468 }, std::plus<ResultT>());
469 }
470 else {
471 for (; iter != values.cend(); ++iter) {
472 if (*iter) total += (**iter);
473 }
474 }
475
476 // set total tree only if a new value was set - if the value
477 // hasn't changed, leave it as inactive background (this is
478 // only possible if a point leaf exists with no points or if a
479 // filter is provided but is not hit for a given leaf)
480 if (totalTree) {
481 manager.foreach([totalTree, &values]
482 (const auto& leaf, const size_t idx) {
483 const auto& v = values[idx];
484 if (v == nullptr) return;
485 const Coord& origin = leaf.origin();
486 totalTree->addTile(1, origin, *v, true);
487 }, false);
488 }
489
490 return true;
491}
492
493template <typename ValueT,
494 typename CodecT,
495 typename FilterT,
496 typename PointDataTreeT>
497std::pair<ValueT, ValueT>
498evalMinMax(const PointDataTreeT& points,
499 const std::string& attribute,
500 const FilterT& filter)
501{
502 std::pair<ValueT, ValueT> results {
504 };
505 evalMinMax<ValueT, CodecT, FilterT, PointDataTreeT>
506 (points, attribute, results.first, results.second, filter);
507 return results;
508}
509
510template <typename ValueT,
511 typename CodecT,
512 typename FilterT,
513 typename PointDataTreeT>
515evalAverage(const PointDataTreeT& points,
516 const std::string& attribute,
517 const FilterT& filter)
518{
519 using ConvertedT = typename ConvertElementType<ValueT, double>::Type;
520 ConvertedT result = zeroVal<ConvertedT>();
521 evalAverage<ValueT, CodecT, FilterT, PointDataTreeT>(points, attribute, result, filter);
522 return result;
523}
524
525template <typename ValueT,
526 typename CodecT,
527 typename FilterT,
528 typename PointDataTreeT>
530accumulate(const PointDataTreeT& points,
531 const std::string& attribute,
532 const FilterT& filter)
533{
534 using PromotedT = typename PromoteType<ValueT>::Highest;
535 PromotedT result = zeroVal<PromotedT>();
536 accumulate<ValueT, CodecT, FilterT, PointDataTreeT>(points, attribute, result, filter);
537 return result;
538}
539
540} // namespace points
541} // namespace OPENVDB_VERSION_NAME
542} // namespace openvdb
543
544#endif // OPENVDB_POINTS_STATISTICS_IMPL_HAS_BEEN_INCLUDED
Signed (x, y, z) 32-bit integer coordinates.
Definition Coord.h:25
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition LeafManager.h:85
LeafType & leaf(size_t leafIdx) const
Return a pointer to the leaf node at index leafIdx in the array.
Definition LeafManager.h:318
size_t leafCount() const
Return the number of leaf nodes.
Definition LeafManager.h:287
LeafRange leafRange(size_t grainsize=1) const
Return a TBB-compatible LeafRange.
Definition LeafManager.h:345
PromoteType< ValueT >::Highest accumulate(const PointDataTreeT &points, const std::string &attribute, const FilterT &filter=NullFilter())
Evaluates the total value of a point attribute.
Definition PointStatisticsImpl.h:530
std::pair< ValueT, ValueT > evalMinMax(const PointDataTreeT &points, const std::string &attribute, const FilterT &filter=NullFilter())
Evaluates the minimum and maximum values of a point attribute.
Definition PointStatisticsImpl.h:498
ConvertElementType< ValueT, double >::Type evalAverage(const PointDataTreeT &points, const std::string &attribute, const FilterT &filter=NullFilter())
Evaluates the average value of a point attribute.
Definition PointStatisticsImpl.h:515
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition Composite.h:110
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition Composite.h:106
Index32 Index
Definition Types.h:54
constexpr T zeroVal()
Return the value of type T that corresponds to zero.
Definition Math.h:70
Definition Exceptions.h:13
SubT Type
Definition Types.h:320
typename TypeT< 64ul >::type Highest
Definition Types.h:371
typename T::ValueType ElementType
Definition Types.h:302
#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