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