OpenVDB 10.1.0
Loading...
Searching...
No Matches
GridStats.h
Go to the documentation of this file.
1// Copyright Contributors to the OpenVDB Project
2// SPDX-License-Identifier: MPL-2.0
3
4/*!
5 \file GridStats.h
6
7 \author Ken Museth
8
9 \date August 29, 2020
10
11 \brief Re-computes min/max/avg/var/bbox information for each node in a
12 pre-existing NanoVDB grid.
13*/
14
15#ifndef NANOVDB_GRIDSTATS_H_HAS_BEEN_INCLUDED
16#define NANOVDB_GRIDSTATS_H_HAS_BEEN_INCLUDED
17
18#include "../NanoVDB.h"
19#include "Range.h"
20#include "ForEach.h"
21
22#ifdef NANOVDB_USE_TBB
23#include <tbb/parallel_reduce.h>
24#endif
25
26#include <atomic>
27#include <iostream>
28
29namespace nanovdb {
30
31/// @brief Grid flags which indicate what extra information is present in the grid buffer
32enum class StatsMode : uint32_t {
33 Disable = 0,// disable the computation of any type of statistics (obviously the FASTEST!)
34 BBox = 1,// only compute the bbox of active values per node and total activeVoxelCount
35 MinMax = 2,// additionally compute extrema values
36 All = 3,// compute all of the statics, i.e. bbox, min/max, average and standard deviation
37 Default = 3,// default computational mode for statistics
38 End = 4,
39};
40
41/// @brief Re-computes the min/max, stats and bbox information for an existing NanoVDB Grid
42///
43/// @param grid Grid whose stats to update
44/// @param mode Mode of computation for the statistics.
45template<typename BuildT>
46void gridStats(NanoGrid<BuildT>& grid, StatsMode mode = StatsMode::Default);
47
48//================================================================================================
49
50template<typename ValueT, int Rank = TensorTraits<ValueT>::Rank>
51class Extrema;
52
53/// @brief Template specialization of Extrema on scalar value types, i.e. rank = 0
54template<typename ValueT>
55class Extrema<ValueT, 0>
56{
57protected:
58 ValueT mMin, mMax;
59
60public:
61 using ValueType = ValueT;
63 : mMin(std::numeric_limits<ValueT>::max())
64 , mMax(std::numeric_limits<ValueT>::lowest())
65 {
66 }
67 Extrema(const ValueT& v)
68 : mMin(v)
69 , mMax(v)
70 {
71 }
72 Extrema(const ValueT& a, const ValueT& b)
73 : mMin(a)
74 , mMax(b)
75 {
76 }
77 Extrema& min(const ValueT& v)
78 {
79 if (v < mMin) {
80 mMin = v;
81 }
82 return *this;
83 }
84 Extrema& max(const ValueT& v)
85 {
86 if (v > mMax) {
87 mMax = v;
88 }
89 return *this;
90 }
91 Extrema& add(const ValueT& v)
92 {
93 this->min(v);
94 this->max(v);
95 return *this;
96 }
97 Extrema& add(const ValueT& v, uint64_t) { return this->add(v); }
99 {
100 this->min(other.mMin);
101 this->max(other.mMax);
102 return *this;
103 }
104 const ValueT& min() const { return mMin; }
105 const ValueT& max() const { return mMax; }
106 operator bool() const { return mMin <= mMax; }
107 static constexpr bool hasMinMax() { return !std::is_same<bool, ValueT>::value; }
108 static constexpr bool hasAverage() { return false; }
109 static constexpr bool hasStdDeviation() { return false; }
110 static constexpr size_t size() { return 0; }
111}; // Extrema<T, 0>
112
113/// @brief Template specialization of Extrema on vector value types, i.e. rank = 1
114template<typename VecT>
115class Extrema<VecT, 1>
116{
117protected:
118 using Real = typename VecT::ValueType; // this works with both nanovdb and openvdb vectors
119 struct Pair
120 {
122 VecT vector;
123
124 Pair(Real s)// is only used by Extrema() default c-tor
125 : scalar(s)
126 , vector(s)
127 {
128 }
129 Pair(const VecT& v)
130 : scalar(v.lengthSqr())
131 , vector(v)
132 {
133 }
134 bool operator<(const Pair& rhs) const { return scalar < rhs.scalar; }
135 } mMin, mMax;
136 Extrema& add(const Pair& p)
137 {
138 if (p < mMin) {
139 mMin = p;
140 }
141 if (mMax < p) {
142 mMax = p;
143 }
144 return *this;
145 }
146
147public:
148 using ValueType = VecT;
150 : mMin(std::numeric_limits<Real>::max())
151 , mMax(std::numeric_limits<Real>::lowest())
152 {
153 }
154 Extrema(const VecT& v)
155 : mMin(v)
156 , mMax(v)
157 {
158 }
159 Extrema(const VecT& a, const VecT& b)
160 : mMin(a)
161 , mMax(b)
162 {
163 }
164 Extrema& min(const VecT& v)
165 {
166 Pair tmp(v);
167 if (tmp < mMin) {
168 mMin = tmp;
169 }
170 return *this;
171 }
172 Extrema& max(const VecT& v)
173 {
174 Pair tmp(v);
175 if (mMax < tmp) {
176 mMax = tmp;
177 }
178 return *this;
179 }
180 Extrema& add(const VecT& v) { return this->add(Pair(v)); }
181 Extrema& add(const VecT& v, uint64_t) { return this->add(Pair(v)); }
183 {
184 if (other.mMin < mMin) {
185 mMin = other.mMin;
186 }
187 if (mMax < other.mMax) {
188 mMax = other.mMax;
189 }
190 return *this;
191 }
192 const VecT& min() const { return mMin.vector; }
193 const VecT& max() const { return mMax.vector; }
194 operator bool() const { return !(mMax < mMin); }
195 static constexpr bool hasMinMax() { return !std::is_same<bool, Real>::value; }
196 static constexpr bool hasAverage() { return false; }
197 static constexpr bool hasStdDeviation() { return false; }
198 static constexpr size_t size() { return 0; }
199}; // Extrema<T, 1>
200
201//================================================================================================
202
203template<typename ValueT, int Rank = TensorTraits<ValueT>::Rank>
204class Stats;
205
206/// @brief This class computes statistics (minimum value, maximum
207/// value, mean, variance and standard deviation) of a population
208/// of floating-point values.
209///
210/// @details variance = Mean[ (X-Mean[X])^2 ] = Mean[X^2] - Mean[X]^2,
211/// standard deviation = sqrt(variance)
212///
213/// @note This class employs incremental computation and double precision.
214template<typename ValueT>
215class Stats<ValueT, 0> : public Extrema<ValueT, 0>
216{
217protected:
219 using RealT = double; // for accuracy the internal precission must be 64 bit floats
220 size_t mSize;
221 double mAvg, mAux;
222
223public:
224 using ValueType = ValueT;
226 : BaseT()
227 , mSize(0)
228 , mAvg(0.0)
229 , mAux(0.0)
230 {
231 }
232 Stats(const ValueT& val)
233 : BaseT(val)
234 , mSize(1)
235 , mAvg(RealT(val))
236 , mAux(0.0)
237 {
238 }
239 /// @brief Add a single sample
240 Stats& add(const ValueT& val)
241 {
242 BaseT::add(val);
243 mSize += 1;
244 const double delta = double(val) - mAvg;
245 mAvg += delta / double(mSize);
246 mAux += delta * (double(val) - mAvg);
247 return *this;
248 }
249 /// @brief Add @a n samples with constant value @a val.
250 Stats& add(const ValueT& val, uint64_t n)
251 {
252 const double denom = 1.0 / double(mSize + n);
253 const double delta = double(val) - mAvg;
254 mAvg += denom * delta * double(n);
255 mAux += denom * delta * delta * double(mSize) * double(n);
256 BaseT::add(val);
257 mSize += n;
258 return *this;
259 }
260
261 /// Add the samples from the other Stats instance.
263 {
264 if (other.mSize > 0) {
265 const double denom = 1.0 / double(mSize + other.mSize);
266 const double delta = other.mAvg - mAvg;
267 mAvg += denom * delta * double(other.mSize);
268 mAux += other.mAux + denom * delta * delta * double(mSize) * double(other.mSize);
269 BaseT::add(other);
270 mSize += other.mSize;
271 }
272 return *this;
273 }
274
275 static constexpr bool hasMinMax() { return !std::is_same<bool, ValueT>::value; }
276 static constexpr bool hasAverage() { return !std::is_same<bool, ValueT>::value; }
277 static constexpr bool hasStdDeviation() { return !std::is_same<bool, ValueT>::value; }
278
279 size_t size() const { return mSize; }
280
281 //@{
282 /// Return the arithmetic mean, i.e. average, value.
283 double avg() const { return mAvg; }
284 double mean() const { return mAvg; }
285 //@}
286
287 //@{
288 /// @brief Return the population variance.
289 ///
290 /// @note The unbiased sample variance = population variance * num/(num-1)
291 double var() const { return mSize < 2 ? 0.0 : mAux / double(mSize); }
292 double variance() const { return this->var(); }
293 //@}
294
295 //@{
296 /// @brief Return the standard deviation (=Sqrt(variance)) as
297 /// defined from the (biased) population variance.
298 double std() const { return sqrt(this->var()); }
299 double stdDev() const { return this->std(); }
300 //@}
301}; // end Stats<T, 0>
302
303/// @brief This class computes statistics (minimum value, maximum
304/// value, mean, variance and standard deviation) of a population
305/// of floating-point values.
306///
307/// @details variance = Mean[ (X-Mean[X])^2 ] = Mean[X^2] - Mean[X]^2,
308/// standard deviation = sqrt(variance)
309///
310/// @note This class employs incremental computation and double precision.
311template<typename ValueT>
312class Stats<ValueT, 1> : public Extrema<ValueT, 1>
313{
314protected:
316 using RealT = double; // for accuracy the internal precision must be 64 bit floats
317 size_t mSize;
318 double mAvg, mAux;
319
320public:
321 using ValueType = ValueT;
323 : BaseT()
324 , mSize(0)
325 , mAvg(0.0)
326 , mAux(0.0)
327 {
328 }
329 /// @brief Add a single sample
330 Stats& add(const ValueT& val)
331 {
332 typename BaseT::Pair tmp(val);
333 BaseT::add(tmp);
334 mSize += 1;
335 const double delta = tmp.scalar - mAvg;
336 mAvg += delta / double(mSize);
337 mAux += delta * (tmp.scalar - mAvg);
338 return *this;
339 }
340 /// @brief Add @a n samples with constant value @a val.
341 Stats& add(const ValueT& val, uint64_t n)
342 {
343 typename BaseT::Pair tmp(val);
344 const double denom = 1.0 / double(mSize + n);
345 const double delta = tmp.scalar - mAvg;
346 mAvg += denom * delta * double(n);
347 mAux += denom * delta * delta * double(mSize) * double(n);
348 BaseT::add(tmp);
349 mSize += n;
350 return *this;
351 }
352
353 /// Add the samples from the other Stats instance.
355 {
356 if (other.mSize > 0) {
357 const double denom = 1.0 / double(mSize + other.mSize);
358 const double delta = other.mAvg - mAvg;
359 mAvg += denom * delta * double(other.mSize);
360 mAux += other.mAux + denom * delta * delta * double(mSize) * double(other.mSize);
361 BaseT::add(other);
362 mSize += other.mSize;
363 }
364 return *this;
365 }
366
367 static constexpr bool hasMinMax() { return !std::is_same<bool, ValueT>::value; }
368 static constexpr bool hasAverage() { return !std::is_same<bool, ValueT>::value; }
369 static constexpr bool hasStdDeviation() { return !std::is_same<bool, ValueT>::value; }
370
371 size_t size() const { return mSize; }
372
373 //@{
374 /// Return the arithmetic mean, i.e. average, value.
375 double avg() const { return mAvg; }
376 double mean() const { return mAvg; }
377 //@}
378
379 //@{
380 /// @brief Return the population variance.
381 ///
382 /// @note The unbiased sample variance = population variance * num/(num-1)
383 double var() const { return mSize < 2 ? 0.0 : mAux / double(mSize); }
384 double variance() const { return this->var(); }
385 //@}
386
387 //@{
388 /// @brief Return the standard deviation (=Sqrt(variance)) as
389 /// defined from the (biased) population variance.
390 double std() const { return sqrt(this->var()); }
391 double stdDev() const { return this->std(); }
392 //@}
393}; // end Stats<T, 1>
394
395/// @brief No-op Stats class
396template<typename ValueT>
398{
399 using ValueType = ValueT;
401 NoopStats(const ValueT&) {}
402 NoopStats& add(const ValueT&) { return *this; }
403 NoopStats& add(const ValueT&, uint64_t) { return *this; }
404 NoopStats& add(const NoopStats&) { return *this; }
405 static constexpr size_t size() { return 0; }
406 static constexpr bool hasMinMax() { return false; }
407 static constexpr bool hasAverage() { return false; }
408 static constexpr bool hasStdDeviation() { return false; }
409}; // end NoopStats<T>
410
411//================================================================================================
412
413/// @brief Allows for the construction of NanoVDB grids without any dependecy
414template<typename GridT, typename StatsT = Stats<typename GridT::ValueType>>
416{
417 struct NodeStats;
418 using TreeT = typename GridT::TreeType;
419 using ValueT = typename TreeT::ValueType;
420 using BuildT = typename TreeT::BuildType;
421 using Node0 = typename TreeT::Node0; // leaf
422 using Node1 = typename TreeT::Node1; // lower
423 using Node2 = typename TreeT::Node2; // upper
424 using RootT = typename TreeT::Node3; // root
425 static_assert(std::is_same<ValueT, typename StatsT::ValueType>::value, "Mismatching type");
426 static constexpr bool DO_STATS = StatsT::hasMinMax() || StatsT::hasAverage() || StatsT::hasStdDeviation();
427
428 ValueT mDelta; // skip rendering of node if: node.max < -mDelta || node.min > mDelta
429
430 void process( GridT& );// process grid and all tree nodes
431 void process( TreeT& );// process Tree, root node and child nodes
432 void process( RootT& );// process root node and child nodes
433 NodeStats process( Node0& );// process leaf node
434
435 template<typename NodeT>
436 NodeStats process( NodeT& );// process internal node and child nodes
437
438 template<typename DataT, int Rank>
439 void setStats(DataT*, const Extrema<ValueT, Rank>&);
440 template<typename DataT, int Rank>
441 void setStats(DataT*, const Stats<ValueT, Rank>&);
442 template<typename DataT>
443 void setStats(DataT*, const NoopStats<ValueT>&) {}
444
445 template<typename T, typename FlagT>
446 typename std::enable_if<!std::is_floating_point<T>::value>::type
447 setFlag(const T&, const T&, FlagT& flag) const { flag &= ~FlagT(1); } // unset 1st bit to enable rendering
448
449 template<typename T, typename FlagT>
450 typename std::enable_if<std::is_floating_point<T>::value>::type
451 setFlag(const T& min, const T& max, FlagT& flag) const;
452
453public:
454 GridStats() = default;
455
456 void operator()(GridT& grid, ValueT delta = ValueT(0));
457
458}; // GridStats
459
460template<typename GridT, typename StatsT>
461struct GridStats<GridT, StatsT>::NodeStats
462{
463 StatsT stats;
464 //uint64_t activeCount;
466
467 NodeStats(): stats(), bbox() {}//activeCount(0), bbox() {};
468
470 {
471 stats.add( other.stats );// no-op for NoopStats?!
472 //activeCount += other.activeCount;
473 bbox[0].minComponent(other.bbox[0]);
474 bbox[1].maxComponent(other.bbox[1]);
475 return *this;
476 }
477};// GridStats::NodeStats
478
479//================================================================================================
480
481template<typename GridT, typename StatsT>
483{
484 mDelta = delta; // delta = voxel size for level sets, else 0
485 this->process( grid );
486}
487
488//================================================================================================
489
490template<typename GridT, typename StatsT>
491template<typename DataT, int Rank>
493 setStats(DataT* data, const Extrema<ValueT, Rank>& e)
494{
495 data->setMin(e.min());
496 data->setMax(e.max());
497}
498
499template<typename GridT, typename StatsT>
500template<typename DataT, int Rank>
501inline void GridStats<GridT, StatsT>::
502 setStats(DataT* data, const Stats<ValueT, Rank>& s)
503{
504 data->setMin(s.min());
505 data->setMax(s.max());
506 data->setAvg(s.avg());
507 data->setDev(s.std());
508}
509
510//================================================================================================
511
512template<typename GridT, typename StatsT>
513template<typename T, typename FlagT>
514inline typename std::enable_if<std::is_floating_point<T>::value>::type
515GridStats<GridT, StatsT>::
516 setFlag(const T& min, const T& max, FlagT& flag) const
517{
518 if (mDelta > 0 && (min > mDelta || max < -mDelta)) {// LS: min > dx || max < -dx
519 flag |= FlagT(1u);// set 1st bit to disable rendering
520 } else {
521 flag &= ~FlagT(1u);// unset 1st bit to enable rendering
522 }
523}
524
525//================================================================================================
526
527template<typename GridT, typename StatsT>
528void GridStats<GridT, StatsT>::process( GridT &grid )
529{
530 this->process( grid.tree() );// this processes tree, root and all nodes
531
532 // set world space AABB
533 auto& data = *grid.data();
534 const auto& indexBBox = grid.tree().root().bbox();
535 if (indexBBox.empty()) {
536 data.mWorldBBox = BBox<Vec3R>();
537 data.setBBoxOn(false);
538 } else {
539 // Note that below max is offset by one since CoordBBox.max is inclusive
540 // while bbox<Vec3R>.max is exclusive. However, min is inclusive in both
541 // CoordBBox and BBox<Vec3R>. This also guarantees that a grid with a single
542 // active voxel, does not have an empty world bbox! E.g. if a grid with a
543 // unit index-to-world transformation only contains the active voxel (0,0,0)
544 // then indeBBox = (0,0,0) -> (0,0,0) and then worldBBox = (0.0, 0.0, 0.0)
545 // -> (1.0, 1.0, 1.0). This is a consequence of the different definitions
546 // of index and world bounding boxes inherited from OpenVDB!
547 const Coord min = indexBBox[0];
548 const Coord max = indexBBox[1] + Coord(1);
549
550 auto& worldBBox = data.mWorldBBox;
551 const auto& map = grid.map();
552 worldBBox[0] = worldBBox[1] = map.applyMap(Vec3d(min[0], min[1], min[2]));
553 worldBBox.expand(map.applyMap(Vec3d(min[0], min[1], max[2])));
554 worldBBox.expand(map.applyMap(Vec3d(min[0], max[1], min[2])));
555 worldBBox.expand(map.applyMap(Vec3d(max[0], min[1], min[2])));
556 worldBBox.expand(map.applyMap(Vec3d(max[0], max[1], min[2])));
557 worldBBox.expand(map.applyMap(Vec3d(max[0], min[1], max[2])));
558 worldBBox.expand(map.applyMap(Vec3d(min[0], max[1], max[2])));
559 worldBBox.expand(map.applyMap(Vec3d(max[0], max[1], max[2])));
560 data.setBBoxOn(true);
561 }
562
563 // set bit flags
564 data.setMinMaxOn(StatsT::hasMinMax());
565 data.setAverageOn(StatsT::hasAverage());
566 data.setStdDeviationOn(StatsT::hasStdDeviation());
567} // GridStats::process( Grid )
568
569//================================================================================================
570
571template<typename GridT, typename StatsT>
572inline void GridStats<GridT, StatsT>::process( typename GridT::TreeType &tree )
573{
574 this->process( tree.root() );
575}
576
577//================================================================================================
578
579template<typename GridT, typename StatsT>
580void GridStats<GridT, StatsT>::process(RootT &root)
581{
582 using ChildT = Node2;
583 auto &data = *root.data();
584 if (data.mTableSize == 0) { // empty root node
585 data.mMinimum = data.mMaximum = data.mBackground;
586 data.mAverage = data.mStdDevi = 0;
587 //data.mActiveVoxelCount = 0;
588 data.mBBox = CoordBBox();
589 } else {
590 NodeStats total;
591 for (uint32_t i = 0; i < data.mTableSize; ++i) {
592 auto* tile = data.tile(i);
593 if (tile->isChild()) { // process child node
594 total.add( this->process( *data.getChild(tile) ) );
595 } else if (tile->state) { // active tile
596 //total.activeCount += ChildT::NUM_VALUES;
597 const Coord ijk = tile->origin();
598 total.bbox[0].minComponent(ijk);
599 total.bbox[1].maxComponent(ijk + Coord(ChildT::DIM - 1));
600 if (DO_STATS) { // resolved at compile time
601 total.stats.add(tile->value, ChildT::NUM_VALUES);
602 }
603 }
604 }
605 this->setStats(&data, total.stats);
606 if (total.bbox.empty()) {
607 std::cerr << "\nWarning: input tree only contained inactive root tiles!"
608 << "\nWhile not strictly an error it's rather suspicious!\n";
609 }
610 //data.mActiveVoxelCount = total.activeCount;
611 data.mBBox = total.bbox;
612 }
613} // GridStats::process( RootNode )
614
615//================================================================================================
616
617template<typename GridT, typename StatsT>
618template<typename NodeT>
619typename GridStats<GridT, StatsT>::NodeStats
620GridStats<GridT, StatsT>::process(NodeT &node)
621{
622 static_assert(is_same<NodeT,Node1>::value || is_same<NodeT,Node2>::value, "Incorrect node type");
623 using ChildT = typename NodeT::ChildNodeType;
624
625 NodeStats total;
626 auto* data = node.data();
627
628 // Serial processing of active tiles
629 if (const auto tileCount = data->mValueMask.countOn()) {
630 //total.activeCount = tileCount * ChildT::NUM_VALUES; // active tiles
631 for (auto it = data->mValueMask.beginOn(); it; ++it) {
632 if (DO_STATS) { // resolved at compile time
633 total.stats.add( data->mTable[*it].value, ChildT::NUM_VALUES );
634 }
635 const Coord ijk = node.offsetToGlobalCoord(*it);
636 total.bbox[0].minComponent(ijk);
637 total.bbox[1].maxComponent(ijk + Coord(int32_t(ChildT::DIM) - 1));
638 }
639 }
640
641 // Serial or parallel processing of child nodes
642 if (const size_t childCount = data->mChildMask.countOn()) {
643#ifndef NANOVDB_USE_TBB
644 for (auto it = data->mChildMask.beginOn(); it; ++it) {
645 total.add( this->process( *data->getChild(*it) ) );
646 }
647#else
648 std::unique_ptr<ChildT*[]> childNodes(new ChildT*[childCount]);
649 ChildT **ptr = childNodes.get();
650 for (auto it = data->mChildMask.beginOn(); it; ++it) {
651 *ptr++ = data->getChild( *it );
652 }
653 using RangeT = tbb::blocked_range<size_t>;
654 total.add( tbb::parallel_reduce(RangeT(0, childCount), NodeStats(),
655 [&](const RangeT &r, NodeStats local)->NodeStats {
656 for(size_t i=r.begin(); i!=r.end(); ++i){
657 local.add( this->process( *childNodes[i] ) );
658 }
659 return local;},
660 [](NodeStats a, const NodeStats &b)->NodeStats { return a.add( b ); }
661 ));
662#endif
663 }
664
665 data->mBBox = total.bbox;
666 if (total.bbox.empty()) {
667 data->mFlags |= uint32_t(1); // set 1st bit on to disable rendering of node
668 data->mFlags &= ~uint32_t(2); // set 2nd bit off since node does not contain active values
669 } else {
670 data->mFlags |= uint32_t(2); // set 2nd bit on since node contains active values
671 if (DO_STATS) { // resolved at compile time
672 this->setStats(data, total.stats);
673 this->setFlag(data->mMinimum, data->mMaximum, data->mFlags);
674 }
675 }
676 return total;
677} // GridStats::process( InternalNode )
678
679//================================================================================================
680
681template<typename GridT, typename StatsT>
682typename GridStats<GridT, StatsT>::NodeStats
683GridStats<GridT, StatsT>::process(Node0 &leaf)
684{
685 static_assert(Node0::SIZE == 512u, "Invalid size of leaf nodes");
686 NodeStats local;
687 auto *data = leaf.data();
688 if (auto activeCount = data->mValueMask.countOn()) {
689 //data->mFlags |= uint8_t(2); // sets 2nd bit on since leaf contains active voxel
690 //local.activeCount += activeCount;
691 leaf.updateBBox(); // optionally update active bounding box (updates data->mFlags)
692 local.bbox[0] = local.bbox[1] = data->mBBoxMin;
693 local.bbox[1] += Coord(data->mBBoxDif[0], data->mBBoxDif[1], data->mBBoxDif[2]);
694 if (DO_STATS) { // resolved at compile time
695 for (auto it = data->mValueMask.beginOn(); it; ++it) {
696 local.stats.add(data->getValue(*it));
697 }
698 this->setStats(data, local.stats);
699 this->setFlag(data->getMin(), data->getMax(), data->mFlags);
700 }
701 } else {
702 data->mFlags &= ~uint8_t(2); // sets 2nd bit off since leaf has no bbox of active active values
703 }
704 return local;
705} // GridStats::process( LeafNode )
706
707//================================================================================================
708
709template<typename BuildT>
711{
712 using GridT = NanoGrid<BuildT>;
713 using ValueT = typename GridT::ValueType;
714 if (mode == StatsMode::Disable) {
715 return;
716 } else if (mode == StatsMode::BBox || std::is_same<bool, ValueT>::value) {
718 stats(grid);
719 } else if (mode == StatsMode::MinMax) {
721 stats(grid);
722 } else if (mode == StatsMode::All) {
724 stats(grid);
725 } else {
726 throw std::runtime_error("gridStats: Unsupported statistics mode.");
727 }
728}
729
730//================================================================================================
731
732namespace {
733
734// returns a bitmask (of size 32^3 or 16^3) that marks all the entries
735// in a node table that intersects with the specified bounding box.
736template<typename NodeT>
737Mask<NodeT::LOG2DIM> getBBoxMask(const CoordBBox &bbox, const NodeT* node)
738{
739 Mask<NodeT::LOG2DIM> mask;// typically 32^3 or 16^3 bit mask
740 auto b = CoordBBox::createCube(node->origin(), node->dim());
741 assert( bbox.hasOverlap(b) );
742 if ( bbox.isInside(b) ) {
743 mask.setOn();//node is completely inside the bbox so early out
744 } else {
745 b.intersect(bbox);// trim bounding box
746 // transform bounding box from global to local coordinates
747 b.min() &= NodeT::DIM-1u;
748 b.min() >>= NodeT::ChildNodeType::TOTAL;
749 b.max() &= NodeT::DIM-1u;
750 b.max() >>= NodeT::ChildNodeType::TOTAL;
751 assert( !b.empty() );
752 auto it = b.begin();// iterates over all the child nodes or tiles that intersects bbox
753 for (const Coord& ijk = *it; it; ++it) {
754 mask.setOn(ijk[2] + (ijk[1] << NodeT::LOG2DIM) + (ijk[0] << 2*NodeT::LOG2DIM));
755 }
756 }
757 return mask;
758}
759}// end of unnamed namespace
760
761/// @brief return the extrema of all the values in a grid that
762/// intersects the specified bounding box.
763template<typename BuildT>
764Extrema<typename NanoGrid<BuildT>::ValueType>
765getExtrema(const NanoGrid<BuildT>& grid, const CoordBBox &bbox)
766{
767 using GridT = NanoGrid<BuildT>;
768 using ValueT = typename GridT::ValueType;
769 using TreeT = typename GridTree<GridT>::type;
770 using RootT = typename NodeTrait<TreeT, 3>::type;// root node
771 using Node2 = typename NodeTrait<TreeT, 2>::type;// upper internal node
772 using Node1 = typename NodeTrait<TreeT, 1>::type;// lower internal node
773 using Node0 = typename NodeTrait<TreeT, 0>::type;// leaf node
774
775 Extrema<ValueT> extrema;
776 const RootT &root = grid.tree().root();
777 const auto &bbox3 = root.bbox();
778 if (bbox.isInside(bbox3)) {// bbox3 is contained inside bbox
779 extrema.min(root.minimum());
780 extrema.max(root.maximum());
781 extrema.add(root.background());
782 } else if (bbox.hasOverlap(bbox3)) {
783 const auto *data3 = root.data();
784 for (uint32_t i=0; i<data3->mTableSize; ++i) {
785 const auto *tile = data3->tile(i);
786 CoordBBox bbox2 = CoordBBox::createCube(tile->origin(), Node2::dim());
787 if (!bbox.hasOverlap(bbox2)) continue;
788 if (tile->isChild()) {
789 const Node2 *node2 = data3->getChild(tile);
790 if (bbox.isInside(bbox2)) {
791 extrema.min(node2->minimum());
792 extrema.max(node2->maximum());
793 } else {// partial intersections at level 2
794 auto *data2 = node2->data();
795 const auto bboxMask2 = getBBoxMask(bbox, node2);
796 for (auto it2 = bboxMask2.beginOn(); it2; ++it2) {
797 if (data2->mChildMask.isOn(*it2)) {
798 const Node1* node1 = data2->getChild(*it2);
799 CoordBBox bbox1 = CoordBBox::createCube(node1->origin(), Node1::dim());
800 if (bbox.isInside(bbox1)) {
801 extrema.min(node1->minimum());
802 extrema.max(node1->maximum());
803 } else {// partial intersection at level 1
804 auto *data1 = node1->data();
805 const auto bboxMask1 = getBBoxMask(bbox, node1);
806 for (auto it1 = bboxMask1.beginOn(); it1; ++it1) {
807 if (data1->mChildMask.isOn(*it1)) {
808 const Node0* node0 = data1->getChild(*it1);
809 CoordBBox bbox0 = CoordBBox::createCube(node0->origin(), Node0::dim());
810 if (bbox.isInside(bbox0)) {
811 extrema.min(node0->minimum());
812 extrema.max(node0->maximum());
813 } else {// partial intersection at level 0
814 auto *data0 = node0->data();
815 const auto bboxMask0 = getBBoxMask(bbox, node0);
816 for (auto it0 = bboxMask0.beginOn(); it0; ++it0) {
817 extrema.add(data0->getValue(*it0));
818 }
819 }// end partial intersection at level 0
820 } else {// tile at level 1
821 extrema.add(data1->mTable[*it1].value);
822 }
823 }
824 }// end of partial intersection at level 1
825 } else {// tile at level 2
826 extrema.add(data2->mTable[*it2].value);
827 }
828 }// loop over tiles and nodes at level 2
829 }// end of partial intersection at level 1
830 } else {// tile at root level
831 extrema.add(tile->value);
832 }
833 }// loop over root table
834 } else {// bbox does not overlap the grid
835 extrema.add(root.background());
836 }
837 return extrema;
838}// getExtrema
839
840
841} // namespace nanovdb
842
843#endif // NANOVDB_GRIDSTATS_H_HAS_BEEN_INCLUDED
A unified wrapper for tbb::parallel_for and a naive std::thread fallback.
Implements a light-weight self-contained VDB data-structure in a single file! In other words,...
Custom Range class that is compatible with the tbb::blocked_range classes.
Definition DenseGrid.h:402
uint8_t * data()
Returns a non-const pointer to the data.
Definition DenseGrid.h:432
const DenseGrid< ValueT > * grid() const
Returns a const pointer to the NanoVDB grid encoded in the DenseGridHandle.
Definition DenseGrid.h:447
Template specialization of Extrema on scalar value types, i.e. rank = 0.
Definition GridStats.h:56
static constexpr size_t size()
Definition GridStats.h:110
Extrema & min(const ValueT &v)
Definition GridStats.h:77
Extrema & add(const ValueT &v, uint64_t)
Definition GridStats.h:97
Extrema(const ValueT &v)
Definition GridStats.h:67
Extrema()
Definition GridStats.h:62
Extrema & add(const Extrema &other)
Definition GridStats.h:98
ValueT ValueType
Definition GridStats.h:61
const ValueT & max() const
Definition GridStats.h:105
Extrema & max(const ValueT &v)
Definition GridStats.h:84
ValueT mMax
Definition GridStats.h:58
Extrema & add(const ValueT &v)
Definition GridStats.h:91
static constexpr bool hasAverage()
Definition GridStats.h:108
static constexpr bool hasMinMax()
Definition GridStats.h:107
Extrema(const ValueT &a, const ValueT &b)
Definition GridStats.h:72
const ValueT & min() const
Definition GridStats.h:104
static constexpr bool hasStdDeviation()
Definition GridStats.h:109
Extrema & min(const VecT &v)
Definition GridStats.h:164
Extrema & add(const VecT &v, uint64_t)
Definition GridStats.h:181
static constexpr size_t size()
Definition GridStats.h:198
Extrema(const VecT &a, const VecT &b)
Definition GridStats.h:159
Extrema()
Definition GridStats.h:149
Extrema & add(const Pair &p)
Definition GridStats.h:136
Extrema & max(const VecT &v)
Definition GridStats.h:172
const VecT & min() const
Definition GridStats.h:192
Extrema & add(const Extrema &other)
Definition GridStats.h:182
Extrema & add(const VecT &v)
Definition GridStats.h:180
typename VecT::ValueType Real
Definition GridStats.h:118
static constexpr bool hasAverage()
Definition GridStats.h:196
static constexpr bool hasMinMax()
Definition GridStats.h:195
Extrema(const VecT &v)
Definition GridStats.h:154
static constexpr bool hasStdDeviation()
Definition GridStats.h:197
const VecT & max() const
Definition GridStats.h:193
VecT ValueType
Definition GridStats.h:148
Definition GridStats.h:51
Allows for the construction of NanoVDB grids without any dependecy.
Definition GridStats.h:416
void operator()(GridT &grid, ValueT delta=ValueT(0))
Definition GridStats.h:482
double var() const
Return the population variance.
Definition GridStats.h:291
size_t size() const
Definition GridStats.h:279
Stats(const ValueT &val)
Definition GridStats.h:232
double mean() const
Definition GridStats.h:284
ValueT ValueType
Definition GridStats.h:224
size_t mSize
Definition GridStats.h:220
Stats & add(const ValueT &val, uint64_t n)
Add n samples with constant value val.
Definition GridStats.h:250
double variance() const
Definition GridStats.h:292
Stats & add(const ValueT &val)
Add a single sample.
Definition GridStats.h:240
double stdDev() const
Definition GridStats.h:299
static constexpr bool hasAverage()
Definition GridStats.h:276
static constexpr bool hasMinMax()
Definition GridStats.h:275
Stats & add(const Stats &other)
Add the samples from the other Stats instance.
Definition GridStats.h:262
double mAux
Definition GridStats.h:221
double std() const
Return the standard deviation (=Sqrt(variance)) as defined from the (biased) population variance.
Definition GridStats.h:298
static constexpr bool hasStdDeviation()
Definition GridStats.h:277
double avg() const
Return the arithmetic mean, i.e. average, value.
Definition GridStats.h:283
Stats()
Definition GridStats.h:225
double var() const
Return the population variance.
Definition GridStats.h:383
size_t size() const
Definition GridStats.h:371
double mean() const
Definition GridStats.h:376
ValueT ValueType
Definition GridStats.h:321
size_t mSize
Definition GridStats.h:317
Stats & add(const ValueT &val, uint64_t n)
Add n samples with constant value val.
Definition GridStats.h:341
double variance() const
Definition GridStats.h:384
Stats & add(const ValueT &val)
Add a single sample.
Definition GridStats.h:330
double stdDev() const
Definition GridStats.h:391
static constexpr bool hasAverage()
Definition GridStats.h:368
static constexpr bool hasMinMax()
Definition GridStats.h:367
Stats & add(const Stats &other)
Add the samples from the other Stats instance.
Definition GridStats.h:354
double mAux
Definition GridStats.h:318
double std() const
Return the standard deviation (=Sqrt(variance)) as defined from the (biased) population variance.
Definition GridStats.h:390
static constexpr bool hasStdDeviation()
Definition GridStats.h:369
double avg() const
Return the arithmetic mean, i.e. average, value.
Definition GridStats.h:375
Stats()
Definition GridStats.h:322
Definition GridStats.h:204
void add(double val)
Add a single sample.
Definition Stats.h:102
double min() const
Return the minimum value.
Definition Stats.h:121
double max() const
Return the maximum value.
Definition Stats.h:124
Definition NanoVDB.h:208
Extrema< typename NanoGrid< BuildT >::ValueType > getExtrema(const NanoGrid< BuildT > &grid, const CoordBBox &bbox)
return the extrema of all the values in a grid that intersects the specified bounding box.
Definition GridStats.h:765
BBox< Coord > CoordBBox
Definition NanoVDB.h:1811
StatsMode
Grid flags which indicate what extra information is present in the grid buffer.
Definition GridStats.h:32
Vec3< double > Vec3d
Definition NanoVDB.h:1290
void gridStats(NanoGrid< BuildT > &grid, StatsMode mode=StatsMode::Default)
Re-computes the min/max, stats and bbox information for an existing NanoVDB Grid.
Definition GridStats.h:710
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
Definition Coord.h:589
Definition NanoVDB.h:1657
VecT vector
Definition GridStats.h:122
Pair(Real s)
Definition GridStats.h:124
Pair(const VecT &v)
Definition GridStats.h:129
bool operator<(const Pair &rhs) const
Definition GridStats.h:134
Real scalar
Definition GridStats.h:121
Definition GridStats.h:462
NodeStats()
Definition GridStats.h:467
StatsT stats
Definition GridStats.h:463
CoordBBox bbox
Definition GridStats.h:465
NodeStats & add(const NodeStats &other)
Definition GridStats.h:469
typename GridT::TreeType type
Definition NanoVDB.h:2787
Struct to derive node type from its level in a given grid, tree or root while preserving constness.
Definition NanoVDB.h:2345
No-op Stats class.
Definition GridStats.h:398
NoopStats(const ValueT &)
Definition GridStats.h:401
static constexpr size_t size()
Definition GridStats.h:405
NoopStats & add(const ValueT &, uint64_t)
Definition GridStats.h:403
NoopStats & add(const NoopStats &)
Definition GridStats.h:404
ValueT ValueType
Definition GridStats.h:399
NoopStats & add(const ValueT &)
Definition GridStats.h:402
static constexpr bool hasAverage()
Definition GridStats.h:407
static constexpr bool hasMinMax()
Definition GridStats.h:406
NoopStats()
Definition GridStats.h:400
static constexpr bool hasStdDeviation()
Definition GridStats.h:408
static constexpr bool value
Definition NanoVDB.h:358