OpenVDB  8.1.1
Morphology.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 //
15 
16 #ifndef OPENVDB_TOOLS_MORPHOLOGY_HAS_BEEN_INCLUDED
17 #define OPENVDB_TOOLS_MORPHOLOGY_HAS_BEEN_INCLUDED
18 
19 #include "Activate.h" // backwards compatibility
20 #include "Prune.h"
21 #include "ValueTransformer.h"
22 #include "openvdb/Types.h"
23 #include "openvdb/Grid.h"
26 
27 #include <tbb/task_scheduler_init.h>
28 #include <tbb/enumerable_thread_specific.h>
29 #include <tbb/parallel_for.h>
30 
31 #include <type_traits>
32 #include <vector>
33 
34 
35 namespace openvdb {
37 namespace OPENVDB_VERSION_NAME {
38 namespace tools {
39 
57 
79 
105 template<typename TreeOrLeafManagerT>
106 inline void dilateActiveValues(TreeOrLeafManagerT& tree,
107  const int iterations = 1,
108  const NearestNeighbors nn = NN_FACE,
109  const TilePolicy mode = PRESERVE_TILES,
110  const bool threaded = true);
111 
136 template<typename TreeOrLeafManagerT>
137 inline void erodeActiveValues(TreeOrLeafManagerT& tree,
138  const int iterations = 1,
139  const NearestNeighbors nn = NN_FACE,
140  const TilePolicy mode = PRESERVE_TILES,
141  const bool threaded = true);
142 
143 
145 
146 
147 namespace morphology {
148 
150 template<typename TreeType>
152 {
153 public:
154  using LeafType = typename TreeType::LeafNodeType;
155  using MaskType = typename LeafType::NodeMaskType;
156  using ValueType = typename TreeType::ValueType;
157  using MaskTreeT = typename TreeType::template ValueConverter<ValueMask>::Type;
158  using MaskLeafT = typename MaskTreeT::LeafNodeType;
160 
161  Morphology(TreeType& tree)
162  : mManagerPtr(new tree::LeafManager<TreeType>(tree))
163  , mManager(*mManagerPtr)
164  , mThreaded(true) {}
165 
167  : mManagerPtr(nullptr)
168  , mManager(tree)
169  , mThreaded(true) {}
170 
172  bool getThreaded() const { return mThreaded; }
175  inline void setThreaded(const bool threaded) { mThreaded = threaded; }
176 
178  inline const tree::LeafManager<TreeType>& leafManager() const { return mManager; }
179 
187  void erodeVoxels(const size_t iter,
188  const NearestNeighbors nn,
189  const bool prune = false);
190 
202  void dilateVoxels(const size_t iter,
203  const NearestNeighbors nn,
204  const bool prune = false,
205  const bool preserveMaskLeafNodes = false);
206 
207 
211  void copyMasks(std::vector<MaskType>& masks) const
212  {
213  if (masks.size() < mManager.leafCount()) {
214  masks.resize(mManager.leafCount());
215  }
216 
217  if (this->getThreaded()) {
218  // @note this is marginally faster than using leafRange or foreach
219  tbb::parallel_for(mManager.getRange(),
220  [&](const tbb::blocked_range<size_t>& r){
221  for (size_t idx = r.begin(); idx < r.end(); ++idx)
222  masks[idx] = mManager.leaf(idx).getValueMask();
223  });
224  }
225  else {
226  for (size_t idx = 0; idx < mManager.leafCount(); ++idx) {
227  masks[idx] = mManager.leaf(idx).getValueMask();
228  }
229  }
230  }
231 
232 public:
238  struct NodeMaskOp
239  {
240  static const Int32 DIM = static_cast<Int32>(LeafType::DIM);
241  static const Int32 LOG2DIM = static_cast<Int32>(LeafType::LOG2DIM);
242 
243  // Select the storage size based off the dimensions of the leaf node
244  using Word = typename std::conditional<LOG2DIM == 3, uint8_t,
245  typename std::conditional<LOG2DIM == 4, uint16_t,
246  typename std::conditional<LOG2DIM == 5, uint32_t,
247  typename std::conditional<LOG2DIM == 6, uint64_t,
248  void>::type>::type>::type>::type;
249 
250  static_assert(!std::is_same<Word, void>::value,
251  "Unsupported Node Dimension for node mask dilation/erosion");
252 
254  const NearestNeighbors op)
255  : mOrigin(nullptr)
256  , mNeighbours(NodeMaskOp::ksize(op), nullptr)
257  , mAccessor(&accessor)
258  , mOnTile(true)
259  , mOffTile(false)
260  , mOp(op) {}
261 
270  inline void dilate(LeafType& leaf)
271  {
272  // copy the mask
273  const MaskType mask = leaf.getValueMask();
274  this->dilate(leaf, mask);
275  }
276 
289  inline void dilate(LeafType& leaf, const MaskType& mask)
290  {
291  this->clear();
292  mNeighbours[0] = &(leaf.getValueMask());
293  this->setOrigin(leaf.origin());
294  switch (mOp) {
295  case NN_FACE_EDGE : { this->dilate18(mask); return; }
296  case NN_FACE_EDGE_VERTEX : { this->dilate26(mask); return; }
297  case NN_FACE : { this->dilate6(mask); return; }
298  default : {
299  assert(false && "Unknown op during dilation."); return;
300  }
301  }
302  }
303 
314  inline MaskType erode(const LeafType& leaf)
315  {
316  // copy the mask
317  MaskType mask = leaf.getValueMask();
318  this->erode(leaf, mask);
319  return mask;
320  }
321 
334  inline void erode(const LeafType& leaf, MaskType& mask)
335  {
336  this->clear();
337  // @note leaf mask will not be modified through gather methods
338  mNeighbours[0] = const_cast<MaskType*>(&leaf.getValueMask());
339  this->setOrigin(leaf.origin());
340  switch (mOp) {
341  case NN_FACE_EDGE : { this->erode18(mask); return; }
342  case NN_FACE_EDGE_VERTEX : { this->erode26(mask); return; }
343  case NN_FACE : { this->erode6(mask); return; }
344  default : {
345  assert(false && "Unknown op during erosion."); return;
346  }
347  }
348  }
349 
350  private:
351  static size_t ksize(const NearestNeighbors op) {
352  switch (op) {
353  case NN_FACE_EDGE : return 19;
354  case NN_FACE_EDGE_VERTEX : return 27;
355  case NN_FACE : return 7;
356  default : return 7;
357  }
358  }
359 
360  void dilate6(const MaskType& mask);
361  void dilate18(const MaskType& mask);
362  void dilate26(const MaskType& mask);
363  void erode6(MaskType& mask);
364 
369  inline void erode18(MaskType&) { OPENVDB_THROW(NotImplementedError, "erode18 is not implemented yet!"); }
370  inline void erode26(MaskType&) { OPENVDB_THROW(NotImplementedError, "erode26 is not implemented yet!"); }
371 
372  inline void setOrigin(const Coord& origin) { mOrigin = &origin; }
373  inline const Coord& getOrigin() const { return *mOrigin; }
374  inline void clear() { std::fill(mNeighbours.begin(), mNeighbours.end(), nullptr); }
375 
376  inline void scatter(size_t n, int indx)
377  {
378  assert(n < mNeighbours.size());
379  assert(mNeighbours[n]);
380  mNeighbours[n]->template getWord<Word>(indx) |= mWord;
381 
382  }
383  template<int DX, int DY, int DZ>
384  inline void scatter(size_t n, int indx)
385  {
386  assert(n < mNeighbours.size());
387  if (!mNeighbours[n]) {
388  mNeighbours[n] = this->getNeighbour<DX,DY,DZ,true>();
389  }
390  assert(mNeighbours[n]);
391  this->scatter(n, indx - (DIM - 1)*(DY + DX*DIM));
392  }
393  inline Word gather(size_t n, int indx)
394  {
395  assert(n < mNeighbours.size());
396  return mNeighbours[n]->template getWord<Word>(indx);
397  }
398  template<int DX, int DY, int DZ>
399  inline Word gather(size_t n, int indx)
400  {
401  assert(n < mNeighbours.size());
402  if (!mNeighbours[n]) {
403  mNeighbours[n] = this->getNeighbour<DX,DY,DZ,false>();
404  }
405  return this->gather(n, indx - (DIM -1)*(DY + DX*DIM));
406  }
407 
408  void scatterFacesXY(int x, int y, int i1, int n, int i2);
409  void scatterEdgesXY(int x, int y, int i1, int n, int i2);
410  Word gatherFacesXY(int x, int y, int i1, int n, int i2);
412  Word gatherEdgesXY(int x, int y, int i1, int n, int i2);
413 
414  template<int DX, int DY, int DZ, bool Create>
415  inline MaskType* getNeighbour()
416  {
417  const Coord xyz = mOrigin->offsetBy(DX*DIM, DY*DIM, DZ*DIM);
418  auto* leaf = mAccessor->probeLeaf(xyz);
419  if (leaf) return &(leaf->getValueMask());
420  if (mAccessor->isValueOn(xyz)) return &mOnTile;
421  if (!Create) return &mOffTile;
422  leaf = mAccessor->touchLeaf(xyz);
423  return &(leaf->getValueMask());
424  }
425 
426  private:
427  const Coord* mOrigin;
428  std::vector<MaskType*> mNeighbours;
429  AccessorType* const mAccessor;
430  Word mWord;
431  MaskType mOnTile, mOffTile;
432  const NearestNeighbors mOp;
433  };// NodeMaskOp
434 
435 private:
436  std::unique_ptr<tree::LeafManager<TreeType>> mManagerPtr;
437  tree::LeafManager<TreeType>& mManager;
438  bool mThreaded;
439 };// Morphology
440 
441 
442 template <typename TreeT>
443 typename std::enable_if<std::is_same<TreeT, typename TreeT::template ValueConverter<ValueMask>::Type>::value,
444  typename TreeT::template ValueConverter<ValueMask>::Type*>::type
445 getMaskTree(TreeT& tree) { return &tree; }
446 
447 template <typename TreeT>
448 typename std::enable_if<!std::is_same<TreeT, typename TreeT::template ValueConverter<ValueMask>::Type>::value,
449  typename TreeT::template ValueConverter<ValueMask>::Type*>::type
450 getMaskTree(TreeT&) { return nullptr; }
451 
452 
453 template <typename TreeType>
454 void Morphology<TreeType>::erodeVoxels(const size_t iter,
455  const NearestNeighbors nn,
456  const bool prune)
457 {
458  if (iter == 0) return;
459  const size_t leafCount = mManager.leafCount();
460  if (leafCount == 0) return;
461  auto& tree = mManager.tree();
462 
463  // If the nearest neighbour mode is not FACE, fall back to an
464  // inverse dilation scheme which executes over a mask topology
465  if (nn != NN_FACE) {
466  // This method 1) dilates the input topology, 2) reverse the node masks,
467  // 3) performs a final dilation and 4) subtracts the result from the original
468  // topology. A cache of the original leaf pointers is required which tracks
469  // the original leaf nodes in a mask topology. These will need their
470  // masks updated in the original tree. The first dilation may create new leaf
471  // nodes in two instances. The first is where no topology existed before. The
472  // second is where an active tile overlaps with dilated topology. These
473  // tiles will be expanded to a dense leaf nodes by topologyUnion. We need
474  // to make sure these tiles are properly turned off.
475 
476  MaskTreeT mask(tree, false, TopologyCopy());
477 
478  // Create a new morphology class to perform dilation over the mask
479  tree::LeafManager<MaskTreeT> manager(mask);
480  Morphology<MaskTreeT> m(manager);
481  m.setThreaded(this->getThreaded());
482 
483  // perform a single dilation using the current scheme. Necessary to
484  // create edge leaf nodes and compute the active wavefront. Note that
485  // the cached array pointers will continue to be valid
486  m.dilateVoxels(1, nn, /*prune=*/false);
487 
488  // compute the wavefront. If the leaf previously existed, compute the
489  // xor activity result which is guaranteed to be equal to but slightly
490  // faster than a subtraction
491  auto computeWavefront = [&](const size_t idx) {
492  auto& leaf = manager.leaf(idx);
493  auto& nodemask = leaf.getValueMask();
494  if (const auto* original = tree.probeConstLeaf(leaf.origin())) {
495  nodemask ^= original->getValueMask();
496  }
497  else {
498  // should never have a dense leaf if it didn't exist in the
499  // original tree (it was previous possible when dilateVoxels()
500  // called topologyUnion without the preservation of active
501  // tiles)
502  assert(!nodemask.isOn());
503  }
504  };
505 
506  if (this->getThreaded()) {
507  tbb::parallel_for(manager.getRange(),
508  [&](const tbb::blocked_range<size_t>& r){
509  for (size_t idx = r.begin(); idx < r.end(); ++idx) {
510  computeWavefront(idx);
511  }
512  });
513  }
514  else {
515  for (size_t idx = 0; idx < manager.leafCount(); ++idx) {
516  computeWavefront(idx);
517  }
518  }
519 
520  // perform the inverse dilation
521  m.dilateVoxels(iter, nn, /*prune=*/false);
522 
523  // subtract the inverse dilation from the original node masks
524  auto subtractTopology = [&](const size_t idx) {
525  auto& leaf = mManager.leaf(idx);
526  const auto* maskleaf = mask.probeConstLeaf(leaf.origin());
527  assert(maskleaf);
528  leaf.getValueMask() -= maskleaf->getValueMask();
529  };
530 
531  if (this->getThreaded()) {
532  tbb::parallel_for(mManager.getRange(),
533  [&](const tbb::blocked_range<size_t>& r){
534  for (size_t idx = r.begin(); idx < r.end(); ++idx) {
535  subtractTopology(idx);
536  }
537  });
538  }
539  else {
540  for (size_t idx = 0; idx < leafCount; ++idx) {
541  subtractTopology(idx);
542  }
543  }
544  }
545  else {
546  // NN_FACE erosion scheme
547 
548  // Save the value masks of all leaf nodes.
549  std::vector<MaskType> nodeMasks;
550  this->copyMasks(nodeMasks);
551 
552  if (this->getThreaded()) {
553  const auto range = mManager.getRange();
554  for (size_t i = 0; i < iter; ++i) {
555  // For each leaf, in parallel, gather neighboring off values
556  // and update the cached value mask
557  tbb::parallel_for(range,
558  [&](const tbb::blocked_range<size_t>& r) {
559  AccessorType accessor(tree);
560  NodeMaskOp cache(accessor, nn);
561  for (size_t idx = r.begin(); idx < r.end(); ++idx) {
562  const auto& leaf = mManager.leaf(idx);
563  if (leaf.isEmpty()) continue;
564  // original bit-mask of current leaf node
565  MaskType& newMask = nodeMasks[idx];
566  cache.erode(leaf, newMask);
567  }
568  });
569 
570  // update the masks after all nodes have been eroded
571  tbb::parallel_for(range,
572  [&](const tbb::blocked_range<size_t>& r){
573  for (size_t idx = r.begin(); idx < r.end(); ++idx)
574  mManager.leaf(idx).setValueMask(nodeMasks[idx]);
575  });
576  }
577  }
578  else {
579  AccessorType accessor(tree);
580  NodeMaskOp cache(accessor, nn);
581  for (size_t i = 0; i < iter; ++i) {
582  // For each leaf, in parallel, gather neighboring off values
583  // and update the cached value mask
584  for (size_t idx = 0; idx < leafCount; ++idx) {
585  const auto& leaf = mManager.leaf(idx);
586  if (leaf.isEmpty()) continue;
587  // original bit-mask of current leaf node
588  MaskType& newMask = nodeMasks[idx];
589  cache.erode(leaf, newMask);
590  }
591 
592  for (size_t idx = 0; idx < leafCount; ++idx) {
593  mManager.leaf(idx).setValueMask(nodeMasks[idx]);
594  }
595  }
596  }
597  }
598 
599  // if prune, replace any inactive nodes
600  if (prune) {
601  tools::prune(mManager.tree(),
602  zeroVal<typename TreeType::ValueType>(),
603  this->getThreaded());
604  mManager.rebuild(!this->getThreaded());
605  }
606 }
607 
608 template <typename TreeType>
609 void Morphology<TreeType>::dilateVoxels(const size_t iter,
610  const NearestNeighbors nn,
611  const bool prune,
612  const bool preserveMaskLeafNodes)
613 {
614  if (iter == 0) return;
615 
616  const bool threaded = this->getThreaded();
617 
618  // Actual dilation op. main implementation is single threaded. Note that this
619  // is templated (auto-ed) as the threaded implemenation may call this with a
620  // different value type to the source morphology class
621  // @note GCC 6.4.0 crashes trying to compile this lambda with [&] capture
622  auto dilate = [iter, nn, threaded](auto& manager, const bool collapse) {
623 
624  using LeafManagerT = typename std::decay<decltype(manager)>::type;
625  using TreeT = typename LeafManagerT::TreeType;
626  using ValueT = typename TreeT::ValueType;
627  using LeafT = typename TreeT::LeafNodeType;
628 
629  // this is only used for the impl of copyMasks
630  Morphology<TreeT> m(manager);
631  m.setThreaded(threaded);
632 
633  TreeT& tree = manager.tree();
634  tree::ValueAccessor<TreeT> accessor(tree);
635 
636  // build cache objects
637  typename Morphology<TreeT>::NodeMaskOp cache(accessor, nn);
638  std::vector<MaskType> nodeMasks;
639  std::vector<std::unique_ptr<LeafT>> nodes;
640  const ValueT& bg = tree.background();
641  const bool steal = iter > 1;
642 
643  for (size_t i = 0; i < iter; ++i) {
644  if (i > 0) manager.rebuild(!threaded);
645  // If the leaf count is zero, we can stop dilation
646  const size_t leafCount = manager.leafCount();
647  if (leafCount == 0) return;
648 
649  // Copy the masks. This only resizes if necessary. As we're stealing/replacing
650  // dense nodes, it's possible we don't need to re-allocate the cache.
651  m.copyMasks(nodeMasks);
652 
653  // For each node, dilate the mask into itself and neighbouring leaf nodes.
654  // If the node was originally dense (all active), steal/replace it so
655  // subsequent iterations are faster
656  manager.foreach([&](LeafT& leaf, const size_t idx) {
657  // original bit-mask of current leaf node
658  const MaskType& oldMask = nodeMasks[idx];
659  const bool dense = oldMask.isOn();
660  cache.dilate(leaf, oldMask);
661  if (!dense) return;
662  // This node does not need to be visited again - replace or steal
663  if (collapse) {
664  // if collapse, replace this dense leaf with an active background tile
665  accessor.addTile(1, leaf.origin(), bg, true);
666  }
667  else if (steal) {
668  // otherwise, temporarily steal this node
669  nodes.emplace_back(
670  tree.template stealNode<LeafT>(leaf.origin(),
671  zeroVal<ValueT>(), true));
672  }
673  }, false);
674  }
675 
676  if (nodes.empty()) return;
677  // Add back all dense nodes
678  for (auto& node : nodes) {
679  accessor.addLeaf(node.release());
680  }
681  };
682 
683  //
684 
685  if (!threaded) {
686  // single threaded dilation. If it's a mask tree we can collapse
687  // nodes during the dilation, otherwise we must call prune afterwards
688  constexpr bool isMask = std::is_same<TreeType, MaskTreeT>::value;
689  dilate(mManager, isMask && prune);
690  if (!isMask && prune) {
691  tools::prune(mManager.tree(),
692  zeroVal<typename TreeType::ValueType>(),
693  threaded);
694  }
695  }
696  else {
697  // multi-threaded dilation
698 
699  // Steal or create mask nodes that represent the current leaf nodes.
700  // If the input is a mask tree, optionally re-allocate the nodes if
701  // preserveMaskLeafNodes is true. This ensures that leaf node
702  // pointers are not changed in the source tree. Stealing the mask
703  // nodes is significantly faster as it also avoids a post union.
704  std::vector<MaskLeafT*> array;
705  MaskTreeT* mask = getMaskTree(mManager.tree());
706 
707  if (!mask) {
708  MaskTreeT topology;
709  topology.topologyUnion(mManager.tree());
710  array.reserve(mManager.leafCount());
711  topology.stealNodes(array);
712  }
713  else if (preserveMaskLeafNodes) {
714  mask = nullptr; // act as if theres no mask tree
715  array.resize(mManager.leafCount());
716  tbb::parallel_for(mManager.getRange(),
717  [&](const tbb::blocked_range<size_t>& r){
718  for (size_t idx = r.begin(); idx < r.end(); ++idx) {
719  array[idx] = new MaskLeafT(mManager.leaf(idx));
720  }
721  });
722  }
723  else {
724  array.reserve(mManager.leafCount());
725  mask->stealNodes(array);
726  }
727 
728  // @note this grain size is used for optimal threading
729  const size_t numThreads = size_t(tbb::task_scheduler_init::default_num_threads());
730  const size_t subTreeSize = math::Max(size_t(1), array.size()/(2*numThreads));
731 
732  // perform recursive dilation to sub trees
733  tbb::enumerable_thread_specific<std::unique_ptr<MaskTreeT>> pool;
734  MaskLeafT** start = array.data();
735  tbb::parallel_for(tbb::blocked_range<MaskLeafT**>(start, start + array.size(), subTreeSize),
736  [&](const tbb::blocked_range<MaskLeafT**>& range) {
737  std::unique_ptr<MaskTreeT> mask(new MaskTreeT);
738  for (MaskLeafT** it = range.begin(); it != range.end(); ++it) mask->addLeaf(*it);
739  tree::LeafManager<MaskTreeT> manager(*mask, range.begin(), range.end());
740  dilate(manager, prune);
741  auto& subtree = pool.local();
742  if (!subtree) subtree = std::move(mask);
743  else subtree->merge(*mask, MERGE_ACTIVE_STATES);
744  });
745 
746  if (!pool.empty()) {
747  auto piter = pool.begin();
748  MaskTreeT& subtree = mask ? *mask : **piter++;
749  for (; piter != pool.end(); ++piter) subtree.merge(**piter);
750  // prune, ensures partially merged nodes that may have become
751  // dense are converted to tiles
752  if (prune) tools::prune(subtree, zeroVal<typename MaskTreeT::ValueType>(), threaded);
753  // copy final topology onto dest. If mask exists, then this
754  // has already been handled by the above subtree merges
755  if (!mask) mManager.tree().topologyUnion(subtree, /*preserve-active-tiles*/true);
756  }
757  }
758 
759  // sync
760  mManager.rebuild(!threaded);
761 }
762 
763 
764 template <typename TreeType>
765 inline void
767 {
768  for (int x = 0; x < DIM; ++x) {
769  for (int y = 0, n = (x << LOG2DIM); y < DIM; ++y, ++n) {
770  // Extract the portion of the original mask that corresponds to a row in z.
771  if (Word& w = mask.template getWord<Word>(n)) {
772  // erode in two z directions (this is first since it uses the original w)
773  w = Word(w &
774  (Word(w<<1 | (this->template gather<0,0,-1>(1, n)>>(DIM-1))) &
775  Word(w>>1 | (this->template gather<0,0, 1>(2, n)<<(DIM-1)))));
776  w = Word(w & this->gatherFacesXY(x, y, 0, n, 3));
777  }
778  }// loop over y
779  }//loop over x
780 }
781 
782 template <typename TreeType>
783 inline void
785 {
786  for (int x = 0; x < DIM; ++x ) {
787  for (int y = 0, n = (x << LOG2DIM);
788  y < DIM; ++y, ++n) {
789  // Extract the portion of the original mask that corresponds to a row in z.
790  if (const Word w = mask.template getWord<Word>(n)) {
791  // Dilate the current leaf in the +z and -z direction
792  this->mWord = Word(w | (w>>1) | (w<<1));
793  this->scatter(0, n);
794  // Dilate into neighbor leaf in the -z direction
795  if ( (this->mWord = Word(w<<(DIM-1))) ) {
796  this->template scatter< 0, 0,-1>(1, n);
797  }
798  // Dilate into neighbor leaf in the +z direction
799  if ( (this->mWord = Word(w>>(DIM-1))) ) {
800  this->template scatter< 0, 0, 1>(2, n);
801  }
802  // Dilate in the xy-face directions relative to the center leaf
803  this->mWord = w;
804  this->scatterFacesXY(x, y, 0, n, 3);
805  }
806  }// loop over y
807  }//loop over x
808 }
809 
810 template <typename TreeType>
811 inline void
813 {
814  //origins of neighbor leaf nodes in the -z and +z directions
815  const Coord origin = this->getOrigin();
816  const Coord orig_mz = origin.offsetBy(0, 0, -DIM);
817  const Coord orig_pz = origin.offsetBy(0, 0, DIM);
818  for (int x = 0; x < DIM; ++x ) {
819  for (int y = 0, n = (x << LOG2DIM); y < DIM; ++y, ++n) {
820  if (const Word w = mask.template getWord<Word>(n)) {
821  {
822  this->mWord = Word(w | (w>>1) | (w<<1));
823  this->setOrigin(origin);
824  this->scatter(0, n);
825  this->scatterFacesXY(x, y, 0, n, 3);
826  this->mWord = w;
827  this->scatterEdgesXY(x, y, 0, n, 3);
828  }
829  if ( (this->mWord = Word(w<<(DIM-1))) ) {
830  this->setOrigin(origin);
831  this->template scatter< 0, 0,-1>(1, n);
832  this->setOrigin(orig_mz);
833  this->scatterFacesXY(x, y, 1, n, 11);
834  }
835  if ( (this->mWord = Word(w>>(DIM-1))) ) {
836  this->setOrigin(origin);
837  this->template scatter< 0, 0, 1>(2, n);
838  this->setOrigin(orig_pz);
839  this->scatterFacesXY(x, y, 2, n, 15);
840  }
841  }
842  }// loop over y
843  }//loop over x
844 }
845 
846 
847 template <typename TreeType>
848 inline void
850 {
851  //origins of neighbor leaf nodes in the -z and +z directions
852  const Coord origin = this->getOrigin();
853  const Coord orig_mz = origin.offsetBy(0, 0, -DIM);
854  const Coord orig_pz = origin.offsetBy(0, 0, DIM);
855  for (int x = 0; x < DIM; ++x) {
856  for (int y = 0, n = (x << LOG2DIM); y < DIM; ++y, ++n) {
857  if (const Word w = mask.template getWord<Word>(n)) {
858  {
859  this->mWord = Word(w | (w>>1) | (w<<1));
860  this->setOrigin(origin);
861  this->scatter(0, n);
862  this->scatterFacesXY(x, y, 0, n, 3);
863  this->scatterEdgesXY(x, y, 0, n, 3);
864  }
865  if ( (this->mWord = Word(w<<(DIM-1))) ) {
866  this->setOrigin(origin);
867  this->template scatter< 0, 0,-1>(1, n);
868  this->setOrigin(orig_mz);
869  this->scatterFacesXY(x, y, 1, n, 11);
870  this->scatterEdgesXY(x, y, 1, n, 11);
871  }
872  if ( (this->mWord = Word(w>>(DIM-1))) ) {
873  this->setOrigin(origin);
874  this->template scatter< 0, 0, 1>(2, n);
875  this->setOrigin(orig_pz);
876  this->scatterFacesXY(x, y, 2, n, 19);
877  this->scatterEdgesXY(x, y, 2, n, 19);
878  }
879  }
880  }// loop over y
881  }//loop over x
882 }
883 
884 template<typename TreeType>
885 inline void
886 Morphology<TreeType>::NodeMaskOp::scatterFacesXY(int x, int y, int i1, int n, int i2)
887 {
888  // dilate current leaf or neighbor in the -x direction
889  if (x > 0) {
890  this->scatter(i1, n-DIM);
891  } else {
892  this->template scatter<-1, 0, 0>(i2, n);
893  }
894  // dilate current leaf or neighbor in the +x direction
895  if (x < DIM-1) {
896  this->scatter(i1, n+DIM);
897  } else {
898  this->template scatter< 1, 0, 0>(i2+1, n);
899  }
900  // dilate current leaf or neighbor in the -y direction
901  if (y > 0) {
902  this->scatter(i1, n-1);
903  } else {
904  this->template scatter< 0,-1, 0>(i2+2, n);
905  }
906  // dilate current leaf or neighbor in the +y direction
907  if (y < DIM-1) {
908  this->scatter(i1, n+1);
909  } else {
910  this->template scatter< 0, 1, 0>(i2+3, n);
911  }
912 }
913 
914 
915 template<typename TreeType>
916 inline void
917 Morphology<TreeType>::NodeMaskOp::scatterEdgesXY(int x, int y, int i1, int n, int i2)
918 {
919  if (x > 0) {
920  if (y > 0) {
921  this->scatter(i1, n-DIM-1);
922  } else {
923  this->template scatter< 0,-1, 0>(i2+2, n-DIM);
924  }
925  if (y < DIM-1) {
926  this->scatter(i1, n-DIM+1);
927  } else {
928  this->template scatter< 0, 1, 0>(i2+3, n-DIM);
929  }
930  } else {
931  if (y < DIM-1) {
932  this->template scatter<-1, 0, 0>(i2 , n+1);
933  } else {
934  this->template scatter<-1, 1, 0>(i2+7, n );
935  }
936  if (y > 0) {
937  this->template scatter<-1, 0, 0>(i2 , n-1);
938  } else {
939  this->template scatter<-1,-1, 0>(i2+4, n );
940  }
941  }
942  if (x < DIM-1) {
943  if (y > 0) {
944  this->scatter(i1, n+DIM-1);
945  } else {
946  this->template scatter< 0,-1, 0>(i2+2, n+DIM);
947  }
948  if (y < DIM-1) {
949  this->scatter(i1, n+DIM+1);
950  } else {
951  this->template scatter< 0, 1, 0>(i2+3, n+DIM);
952  }
953  } else {
954  if (y > 0) {
955  this->template scatter< 1, 0, 0>(i2+1, n-1);
956  } else {
957  this->template scatter< 1,-1, 0>(i2+6, n );
958  }
959  if (y < DIM-1) {
960  this->template scatter< 1, 0, 0>(i2+1, n+1);
961  } else {
962  this->template scatter< 1, 1, 0>(i2+5, n );
963  }
964  }
965 }
966 
967 
968 template<typename TreeType>
970 Morphology<TreeType>::NodeMaskOp::gatherFacesXY(int x, int y, int i1, int n, int i2)
971 {
972  // erode current leaf or neighbor in negative x-direction
973  Word w = x > 0 ?
974  this->gather(i1, n - DIM) :
975  this->template gather<-1,0,0>(i2, n);
976 
977  // erode current leaf or neighbor in positive x-direction
978  w = Word(w & (x < DIM - 1 ?
979  this->gather(i1, n + DIM) :
980  this->template gather<1,0,0>(i2 + 1, n)));
981 
982  // erode current leaf or neighbor in negative y-direction
983  w = Word(w & (y > 0 ?
984  this->gather(i1, n - 1) :
985  this->template gather<0,-1,0>(i2 + 2, n)));
986 
987  // erode current leaf or neighbor in positive y-direction
988  w = Word(w & (y < DIM - 1 ?
989  this->gather(i1, n + 1) :
990  this->template gather<0,1,0>(i2+3, n)));
991 
992  return w;
993 }
994 
995 
996 template<typename TreeType>
998 Morphology<TreeType>::NodeMaskOp::gatherEdgesXY(int x, int y, int i1, int n, int i2)
999 {
1000  Word w = ~Word(0);
1001 
1002  if (x > 0) {
1003  w &= y > 0 ? this->gather(i1, n-DIM-1) :
1004  this->template gather< 0,-1, 0>(i2+2, n-DIM);
1005  w &= y < DIM-1 ? this->gather(i1, n-DIM+1) :
1006  this->template gather< 0, 1, 0>(i2+3, n-DIM);
1007  } else {
1008  w &= y < DIM-1 ? this->template gather<-1, 0, 0>(i2 , n+1):
1009  this->template gather<-1, 1, 0>(i2+7, n );
1010  w &= y > 0 ? this->template gather<-1, 0, 0>(i2 , n-1):
1011  this->template gather<-1,-1, 0>(i2+4, n );
1012  }
1013  if (x < DIM-1) {
1014  w &= y > 0 ? this->gather(i1, n+DIM-1) :
1015  this->template gather< 0,-1, 0>(i2+2, n+DIM);
1016  w &= y < DIM-1 ? this->gather(i1, n+DIM+1) :
1017  this->template gather< 0, 1, 0>(i2+3, n+DIM);
1018  } else {
1019  w &= y > 0 ? this->template gather< 1, 0, 0>(i2+1, n-1):
1020  this->template gather< 1,-1, 0>(i2+6, n );
1021  w &= y < DIM-1 ? this->template gather< 1, 0, 0>(i2+1, n+1):
1022  this->template gather< 1, 1, 0>(i2+5, n );
1023  }
1024 
1025  return w;
1026 }
1027 
1028 } // namespace morphology
1029 
1030 
1033 
1034 namespace morph_internal {
1035 template <typename T> struct Adapter {
1036  using TreeType = T;
1037  static TreeType& get(T& tree) { return tree; }
1038  static void sync(T&) {} // no-op
1039 };
1040 template <typename T>
1042  using TreeType = T;
1043  static TreeType& get(openvdb::tree::LeafManager<T>& M) { return M.tree(); }
1044  static void sync(openvdb::tree::LeafManager<T>& M) { M.rebuild(); }
1045 };
1046 }
1047 
1048 template<typename TreeOrLeafManagerT>
1049 inline void dilateActiveValues(TreeOrLeafManagerT& treeOrLeafM,
1050  const int iterations,
1051  const NearestNeighbors nn,
1052  const TilePolicy mode,
1053  const bool threaded)
1054 {
1056  using TreeT = typename AdapterT::TreeType;
1057  using MaskT = typename TreeT::template ValueConverter<ValueMask>::Type;
1058 
1059  if (iterations <= 0) return;
1060 
1061  if (mode == IGNORE_TILES) {
1062  morphology::Morphology<TreeT> morph(treeOrLeafM);
1063  morph.setThreaded(threaded);
1064  // This will also sync the leaf manager
1065  morph.dilateVoxels(static_cast<size_t>(iterations), nn, /*prune=*/false);
1066  return;
1067  }
1068 
1069  // The following branching optimises from the different tree types
1070  // and TilePolicy combinations
1071 
1072  auto& tree = AdapterT::get(treeOrLeafM);
1073 
1074  // If the input is a mask tree, don't copy the topology - voxelize
1075  // it directly and let the morphology class directly steal/prune
1076  // its nodes
1077  constexpr bool isMask = std::is_same<TreeT, MaskT>::value;
1078 
1079  if (isMask || mode == EXPAND_TILES) {
1080  tree.voxelizeActiveTiles();
1081  AdapterT::sync(treeOrLeafM);
1082  morphology::Morphology<TreeT> morph(treeOrLeafM);
1083  morph.setThreaded(threaded);
1084 
1085  if (mode == PRESERVE_TILES) {
1086  morph.dilateVoxels(static_cast<size_t>(iterations), nn, /*prune=*/true);
1087  }
1088  else {
1089  assert(mode == EXPAND_TILES);
1090  morph.dilateVoxels(static_cast<size_t>(iterations), nn, /*prune=*/false);
1091  }
1092  return;
1093  }
1094 
1095  // If the tree TreeType being dilated is not a MaskTree, always copy
1096  // the topology over onto a MaskTree, perform the required dilation
1097  // and copy the final topology back. This technique avoids unnecessary
1098  // allocation with tile expansion and correctly preserves the tree
1099  // topology.
1100  //
1101  // Note that we also always use a mask if the tile policy is PRESERVE_TILES
1102  // due to the way the underlying dilation only works on voxels.
1103  // @todo Investigate tile based dilation
1104  assert(mode == PRESERVE_TILES);
1105 
1106  MaskT topology;
1107  topology.topologyUnion(tree);
1108  topology.voxelizeActiveTiles();
1109 
1110  morphology::Morphology<MaskT> morph(topology);
1111  morph.setThreaded(threaded);
1112  morph.dilateVoxels(static_cast<size_t>(iterations), nn, /*prune=*/true);
1113 
1114  tree.topologyUnion(topology, /*preserve-tiles*/true);
1115  topology.clear();
1116 
1117  // @note this is necessary to match the behaviour of mask tree dilation
1118  // where source partial leaf nodes that become dense are also
1119  // converted into tiles, not simply newly created dense nodes
1120  tools::prune(tree, zeroVal<typename TreeT::ValueType>(), threaded);
1121  AdapterT::sync(treeOrLeafM);
1122 }
1123 
1124 
1125 template<typename TreeOrLeafManagerT>
1126 inline void erodeActiveValues(TreeOrLeafManagerT& treeOrLeafM,
1127  const int iterations,
1128  const NearestNeighbors nn,
1129  const TilePolicy mode,
1130  const bool threaded)
1131 {
1133  using TreeT = typename AdapterT::TreeType;
1134  using MaskT = typename TreeT::template ValueConverter<ValueMask>::Type;
1135 
1136  if (iterations <= 0) return;
1137 
1138  // If the tile policiy is PRESERVE_TILES, peform the erosion on a
1139  // voxelized mask grid followed by a topology intersection such that
1140  // the original uneroded topology is preserved.
1141  if (mode == PRESERVE_TILES) {
1142  auto& tree = AdapterT::get(treeOrLeafM);
1143  MaskT topology;
1144  topology.topologyUnion(tree);
1145  topology.voxelizeActiveTiles();
1146 
1147  {
1148  morphology::Morphology<MaskT> morph(topology);
1149  morph.setThreaded(threaded);
1150  morph.erodeVoxels(static_cast<size_t>(iterations), nn, /*prune=*/false);
1151  }
1152 
1153  // prune to ensure topologyIntersection does not expand tiles
1154  // which have not been changed
1155  tools::prune(topology, zeroVal<typename MaskT::ValueType>(), threaded);
1156  tree.topologyIntersection(topology);
1157  AdapterT::sync(treeOrLeafM);
1158  return;
1159  }
1160 
1161  if (mode == EXPAND_TILES) {
1162  // if expanding, voxelize everything first if there are active tiles
1163  // @note check first to avoid any unnecessary rebuilds
1164  auto& tree = AdapterT::get(treeOrLeafM);
1165  if (tree.hasActiveTiles()) {
1166  tree.voxelizeActiveTiles();
1167  AdapterT::sync(treeOrLeafM);
1168  }
1169  }
1170 
1171  // ignoring tiles. They won't be eroded
1172  morphology::Morphology<TreeT> morph(treeOrLeafM);
1173  morph.setThreaded(threaded);
1174  morph.erodeVoxels(static_cast<size_t>(iterations), nn, /*prune=*/false);
1175 }
1176 
1177 
1180 
1181 
1194 template<typename TreeType>
1195 OPENVDB_DEPRECATED_MESSAGE("Switch to tools::dilateActiveValues. Use tools::IGNORE_TILES to maintain same (but perhaps unintended) behaviour")
1196 inline void dilateVoxels(TreeType& tree,
1197  int iterations = 1,
1199 {
1200  if (iterations <= 0) return;
1202  morph.setThreaded(false); // backwards compatible
1203  // This will also sync the leaf manager
1204  morph.dilateVoxels(static_cast<size_t>(iterations), nn, /*prune=*/false);
1205 }
1206 
1221 template<typename TreeType>
1222 OPENVDB_DEPRECATED_MESSAGE("Switch to tools::dilateActiveValues. Use tools::IGNORE_TILES to maintain same (but perhaps unintended) behaviour")
1223 inline void dilateVoxels(tree::LeafManager<TreeType>& manager,
1224  int iterations = 1,
1226 {
1227  if (iterations <= 0) return;
1228  morphology::Morphology<TreeType> morph(manager);
1229  morph.setThreaded(false);
1230  morph.dilateVoxels(static_cast<size_t>(iterations), nn, /*prune=*/false);
1231 }
1232 
1234 template<typename TreeType>
1240 OPENVDB_DEPRECATED_MESSAGE("Switch to tools::erodeActiveValues. Use tools::IGNORE_TILES to maintain same (but perhaps unintended) behaviour")
1241 inline void erodeVoxels(TreeType& tree,
1242  int iterations=1,
1244 {
1245  if (iterations > 0) {
1247  morph.setThreaded(true);
1248  morph.erodeVoxels(static_cast<size_t>(iterations), nn, /*prune=*/false);
1249  }
1250 
1251  tools::pruneLevelSet(tree); // matches old behaviour
1252 }
1253 
1254 template<typename TreeType>
1255 OPENVDB_DEPRECATED_MESSAGE("Switch to tools::erodeActiveValues. Use tools::IGNORE_TILES to maintain same (but perhaps unintended) behaviour")
1256 inline void erodeVoxels(tree::LeafManager<TreeType>& manager,
1257  int iterations = 1,
1259 {
1260  if (iterations <= 0) return;
1261  morphology::Morphology<TreeType> morph(manager);
1262  morph.setThreaded(true);
1263  morph.erodeVoxels(static_cast<size_t>(iterations), nn, /*prune=*/false);
1264  tools::pruneLevelSet(manager.tree()); // matches old behaviour
1265 }
1267 
1268 } // namespace tools
1269 } // namespace OPENVDB_VERSION_NAME
1270 } // namespace openvdb
1271 
1272 #endif // OPENVDB_TOOLS_MORPHOLOGY_HAS_BEEN_INCLUDED
Definition: Morphology.h:56
void dilateVoxels(const size_t iter, const NearestNeighbors nn, const bool prune=false, const bool preserveMaskLeafNodes=false)
Topologically dilate all voxels by the provided nearest neighbour scheme and optionally collapse cons...
Definition: Morphology.h:609
void copyMasks(std::vector< MaskType > &masks) const
Copy the current node masks onto the provided vector. The vector is resized if necessary.
Definition: Morphology.h:211
const tree::LeafManager< TreeType > & leafManager() const
Return a const reference to the leaf manager.
Definition: Morphology.h:178
#define OPENVDB_THROW(exception, message)
Definition: openvdb/Exceptions.h:74
void addTile(Index level, const Coord &xyz, const ValueType &value, bool state)
Add a tile at the specified tree level that contains voxel (x, y, z), possibly deleting existing node...
Definition: ValueAccessor.h:337
void prune(TreeT &tree, typename TreeT::ValueType tolerance=zeroVal< typename TreeT::ValueType >(), bool threaded=true, size_t grainSize=1)
Reduce the memory footprint of a tree by replacing with tiles any nodes whose values are all the same...
Definition: Prune.h:334
T TreeType
Definition: Morphology.h:1036
size_t leafCount() const
Return the number of leaf nodes.
Definition: LeafManager.h:287
typename std::conditional< LOG2DIM==3, uint8_t, typename std::conditional< LOG2DIM==4, uint16_t, typename std::conditional< LOG2DIM==5, uint32_t, typename std::conditional< LOG2DIM==6, uint64_t, void >::type >::type >::type >::type Word
Definition: Morphology.h:248
RangeType getRange(size_t grainsize=1) const
Return a tbb::blocked_range of leaf array indices.
Definition: LeafManager.h:342
void dilate(LeafType &leaf)
Dilate a single leaf node by the current spatial scheme stored on the instance of this NodeMaskOp...
Definition: Morphology.h:270
Definition: Morphology.h:56
void setThreaded(const bool threaded)
Set whether to use multi-threading.
Definition: Morphology.h:175
Definition: openvdb/Exceptions.h:61
Definition: Morphology.h:78
Morphology(TreeType &tree)
Definition: Morphology.h:161
Implementation of topological activation/deactivation.
void pruneLevelSet(TreeT &tree, bool threaded=true, size_t grainSize=1)
Reduce the memory footprint of a tree by replacing nodes whose values are all inactive with inactive ...
Definition: Prune.h:389
typename TreeType::LeafNodeType LeafType
Definition: Morphology.h:154
Defined various multi-threaded utility functions for trees.
std::enable_if<!std::is_same< TreeT, typename TreeT::template ValueConverter< ValueMask >::Type >::value, typename TreeT::template ValueConverter< ValueMask >::Type * >::type getMaskTree(TreeT &)
Definition: Morphology.h:450
void dilate(LeafType &leaf, const MaskType &mask)
Dilate a single leaf node by the current spatial scheme stored on the instance of this NodeMaskOp...
Definition: Morphology.h:289
LeafType & leaf(size_t leafIdx) const
Return a pointer to the leaf node at index leafIdx in the array.
Definition: LeafManager.h:318
void addLeaf(LeafNodeT *leaf)
Add the specified leaf to this tree, possibly creating a child branch in the process. If the leaf node already exists, replace it.
Definition: ValueAccessor.h:329
void dilateActiveValues(TreeOrLeafManagerT &tree, const int iterations=1, const NearestNeighbors nn=NN_FACE, const TilePolicy mode=PRESERVE_TILES, const bool threaded=true)
Topologically dilate all active values (i.e. both voxels and tiles) in a tree using one of three near...
Definition: Morphology.h:1049
void erode(const LeafType &leaf, MaskType &mask)
Erode a single leaf node by the current spatial scheme stored on the instance of this NodeMaskOp...
Definition: Morphology.h:334
void erodeVoxels(const size_t iter, const NearestNeighbors nn, const bool prune=false)
Topologically erode all voxels by the provided nearest neighbour scheme and optionally collapse const...
Definition: Morphology.h:454
TilePolicy
Different policies when dilating trees with active tiles.
Definition: Morphology.h:78
typename MaskTreeT::LeafNodeType MaskLeafT
Definition: Morphology.h:158
Definition: openvdb/Exceptions.h:13
NodeMaskOp(AccessorType &accessor, const NearestNeighbors op)
Definition: Morphology.h:253
Definition: Morphology.h:78
Definition: Morphology.h:78
MaskType erode(const LeafType &leaf)
Erode a single leaf node by the current spatial scheme stored on the instance of this NodeMaskOp...
Definition: Morphology.h:314
NearestNeighbors
Voxel topology of nearest neighbors.
Definition: Morphology.h:56
bool getThreaded() const
Return whether this class is using multi-threading.
Definition: Morphology.h:172
typename LeafType::NodeMaskType MaskType
Definition: Morphology.h:155
typename TreeType::ValueType ValueType
Definition: Morphology.h:156
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: openvdb/Types.h:560
Morphology(tree::LeafManager< TreeType > &tree)
Definition: Morphology.h:166
Definition: openvdb/Types.h:385
Node Mask dilation/erosion operations for individual leaf nodes on a given tree. The leaf node may op...
Definition: Morphology.h:238
void dilateVoxels(tree::LeafManager< TreeType > &manager, int iterations=1, NearestNeighbors nn=NN_FACE)
Topologically dilate all leaf-level active voxels in a tree using one of three nearest neighbor conne...
Definition: Morphology.h:1223
void erodeVoxels(tree::LeafManager< TreeType > &manager, int iterations=1, NearestNeighbors nn=NN_FACE)
Topologically erode all leaf-level active voxels in the given tree.
Definition: Morphology.h:1256
int32_t Int32
Definition: openvdb/Types.h:52
void erodeActiveValues(TreeOrLeafManagerT &tree, const int iterations=1, const NearestNeighbors nn=NN_FACE, const TilePolicy mode=PRESERVE_TILES, const bool threaded=true)
Topologically erode all active values (i.e. both voxels and tiles) in a tree using one of three neare...
Definition: Morphology.h:1126
A LeafManager manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional au...
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:598
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:116
typename TreeType::template ValueConverter< ValueMask >::Type MaskTreeT
Definition: Morphology.h:157
Dilation/Erosion operations over a Trees leaf level voxel topology.
Definition: Morphology.h:151
static void sync(T &)
Definition: Morphology.h:1038
#define OPENVDB_DEPRECATED_MESSAGE(msg)
Definition: Platform.h:123
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:178
static void sync(openvdb::tree::LeafManager< T > &M)
Definition: Morphology.h:1044