OpenVDB  12.0.0
PointStatisticsImpl.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: Apache-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 
12 namespace openvdb {
14 namespace OPENVDB_VERSION_NAME {
15 namespace points {
16 
17 /// @cond OPENVDB_DOCS_INTERNAL
18 
19 namespace 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
24 template <typename ValueT>
25 struct 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
47 template <typename VecT, bool MagResult = true>
48 struct 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
62 template <typename VecT>
63 struct 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
111 template <typename VecT>
112 struct 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 
132 template <typename ValueT,
133  typename CodecT,
134  typename FilterT,
135  typename ExtentOp,
136  typename PointDataTreeT>
137 bool 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) {
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 
231 template <typename ValueT,
232  typename CodecT,
233  typename FilterT,
234  typename PointDataTreeT,
235  typename std::enable_if<ValueTraits<ValueT>::IsVec, int>::type = 0>
236 bool 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 
252 template <typename ValueT,
253  typename CodecT,
254  typename FilterT,
255  typename PointDataTreeT,
256  typename std::enable_if<!ValueTraits<ValueT>::IsVec, int>::type = 0>
257 bool 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 
277 template <typename ValueT,
278  typename CodecT,
279  typename FilterT,
280  typename PointDataTreeT>
281 bool 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 
293 template <typename ValueT,
294  typename CodecT,
295  typename FilterT,
296  typename PointDataTreeT,
297  typename ResultTreeT>
298 bool 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  OPENVDB_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) {
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  OPENVDB_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 
404 template <typename ValueT,
405  typename CodecT,
406  typename FilterT,
407  typename PointDataTreeT,
408  typename ResultTreeT>
409 bool 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) {
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  OPENVDB_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 
493 template <typename ValueT,
494  typename CodecT,
495  typename FilterT,
496  typename PointDataTreeT>
497 std::pair<ValueT, ValueT>
498 evalMinMax(const PointDataTreeT& points,
499  const std::string& attribute,
500  const FilterT& filter)
501 {
502  std::pair<ValueT, ValueT> results {
503  zeroVal<ValueT>(), zeroVal<ValueT>()
504  };
505  evalMinMax<ValueT, CodecT, FilterT, PointDataTreeT>
506  (points, attribute, results.first, results.second, filter);
507  return results;
508 }
509 
510 template <typename ValueT,
511  typename CodecT,
512  typename FilterT,
513  typename PointDataTreeT>
515 evalAverage(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 
525 template <typename ValueT,
526  typename CodecT,
527  typename FilterT,
528  typename PointDataTreeT>
530 accumulate(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
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
typename TypeT< 64ul >::type Highest
Definition: Types.h:371
bool isUniform() const
Definition: AttributeArray.h:2087
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
Definition: AttributeArray.h:763
Index32 Index
Definition: Types.h:54
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:110
OutGridT XformOp & op
Definition: ValueTransformer.h:139
Definition: IndexIterator.h:44
size_t leafCount() const
Return the number of leaf nodes.
Definition: LeafManager.h:288
#define OPENVDB_ASSERT(X)
Definition: Assert.h:41
LeafRange leafRange(size_t grainsize=1) const
Return a TBB-compatible LeafRange.
Definition: LeafManager.h:346
LeafType & leaf(size_t leafIdx) const
Return a pointer to the leaf node at index leafIdx in the array.
Definition: LeafManager.h:319
Definition: Exceptions.h:13
Index size() const
Definition: AttributeArray.h:786
void foreach(const LeafOp &op, bool threaded=true, size_t grainSize=1)
Threaded method that applies a user-supplied functor to each leaf node in the LeafManager.
Definition: LeafManager.h:484
SubT Type
Definition: Types.h:320
This class manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional auxil...
Definition: LeafManager.h:85
Vec2< T > maxComponent(const Vec2< T > &v1, const Vec2< T > &v2)
Return component-wise maximum of the two vectors.
Definition: Vec2.h:513
Vec2< T > minComponent(const Vec2< T > &v1, const Vec2< T > &v2)
Return component-wise minimum of the two vectors.
Definition: Vec2.h:504
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:106
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
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:121
typename T::ValueType ElementType
Definition: Types.h:302
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:218