1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 //
4 /// @file   Morphology.h
5 ///
6 /// @authors Ken Museth, Nick Avramoussis
7 ///
8 /// @brief  Implementation of morphological dilation and erosion.
9 ///
10 /// @note   By design the morphological operations only change the
11 ///         state of voxels, not their values. If one desires to
12 ///         change the values of voxels that change state an efficient
13 ///         technique is to construct a boolean mask by performing a
14 ///         topology difference between the original and final grids.
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"
24 #include "openvdb/tree/ValueAccessor.h"
25 #include "openvdb/tree/LeafManager.h"
26 #include <openvdb/openvdb.h>
27 #include <openvdb/points/PointDataGrid.h>
28 
29 #include <tbb/task_arena.h>
30 #include <tbb/enumerable_thread_specific.h>
31 #include <tbb/parallel_for.h>
32 
33 #include <type_traits>
34 #include <vector>
35 
36 
37 namespace openvdb {
38 OPENVDB_USE_VERSION_NAMESPACE
39 namespace OPENVDB_VERSION_NAME {
40 namespace tools {
41 
42 /// @brief Voxel topology of nearest neighbors
43 /// @details
44 /// <dl>
45 /// <dt><b>NN_FACE</b>
46 /// <dd>face adjacency (6 nearest neighbors, defined as all neighbor
47 /// voxels connected along one of the primary axes)
48 ///
49 /// <dt><b>NN_FACE_EDGE</b>
50 /// <dd>face and edge adjacency (18 nearest neighbors, defined as all
51 /// neighbor voxels connected along either one or two of the primary axes)
52 ///
53 /// <dt><b>NN_FACE_EDGE_VERTEX</b>
54 /// <dd>face, edge and vertex adjacency (26 nearest neighbors, defined
55 /// as all neighbor voxels connected along either one, two or all
56 /// three of the primary axes)
57 /// </dl>
58 enum NearestNeighbors { NN_FACE = 6, NN_FACE_EDGE = 18, NN_FACE_EDGE_VERTEX = 26 };
59 
60 /// @brief Different policies when dilating trees with active tiles
61 /// @details
62 /// <dl>
63 /// <dt><b>IGNORE_TILES</b>
64 /// <dd>Active tiles are ignores. For dilation, only active voxels are
65 /// dilated. For erosion, active tiles still appear as neighboring
66 /// activity however will themselves not be eroded.
67 ///
68 /// <dt><b>EXPAND_TILES</b>
69 /// <dd>For dilation and erosion, active tiles are voxelized (expanded),
70 /// dilated or eroded and left in their voxelized state irrespective of
71 /// their final state.
72 ///
73 /// <dt><b>PRESERVE_TILES</b>
74 /// <dd>For dilation, active tiles remain unchanged but they still
75 /// contribute to the dilation as if they were active voxels. For
76 /// erosion, active tiles are only eroded should the erosion wavefront
77 /// reach them, otherwise they are left unchanged. Additionally, dense
78 /// or empty nodes with constant values are pruned.
79 /// </dl>
80 enum TilePolicy { IGNORE_TILES, EXPAND_TILES, PRESERVE_TILES };
81 
82 /// @brief Topologically dilate all active values (i.e. both voxels
83 ///   and tiles) in a tree using one of three nearest neighbor
84 ///   connectivity patterns.
85 /// @details If the input is *not* a MaskTree OR if tiles are being
86 ///   preserved, this algorithm will copy the input tree topology onto a
87 ///   MaskTree, performs the dilation on the mask and copies the resulting
88 ///   topology back. This algorithm guarantees topology preservation
89 ///   (non-pruned leaf nodes will persists) EXCEPT for direct MaskTree
90 ///   dilation. MaskTree dilation is optimised for performance and may
91 ///   replace existing leaf nodes i.e. any held leaf node pointers may
92 ///   become invalid. See the Morphology class for more granular control.
93 /// @note This method is fully multi-threaded and support active tiles,
94 ///   however only the PRESERVE_TILES policy ensures a pruned topology.
95 ///   The values of any voxels are unchanged.
96 ///
97 /// @param tree          tree or leaf manager to be dilated. The leaf
98 ///                      manager will be synchronized with the result.
99 /// @param iterations    number of iterations to apply the dilation
100 /// @param nn            connectivity pattern of the dilation: either
101 ///     face-adjacent (6 nearest neighbors), face- and edge-adjacent
102 ///     (18 nearest neighbors) or face-, edge- and vertex-adjacent (26
103 ///     nearest neighbors).
104 /// @param mode          Defined the policy for handling active tiles
105 ///                      (see above for details)
106 /// @param threaded      Whether to multi-thread execution
107 template<typename TreeOrLeafManagerT>
108 void dilateActiveValues(TreeOrLeafManagerT& tree,
109     const int iterations = 1,
110     const NearestNeighbors nn = NN_FACE,
111     const TilePolicy mode = PRESERVE_TILES,
112     const bool threaded = true);
113 
114 /// @brief Topologically erode all active values (i.e. both voxels
115 ///   and tiles) in a tree using one of three nearest neighbor
116 ///   connectivity patterns.
117 /// @details If tiles are being preserve, this algorithm will copy the input
118 ///   tree topology onto a MaskTree, performs the erosion on the mask and
119 ///   intersects the resulting topology back. This algorithm guarantees
120 ///   topology preservation (non-pruned leaf nodes will persists). See the
121 ///   Morphology class for more granular control.
122 /// @note This method is fully multi-threaded and support active tiles,
123 ///   however only the PRESERVE_TILES policy ensures a pruned topology.
124 ///   The values of any voxels are unchanged. Erosion by NN_FACE neighbors
125 ///   is usually faster than other neighbor schemes. NN_FACE_EDGE and
126 ///   NN_FACE_EDGE_VERTEX operate at comparable dilation speeds.
127 ///
128 /// @param tree          tree or leaf manager to be eroded. The leaf
129 ///                      manager will be synchronized with the result.
130 /// @param iterations    number of iterations to apply the erosion
131 /// @param nn            connectivity pattern of the erosion: either
132 ///     face-adjacent (6 nearest neighbors), face- and edge-adjacent
133 ///     (18 nearest neighbors) or face-, edge- and vertex-adjacent (26
134 ///     nearest neighbors).
135 /// @param mode          Defined the policy for handling active tiles
136 ///                      (see above for details)
137 /// @param threaded      Whether to multi-thread execution
138 template<typename TreeOrLeafManagerT>
139 void erodeActiveValues(TreeOrLeafManagerT& tree,
140     const int iterations = 1,
141     const NearestNeighbors nn = NN_FACE,
142     const TilePolicy mode = PRESERVE_TILES,
143     const bool threaded = true);
144 
145 
146 ////////////////////////////////////////
147 
148 
149 namespace morphology {
150 
151 /// @brief  Dilation/Erosion operations over a Trees leaf level voxel topology.
152 template<typename TreeType>
153 class Morphology
154 {
155 public:
156     using LeafType = typename TreeType::LeafNodeType;
157     using MaskType = typename LeafType::NodeMaskType;
158     using ValueType = typename TreeType::ValueType;
159     using MaskTreeT = typename TreeType::template ValueConverter<ValueMask>::Type;
160     using MaskLeafT = typename MaskTreeT::LeafNodeType;
161     using AccessorType = tree::ValueAccessor<TreeType>;
162 
Morphology(TreeType & tree)163     Morphology(TreeType& tree)
164         : mManagerPtr(new tree::LeafManager<TreeType>(tree))
165         , mManager(*mManagerPtr)
166         , mThreaded(true) {}
167 
Morphology(tree::LeafManager<TreeType> & tree)168     Morphology(tree::LeafManager<TreeType>& tree)
169         : mManagerPtr(nullptr)
170         , mManager(tree)
171         , mThreaded(true) {}
172 
173     /// @brief  Return whether this class is using multi-threading.
getThreaded()174     bool getThreaded() const { return mThreaded; }
175     /// @brief Set whether to use multi-threading.
176     /// @note  The grain size is not exposed
setThreaded(const bool threaded)177     inline void setThreaded(const bool threaded) { mThreaded = threaded; }
178 
179     /// @brief  Return a const reference to the leaf manager
leafManager()180     inline const tree::LeafManager<TreeType>& leafManager() const { return mManager; }
181 
182     /// @brief Topologically erode all voxels by the provided nearest neighbor
183     ///    scheme and optionally collapse constant leaf nodes
184     /// @details Inactive Tiles contribute to the erosion but active tiles are
185     ///    not modified.
186     /// @param iter Number of erosion iterations
187     /// @param nn Connectivity pattern of the erosion
188     /// @param prune Whether to collapse constant leaf nodes after the erosion
189     void erodeVoxels(const size_t iter,
190         const NearestNeighbors nn,
191         const bool prune = false);
192 
193     /// @brief Topologically dilate all voxels by the provided nearest neighbor
194     ///    scheme and optionally collapse constant leaf nodes
195     /// @details Voxel values are unchanged and only leaf nodes are used to
196     ///    propagate the dilation.
197     /// @param iter Number of dilation iterations
198     /// @param nn Connectivity pattern of the dilation
199     /// @param prune Whether to collapse constant leaf nodes after the dilation
200     /// @param preserveMaskLeafNodes When dilating mask trees, the default behaviour
201     ///    chooses to steal the mask nodes rather than copy them. Although faster,
202     ///    this means that leaf nodes may be re-allocated. Set this to true if you
203     ///    need the original topology pointers to be preserved.
204     void dilateVoxels(const size_t iter,
205         const NearestNeighbors nn,
206         const bool prune = false,
207         const bool preserveMaskLeafNodes = false);
208 
209 
210     /// @brief Copy the current node masks onto the provided vector. The vector
211     ///    is resized if necessary.
212     /// @param masks The vector of NodeMasks to copy onto
copyMasks(std::vector<MaskType> & masks)213     void copyMasks(std::vector<MaskType>& masks) const
214     {
215         if (masks.size() < mManager.leafCount()) {
216             masks.resize(mManager.leafCount());
217         }
218 
219         if (this->getThreaded()) {
220             // @note this is marginally faster than using leafRange or foreach
221             tbb::parallel_for(mManager.getRange(),
222                 [&](const tbb::blocked_range<size_t>& r){
223                 for (size_t idx = r.begin(); idx < r.end(); ++idx)
224                     masks[idx] = mManager.leaf(idx).getValueMask();
225             });
226         }
227         else {
228             for (size_t idx = 0; idx < mManager.leafCount(); ++idx) {
229                 masks[idx] = mManager.leaf(idx).getValueMask();
230             }
231         }
232     }
233 
234 public:
235     /// @brief  Node Mask dilation/erosion operations for individual leaf nodes on
236     ///   a given tree. The leaf node may optionally belong to a different tree
237     ///   than the provided accessor, which will have the effect of dilating the
238     ///   leaf node mask into a different tree, or eroding the node mask based
239     ///   on corresponding neighbors in a different tree.
240     struct NodeMaskOp
241     {
242         static const Int32 DIM = static_cast<Int32>(LeafType::DIM);
243         static const Int32 LOG2DIM = static_cast<Int32>(LeafType::LOG2DIM);
244 
245         // Select the storage size based off the dimensions of the leaf node
246         using Word = typename std::conditional<LOG2DIM == 3, uint8_t,
247             typename std::conditional<LOG2DIM == 4, uint16_t,
248             typename std::conditional<LOG2DIM == 5, uint32_t,
249             typename std::conditional<LOG2DIM == 6, uint64_t,
250                 void>::type>::type>::type>::type;
251 
252         static_assert(!std::is_same<Word, void>::value,
253             "Unsupported Node Dimension for node mask dilation/erosion");
254 
NodeMaskOpNodeMaskOp255         NodeMaskOp(AccessorType& accessor,
256             const NearestNeighbors op)
257             : mOrigin(nullptr)
258             , mNeighbors(NodeMaskOp::ksize(op), nullptr)
259             , mAccessor(&accessor)
260             , mOnTile(true)
261             , mOffTile(false)
262             , mOp(op) {}
263 
264         /// @brief Dilate a single leaf node by the current spatial scheme
265         ///        stored on the instance of this NodeMaskOp. Neighbor leaf
266         ///        nodes are also updated.
267         /// @details  Unlike erode, dilate is expected to be called in a
268         ///           single threaded context as it will update the node masks
269         ///           of neighboring leaf nodes as well as the provided leaf.
270         /// @param  leaf  The leaf to dilate. The leaf's origin and value mask
271         ///               are used to calculate the result of the dilation.
dilateNodeMaskOp272         inline void dilate(LeafType& leaf)
273         {
274             // copy the mask
275             const MaskType mask = leaf.getValueMask();
276             this->dilate(leaf, mask);
277         }
278 
279         /// @brief Dilate a single leaf node by the current spatial scheme
280         ///        stored on the instance of this NodeMaskOp. The provided
281         ///        mask is used in place of the actual leaf's node mask and
282         ///        applied to the leaf afterwards. Neighbor leaf nodes are
283         ///        also updated.
284         /// @details  Unlike erode, dilate is expected to be called in a
285         ///           single threaded context as it will update the node masks
286         ///           of neighboring leaf nodes as well as the provided leaf.
287         /// @param  leaf  The leaf to dilate. The leaf's origin is used to
288         ///               calculate the result of the dilation.
289         /// @param  mask  The node mask to use in place of the current leaf
290         ///               node mask.
dilateNodeMaskOp291         inline void dilate(LeafType& leaf, const MaskType& mask)
292         {
293             this->clear();
294             mNeighbors[0] = &(leaf.getValueMask());
295             this->setOrigin(leaf.origin());
296             switch (mOp) {
297                 case NN_FACE_EDGE        : { this->dilate18(mask); return; }
298                 case NN_FACE_EDGE_VERTEX : { this->dilate26(mask); return; }
299                 case NN_FACE             : { this->dilate6(mask);  return; }
300                 default                  : {
301                     assert(false && "Unknown op during dilation."); return;
302                 }
303             }
304         }
305 
306         /// @brief Erode a single leaf node by the current spatial scheme
307         ///        stored on the instance of this NodeMaskOp.
308         /// @details  Unlike dilate, this method updates the provided mask
309         ///           and does not apply the result to the leaf node. The
310         ///           leaf node is simply used to infer the position in the
311         ///           tree to find it's neighbors. This allows erode to be
312         ///           called from multiple threads
313         /// @param  leaf  The leaf to erode. The leaf's origin is used to
314         ///               calculate the result of the erosion.
315         /// @return The eroded mask
erodeNodeMaskOp316         inline MaskType erode(const LeafType& leaf)
317         {
318             // copy the mask
319             MaskType mask = leaf.getValueMask();
320             this->erode(leaf, mask);
321             return mask;
322         }
323 
324         /// @brief Erode a single leaf node by the current spatial scheme
325         ///        stored on the instance of this NodeMaskOp. The provided
326         ///        mask is used in place of the actual leaf's node mask and
327         ///        stores the erosion result.
328         /// @details  Unlike dilate, this method updates the provided mask
329         ///           and does not apply the result to the leaf node. The
330         ///           leaf node is simply used to infer the position in the
331         ///           tree to find it's neighbors.
332         /// @param  leaf  The leaf to erode. The leaf's origin is used to
333         ///               calculate the result of the erosion.
334         /// @param  mask  The node mask to use in place of the current leaf
335         ///               node mask.
erodeNodeMaskOp336         inline void erode(const LeafType& leaf, MaskType& mask)
337         {
338             this->clear();
339             // @note leaf mask will not be modified through gather methods
340             mNeighbors[0] = const_cast<MaskType*>(&leaf.getValueMask());
341             this->setOrigin(leaf.origin());
342             switch (mOp) {
343                 case NN_FACE_EDGE        : { this->erode18(mask); return; }
344                 case NN_FACE_EDGE_VERTEX : { this->erode26(mask); return; }
345                 case NN_FACE             : { this->erode6(mask);  return; }
346                 default                  : {
347                     assert(false && "Unknown op during erosion."); return;
348                 }
349             }
350         }
351 
352     private:
ksizeNodeMaskOp353         static size_t ksize(const NearestNeighbors op) {
354             switch (op) {
355                 case NN_FACE_EDGE        : return 19;
356                 case NN_FACE_EDGE_VERTEX : return 27;
357                 case NN_FACE             : return 7;
358                 default                  : return 7;
359             }
360         }
361 
362         void dilate6(const MaskType& mask);
363         void dilate18(const MaskType& mask);
364         void dilate26(const MaskType& mask);
365         void erode6(MaskType& mask);
366 
367         /// @note  Forward API for erosion of 18/26 trees is to use erodeActiveValues
368         ///        which falls back to an inverse dilation
369         /// @todo  It may still be worth investigating more optimal gathering
370         ///        techniques
erode18NodeMaskOp371         inline void erode18(MaskType&) { OPENVDB_THROW(NotImplementedError, "erode18 is not implemented yet!"); }
erode26NodeMaskOp372         inline void erode26(MaskType&) { OPENVDB_THROW(NotImplementedError, "erode26 is not implemented yet!"); }
373 
setOriginNodeMaskOp374         inline void setOrigin(const Coord& origin) { mOrigin = &origin; }
getOriginNodeMaskOp375         inline const Coord& getOrigin() const { return *mOrigin; }
clearNodeMaskOp376         inline void clear() { std::fill(mNeighbors.begin(), mNeighbors.end(), nullptr); }
377 
scatterNodeMaskOp378         inline void scatter(size_t n, int indx)
379         {
380             assert(n < mNeighbors.size());
381             assert(mNeighbors[n]);
382             mNeighbors[n]->template getWord<Word>(indx) |= mWord;
383 
384         }
385         template<int DX, int DY, int DZ>
scatterNodeMaskOp386         inline void scatter(size_t n, int indx)
387         {
388             assert(n < mNeighbors.size());
389             if (!mNeighbors[n]) {
390                 mNeighbors[n] = this->getNeighbor<DX,DY,DZ,true>();
391             }
392             assert(mNeighbors[n]);
393             this->scatter(n, indx - (DIM - 1)*(DY + DX*DIM));
394         }
gatherNodeMaskOp395         inline Word gather(size_t n, int indx)
396         {
397             assert(n < mNeighbors.size());
398             return mNeighbors[n]->template getWord<Word>(indx);
399         }
400         template<int DX, int DY, int DZ>
gatherNodeMaskOp401         inline Word gather(size_t n, int indx)
402         {
403             assert(n < mNeighbors.size());
404             if (!mNeighbors[n]) {
405                 mNeighbors[n] = this->getNeighbor<DX,DY,DZ,false>();
406             }
407             return this->gather(n, indx - (DIM -1)*(DY + DX*DIM));
408         }
409 
410         void scatterFacesXY(int x, int y, int i1, int n, int i2);
411         void scatterEdgesXY(int x, int y, int i1, int n, int i2);
412         Word gatherFacesXY(int x, int y, int i1, int n, int i2);
413         /// @note Currently unused
414         Word gatherEdgesXY(int x, int y, int i1, int n, int i2);
415 
416         template<int DX, int DY, int DZ, bool Create>
getNeighborNodeMaskOp417         inline MaskType* getNeighbor()
418         {
419             const Coord xyz = mOrigin->offsetBy(DX*DIM, DY*DIM, DZ*DIM);
420             auto* leaf = mAccessor->probeLeaf(xyz);
421             if (leaf) return &(leaf->getValueMask());
422             if (mAccessor->isValueOn(xyz)) return &mOnTile;
423             if (!Create)                   return &mOffTile;
424             leaf = mAccessor->touchLeaf(xyz);
425             return &(leaf->getValueMask());
426         }
427 
428     private:
429         const Coord* mOrigin;
430         std::vector<MaskType*> mNeighbors;
431         AccessorType* const mAccessor;
432         Word mWord;
433         MaskType mOnTile, mOffTile;
434         const NearestNeighbors mOp;
435     };// NodeMaskOp
436 
437 private:
438     std::unique_ptr<tree::LeafManager<TreeType>> mManagerPtr;
439     tree::LeafManager<TreeType>& mManager;
440     bool mThreaded;
441 };// Morphology
442 
443 
444 template <typename TreeT>
445 typename std::enable_if<std::is_same<TreeT, typename TreeT::template ValueConverter<ValueMask>::Type>::value,
446     typename TreeT::template ValueConverter<ValueMask>::Type*>::type
getMaskTree(TreeT & tree)447 getMaskTree(TreeT& tree) { return &tree; }
448 
449 template <typename TreeT>
450 typename std::enable_if<!std::is_same<TreeT, typename TreeT::template ValueConverter<ValueMask>::Type>::value,
451     typename TreeT::template ValueConverter<ValueMask>::Type*>::type
getMaskTree(TreeT &)452 getMaskTree(TreeT&) { return nullptr; }
453 
454 
455 template <typename TreeType>
erodeVoxels(const size_t iter,const NearestNeighbors nn,const bool prune)456 void Morphology<TreeType>::erodeVoxels(const size_t iter,
457     const NearestNeighbors nn,
458     const bool prune)
459 {
460     if (iter == 0) return;
461     const size_t leafCount = mManager.leafCount();
462     if (leafCount == 0) return;
463     auto& tree = mManager.tree();
464 
465     // If the nearest neighbor mode is not FACE, fall back to an
466     // inverse dilation scheme which executes over a mask topology
467     if (nn != NN_FACE) {
468         // This method 1) dilates the input topology, 2) reverse the node masks,
469         // 3) performs a final dilation and 4) subtracts the result from the original
470         // topology. A cache of the original leaf pointers is required which tracks
471         // the original leaf nodes in a mask topology. These will need their
472         // masks updated in the original tree. The first dilation may create new leaf
473         // nodes in two instances. The first is where no topology existed before. The
474         // second is where an active tile overlaps with dilated topology. These
475         // tiles will be expanded to a dense leaf nodes by topologyUnion. We need
476         // to make sure these tiles are properly turned off.
477 
478         MaskTreeT mask(tree, false, TopologyCopy());
479 
480         // Create a new morphology class to perform dilation over the mask
481         tree::LeafManager<MaskTreeT> manager(mask);
482         Morphology<MaskTreeT> m(manager);
483         m.setThreaded(this->getThreaded());
484 
485         // perform a single dilation using the current scheme. Necessary to
486         // create edge leaf nodes and compute the active wavefront. Note that
487         // the cached array pointers will continue to be valid
488         m.dilateVoxels(1, nn, /*prune=*/false);
489 
490         // compute the wavefront. If the leaf previously existed, compute the
491         // xor activity result which is guaranteed to be equal to but slightly
492         // faster than a subtraction
493         auto computeWavefront = [&](const size_t idx) {
494             auto& leaf = manager.leaf(idx);
495             auto& nodemask = leaf.getValueMask();
496             if (const auto* original = tree.probeConstLeaf(leaf.origin())) {
497                 nodemask ^= original->getValueMask();
498             }
499             else {
500                 // should never have a dense leaf if it didn't exist in the
501                 // original tree (it was previous possible when dilateVoxels()
502                 // called topologyUnion without the preservation of active
503                 // tiles)
504                 assert(!nodemask.isOn());
505             }
506         };
507 
508         if (this->getThreaded()) {
509             tbb::parallel_for(manager.getRange(),
510                 [&](const tbb::blocked_range<size_t>& r){
511                 for (size_t idx = r.begin(); idx < r.end(); ++idx) {
512                     computeWavefront(idx);
513                 }
514             });
515         }
516         else {
517             for (size_t idx = 0; idx < manager.leafCount(); ++idx) {
518                 computeWavefront(idx);
519             }
520         }
521 
522         // perform the inverse dilation
523         m.dilateVoxels(iter, nn, /*prune=*/false);
524 
525         // subtract the inverse dilation from the original node masks
526         auto subtractTopology = [&](const size_t idx) {
527             auto& leaf = mManager.leaf(idx);
528             const auto* maskleaf = mask.probeConstLeaf(leaf.origin());
529             assert(maskleaf);
530             leaf.getValueMask() -= maskleaf->getValueMask();
531         };
532 
533         if (this->getThreaded()) {
534             tbb::parallel_for(mManager.getRange(),
535                 [&](const tbb::blocked_range<size_t>& r){
536                 for (size_t idx = r.begin(); idx < r.end(); ++idx) {
537                     subtractTopology(idx);
538                 }
539             });
540         }
541         else {
542             for (size_t idx = 0; idx < leafCount; ++idx) {
543                 subtractTopology(idx);
544             }
545         }
546     }
547     else {
548         // NN_FACE erosion scheme
549 
550         // Save the value masks of all leaf nodes.
551         std::vector<MaskType> nodeMasks;
552         this->copyMasks(nodeMasks);
553 
554         if (this->getThreaded()) {
555             const auto range = mManager.getRange();
556             for (size_t i = 0; i < iter; ++i) {
557                 // For each leaf, in parallel, gather neighboring off values
558                 // and update the cached value mask
559                 tbb::parallel_for(range,
560                     [&](const tbb::blocked_range<size_t>& r) {
561                     AccessorType accessor(tree);
562                     NodeMaskOp cache(accessor, nn);
563                     for (size_t idx = r.begin(); idx < r.end(); ++idx) {
564                         const auto& leaf = mManager.leaf(idx);
565                         if (leaf.isEmpty()) continue;
566                         // original bit-mask of current leaf node
567                         MaskType& newMask = nodeMasks[idx];
568                         cache.erode(leaf, newMask);
569                     }
570                 });
571 
572                 // update the masks after all nodes have been eroded
573                 tbb::parallel_for(range,
574                     [&](const tbb::blocked_range<size_t>& r){
575                     for (size_t idx = r.begin(); idx < r.end(); ++idx)
576                         mManager.leaf(idx).setValueMask(nodeMasks[idx]);
577                 });
578             }
579         }
580         else {
581             AccessorType accessor(tree);
582             NodeMaskOp cache(accessor, nn);
583             for (size_t i = 0; i < iter; ++i) {
584                 // For each leaf, in parallel, gather neighboring off values
585                 // and update the cached value mask
586                 for (size_t idx = 0; idx < leafCount; ++idx) {
587                     const auto& leaf = mManager.leaf(idx);
588                     if (leaf.isEmpty()) continue;
589                     // original bit-mask of current leaf node
590                     MaskType& newMask = nodeMasks[idx];
591                     cache.erode(leaf, newMask);
592                 }
593 
594                 for (size_t idx = 0; idx < leafCount; ++idx) {
595                     mManager.leaf(idx).setValueMask(nodeMasks[idx]);
596                 }
597             }
598         }
599     }
600 
601     // if prune, replace any inactive nodes
602     if (prune) {
603         tools::prune(mManager.tree(),
604             zeroVal<typename TreeType::ValueType>(),
605             this->getThreaded());
606         mManager.rebuild(!this->getThreaded());
607     }
608 }
609 
610 template <typename TreeType>
dilateVoxels(const size_t iter,const NearestNeighbors nn,const bool prune,const bool preserveMaskLeafNodes)611 void Morphology<TreeType>::dilateVoxels(const size_t iter,
612     const NearestNeighbors nn,
613     const bool prune,
614     const bool preserveMaskLeafNodes)
615 {
616     if (iter == 0) return;
617 
618     const bool threaded = this->getThreaded();
619 
620     // Actual dilation op. main implementation is single threaded. Note that this
621     // is templated (auto-ed) as the threaded implemenation may call this with a
622     // different value type to the source morphology class
623     // @note  GCC 6.4.0 crashes trying to compile this lambda with [&] capture
624     auto dilate = [iter, nn, threaded](auto& manager, const bool collapse) {
625 
626         using LeafManagerT = typename std::decay<decltype(manager)>::type;
627         using TreeT = typename LeafManagerT::TreeType;
628         using ValueT = typename TreeT::ValueType;
629         using LeafT = typename TreeT::LeafNodeType;
630 
631         // this is only used for the impl of copyMasks
632         Morphology<TreeT> m(manager);
633         m.setThreaded(threaded);
634 
635         TreeT& tree = manager.tree();
636         tree::ValueAccessor<TreeT> accessor(tree);
637 
638         // build cache objects
639         typename Morphology<TreeT>::NodeMaskOp cache(accessor, nn);
640         std::vector<MaskType> nodeMasks;
641         std::vector<std::unique_ptr<LeafT>> nodes;
642         const ValueT& bg = tree.background();
643         const bool steal = iter > 1;
644 
645         for (size_t i = 0; i < iter; ++i) {
646             if (i > 0) manager.rebuild(!threaded);
647             // If the leaf count is zero, we can stop dilation
648             const size_t leafCount = manager.leafCount();
649             if (leafCount == 0) return;
650 
651             // Copy the masks. This only resizes if necessary. As we're stealing/replacing
652             // dense nodes, it's possible we don't need to re-allocate the cache.
653             m.copyMasks(nodeMasks);
654 
655             // For each node, dilate the mask into itself and neighboring leaf nodes.
656             // If the node was originally dense (all active), steal/replace it so
657             // subsequent iterations are faster
658             manager.foreach([&](LeafT& leaf, const size_t idx) {
659                 // original bit-mask of current leaf node
660                 const MaskType& oldMask = nodeMasks[idx];
661                 const bool dense = oldMask.isOn();
662                 cache.dilate(leaf, oldMask);
663                 if (!dense) return;
664                 // This node does not need to be visited again - replace or steal
665                 if (collapse) {
666                     // if collapse, replace this dense leaf with an active background tile
667                     accessor.addTile(1, leaf.origin(), bg, true);
668                 }
669                 else if (steal) {
670                     // otherwise, temporarily steal this node
671                     nodes.emplace_back(
672                         tree.template stealNode<LeafT>(leaf.origin(),
673                             zeroVal<ValueT>(), true));
674                 }
675             }, false);
676         }
677 
678         if (nodes.empty()) return;
679         // Add back all dense nodes
680         for (auto& node : nodes) {
681             accessor.addLeaf(node.release());
682         }
683     };
684 
685     //
686 
687     if (!threaded) {
688         // single threaded dilation. If it's a mask tree we can collapse
689         // nodes during the dilation, otherwise we must call prune afterwards
690         constexpr bool isMask = std::is_same<TreeType, MaskTreeT>::value;
691         dilate(mManager, isMask && prune);
692         if (!isMask && prune) {
693             tools::prune(mManager.tree(),
694                 zeroVal<typename TreeType::ValueType>(),
695                 threaded);
696         }
697     }
698     else {
699         // multi-threaded dilation
700 
701         // Steal or create mask nodes that represent the current leaf nodes.
702         // If the input is a mask tree, optionally re-allocate the nodes if
703         // preserveMaskLeafNodes is true. This ensures that leaf node
704         // pointers are not changed in the source tree. Stealing the mask
705         // nodes is significantly faster as it also avoids a post union.
706         std::vector<MaskLeafT*> array;
707         MaskTreeT* mask = getMaskTree(mManager.tree());
708 
709         if (!mask) {
710             MaskTreeT topology;
711             topology.topologyUnion(mManager.tree());
712             array.reserve(mManager.leafCount());
713             topology.stealNodes(array);
714         }
715         else if (preserveMaskLeafNodes) {
716             mask = nullptr; // act as if theres no mask tree
717             array.resize(mManager.leafCount());
718             tbb::parallel_for(mManager.getRange(),
719                 [&](const tbb::blocked_range<size_t>& r){
720                 for (size_t idx = r.begin(); idx < r.end(); ++idx) {
721                     array[idx] = new MaskLeafT(mManager.leaf(idx));
722                 }
723             });
724         }
725         else {
726             array.reserve(mManager.leafCount());
727             mask->stealNodes(array);
728         }
729 
730         // @note this grain size is used for optimal threading
731         const size_t numThreads = size_t(tbb::this_task_arena::max_concurrency());
732         const size_t subTreeSize = math::Max(size_t(1), array.size()/(2*numThreads));
733 
734         // perform recursive dilation to sub trees
735         tbb::enumerable_thread_specific<std::unique_ptr<MaskTreeT>> pool;
736         MaskLeafT** start = array.data();
737         tbb::parallel_for(tbb::blocked_range<MaskLeafT**>(start, start + array.size(), subTreeSize),
738             [&](const tbb::blocked_range<MaskLeafT**>& range) {
739                 std::unique_ptr<MaskTreeT> mask(new MaskTreeT);
740                 for (MaskLeafT** it = range.begin(); it != range.end(); ++it) mask->addLeaf(*it);
741                 tree::LeafManager<MaskTreeT> manager(*mask, range.begin(), range.end());
742                 dilate(manager, prune);
743                 auto& subtree = pool.local();
744                 if (!subtree) subtree = std::move(mask);
745                 else          subtree->merge(*mask, MERGE_ACTIVE_STATES);
746             });
747 
748         if (!pool.empty()) {
749             auto piter = pool.begin();
750             MaskTreeT& subtree = mask ? *mask : **piter++;
751             for (; piter != pool.end(); ++piter) subtree.merge(**piter);
752             // prune, ensures partially merged nodes that may have become
753             // dense are converted to tiles
754             if (prune) tools::prune(subtree, zeroVal<typename MaskTreeT::ValueType>(), threaded);
755             // copy final topology onto dest. If mask exists, then this
756             // has already been handled by the above subtree merges
757             if (!mask) mManager.tree().topologyUnion(subtree, /*preserve-active-tiles*/true);
758         }
759     }
760 
761     // sync
762     mManager.rebuild(!threaded);
763 }
764 
765 
766 template <typename TreeType>
767 inline void
erode6(MaskType & mask)768 Morphology<TreeType>::NodeMaskOp::erode6(MaskType& mask)
769 {
770     for (int x = 0; x < DIM; ++x) {
771         for (int y = 0, n = (x << LOG2DIM); y < DIM; ++y, ++n) {
772             // Extract the portion of the original mask that corresponds to a row in z.
773             if (Word& w = mask.template getWord<Word>(n)) {
774                 // erode in two z directions (this is first since it uses the original w)
775                 w = Word(w &
776                     (Word(w<<1 | (this->template gather<0,0,-1>(1, n)>>(DIM-1))) &
777                      Word(w>>1 | (this->template gather<0,0, 1>(2, n)<<(DIM-1)))));
778                 w = Word(w & this->gatherFacesXY(x, y, 0, n, 3));
779             }
780         }// loop over y
781     }//loop over x
782 }
783 
784 template <typename TreeType>
785 inline void
dilate6(const MaskType & mask)786 Morphology<TreeType>::NodeMaskOp::dilate6(const MaskType& mask)
787 {
788     for (int x = 0; x < DIM; ++x ) {
789         for (int y = 0, n = (x << LOG2DIM);
790                  y < DIM; ++y, ++n) {
791             // Extract the portion of the original mask that corresponds to a row in z.
792             if (const Word w = mask.template getWord<Word>(n)) {
793                 // Dilate the current leaf in the +z and -z direction
794                 this->mWord = Word(w | (w>>1) | (w<<1));
795                 this->scatter(0, n);
796                 // Dilate into neighbor leaf in the -z direction
797                 if ( (this->mWord = Word(w<<(DIM-1))) ) {
798                     this->template scatter< 0, 0,-1>(1, n);
799                 }
800                 // Dilate into neighbor leaf in the +z direction
801                 if ( (this->mWord = Word(w>>(DIM-1))) ) {
802                     this->template scatter< 0, 0, 1>(2, n);
803                 }
804                 // Dilate in the xy-face directions relative to the center leaf
805                 this->mWord = w;
806                 this->scatterFacesXY(x, y, 0, n, 3);
807             }
808         }// loop over y
809     }//loop over x
810 }
811 
812 template <typename TreeType>
813 inline void
dilate18(const MaskType & mask)814 Morphology<TreeType>::NodeMaskOp::dilate18(const MaskType& mask)
815 {
816     //origins of neighbor leaf nodes in the -z and +z directions
817     const Coord origin = this->getOrigin();
818     const Coord orig_mz = origin.offsetBy(0, 0, -DIM);
819     const Coord orig_pz = origin.offsetBy(0, 0,  DIM);
820     for (int x = 0; x < DIM; ++x ) {
821         for (int y = 0, n = (x << LOG2DIM); y < DIM; ++y, ++n) {
822             if (const Word w = mask.template getWord<Word>(n)) {
823                 {
824                     this->mWord = Word(w | (w>>1) | (w<<1));
825                     this->setOrigin(origin);
826                     this->scatter(0, n);
827                     this->scatterFacesXY(x, y, 0, n, 3);
828                     this->mWord = w;
829                     this->scatterEdgesXY(x, y, 0, n, 3);
830                 }
831                 if ( (this->mWord = Word(w<<(DIM-1))) ) {
832                     this->setOrigin(origin);
833                     this->template scatter< 0, 0,-1>(1, n);
834                     this->setOrigin(orig_mz);
835                     this->scatterFacesXY(x, y, 1, n, 11);
836                 }
837                 if ( (this->mWord = Word(w>>(DIM-1))) ) {
838                     this->setOrigin(origin);
839                     this->template scatter< 0, 0, 1>(2, n);
840                     this->setOrigin(orig_pz);
841                     this->scatterFacesXY(x, y, 2, n, 15);
842                 }
843             }
844         }// loop over y
845     }//loop over x
846 }
847 
848 
849 template <typename TreeType>
850 inline void
dilate26(const MaskType & mask)851 Morphology<TreeType>::NodeMaskOp::dilate26(const MaskType& mask)
852 {
853     //origins of neighbor leaf nodes in the -z and +z directions
854     const Coord origin = this->getOrigin();
855     const Coord orig_mz = origin.offsetBy(0, 0, -DIM);
856     const Coord orig_pz = origin.offsetBy(0, 0,  DIM);
857     for (int x = 0; x < DIM; ++x) {
858         for (int y = 0, n = (x << LOG2DIM); y < DIM; ++y, ++n) {
859             if (const Word w = mask.template getWord<Word>(n)) {
860                 {
861                     this->mWord = Word(w | (w>>1) | (w<<1));
862                     this->setOrigin(origin);
863                     this->scatter(0, n);
864                     this->scatterFacesXY(x, y, 0, n, 3);
865                     this->scatterEdgesXY(x, y, 0, n, 3);
866                 }
867                 if ( (this->mWord = Word(w<<(DIM-1))) ) {
868                     this->setOrigin(origin);
869                     this->template scatter< 0, 0,-1>(1, n);
870                     this->setOrigin(orig_mz);
871                     this->scatterFacesXY(x, y, 1, n, 11);
872                     this->scatterEdgesXY(x, y, 1, n, 11);
873                 }
874                 if ( (this->mWord = Word(w>>(DIM-1))) ) {
875                     this->setOrigin(origin);
876                     this->template scatter< 0, 0, 1>(2, n);
877                     this->setOrigin(orig_pz);
878                     this->scatterFacesXY(x, y, 2, n, 19);
879                     this->scatterEdgesXY(x, y, 2, n, 19);
880                 }
881             }
882         }// loop over y
883     }//loop over x
884 }
885 
886 template<typename TreeType>
887 inline void
scatterFacesXY(int x,int y,int i1,int n,int i2)888 Morphology<TreeType>::NodeMaskOp::scatterFacesXY(int x, int y, int i1, int n, int i2)
889 {
890     // dilate current leaf or neighbor in the -x direction
891     if (x > 0) {
892         this->scatter(i1, n-DIM);
893     } else {
894         this->template scatter<-1, 0, 0>(i2, n);
895     }
896     // dilate current leaf or neighbor in the +x direction
897     if (x < DIM-1) {
898         this->scatter(i1, n+DIM);
899     } else {
900         this->template scatter< 1, 0, 0>(i2+1, n);
901     }
902     // dilate current leaf or neighbor in the -y direction
903     if (y > 0) {
904         this->scatter(i1, n-1);
905     } else {
906         this->template scatter< 0,-1, 0>(i2+2, n);
907     }
908     // dilate current leaf or neighbor in the +y direction
909     if (y < DIM-1) {
910         this->scatter(i1, n+1);
911     } else {
912         this->template scatter< 0, 1, 0>(i2+3, n);
913     }
914 }
915 
916 
917 template<typename TreeType>
918 inline void
scatterEdgesXY(int x,int y,int i1,int n,int i2)919 Morphology<TreeType>::NodeMaskOp::scatterEdgesXY(int x, int y, int i1, int n, int i2)
920 {
921     if (x > 0) {
922         if (y > 0) {
923             this->scatter(i1, n-DIM-1);
924         } else {
925             this->template scatter< 0,-1, 0>(i2+2, n-DIM);
926         }
927         if (y < DIM-1) {
928             this->scatter(i1, n-DIM+1);
929         } else {
930             this->template scatter< 0, 1, 0>(i2+3, n-DIM);
931         }
932     } else {
933         if (y < DIM-1) {
934             this->template scatter<-1, 0, 0>(i2  , n+1);
935         } else {
936             this->template scatter<-1, 1, 0>(i2+7, n  );
937         }
938         if (y > 0) {
939             this->template scatter<-1, 0, 0>(i2  , n-1);
940         } else {
941             this->template scatter<-1,-1, 0>(i2+4, n  );
942         }
943     }
944     if (x < DIM-1) {
945         if (y > 0) {
946             this->scatter(i1, n+DIM-1);
947         } else {
948             this->template scatter< 0,-1, 0>(i2+2, n+DIM);
949         }
950         if (y < DIM-1) {
951             this->scatter(i1, n+DIM+1);
952         } else {
953             this->template scatter< 0, 1, 0>(i2+3, n+DIM);
954         }
955     } else {
956         if (y > 0) {
957             this->template scatter< 1, 0, 0>(i2+1, n-1);
958         } else {
959             this->template scatter< 1,-1, 0>(i2+6, n  );
960         }
961         if (y < DIM-1) {
962             this->template scatter< 1, 0, 0>(i2+1, n+1);
963         } else {
964             this->template scatter< 1, 1, 0>(i2+5, n  );
965         }
966     }
967 }
968 
969 
970 template<typename TreeType>
971 inline typename Morphology<TreeType>::NodeMaskOp::Word
gatherFacesXY(int x,int y,int i1,int n,int i2)972 Morphology<TreeType>::NodeMaskOp::gatherFacesXY(int x, int y, int i1, int n, int i2)
973 {
974     // erode current leaf or neighbor in negative x-direction
975     Word w = x > 0 ?
976         this->gather(i1, n - DIM) :
977         this->template gather<-1,0,0>(i2, n);
978 
979     // erode current leaf or neighbor in positive x-direction
980     w = Word(w & (x < DIM - 1 ?
981         this->gather(i1, n + DIM) :
982         this->template gather<1,0,0>(i2 + 1, n)));
983 
984     // erode current leaf or neighbor in negative y-direction
985     w = Word(w & (y > 0 ?
986         this->gather(i1, n - 1) :
987         this->template gather<0,-1,0>(i2 + 2, n)));
988 
989     // erode current leaf or neighbor in positive y-direction
990     w = Word(w & (y < DIM - 1 ?
991         this->gather(i1, n + 1) :
992         this->template gather<0,1,0>(i2+3, n)));
993 
994     return w;
995 }
996 
997 
998 template<typename TreeType>
999 inline typename Morphology<TreeType>::NodeMaskOp::Word
gatherEdgesXY(int x,int y,int i1,int n,int i2)1000 Morphology<TreeType>::NodeMaskOp::gatherEdgesXY(int x, int y, int i1, int n, int i2)
1001 {
1002     Word w = ~Word(0);
1003 
1004     if (x > 0) {
1005         w &= y > 0 ?          this->gather(i1, n-DIM-1) :
1006                               this->template gather< 0,-1, 0>(i2+2, n-DIM);
1007         w &= y < DIM-1 ? this->gather(i1, n-DIM+1) :
1008                               this->template gather< 0, 1, 0>(i2+3, n-DIM);
1009     } else {
1010         w &= y < DIM-1 ? this->template gather<-1, 0, 0>(i2  , n+1):
1011                               this->template gather<-1, 1, 0>(i2+7, n  );
1012         w &= y > 0 ?          this->template gather<-1, 0, 0>(i2  , n-1):
1013                               this->template gather<-1,-1, 0>(i2+4, n  );
1014     }
1015     if (x < DIM-1) {
1016         w &= y > 0 ?          this->gather(i1, n+DIM-1) :
1017                               this->template gather< 0,-1, 0>(i2+2, n+DIM);
1018         w &= y < DIM-1 ? this->gather(i1, n+DIM+1) :
1019                               this->template gather< 0, 1, 0>(i2+3, n+DIM);
1020     } else {
1021         w &= y > 0          ? this->template gather< 1, 0, 0>(i2+1, n-1):
1022                               this->template gather< 1,-1, 0>(i2+6, n  );
1023         w &= y < DIM-1 ? this->template gather< 1, 0, 0>(i2+1, n+1):
1024                               this->template gather< 1, 1, 0>(i2+5, n  );
1025     }
1026 
1027     return w;
1028 }
1029 
1030 } // namespace morphology
1031 
1032 
1033 /////////////////////////////////////////////////////////////////////
1034 /////////////////////////////////////////////////////////////////////
1035 
1036 /// @cond OPENVDB_DOCS_INTERNAL
1037 
1038 namespace morph_internal {
1039 template <typename T> struct Adapter {
1040     using TreeType = T;
getAdapter1041     static TreeType& get(T& tree) { return tree; }
syncAdapter1042     static void sync(T&) {} // no-op
1043 };
1044 template <typename T>
1045 struct Adapter<openvdb::tree::LeafManager<T>> {
1046     using TreeType = T;
1047     static TreeType& get(openvdb::tree::LeafManager<T>& M) { return M.tree(); }
1048     static void sync(openvdb::tree::LeafManager<T>& M) { M.rebuild(); }
1049 };
1050 }
1051 
1052 /// @endcond
1053 
1054 template<typename TreeOrLeafManagerT>
1055 void dilateActiveValues(TreeOrLeafManagerT& treeOrLeafM,
1056                    const int iterations,
1057                    const NearestNeighbors nn,
1058                    const TilePolicy mode,
1059                    const bool threaded)
1060 {
1061     using AdapterT = morph_internal::Adapter<TreeOrLeafManagerT>;
1062     using TreeT = typename AdapterT::TreeType;
1063     using MaskT = typename TreeT::template ValueConverter<ValueMask>::Type;
1064 
1065     if (iterations <= 0) return;
1066 
1067     if (mode == IGNORE_TILES) {
1068         morphology::Morphology<TreeT> morph(treeOrLeafM);
1069         morph.setThreaded(threaded);
1070         // This will also sync the leaf manager
1071         morph.dilateVoxels(static_cast<size_t>(iterations), nn, /*prune=*/false);
1072         return;
1073     }
1074 
1075     // The following branching optimises from the different tree types
1076     // and TilePolicy combinations
1077 
1078     auto& tree = AdapterT::get(treeOrLeafM);
1079 
1080     // If the input is a mask tree, don't copy the topology - voxelize
1081     // it directly and let the morphology class directly steal/prune
1082     // its nodes
1083     constexpr bool isMask = std::is_same<TreeT, MaskT>::value;
1084 
1085     if (isMask || mode == EXPAND_TILES) {
1086         tree.voxelizeActiveTiles();
1087         AdapterT::sync(treeOrLeafM);
1088         morphology::Morphology<TreeT> morph(treeOrLeafM);
1089         morph.setThreaded(threaded);
1090 
1091         if (mode == PRESERVE_TILES) {
1092             morph.dilateVoxels(static_cast<size_t>(iterations), nn, /*prune=*/true);
1093         }
1094         else {
1095             assert(mode == EXPAND_TILES);
1096             morph.dilateVoxels(static_cast<size_t>(iterations), nn, /*prune=*/false);
1097         }
1098         return;
1099     }
1100 
1101     // If the tree TreeType being dilated is not a MaskTree, always copy
1102     // the topology over onto a MaskTree, perform the required dilation
1103     // and copy the final topology back. This technique avoids unnecessary
1104     // allocation with tile expansion and correctly preserves the tree
1105     // topology.
1106     //
1107     // Note that we also always use a mask if the tile policy is PRESERVE_TILES
1108     // due to the way the underlying dilation only works on voxels.
1109     // @todo Investigate tile based dilation
1110     assert(mode == PRESERVE_TILES);
1111 
1112     MaskT topology;
1113     topology.topologyUnion(tree);
1114     topology.voxelizeActiveTiles();
1115 
1116     morphology::Morphology<MaskT> morph(topology);
1117     morph.setThreaded(threaded);
1118     morph.dilateVoxels(static_cast<size_t>(iterations), nn, /*prune=*/true);
1119 
1120     tree.topologyUnion(topology, /*preserve-tiles*/true);
1121     topology.clear();
1122 
1123     // @note  this is necessary to match the behaviour of mask tree dilation
1124     //        where  source partial leaf nodes that become dense are also
1125     //        converted into tiles, not simply newly created dense nodes
1126     tools::prune(tree, zeroVal<typename TreeT::ValueType>(), threaded);
1127     AdapterT::sync(treeOrLeafM);
1128 }
1129 
1130 
1131 template<typename TreeOrLeafManagerT>
1132 void erodeActiveValues(TreeOrLeafManagerT& treeOrLeafM,
1133                       const int iterations,
1134                       const NearestNeighbors nn,
1135                       const TilePolicy mode,
1136                       const bool threaded)
1137 {
1138     using AdapterT = morph_internal::Adapter<TreeOrLeafManagerT>;
1139     using TreeT = typename AdapterT::TreeType;
1140     using MaskT = typename TreeT::template ValueConverter<ValueMask>::Type;
1141 
1142     if (iterations <= 0) return;
1143 
1144     // If the tile policiy is PRESERVE_TILES, peform the erosion on a
1145     // voxelized mask grid followed by a topology intersection such that
1146     // the original uneroded topology is preserved.
1147     if (mode == PRESERVE_TILES) {
1148         auto& tree = AdapterT::get(treeOrLeafM);
1149         MaskT topology;
1150         topology.topologyUnion(tree);
1151         topology.voxelizeActiveTiles();
1152 
1153         {
1154             morphology::Morphology<MaskT> morph(topology);
1155             morph.setThreaded(threaded);
1156             morph.erodeVoxels(static_cast<size_t>(iterations), nn, /*prune=*/false);
1157         }
1158 
1159         // prune to ensure topologyIntersection does not expand tiles
1160         // which have not been changed
1161         tools::prune(topology, zeroVal<typename MaskT::ValueType>(), threaded);
1162         tree.topologyIntersection(topology);
1163         AdapterT::sync(treeOrLeafM);
1164         return;
1165     }
1166 
1167     if (mode == EXPAND_TILES) {
1168         // if expanding, voxelize everything first if there are active tiles
1169         // @note  check first to avoid any unnecessary rebuilds
1170         auto& tree = AdapterT::get(treeOrLeafM);
1171         if (tree.hasActiveTiles()) {
1172             tree.voxelizeActiveTiles();
1173             AdapterT::sync(treeOrLeafM);
1174         }
1175     }
1176 
1177     // ignoring tiles. They won't be eroded
1178     morphology::Morphology<TreeT> morph(treeOrLeafM);
1179     morph.setThreaded(threaded);
1180     morph.erodeVoxels(static_cast<size_t>(iterations), nn, /*prune=*/false);
1181 }
1182 
1183 
1184 /////////////////////////////////////////////////////////////////////
1185 /////////////////////////////////////////////////////////////////////
1186 
1187 
1188 /// @brief Topologically dilate all leaf-level active voxels in a tree
1189 /// using one of three nearest neighbor connectivity patterns.
1190 /// @warning This method is NOT multi-threaded and ignores active tiles!
1191 ///
1192 /// @param tree          tree to be dilated
1193 /// @param iterations    number of iterations to apply the dilation
1194 /// @param nn            connectivity pattern of the dilation: either
1195 ///     face-adjacent (6 nearest neighbors), face- and edge-adjacent
1196 ///     (18 nearest neighbors) or face-, edge- and vertex-adjacent (26
1197 ///     nearest neighbors).
1198 ///
1199 /// @note The values of any voxels are unchanged.
1200 template<typename TreeType>
1201 OPENVDB_DEPRECATED_MESSAGE("Switch to tools::dilateActiveValues. Use tools::IGNORE_TILES to maintain same (but perhaps unintended) behaviour")
1202 inline void dilateVoxels(TreeType& tree,
1203                          int iterations = 1,
1204                          NearestNeighbors nn = NN_FACE)
1205 {
1206     if (iterations <= 0) return;
1207     morphology::Morphology<TreeType> morph(tree);
1208     morph.setThreaded(false); // backwards compatible
1209     // This will also sync the leaf manager
1210     morph.dilateVoxels(static_cast<size_t>(iterations), nn, /*prune=*/false);
1211 }
1212 
1213 /// @brief Topologically dilate all leaf-level active voxels in a tree
1214 /// using one of three nearest neighbor connectivity patterns.
1215 /// @warning This method is NOT multi-threaded and ignores active tiles!
1216 ///
1217 /// @param manager       LeafManager containing the tree to be dilated.
1218 ///                      On exit it is updated to include all the leaf
1219 ///                      nodes of the dilated tree.
1220 /// @param iterations    number of iterations to apply the dilation
1221 /// @param nn           connectivity pattern of the dilation: either
1222 ///     face-adjacent (6 nearest neighbors), face- and edge-adjacent
1223 ///     (18 nearest neighbors) or face-, edge- and vertex-adjacent (26
1224 ///     nearest neighbors).
1225 ///
1226 /// @note The values of any voxels are unchanged.
1227 template<typename TreeType>
1228 OPENVDB_DEPRECATED_MESSAGE("Switch to tools::dilateActiveValues. Use tools::IGNORE_TILES to maintain same (but perhaps unintended) behaviour")
1229 inline void dilateVoxels(tree::LeafManager<TreeType>& manager,
1230                          int iterations = 1,
1231                          NearestNeighbors nn = NN_FACE)
1232 {
1233     if (iterations <= 0) return;
1234     morphology::Morphology<TreeType> morph(manager);
1235     morph.setThreaded(false);
1236     morph.dilateVoxels(static_cast<size_t>(iterations), nn, /*prune=*/false);
1237 }
1238 
1239 //@{
1240 /// @brief Topologically erode all leaf-level active voxels in the given tree.
1241 /// @details That is, shrink the set of active voxels by @a iterations voxels
1242 /// in the +x, -x, +y, -y, +z and -z directions, but don't change the values
1243 /// of any voxels, only their active states.
1244 /// @todo Currently operates only on leaf voxels; need to extend to tiles.
1245 template<typename TreeType>
1246 OPENVDB_DEPRECATED_MESSAGE("Switch to tools::erodeActiveValues. Use tools::IGNORE_TILES to maintain same (but perhaps unintended) behaviour")
1247 inline void erodeVoxels(TreeType& tree,
1248                         int iterations=1,
1249                         NearestNeighbors nn = NN_FACE)
1250 {
1251     if (iterations > 0) {
1252         morphology::Morphology<TreeType> morph(tree);
1253         morph.setThreaded(true);
1254         morph.erodeVoxels(static_cast<size_t>(iterations), nn, /*prune=*/false);
1255     }
1256 
1257     tools::pruneLevelSet(tree); // matches old behaviour
1258 }
1259 
1260 template<typename TreeType>
1261 OPENVDB_DEPRECATED_MESSAGE("Switch to tools::erodeActiveValues. Use tools::IGNORE_TILES to maintain same (but perhaps unintended) behaviour")
1262 inline void erodeVoxels(tree::LeafManager<TreeType>& manager,
1263                         int iterations = 1,
1264                         NearestNeighbors nn = NN_FACE)
1265 {
1266     if (iterations <= 0) return;
1267     morphology::Morphology<TreeType> morph(manager);
1268     morph.setThreaded(true);
1269     morph.erodeVoxels(static_cast<size_t>(iterations), nn, /*prune=*/false);
1270     tools::pruneLevelSet(manager.tree()); // matches old behaviour
1271 }
1272 //@}
1273 
1274 
1275 ////////////////////////////////////////
1276 
1277 
1278 // Explicit Template Instantiation
1279 
1280 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
1281 
1282 #ifdef OPENVDB_INSTANTIATE_MORPHOLOGY
1283 #include <openvdb/util/ExplicitInstantiation.h>
1284 #endif
1285 
1286 #define _FUNCTION(TreeT) \
1287     void dilateActiveValues(TreeT&, \
1288         const int, const NearestNeighbors, const TilePolicy, const bool)
1289 OPENVDB_ALL_TREE_INSTANTIATE(_FUNCTION)
1290 #undef _FUNCTION
1291 
1292 #define _FUNCTION(TreeT) \
1293     void dilateActiveValues(tree::LeafManager<TreeT>&, \
1294         const int, const NearestNeighbors, const TilePolicy, const bool)
1295 OPENVDB_ALL_TREE_INSTANTIATE(_FUNCTION)
1296 #undef _FUNCTION
1297 
1298 #define _FUNCTION(TreeT) \
1299     void erodeActiveValues(TreeT&, \
1300         const int, const NearestNeighbors, const TilePolicy, const bool)
1301 OPENVDB_ALL_TREE_INSTANTIATE(_FUNCTION)
1302 #undef _FUNCTION
1303 
1304 #define _FUNCTION(TreeT) \
1305     void erodeActiveValues(tree::LeafManager<TreeT>&, \
1306         const int, const NearestNeighbors, const TilePolicy, const bool)
1307 OPENVDB_ALL_TREE_INSTANTIATE(_FUNCTION)
1308 #undef _FUNCTION
1309 
1310 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
1311 
1312 
1313 } // namespace tools
1314 } // namespace OPENVDB_VERSION_NAME
1315 } // namespace openvdb
1316 
1317 #endif // OPENVDB_TOOLS_MORPHOLOGY_HAS_BEEN_INCLUDED
1318