HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Filter.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 //
4 /// @author Ken Museth
5 ///
6 /// @file tools/Filter.h
7 ///
8 /// @brief Filtering of VDB volumes. All operations can optionally be masked
9 /// with another grid that acts as an alpha-mask. By default, filtering
10 /// operations do not modify the topology of the input tree and thus do
11 /// not process active tiles. However Filter::setProcessTiles can be
12 /// used to process active tiles, densifying them on demand when necessary.
13 
14 #ifndef OPENVDB_TOOLS_FILTER_HAS_BEEN_INCLUDED
15 #define OPENVDB_TOOLS_FILTER_HAS_BEEN_INCLUDED
16 
17 #include <openvdb/Types.h>
18 #include <openvdb/Grid.h>
19 #include <openvdb/math/Math.h>
20 #include <openvdb/math/Stencils.h>
21 #include <openvdb/math/Transform.h>
25 #include <openvdb/util/Util.h>
27 #include "Interpolation.h"
28 
29 #include <tbb/parallel_for.h>
30 #include <tbb/concurrent_vector.h>
31 
32 #include <algorithm> // for std::max()
33 #include <functional>
34 #include <type_traits>
35 
36 namespace openvdb {
38 namespace OPENVDB_VERSION_NAME {
39 namespace tools {
40 
41 /// @brief Volume filtering (e.g., diffusion) with optional alpha masking
42 template<typename GridT,
43  typename MaskT = typename GridT::template ValueConverter<float>::Type,
44  typename InterruptT = util::NullInterrupter>
45 class Filter
46 {
47 public:
48  using GridType = GridT;
49  using MaskType = MaskT;
50  using TreeType = typename GridType::TreeType;
51  using LeafType = typename TreeType::LeafNodeType;
52  using ValueType = typename GridType::ValueType;
53  using AlphaType = typename MaskType::ValueType;
55  using RangeType = typename LeafManagerType::LeafRange;
56  using BufferType = typename LeafManagerType::BufferType;
58  "openvdb::tools::Filter requires a mask grid with floating-point values");
59 
60  /// Constructor
61  /// @param grid Grid to be filtered.
62  /// @param interrupt Optional interrupter.
63  Filter(GridT& grid, InterruptT* interrupt = nullptr)
64  : mGrid(&grid)
65  , mTask(nullptr)
66  , mInterrupter(interrupt)
67  , mMask(nullptr)
68  , mGrainSize(1)
69  , mMinMask(0)
70  , mMaxMask(1)
71  , mInvertMask(false)
72  , mTiles(false) {}
73 
74  /// @brief Shallow copy constructor called by tbb::parallel_for()
75  /// threads during filtering.
76  /// @param other The other Filter from which to copy.
77  Filter(const Filter& other)
78  : mGrid(other.mGrid)
79  , mTask(other.mTask)
80  , mInterrupter(other.mInterrupter)
81  , mMask(other.mMask)
82  , mGrainSize(other.mGrainSize)
83  , mMinMask(other.mMinMask)
84  , mMaxMask(other.mMaxMask)
85  , mInvertMask(other.mInvertMask)
86  , mTiles(other.mTiles) {}
87 
88  /// @return the grain-size used for multi-threading
89  int getGrainSize() const { return mGrainSize; }
90  /// @brief Set the grain-size used for multi-threading.
91  /// @note A grain size of 0 or less disables multi-threading!
92  void setGrainSize(int grainsize) { mGrainSize = grainsize; }
93 
94  /// @return whether active tiles are being processed
95  bool getProcessTiles() const { return mTiles; }
96  /// @brief Set whether active tiles should also be processed.
97  /// @note If true, some tiles may become voxelized
98  /// @warning If using with a mask, ensure that the mask topology matches the
99  /// tile topology of the filter grid as tiles will not respect overlapping
100  /// mask values at tree levels finer than themselves e.g. a leaf level tile
101  /// will only use the corresponding tile ijk value in the mask grid
102  void setProcessTiles(bool flag) { mTiles = flag; }
103 
104  /// @brief Return the minimum value of the mask to be used for the
105  /// derivation of a smooth alpha value.
106  AlphaType minMask() const { return mMinMask; }
107  /// @brief Return the maximum value of the mask to be used for the
108  /// derivation of a smooth alpha value.
109  AlphaType maxMask() const { return mMaxMask; }
110  /// @brief Define the range for the (optional) scalar mask.
111  /// @param min Minimum value of the range.
112  /// @param max Maximum value of the range.
113  /// @details Mask values outside the range are clamped to zero or one, and
114  /// values inside the range map smoothly to 0->1 (unless the mask is inverted).
115  /// @throw ValueError if @a min is not smaller than @a max.
117  {
118  if (!(min < max)) OPENVDB_THROW(ValueError, "Invalid mask range (expects min < max)");
119  mMinMask = min;
120  mMaxMask = max;
121  }
122 
123  /// @brief Return true if the mask is inverted, i.e. min->max in the
124  /// original mask maps to 1->0 in the inverted alpha mask.
125  bool isMaskInverted() const { return mInvertMask; }
126  /// @brief Invert the optional mask, i.e. min->max in the original
127  /// mask maps to 1->0 in the inverted alpha mask.
128  void invertMask(bool invert=true) { mInvertMask = invert; }
129 
130  /// @brief One iteration of a fast separable mean-value (i.e. box) filter.
131  /// @param width The width of the mean-value filter is 2*width+1 voxels.
132  /// @param iterations Number of times the mean-value filter is applied.
133  /// @param mask Optional alpha mask.
134  void mean(int width = 1, int iterations = 1, const MaskType* mask = nullptr);
135 
136  /// @brief One iteration of a fast separable Gaussian filter.
137  ///
138  /// @note This is approximated as 4 iterations of a separable mean filter
139  /// which typically leads an approximation that's better than 95%!
140  /// @param width The width of the mean-value filter is 2*width+1 voxels.
141  /// @param iterations Number of times the mean-value filter is applied.
142  /// @param mask Optional alpha mask.
143  void gaussian(int width = 1, int iterations = 1, const MaskType* mask = nullptr);
144 
145  /// @brief One iteration of a median-value filter
146  ///
147  /// @note This filter is not separable and is hence relatively slow!
148  /// @param width The width of the mean-value filter is 2*width+1 voxels.
149  /// @param iterations Number of times the mean-value filter is applied.
150  /// @param mask Optional alpha mask.
151  void median(int width = 1, int iterations = 1, const MaskType* mask = nullptr);
152 
153  /// Offsets (i.e. adds) a constant value to all active voxels.
154  /// @param offset Offset in the same units as the grid.
155  /// @param mask Optional alpha mask.
156  void offset(ValueType offset, const MaskType* mask = nullptr);
157 
158  /// @brief Used internally by tbb::parallel_for()
159  /// @param range Range of LeafNodes over which to multi-thread.
160  ///
161  /// @warning Never call this method directly!
162  void operator()(const RangeType& range) const
163  {
164  if (mTask) mTask(const_cast<Filter*>(this), range);
165  else OPENVDB_THROW(ValueError, "task is undefined - call median(), mean(), etc.");
166  }
167 
168 private:
169  using LeafT = typename TreeType::LeafNodeType;
170  using VoxelIterT = typename LeafT::ValueOnIter;
171  using VoxelCIterT = typename LeafT::ValueOnCIter;
172  using BufferT = typename tree::LeafManager<TreeType>::BufferType;
173  using LeafIterT = typename RangeType::Iterator;
174  using AlphaMaskT = tools::AlphaMask<GridT, MaskT>;
175 
176  void cook(LeafManagerType& leafs);
177 
178  template<size_t Axis>
179  struct Avg {
180  Avg(const GridT* grid, Int32 w): acc(grid->tree()), width(w), frac(1.f/float(2*w+1)) {}
181  inline ValueType operator()(Coord xyz);
182  typename GridT::ConstAccessor acc;
183  const Int32 width;
184  const float frac;
185  };
186 
187  // Private filter methods called by tbb::parallel_for threads
188  template <typename AvgT>
189  void doBox(const RangeType& r, Int32 w);
190  void doBoxX(const RangeType& r, Int32 w) { this->doBox<Avg<0> >(r,w); }
191  void doBoxY(const RangeType& r, Int32 w) { this->doBox<Avg<1> >(r,w); }
192  void doBoxZ(const RangeType& r, Int32 w) { this->doBox<Avg<2> >(r,w); }
193  void doMedian(const RangeType&, int);
194  void doOffset(const RangeType&, ValueType);
195  /// @return true if the process was interrupted
196  bool wasInterrupted();
197 
198  GridType* mGrid;
199  typename std::function<void (Filter*, const RangeType&)> mTask;
200  InterruptT* mInterrupter;
201  const MaskType* mMask;
202  int mGrainSize;
203  AlphaType mMinMask, mMaxMask;
204  bool mInvertMask;
205  bool mTiles;
206 }; // end of Filter class
207 
208 
209 ////////////////////////////////////////
210 
211 /// @cond OPENVDB_DOCS_INTERNAL
212 
213 namespace filter_internal {
214 
215 template<typename TreeT>
216 struct Voxelizer
217 {
218  // NodeManager for processing internal/root node values
219  // @note Should not cache leaf nodes
220  using NodeManagerT = tree::NodeManager<TreeT, TreeT::RootNodeType::LEVEL-1>;
221  using MaskT = typename TreeT::template ValueConverter<ValueMask>::Type;
222 
223  Voxelizer(TreeT& tree, const bool allNeighbors, const size_t grainSize)
224  : mVoxelTopology()
225  , mManager(nullptr)
226  , mGrainSize(grainSize)
227  , mOp(tree, mVoxelTopology, allNeighbors ? 26 : 6) {}
228 
229  /// @brief Convert tiles to leaf nodes that exist at a particular
230  /// voxel distance away
231  /// @param width distance in voxels to seach for tiles from each leaf
232  /// @return Returns how many search iterations were performed, which
233  /// also represents how many leaf node neighbors may have been
234  /// created. Returns 0 if the tree is already entirely voxelized
235  int run(const int width)
236  {
237  if (!mOp.tree().hasActiveTiles()) return 0;
238  this->init();
239  int count = 0;
240  for (int i = 0; i < width; i += int(TreeT::LeafNodeType::DIM), ++count) {
241  if (i > 0) mManager->rebuild();
242  mManager->foreachBottomUp(mOp, mGrainSize > 0, mGrainSize);
243  mOp.tree().topologyUnion(mVoxelTopology);
244  }
245  return count;
246  }
247 
248 private:
249  void init()
250  {
251  if (mManager) {
252  mManager->rebuild();
253  }
254  else {
255  // @note We don't actually need the leaf topology here, just the
256  // internal node structure so that we can generate leaf nodes in parallel
257  mVoxelTopology.topologyUnion(mOp.tree());
258  mManager.reset(new NodeManagerT(mOp.tree()));
259  }
260  }
261 
262  struct CreateVoxelMask
263  {
264  using LeafT = typename TreeT::LeafNodeType;
265  using RootT = typename TreeT::RootNodeType;
266 
267  CreateVoxelMask(TreeT& tree, MaskT& mask, const size_t NN)
268  : mTree(tree), mVoxelTopology(mask), mNeighbors(NN) {}
269 
270  TreeT& tree() { return mTree; }
271 
272  // do nothing for leaf nodes. They shouldn't even be cached as
273  // part of the NodeManager used with this method.
274  void operator()(const LeafT&) const { assert(false); }
275 
276  void operator()(const RootT& node) const
277  {
278  using ChildT = typename RootT::ChildNodeType;
279  static constexpr Int32 CHILDDIM = Int32(ChildT::DIM);
280  static constexpr Int32 LEAFDIM = Int32(LeafT::DIM);
281  const Tester op(mTree, mNeighbors);
282 
283  auto step =
284  [&](const Coord& ijk,
285  const size_t axis1,
286  const size_t axis2,
287  const auto& val)
288  {
289  Coord offset(0);
290  Int32& a = offset[axis1];
291  Int32& b = offset[axis2];
292  for (a = 0; a < CHILDDIM; a+=LEAFDIM) {
293  for (b = 0; b < CHILDDIM; b+=LEAFDIM) {
294  const Coord childijk = ijk + offset;
295  if (op.test(childijk, val)) {
296  mVoxelTopology.touchLeaf(childijk);
297  }
298  }
299  }
300 
301  offset.reset(CHILDDIM-1);
302  for (a = 0; a < CHILDDIM; a+=LEAFDIM) {
303  for (b = 0; b < CHILDDIM; b+=LEAFDIM) {
304  const Coord childijk = ijk + offset;
305  if (op.test(childijk, val)) {
306  mVoxelTopology.touchLeaf(childijk);
307  }
308  }
309  }
310  };
311 
312  for (auto iter = node.cbeginValueOn(); iter; ++iter) {
313  const Coord& ijk = iter.getCoord();
314  // @todo step only needs to search if a given direction
315  // depending on the face
316  step(ijk, 0, 1, *iter);
317  step(ijk, 0, 2, *iter);
318  step(ijk, 1, 2, *iter);
319  }
320  }
321 
322  template<typename NodeT>
323  void operator()(const NodeT& node) const
324  {
325  using ChildT = typename NodeT::ChildNodeType;
326  static constexpr Int32 CHILDDIM = Int32(ChildT::DIM);
327  static constexpr Int32 LEAFDIM = Int32(LeafT::DIM);
328 
329  static auto step =
330  [](const Tester& op,
331  const Coord& ijk,
332  const size_t axis1,
333  const size_t axis2,
334  const auto& val,
335  std::vector<Coord>& coords)
336  {
337  Coord offset(0);
338  Int32& a = offset[axis1];
339  Int32& b = offset[axis2];
340  for (a = 0; a < CHILDDIM; a+=LEAFDIM) {
341  for (b = 0; b < CHILDDIM; b+=LEAFDIM) {
342  const Coord childijk = ijk + offset;
343  if (op.test(childijk, val)) {
344  coords.emplace_back(childijk);
345  }
346  }
347  }
348 
349  offset.reset(CHILDDIM-1);
350  for (a = 0; a < CHILDDIM; a+=LEAFDIM) {
351  for (b = 0; b < CHILDDIM; b+=LEAFDIM) {
352  const Coord childijk = ijk + offset;
353  if (op.test(childijk, val)) {
354  coords.emplace_back(childijk);
355  }
356  }
357  }
358  };
359 
360  /// Two types of algorithms here
361  /// 1) For the case where this node is the direct parent of leaf nodes
362  /// 2) For all other node types
363  ///
364  /// In general, given a tile's ijk, search its faces/edges/corners for
365  /// values which differ from its own or leaf level topology. When a
366  /// difference is detected, mask topology is generated which can be used
367  /// with topologyUnion to ensure valid voxel values exist in the source
368  /// grid.
369  ///
370  /// This operator handles all internal node types. For example, for the
371  /// lowest level internal node (which contains leaf nodes as children)
372  /// each tile is at the leaf level (a single tile represents an 8x8x8
373  /// node). CHILDDIM is this case will match the valid of LEAFDIM, as we
374  /// only need to check each tiles immediate neighbors. For higher level
375  /// internal nodes (and the root node) each child tile will have a
376  /// significantly larger CHILDDIM than the grid's LEAFDIM. We
377  /// consistently probe values along the LEAFDIM stride to ensure no
378  /// changes are missed.
379 
380  if (CHILDDIM == LEAFDIM) {
381  // If the current node is the parent of leaf nodes, search each
382  // neighbor directly and use a flag buffer to test offsets in
383  // this node which need converting to leaf level topology.
384  // This is faster than the more general method which steps across
385  // faces (unecessary due to CHILDDIM == LEAFDIM) and provides
386  // a simpler way of tracking new topology
387 
388  std::vector<char> flags(NodeT::NUM_VALUES, char(0));
389  tbb::parallel_for(tbb::blocked_range<size_t>(0, NodeT::NUM_VALUES),
390  [&](const tbb::blocked_range<size_t>& range) {
391  const Tester op(mTree, mNeighbors);
392  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
393  if (node.isValueMaskOn(Index(n))) {
394  // if index is a tile, search its neighbors
395  const Coord ijk = node.offsetToGlobalCoord(Index(n));
396  flags[n] = op.test(ijk, node.getValue(ijk));
397  }
398  }
399  });
400 
401  // create leaf level topology in this internal node
402  Index idx = 0;
403  for (auto iter = flags.begin(); iter != flags.end(); ++iter, ++idx) {
404  if (*iter) mVoxelTopology.touchLeaf(node.offsetToGlobalCoord(idx));
405  }
406  }
407  else {
408  // If this is a higher level internal node, we only need to search its
409  // face/edge/vertex neighbors for values which differ or leaf level
410  // topology. When a difference is detected, store the coordinate which
411  // needs to be voxelized.
412  // @todo investigate better threaded impl
413 
414  tbb::concurrent_vector<Coord> nodes;
415  tbb::parallel_for(tbb::blocked_range<size_t>(0, NodeT::NUM_VALUES),
416  [&](const tbb::blocked_range<size_t>& range)
417  {
418  const Tester op(mTree, mNeighbors);
419  std::vector<Coord> coords;
420 
421  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
422  if (!node.isValueMaskOn(Index(n))) continue;
423 
424  const Coord ijk = node.offsetToGlobalCoord(Index(n));
425  const auto& val = node.getValue(ijk);
426  // @todo step only needs to search if a given direction
427  // depending on the face
428  step(op, ijk, 0, 1, val, coords);
429  step(op, ijk, 0, 2, val, coords);
430  step(op, ijk, 1, 2, val, coords);
431  }
432 
433  if (!coords.empty()) {
434  std::copy(coords.begin(), coords.end(),
435  nodes.grow_by(coords.size()));
436  }
437  });
438 
439  // create leaf level topology in this internal node
440  // @note nodes may contain duplicate coords
441  for (const auto& coord : nodes) {
442  mVoxelTopology.touchLeaf(coord);
443  }
444  }
445  }
446 
447  private:
448  struct Tester
449  {
450  Tester(const TreeT& tree, const size_t NN)
451  : mAcc(tree), mNeighbors(NN) {}
452 
453  inline bool test(const Coord& ijk,
454  const typename TreeT::ValueType& val) const
455  {
456  static constexpr Int32 LEAFDIM = Int32(LeafT::DIM);
457  const Coord* NN = util::COORD_OFFSETS;
458  for (size_t i = 0; i < mNeighbors; ++i, ++NN) {
459  Coord neighbor(*NN);
460  neighbor.x() *= LEAFDIM;
461  neighbor.y() *= LEAFDIM;
462  neighbor.z() *= LEAFDIM;
463  neighbor += ijk;
464  // if a leaf exists, assume its buffer is not constant
465  if (mAcc.getValue(neighbor) != val ||
466  mAcc.probeConstLeaf(neighbor)) {
467  return true;
468  }
469  }
470  return false;
471  }
472  private:
473  const tree::ValueAccessor<const TreeT> mAcc;
474  const size_t mNeighbors;
475  };
476 
477  private:
478  TreeT& mTree;
479  MaskT& mVoxelTopology;
480  const size_t mNeighbors;
481  };// CreateVoxelMask
482 
483 private:
484  MaskT mVoxelTopology;
485  std::unique_ptr<NodeManagerT> mManager;
486  const size_t mGrainSize;
487  CreateVoxelMask mOp;
488 };
489 
490 // Helper function for Filter::Avg::operator()
491 template<typename T> static inline void accum(T& sum, T addend) { sum += addend; }
492 // Overload for bool ValueType
493 inline void accum(bool& sum, bool addend) { sum = sum || addend; }
494 
495 } // namespace filter_internal
496 
497 /// @endcond
498 
499 ////////////////////////////////////////
500 
501 
502 template<typename GridT, typename MaskT, typename InterruptT>
503 template<size_t Axis>
504 inline typename GridT::ValueType
506 {
507  ValueType sum = zeroVal<ValueType>();
508  Int32 &i = xyz[Axis], j = i + width;
509  for (i -= width; i <= j; ++i) filter_internal::accum(sum, acc.getValue(xyz));
511  return sum && frac > 0.0f;
512  } else {
514  ValueType value = static_cast<ValueType>(sum * frac);
516  return value;
517  }
518 }
519 
520 
521 ////////////////////////////////////////
522 
523 
524 template<typename GridT, typename MaskT, typename InterruptT>
525 void
526 Filter<GridT, MaskT, InterruptT>::mean(int width, int iterations, const MaskType* mask)
527 {
528  if (iterations <= 0) return;
529  mMask = mask;
530  const int w = std::max(1, width);
531  const bool serial = mGrainSize == 0;
532 
533  if (mInterrupter) mInterrupter->start("Applying mean filter");
534 
535  std::unique_ptr<filter_internal::Voxelizer<TreeType>> voxelizer;
536  if (this->getProcessTiles()) {
537  // if performing multiple iterations, also search edge/vertex
538  // neighbors for difference topology.
539  const bool allNeighbors = iterations > 1;
540  // If processing tiles, create a voxelizer and run a single
541  // width based search for tiles that need to be voxelized
542  voxelizer.reset(new filter_internal::Voxelizer<TreeType>
543  (mGrid->tree(), allNeighbors, mGrainSize));
544  if (!voxelizer->run(w)) voxelizer.reset();
545  }
546 
547  LeafManagerType leafs(mGrid->tree(), 1, serial);
548 
549  int iter = 1; // num of leaf level neighbor based searches performed
550  int dist = w; // kernel distance of the current iteration
551  for (int i=0; i<iterations && !this->wasInterrupted(); ++i, dist+=w) {
552  if (i > 0 && voxelizer) {
553  // the total influence distance in voxels of this iteration
554  // minus how far we've already accounted for
555  const int remain = dist - iter * int(TreeType::LeafNodeType::DIM);
556  if (remain > 0) {
557  const int searches = voxelizer->run(remain);
558  if (searches == 0) voxelizer.reset();
559  else leafs.rebuild(serial);
560  iter += searches;
561  }
562  }
563 
564  mTask = std::bind(&Filter::doBoxX, std::placeholders::_1, std::placeholders::_2, w);
565  this->cook(leafs);
566  // note that the order of the YZ passes are flipped to maintain backwards-compatibility
567  // with an indexing typo in the original logic
568  mTask = std::bind(&Filter::doBoxZ, std::placeholders::_1, std::placeholders::_2, w);
569  this->cook(leafs);
570  mTask = std::bind(&Filter::doBoxY, std::placeholders::_1, std::placeholders::_2, w);
571  this->cook(leafs);
572  }
573 
574  if (mInterrupter) mInterrupter->end();
575 }
576 
577 
578 template<typename GridT, typename MaskT, typename InterruptT>
579 void
580 Filter<GridT, MaskT, InterruptT>::gaussian(int width, int iterations, const MaskType* mask)
581 {
582  if (iterations <= 0) return;
583  mMask = mask;
584  const int w = std::max(1, width);
585  const bool serial = mGrainSize == 0;
586 
587  if (mInterrupter) mInterrupter->start("Applying Gaussian filter");
588 
589  std::unique_ptr<filter_internal::Voxelizer<TreeType>> voxelizer;
590  if (this->getProcessTiles()) {
591  // if performing multiple iterations, also search edge/vertex
592  // neighbors for difference topology.
593  const bool allNeighbors = iterations > 1;
594  // If processing tiles, create a voxelizer and run a single
595  // width based search for tiles that need to be voxelized
596  // @note account for sub iteration due to gaussian filter
597  voxelizer.reset(new filter_internal::Voxelizer<TreeType>
598  (mGrid->tree(), allNeighbors, mGrainSize));
599  if (!voxelizer->run(w*4)) voxelizer.reset();
600  }
601 
602  LeafManagerType leafs(mGrid->tree(), 1, serial);
603 
604  int iter = 1; // num of leaf level neighbor based searches performed
605  int dist = w*4; // kernel distance of the current iteration
606  for (int i=0; i<iterations; ++i, dist+=(w*4)) {
607  if (i > 0 && voxelizer) {
608  // the total influence distance in voxels of this iteration
609  // minus how far we've already accounted for
610  const int remain = dist - iter * int(TreeType::LeafNodeType::DIM);
611  if (remain > 0) {
612  const int searches = voxelizer->run(remain);
613  if (searches == 0) voxelizer.reset();
614  else leafs.rebuild(serial);
615  iter += searches;
616  }
617  }
618 
619  for (int n=0; n<4 && !this->wasInterrupted(); ++n) {
620  mTask = std::bind(&Filter::doBoxX, std::placeholders::_1, std::placeholders::_2, w);
621  this->cook(leafs);
622  // note that the order of the YZ passes are flipped to maintain backwards-compatibility
623  // with an indexing typo in the original logic
624  mTask = std::bind(&Filter::doBoxZ, std::placeholders::_1, std::placeholders::_2, w);
625  this->cook(leafs);
626  mTask = std::bind(&Filter::doBoxY, std::placeholders::_1, std::placeholders::_2, w);
627  this->cook(leafs);
628  }
629  }
630 
631  if (mInterrupter) mInterrupter->end();
632 }
633 
634 
635 template<typename GridT, typename MaskT, typename InterruptT>
636 void
637 Filter<GridT, MaskT, InterruptT>::median(int width, int iterations, const MaskType* mask)
638 {
639  if (iterations <= 0) return;
640  mMask = mask;
641  const int w = std::max(1, width);
642  const bool serial = mGrainSize == 0;
643 
644  if (mInterrupter) mInterrupter->start("Applying median filter");
645 
646  std::unique_ptr<filter_internal::Voxelizer<TreeType>> voxelizer;
647  if (this->getProcessTiles()) {
648  // If processing tiles, create a voxelizer and run a single
649  // width based search for tiles that need to be voxelized
650  voxelizer.reset(new filter_internal::Voxelizer<TreeType>
651  (mGrid->tree(), /*allNeighbors*/true, mGrainSize));
652  if (!voxelizer->run(w)) voxelizer.reset();
653  }
654 
655  LeafManagerType leafs(mGrid->tree(), 1, serial);
656 
657  mTask = std::bind(&Filter::doMedian, std::placeholders::_1, std::placeholders::_2, w);
658 
659  int iter = 1; // num of leaf level neighbor based searches performed
660  int dist = w; // kernel distance of the current iteration
661  for (int i=0; i<iterations && !this->wasInterrupted(); ++i, dist+=w) {
662  if (i > 0 && voxelizer) {
663  // the total influence distance in voxels of this iteration
664  // minus how far we've already accounted for
665  const int remain = dist - iter * int(TreeType::LeafNodeType::DIM);
666  if (remain > 0) {
667  const int searches = voxelizer->run(remain);
668  if (searches == 0) voxelizer.reset();
669  else leafs.rebuild(serial);
670  iter += searches;
671  }
672  }
673 
674  this->cook(leafs);
675  }
676 
677  if (mInterrupter) mInterrupter->end();
678 }
679 
680 
681 template<typename GridT, typename MaskT, typename InterruptT>
682 void
684 {
685  mMask = mask;
686 
687  if (mInterrupter) mInterrupter->start("Applying offset");
688 
689  if (this->getProcessTiles()) {
690  // Don't process leaf nodes with the node manager - we'll do them
691  // separately to allow for cleaner branching
692  using NodeManagerT = tree::NodeManager<TreeType, TreeType::RootNodeType::LEVEL-1>;
693  NodeManagerT manager(mGrid->tree());
694 
695  if (mask) {
696  manager.foreachBottomUp([&](auto& node) {
697  this->wasInterrupted();
698  AlphaMaskT alpha(*mGrid, *mMask, mMinMask, mMaxMask, mInvertMask);
699  typename AlphaMaskT::FloatType a, b;
700  for (auto iter = node.beginValueOn(); iter; ++iter) {
701  if (!alpha(iter.getCoord(), a, b)) continue;
703  iter.modifyValue([&](ValueType& v) { v += a*value; });
705  }
706  });
707  }
708  else {
709  manager.foreachBottomUp([&](auto& node) {
710  this->wasInterrupted();
711  for (auto iter = node.beginValueOn(); iter; ++iter) {
712  iter.modifyValue([&](ValueType& v) { v += value; });
713  }
714  });
715  }
716  }
717 
718  LeafManagerType leafs(mGrid->tree(), 0, mGrainSize==0);
719  mTask = std::bind(&Filter::doOffset, std::placeholders::_1, std::placeholders::_2, value);
720  this->cook(leafs);
721 
722  if (mInterrupter) mInterrupter->end();
723 }
724 
725 
726 ////////////////////////////////////////
727 
728 
729 /// Private method to perform the task (serial or threaded) and
730 /// subsequently swap the leaf buffers.
731 template<typename GridT, typename MaskT, typename InterruptT>
732 void
733 Filter<GridT, MaskT, InterruptT>::cook(LeafManagerType& leafs)
734 {
735  if (mGrainSize>0) {
736  tbb::parallel_for(leafs.leafRange(mGrainSize), *this);
737  } else {
738  (*this)(leafs.leafRange());
739  }
740  leafs.swapLeafBuffer(1, mGrainSize==0);
741 }
742 
743 
744 /// One dimensional convolution of a separable box filter
745 template<typename GridT, typename MaskT, typename InterruptT>
746 template <typename AvgT>
747 void
748 Filter<GridT, MaskT, InterruptT>::doBox(const RangeType& range, Int32 w)
749 {
750  this->wasInterrupted();
751  AvgT avg(mGrid, w);
752  if (mMask) {
753  typename AlphaMaskT::FloatType a, b;
754  AlphaMaskT alpha(*mGrid, *mMask, mMinMask, mMaxMask, mInvertMask);
755  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
756  BufferT& buffer = leafIter.buffer(1);
757  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
758  const Coord xyz = iter.getCoord();
759  if (alpha(xyz, a, b)) {
761  const ValueType value(b*(*iter) + a*avg(xyz));
763  buffer.setValue(iter.pos(), value);
764  }
765  }
766  }
767  } else {
768  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
769  BufferT& buffer = leafIter.buffer(1);
770  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
771  buffer.setValue(iter.pos(), avg(iter.getCoord()));
772  }
773  }
774  }
775 }
776 
777 
778 /// Performs simple but slow median-value diffusion
779 template<typename GridT, typename MaskT, typename InterruptT>
780 void
781 Filter<GridT, MaskT, InterruptT>::doMedian(const RangeType& range, int width)
782 {
783  this->wasInterrupted();
784  typename math::DenseStencil<GridType> stencil(*mGrid, width);//creates local cache!
785  if (mMask) {
786  typename AlphaMaskT::FloatType a, b;
787  AlphaMaskT alpha(*mGrid, *mMask, mMinMask, mMaxMask, mInvertMask);
788  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
789  BufferT& buffer = leafIter.buffer(1);
790  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
791  if (alpha(iter.getCoord(), a, b)) {
792  stencil.moveTo(iter);
794  ValueType value(b*(*iter) + a*stencil.median());
796  buffer.setValue(iter.pos(), value);
797  }
798  }
799  }
800  } else {
801  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
802  BufferT& buffer = leafIter.buffer(1);
803  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
804  stencil.moveTo(iter);
805  buffer.setValue(iter.pos(), stencil.median());
806  }
807  }
808  }
809 }
810 
811 
812 /// Offsets the values by a constant
813 template<typename GridT, typename MaskT, typename InterruptT>
814 void
815 Filter<GridT, MaskT, InterruptT>::doOffset(const RangeType& range, ValueType offset)
816 {
817  this->wasInterrupted();
818  if (mMask) {
819  typename AlphaMaskT::FloatType a, b;
820  AlphaMaskT alpha(*mGrid, *mMask, mMinMask, mMaxMask, mInvertMask);
821  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
822  for (VoxelIterT iter = leafIter->beginValueOn(); iter; ++iter) {
823  if (alpha(iter.getCoord(), a, b)) {
825  ValueType value(*iter + a*offset);
827  iter.setValue(value);
828  }
829  }
830  }
831  } else {
832  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
833  for (VoxelIterT iter = leafIter->beginValueOn(); iter; ++iter) {
834  iter.setValue(*iter + offset);
835  }
836  }
837  }
838 }
839 
840 
841 template<typename GridT, typename MaskT, typename InterruptT>
842 inline bool
844 {
845  if (util::wasInterrupted(mInterrupter)) {
846  thread::cancelGroupExecution();
847  return true;
848  }
849  return false;
850 }
851 
852 
853 ////////////////////////////////////////
854 
855 
856 // Explicit Template Instantiation
857 
858 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
859 
860 #ifdef OPENVDB_INSTANTIATE_FILTER
862 #endif
863 
864 OPENVDB_INSTANTIATE_CLASS Filter<FloatGrid, FloatGrid, util::NullInterrupter>;
865 OPENVDB_INSTANTIATE_CLASS Filter<DoubleGrid, FloatGrid, util::NullInterrupter>;
866 
867 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
868 
869 
870 } // namespace tools
871 } // namespace OPENVDB_VERSION_NAME
872 } // namespace openvdb
873 
874 #endif // OPENVDB_TOOLS_FILTER_HAS_BEEN_INCLUDED
#define OPENVDB_NO_TYPE_CONVERSION_WARNING_END
Definition: Platform.h:234
void parallel_for(int64_t start, int64_t end, std::function< void(int64_t index)> &&task, parallel_options opt=parallel_options(0, Split_Y, 1))
Definition: parallel.h:127
typedef int(APIENTRYP RE_PFNGLXSWAPINTERVALSGIPROC)(int)
GA_API const UT_StringHolder dist
GLbitfield flags
Definition: glcorearb.h:1596
cvex test(vector P=0;int unbound=3;export float s=0;export vector Cf=0;)
Definition: test.vfl:11
void mean(int width=1, int iterations=1, const MaskType *mask=nullptr)
One iteration of a fast separable mean-value (i.e. box) filter.
Definition: Filter.h:526
Volume filtering (e.g., diffusion) with optional alpha masking.
Definition: Filter.h:45
GLenum GLint * range
Definition: glcorearb.h:1925
typename CopyConstness< TreeType, NonConstBufferType >::Type BufferType
Definition: LeafManager.h:95
AlphaType minMask() const
Return the minimum value of the mask to be used for the derivation of a smooth alpha value...
Definition: Filter.h:106
GLboolean invert
Definition: glcorearb.h:549
constexpr Coord COORD_OFFSETS[26]
coordinate offset table for neighboring voxels
Definition: Util.h:22
OIIO_UTIL_API bool copy(string_view from, string_view to, std::string &err)
const GLdouble * v
Definition: glcorearb.h:837
GLsizei const GLfloat * value
Definition: glcorearb.h:824
void invertMask(bool invert=true)
Invert the optional mask, i.e. min->max in the original mask maps to 1->0 in the inverted alpha mask...
Definition: Filter.h:128
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1222
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:239
ImageBuf OIIO_API min(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
void offset(ValueType offset, const MaskType *mask=nullptr)
Definition: Filter.h:683
SYS_FORCE_INLINE const_iterator end() const
To facilitate threading over the nodes of a tree, cache node pointers in linear arrays, one for each level of the tree.
Definition: NodeManager.h:30
GLdouble n
Definition: glcorearb.h:2008
GLfloat f
Definition: glcorearb.h:1926
typename tree::LeafManager< TreeType > LeafManagerType
Definition: Filter.h:54
void setMaskRange(AlphaType min, AlphaType max)
Define the range for the (optional) scalar mask.
Definition: Filter.h:116
GLintptr offset
Definition: glcorearb.h:665
Definition: core.h:760
#define OPENVDB_INSTANTIATE_CLASS
typename MaskType::ValueType AlphaType
Definition: Filter.h:53
Filter(GridT &grid, InterruptT *interrupt=nullptr)
Definition: Filter.h:63
void foreachBottomUp(const NodeOp &op, bool threaded=true, size_t grainSize=1)
Threaded method that applies a user-supplied functor to all the nodes in the tree.
Definition: NodeManager.h:624
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
Filter(const Filter &other)
Shallow copy constructor called by tbb::parallel_for() threads during filtering.
Definition: Filter.h:77
bool isMaskInverted() const
Return true if the mask is inverted, i.e. min->max in the original mask maps to 1->0 in the inverted ...
Definition: Filter.h:125
GLint GLuint mask
Definition: glcorearb.h:124
GLuint coords
Definition: glad.h:4091
GLfloat GLfloat GLfloat alpha
Definition: glcorearb.h:112
AlphaType maxMask() const
Return the maximum value of the mask to be used for the derivation of a smooth alpha value...
Definition: Filter.h:109
void setGrainSize(int grainsize)
Set the grain-size used for multi-threading.
Definition: Filter.h:92
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1222
typename GridType::TreeType TreeType
Definition: Filter.h:50
buffer(size_t sz) FMT_NOEXCEPT
Definition: core.h:769
void operator()(const RangeType &range) const
Used internally by tbb::parallel_for()
Definition: Filter.h:162
void gaussian(int width=1, int iterations=1, const MaskType *mask=nullptr)
One iteration of a fast separable Gaussian filter.
Definition: Filter.h:580
GLint j
Definition: glad.h:2733
typename TreeType::LeafNodeType LeafType
Definition: Filter.h:51
typename LeafManagerType::BufferType BufferType
Definition: Filter.h:56
#define OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN
Bracket code with OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN/_END, to inhibit warnings about type conve...
Definition: Platform.h:233
NodeManager produces linear arrays of all tree nodes allowing for efficient threading and bottom-up p...
GLuint GLfloat * val
Definition: glcorearb.h:1608
ImageBuf OIIO_API max(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
GA_API const UT_StringHolder N
void setProcessTiles(bool flag)
Set whether active tiles should also be processed.
Definition: Filter.h:102
Defines various finite difference stencils by means of the "curiously recurring template pattern" on ...
typename GridType::ValueType ValueType
Definition: Filter.h:52
GLint GLsizei width
Definition: glcorearb.h:103
void median(int width=1, int iterations=1, const MaskType *mask=nullptr)
One iteration of a median-value filter.
Definition: Filter.h:637
GLubyte GLubyte GLubyte GLubyte w
Definition: glcorearb.h:857
Definition: core.h:1131
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...
GLboolean r
Definition: glcorearb.h:1222
GLint GLfloat GLint stencil
Definition: glcorearb.h:1278
bool wasInterrupted(T *i, int percent=-1)
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:119
GLint GLsizei count
Definition: glcorearb.h:405
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
typename LeafManagerType::LeafRange RangeType
Definition: Filter.h:55