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 static inline double get(const T& val) { return double(val); } 222 }; 223 224 template<typename T> 225 struct GetValImpl<T, /*IsVector=*/true> { 226 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 inline void operator()(const IterT& iter, StatsT& stats) const { 240 if (iter.isVoxelValue()) stats.add(ImplT::get(*iter)); 241 else stats.add(ImplT::get(*iter), iter.getVoxelCount()); 242 } 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 struct StatsOp 249 { 250 StatsOp(const ValueOp& op): getValue(op) {} 251 252 // Accumulate voxel and tile values into this functor's Stats object. 253 inline void operator()(const IterT& iter) { getValue(iter, stats); } 254 255 // Accumulate another functor's Stats object into this functor's. 256 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 struct HistOp 267 { 268 HistOp(const ValueOp& op, double vmin, double vmax, size_t bins): 269 hist(vmin, vmax, bins), getValue(op) 270 {} 271 272 // Accumulate voxel and tile values into this functor's Histogram object. 273 inline void operator()(const IterT& iter) { getValue(iter, hist); } 274 275 // Accumulate another functor's Histogram object into this functor's. 276 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 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 static inline TreeT* THROW_IF_NULL(TreeT* ptr) { 300 if (ptr == nullptr) OPENVDB_THROW(ValueError, "iterator references a null tree"); 301 return ptr; 302 } 303 304 MathOp(const IterT& iter, const OpT& op): 305 mAcc(*THROW_IF_NULL(iter.getTree())), mOp(op) 306 {} 307 308 // Accumulate voxel and tile values into this functor's Stats object. 309 void operator()(const IterT& it) 310 { 311 if (it.isVoxelValue()) { 312 // Add the magnitude of the gradient at a single voxel. 313 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 for (x = bbox.min().x(); x <= bbox.max().x(); ++x) { 324 for (y = bbox.min().y(); y <= bbox.max().y(); ++y) { 325 for (z = bbox.min().z(); z <= bbox.max().z(); ++z) { 326 mStats.add(mOp.result(mAcc, it.getCoord())); 327 } 328 } 329 } 330 } 331 } 332 333 // Accumulate another functor's Stats object into this functor's. 334 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 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 tools::accumulate(iter, op, threaded); 349 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 return statistics(iter, valOp, threaded); 366 } 367 368 template<typename IterT, typename ValueOp> 369 inline math::Extrema 370 extrema(const IterT& iter, const ValueOp& valOp, bool threaded) 371 { 372 stats_internal::StatsOp<IterT, const ValueOp, math::Extrema> op(valOp); 373 tools::accumulate(iter, op, threaded); 374 return op.stats; 375 } 376 377 template<typename IterT, typename ValueOp> 378 inline math::Stats 379 statistics(const IterT& iter, const ValueOp& valOp, bool threaded) 380 { 381 stats_internal::StatsOp<IterT, const ValueOp, math::Stats> op(valOp); 382 tools::accumulate(iter, op, threaded); 383 return op.stats; 384 } 385 386 387 template<typename OperatorT, typename IterT> 388 inline math::Extrema 389 opExtrema(const IterT& iter, const OperatorT& op, bool threaded) 390 { 391 stats_internal::MathOp<IterT, OperatorT, math::Extrema> func(iter, op); 392 tools::accumulate(iter, func, threaded); 393 return func.mStats; 394 } 395 396 template<typename OperatorT, typename IterT> 397 inline math::Stats 398 opStatistics(const IterT& iter, const OperatorT& op, bool threaded) 399 { 400 stats_internal::MathOp<IterT, OperatorT, math::Stats> func(iter, op); 401 tools::accumulate(iter, func, threaded); 402 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