GCC Code Coverage Report


Directory: ./
File: openvdb/openvdb/tools/Statistics.h
Date: 2022-07-25 17:40:05
Exec Total Coverage
Lines: 48 48 100.0%
Functions: 39 40 97.5%
Branches: 77 138 55.8%

Line Branch Exec Source
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 //
4 /// @file Statistics.h
5 ///
6 /// @brief Functions to efficiently compute histograms, extrema
7 /// (min/max) and statistics (mean, variance, etc.) of grid values
8
9 #ifndef OPENVDB_TOOLS_STATISTICS_HAS_BEEN_INCLUDED
10 #define OPENVDB_TOOLS_STATISTICS_HAS_BEEN_INCLUDED
11
12 #include <openvdb/Types.h>
13 #include <openvdb/Exceptions.h>
14 #include <openvdb/math/Stats.h>
15 #include "ValueTransformer.h"
16
17
18 namespace openvdb {
19 OPENVDB_USE_VERSION_NAMESPACE
20 namespace OPENVDB_VERSION_NAME {
21 namespace tools {
22
23 /// @brief Iterate over a scalar grid and compute a histogram of the values
24 /// of the voxels that are visited, or iterate over a vector-valued grid
25 /// and compute a histogram of the magnitudes of the vectors.
26 /// @param iter an iterator over the values of a grid or its tree
27 /// (@c Grid::ValueOnCIter, @c Tree::ValueOffIter, etc.)
28 /// @param minVal the smallest value that can be added to the histogram
29 /// @param maxVal the largest value that can be added to the histogram
30 /// @param numBins the number of histogram bins
31 /// @param threaded if true, iterate over the grid in parallel
32 template<typename IterT>
33 inline math::Histogram
34 histogram(const IterT& iter, double minVal, double maxVal,
35 size_t numBins = 10, bool threaded = true);
36
37 /// @brief Iterate over a scalar grid and compute extrema (min/max) of the
38 /// values of the voxels that are visited, or iterate over a vector-valued grid
39 /// and compute extrema of the magnitudes of the vectors.
40 /// @param iter an iterator over the values of a grid or its tree
41 /// (@c Grid::ValueOnCIter, @c Tree::ValueOffIter, etc.)
42 /// @param threaded if true, iterate over the grid in parallel
43 template<typename IterT>
44 inline math::Extrema
45 extrema(const IterT& iter, bool threaded = true);
46
47 /// @brief Iterate over a scalar grid and compute statistics (mean, variance, etc.)
48 /// of the values of the voxels that are visited, or iterate over a vector-valued grid
49 /// and compute statistics of the magnitudes of the vectors.
50 /// @param iter an iterator over the values of a grid or its tree
51 /// (@c Grid::ValueOnCIter, @c Tree::ValueOffIter, etc.)
52 /// @param threaded if true, iterate over the grid in parallel
53 template<typename IterT>
54 inline math::Stats
55 statistics(const IterT& iter, bool threaded = true);
56
57 /// @brief Iterate over a grid and compute extrema (min/max) of
58 /// the values produced by applying the given functor at each voxel that is visited.
59 /// @param iter an iterator over the values of a grid or its tree
60 /// (@c Grid::ValueOnCIter, @c Tree::ValueOffIter, etc.)
61 /// @param op a functor of the form <tt>void op(const IterT&, math::Stats&)</tt>,
62 /// where @c IterT is the type of @a iter, that inserts zero or more
63 /// floating-point values into the provided @c math::Stats object
64 /// @param threaded if true, iterate over the grid in parallel
65 /// @note When @a threaded is true, each thread gets its own copy of the functor.
66 ///
67 /// @par Example:
68 /// Compute statistics of just the active and positive-valued voxels of a scalar,
69 /// floating-point grid.
70 /// @code
71 /// struct Local {
72 /// static inline
73 /// void addIfPositive(const FloatGrid::ValueOnCIter& iter, math::Extrema& ex)
74 /// {
75 /// const float f = *iter;
76 /// if (f > 0.0) {
77 /// if (iter.isVoxelValue()) ex.add(f);
78 /// else ex.add(f, iter.getVoxelCount());
79 /// }
80 /// }
81 /// };
82 /// FloatGrid grid = ...;
83 /// math::Extrema stats =
84 /// tools::extrema(grid.cbeginValueOn(), Local::addIfPositive, /*threaded=*/true);
85 /// @endcode
86 template<typename IterT, typename ValueOp>
87 inline math::Extrema
88 extrema(const IterT& iter, const ValueOp& op, bool threaded);
89
90 /// @brief Iterate over a grid and compute statistics (mean, variance, etc.) of
91 /// the values produced by applying the given functor at each voxel that is visited.
92 /// @param iter an iterator over the values of a grid or its tree
93 /// (@c Grid::ValueOnCIter, @c Tree::ValueOffIter, etc.)
94 /// @param op a functor of the form <tt>void op(const IterT&, math::Stats&)</tt>,
95 /// where @c IterT is the type of @a iter, that inserts zero or more
96 /// floating-point values into the provided @c math::Stats object
97 /// @param threaded if true, iterate over the grid in parallel
98 /// @note When @a threaded is true, each thread gets its own copy of the functor.
99 ///
100 /// @par Example:
101 /// Compute statistics of just the active and positive-valued voxels of a scalar,
102 /// floating-point grid.
103 /// @code
104 /// struct Local {
105 /// static inline
106 /// void addIfPositive(const FloatGrid::ValueOnCIter& iter, math::Stats& stats)
107 /// {
108 /// const float f = *iter;
109 /// if (f > 0.0) {
110 /// if (iter.isVoxelValue()) stats.add(f);
111 /// else stats.add(f, iter.getVoxelCount());
112 /// }
113 /// }
114 /// };
115 /// FloatGrid grid = ...;
116 /// math::Stats stats =
117 /// tools::statistics(grid.cbeginValueOn(), Local::addIfPositive, /*threaded=*/true);
118 /// @endcode
119 template<typename IterT, typename ValueOp>
120 inline math::Stats
121 statistics(const IterT& iter, const ValueOp& op, bool threaded);
122
123
124 /// @brief Iterate over a grid and compute statistics (mean, variance, etc.)
125 /// of the values produced by applying a given operator (see math/Operators.h)
126 /// at each voxel that is visited.
127 /// @param iter an iterator over the values of a grid or its tree
128 /// (@c Grid::ValueOnCIter, @c Tree::ValueOffIter, etc.)
129 /// @param op an operator object with a method of the form
130 /// <tt>double result(Accessor&, const Coord&)</tt>
131 /// @param threaded if true, iterate over the grid in parallel
132 /// @note World-space operators, whose @c result() methods are of the form
133 /// <tt>double result(const Map&, Accessor&, const Coord&)</tt>, must be wrapped
134 /// in a math::MapAdapter.
135 /// @note Vector-valued operators like math::Gradient must be wrapped in an adapter
136 /// such as math::OpMagnitude.
137 ///
138 /// @par Example:
139 /// Compute statistics of the magnitude of the gradient at the active voxels of
140 /// a scalar, floating-point grid. (Note the use of the math::MapAdapter and
141 /// math::OpMagnitude adapters.)
142 /// @code
143 /// FloatGrid grid = ...;
144 ///
145 /// // Assume that we know that the grid has a uniform scale map.
146 /// using MapType = math::UniformScaleMap;
147 /// // Specify a world-space gradient operator that uses first-order differencing.
148 /// using GradientOp = math::Gradient<MapType, math::FD_1ST>;
149 /// // Wrap the operator with an adapter that computes the magnitude of the gradient.
150 /// using MagnitudeOp = math::OpMagnitude<GradientOp, MapType>;
151 /// // Wrap the operator with an adapter that associates a map with it.
152 /// using CompoundOp = math::MapAdapter<MapType, GradientOp, double>;
153 ///
154 /// if (MapType::Ptr map = grid.constTransform().constMap<MapType>()) {
155 /// math::Stats stats = tools::opStatistics(grid.cbeginValueOn(), CompoundOp(*map));
156 /// }
157 /// @endcode
158 ///
159 /// @par Example:
160 /// Compute statistics of the divergence at the active voxels of a vector-valued grid.
161 /// @code
162 /// Vec3SGrid grid = ...;
163 ///
164 /// // Assume that we know that the grid has a uniform scale map.
165 /// using MapType = math::UniformScaleMap;
166 /// // Specify a world-space divergence operator that uses first-order differencing.
167 /// using DivergenceOp = math::Divergence<MapType, math::FD_1ST>;
168 /// // Wrap the operator with an adapter that associates a map with it.
169 /// using CompoundOp = math::MapAdapter<MapType, DivergenceOp, double>;
170 ///
171 /// if (MapType::Ptr map = grid.constTransform().constMap<MapType>()) {
172 /// math::Stats stats = tools::opStatistics(grid.cbeginValueOn(), CompoundOp(*map));
173 /// }
174 /// @endcode
175 ///
176 /// @par Example:
177 /// As above, but computing the divergence in index space.
178 /// @code
179 /// Vec3SGrid grid = ...;
180 ///
181 /// // Specify an index-space divergence operator that uses first-order differencing.
182 /// using DivergenceOp = math::ISDivergence<math::FD_1ST>;
183 ///
184 /// math::Stats stats = tools::opStatistics(grid.cbeginValueOn(), DivergenceOp());
185 /// @endcode
186 template<typename OperatorT, typename IterT>
187 inline math::Stats
188 opStatistics(const IterT& iter, const OperatorT& op = OperatorT(), bool threaded = true);
189
190 /// @brief Same as opStatistics except it returns a math::Extrema vs a math::Stats
191 template<typename OperatorT, typename IterT>
192 inline math::Extrema
193 opExtrema(const IterT& iter, const OperatorT& op = OperatorT(), bool threaded = true);
194
195 ////////////////////////////////////////
196
197 /// @cond OPENVDB_DOCS_INTERNAL
198
199 namespace stats_internal {
200
201 /// @todo This traits class is needed because tree::TreeValueIteratorBase uses
202 /// the name ValueT for the type of the value to which the iterator points,
203 /// whereas node-level iterators use the name ValueType.
204 template<typename IterT, typename AuxT = void>
205 struct IterTraits {
206 using ValueType = typename IterT::ValueType;
207 };
208
209 template<typename TreeT, typename ValueIterT>
210 struct IterTraits<tree::TreeValueIteratorBase<TreeT, ValueIterT> > {
211 using ValueType = typename tree::TreeValueIteratorBase<TreeT, ValueIterT>::ValueT;
212 };
213
214
215 // Helper class to compute a scalar value from either a scalar or a vector value
216 // (the latter by computing the vector's magnitude)
217 template<typename T, bool IsVector> struct GetValImpl;
218
219 template<typename T>
220 struct GetValImpl<T, /*IsVector=*/false> {
221
10/12
✓ Branch 0 taken 241453 times.
✓ Branch 1 taken 98515 times.
✓ Branch 2 taken 20413 times.
✓ Branch 3 taken 73124 times.
✓ Branch 4 taken 262771 times.
✓ Branch 5 taken 510936 times.
✓ Branch 6 taken 2197 times.
✓ Branch 7 taken 6591 times.
✓ Branch 8 taken 680661 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 8788 times.
✗ Branch 11 not taken.
1905449 static inline double get(const T& val) { return double(val); }
222 };
223
224 template<typename T>
225 struct GetValImpl<T, /*IsVector=*/true> {
226
6/8
✓ Branch 0 taken 170165 times.
✓ Branch 1 taken 510495 times.
✓ Branch 2 taken 2197 times.
✓ Branch 3 taken 6591 times.
✓ Branch 4 taken 680660 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 8788 times.
✗ Branch 7 not taken.
1378896 static inline double get(const T& val) { return val.length(); }
227 };
228
229
230 // Helper class to compute a scalar value from a tree or node iterator
231 // that points to a value in either a scalar or a vector grid, and to
232 // add that value to a math::Stats object.
233 template<typename IterT, typename StatsT>
234 struct GetVal
235 {
236 using ValueT = typename IterTraits<IterT>::ValueType;
237 using ImplT = GetValImpl<ValueT, VecTraits<ValueT>::IsVec>;
238
239
2/2
✓ Branch 0 taken 4517489 times.
✓ Branch 1 taken 183127 times.
9401232 inline void operator()(const IterT& iter, StatsT& stats) const {
240 6312336 if (iter.isVoxelValue()) stats.add(ImplT::get(*iter));
241 331102 else stats.add(ImplT::get(*iter), iter.getVoxelCount());
242 9401232 }
243 };
244
245 // Helper class to accumulate scalar voxel values or vector voxel magnitudes
246 // into a math::Stats object
247 template<typename IterT, typename ValueOp, typename StatsT>
248
4/8
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
34 struct StatsOp
249 {
250
5/10
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2 times.
✗ Branch 14 not taken.
15 StatsOp(const ValueOp& op): getValue(op) {}
251
252 // Accumulate voxel and tile values into this functor's Stats object.
253 6079511 inline void operator()(const IterT& iter) { getValue(iter, stats); }
254
255 // Accumulate another functor's Stats object into this functor's.
256
6/13
✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 12 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 11 taken 9 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✓ Branch 14 taken 1 times.
59 inline void join(StatsOp& other) { stats.add(other.stats); }
257
258 StatsT stats;
259 ValueOp getValue;
260 };
261
262
263 // Helper class to accumulate scalar voxel values or vector voxel magnitudes
264 // into a math::Histogram object
265 template<typename IterT, typename ValueOp>
266
4/8
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✓ Branch 14 taken 3 times.
✗ Branch 15 not taken.
42 struct HistOp
267 {
268 5 HistOp(const ValueOp& op, double vmin, double vmax, size_t bins):
269 5 hist(vmin, vmax, bins), getValue(op)
270 {}
271
272 // Accumulate voxel and tile values into this functor's Histogram object.
273 1378897 inline void operator()(const IterT& iter) { getValue(iter, hist); }
274
275 // Accumulate another functor's Histogram object into this functor's.
276 16 inline void join(HistOp& other) { hist.add(other.hist); }
277
278 math::Histogram hist;
279 ValueOp getValue;
280 };
281
282
283 // Helper class to apply an operator such as math::Gradient or math::Laplacian
284 // to voxels and accumulate the scalar results or the magnitudes of vector results
285 // into a math::Stats object
286 template<typename IterT, typename OpT, typename StatsT>
287
8/16
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 7 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 6 times.
✗ Branch 11 not taken.
✓ Branch 15 taken 2 times.
✗ Branch 16 not taken.
✓ Branch 24 taken 2 times.
✗ Branch 25 not taken.
✓ Branch 31 taken 4 times.
✗ Branch 32 not taken.
✓ Branch 36 taken 2 times.
✗ Branch 37 not taken.
105 struct MathOp
288 {
289 using TreeT = typename IterT::TreeT;
290 using ValueT = typename TreeT::ValueType;
291 using ConstAccessor = typename tree::ValueAccessor<const TreeT>;
292
293 // Each thread gets its own accessor and its own copy of the operator.
294 ConstAccessor mAcc;
295 OpT mOp;
296 StatsT mStats;
297
298 template<typename TreeT>
299 32 static inline TreeT* THROW_IF_NULL(TreeT* ptr) {
300
1/8
✗ Branch 0 not taken.
✓ Branch 1 taken 16 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
32 if (ptr == nullptr) OPENVDB_THROW(ValueError, "iterator references a null tree");
301 32 return ptr;
302 }
303
304 32 MathOp(const IterT& iter, const OpT& op):
305 32 mAcc(*THROW_IF_NULL(iter.getTree())), mOp(op)
306 32 {}
307
308 // Accumulate voxel and tile values into this functor's Stats object.
309
2/2
✓ Branch 0 taken 5445280 times.
✓ Branch 1 taken 70304 times.
11031168 void operator()(const IterT& it)
310 {
311 if (it.isVoxelValue()) {
312 // Add the magnitude of the gradient at a single voxel.
313 12251880 mStats.add(mOp.result(mAcc, it.getCoord()));
314 } else {
315 // Iterate over the voxels enclosed by a tile and add the results
316 // of applying the operator at each voxel.
317 /// @todo This could be specialized to be done more efficiently for some operators.
318 /// For example, all voxels in the interior of a tile (i.e., not on the borders)
319 /// have gradient zero, so there's no need to apply the operator to every voxel.
320 CoordBBox bbox = it.getBoundingBox();
321 Coord xyz;
322 int &x = xyz.x(), &y = xyz.y(), &z = xyz.z();
323
2/2
✓ Branch 0 taken 562432 times.
✓ Branch 1 taken 70304 times.
1265472 for (x = bbox.min().x(); x <= bbox.max().x(); ++x) {
324
2/2
✓ Branch 0 taken 4499456 times.
✓ Branch 1 taken 562432 times.
10123776 for (y = bbox.min().y(); y <= bbox.max().y(); ++y) {
325
2/2
✓ Branch 0 taken 35995648 times.
✓ Branch 1 taken 4499456 times.
80990208 for (z = bbox.min().z(); z <= bbox.max().z(); ++z) {
326 80990208 mStats.add(mOp.result(mAcc, it.getCoord()));
327 }
328 }
329 }
330 }
331 11031168 }
332
333 // Accumulate another functor's Stats object into this functor's.
334
1/2
✓ Branch 6 taken 9 times.
✗ Branch 7 not taken.
48 inline void join(MathOp& other) { mStats.add(other.mStats); }
335 }; // struct MathOp
336
337 } // namespace stats_internal
338
339 /// @endcond
340
341 template<typename IterT>
342 inline math::Histogram
343 10 histogram(const IterT& iter, double vmin, double vmax, size_t numBins, bool threaded)
344 {
345 using ValueOp = stats_internal::GetVal<IterT, math::Histogram>;
346 ValueOp valOp;
347 stats_internal::HistOp<IterT, ValueOp> op(valOp, vmin, vmax, numBins);
348
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
10 tools::accumulate(iter, op, threaded);
349 10 return op.hist;
350 }
351
352 template<typename IterT>
353 inline math::Extrema
354 extrema(const IterT& iter, bool threaded)
355 {
356 stats_internal::GetVal<IterT, math::Extrema> valOp;
357 return extrema(iter, valOp, threaded);
358 }
359
360 template<typename IterT>
361 inline math::Stats
362 statistics(const IterT& iter, bool threaded)
363 {
364 stats_internal::GetVal<IterT, math::Stats> valOp;
365
4/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
6 return statistics(iter, valOp, threaded);
366 }
367
368 template<typename IterT, typename ValueOp>
369 inline math::Extrema
370 2 extrema(const IterT& iter, const ValueOp& valOp, bool threaded)
371 {
372 stats_internal::StatsOp<IterT, const ValueOp, math::Extrema> op(valOp);
373
6/12
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 2 times.
✗ Branch 17 not taken.
17 tools::accumulate(iter, op, threaded);
374
7/13
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 2 times.
✗ Branch 16 not taken.
19 return op.stats;
375 }
376
377 template<typename IterT, typename ValueOp>
378 inline math::Stats
379 20 statistics(const IterT& iter, const ValueOp& valOp, bool threaded)
380 {
381 stats_internal::StatsOp<IterT, const ValueOp, math::Stats> op(valOp);
382
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
20 tools::accumulate(iter, op, threaded);
383
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
24 return op.stats;
384 }
385
386
387 template<typename OperatorT, typename IterT>
388 inline math::Extrema
389 2 opExtrema(const IterT& iter, const OperatorT& op, bool threaded)
390 {
391 2 stats_internal::MathOp<IterT, OperatorT, math::Extrema> func(iter, op);
392
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
2 tools::accumulate(iter, func, threaded);
393
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
4 return func.mStats;
394 }
395
396 template<typename OperatorT, typename IterT>
397 inline math::Stats
398 14 opStatistics(const IterT& iter, const OperatorT& op, bool threaded)
399 {
400 14 stats_internal::MathOp<IterT, OperatorT, math::Stats> func(iter, op);
401 14 tools::accumulate(iter, func, threaded);
402 28 return func.mStats;
403 }
404
405 } // namespace tools
406 } // namespace OPENVDB_VERSION_NAME
407 } // namespace openvdb
408
409 #endif // OPENVDB_TOOLS_STATISTICS_HAS_BEEN_INCLUDED
410