16 #ifndef OPENVDB_TOOLS_MORPHOLOGY_HAS_BEEN_INCLUDED
17 #define OPENVDB_TOOLS_MORPHOLOGY_HAS_BEEN_INCLUDED
31 #include <tbb/task_arena.h>
32 #include <tbb/enumerable_thread_specific.h>
33 #include <tbb/parallel_for.h>
35 #include <type_traits>
109 template<
typename TreeOrLeafManagerT>
111 const int iterations = 1,
140 template<
typename TreeOrLeafManagerT>
142 const int iterations = 1,
151 namespace morphology {
154 template<
typename TreeType>
166 : mManagerPtr(new tree::LeafManager<TreeType>(tree))
167 , mManager(*mManagerPtr)
171 : mManagerPtr(nullptr)
193 const bool prune =
false);
208 const bool prune =
false,
209 const bool preserveMaskLeafNodes =
false);
217 if (masks.size() < mManager.
leafCount()) {
224 [&](
const tbb::blocked_range<size_t>&
r){
225 for (
size_t idx =
r.begin(); idx <
r.end(); ++idx)
226 masks[idx] = mManager.
leaf(idx).getValueMask();
230 for (
size_t idx = 0; idx < mManager.
leafCount(); ++idx) {
231 masks[idx] = mManager.
leaf(idx).getValueMask();
248 using Word =
typename std::conditional<
LOG2DIM == 3, uint8_t,
249 typename std::conditional<
LOG2DIM == 4, uint16_t,
250 typename std::conditional<
LOG2DIM == 5, uint32_t,
251 typename std::conditional<
LOG2DIM == 6, uint64_t,
255 "Unsupported Node Dimension for node mask dilation/erosion");
261 , mAccessor(&accessor)
296 mNeighbors[0] = &(leaf.getValueMask());
297 this->setOrigin(leaf.origin());
301 case NN_FACE : { this->dilate6(mask);
return; }
322 this->
erode(leaf, mask);
342 mNeighbors[0] =
const_cast<MaskType*
>(&leaf.getValueMask());
343 this->setOrigin(leaf.origin());
347 case NN_FACE : { this->erode6(mask);
return; }
373 inline void erode18(
MaskType&) {
OPENVDB_THROW(NotImplementedError,
"erode18 is not implemented yet!"); }
374 inline void erode26(
MaskType&) {
OPENVDB_THROW(NotImplementedError,
"erode26 is not implemented yet!"); }
376 inline void setOrigin(
const Coord& origin) { mOrigin = &origin; }
377 inline const Coord& getOrigin()
const {
return *mOrigin; }
378 inline void clear() {
std::fill(mNeighbors.begin(), mNeighbors.end(),
nullptr); }
380 inline void scatter(
size_t n,
int indx)
384 mNeighbors[
n]->template getWord<Word>(indx) |= mWord;
387 template<
int DX,
int DY,
int DZ>
388 inline void scatter(
size_t n,
int indx)
391 if (!mNeighbors[n]) {
392 mNeighbors[
n] = this->getNeighbor<DX,DY,DZ,true>();
395 this->scatter(n, indx - (
DIM - 1)*(DY + DX*
DIM));
397 inline Word gather(
size_t n,
int indx)
400 return mNeighbors[
n]->template getWord<Word>(indx);
402 template<
int DX,
int DY,
int DZ>
403 inline Word gather(
size_t n,
int indx)
406 if (!mNeighbors[n]) {
407 mNeighbors[
n] = this->getNeighbor<DX,DY,DZ,false>();
409 return this->gather(n, indx - (
DIM -1)*(DY + DX*
DIM));
412 void scatterFacesXY(
int x,
int y,
int i1,
int n,
int i2);
413 void scatterEdgesXY(
int x,
int y,
int i1,
int n,
int i2);
414 Word gatherFacesXY(
int x,
int y,
int i1,
int n,
int i2);
416 Word gatherEdgesXY(
int x,
int y,
int i1,
int n,
int i2);
418 template<
int DX,
int DY,
int DZ,
bool Create>
421 const Coord xyz = mOrigin->offsetBy(DX*
DIM, DY*DIM, DZ*DIM);
423 if (leaf)
return &(leaf->getValueMask());
424 if (mAccessor->
isValueOn(xyz))
return &mOnTile;
425 if (!Create)
return &mOffTile;
427 return &(leaf->getValueMask());
431 const Coord* mOrigin;
432 std::vector<MaskType*> mNeighbors;
440 std::unique_ptr<tree::LeafManager<TreeType>> mManagerPtr;
441 tree::LeafManager<TreeType>& mManager;
446 template <
typename TreeT>
451 template <
typename TreeT>
457 template <
typename TreeType>
462 if (iter == 0)
return;
463 const size_t leafCount = mManager.leafCount();
464 if (leafCount == 0)
return;
465 auto& tree = mManager.tree();
495 auto computeWavefront = [&](
const size_t idx) {
496 auto& leaf = manager.
leaf(idx);
497 auto& nodemask = leaf.getValueMask();
498 if (
const auto* original = tree.probeConstLeaf(leaf.origin())) {
499 nodemask ^= original->getValueMask();
510 if (this->getThreaded()) {
512 [&](
const tbb::blocked_range<size_t>&
r){
513 for (
size_t idx =
r.begin(); idx <
r.end(); ++idx) {
514 computeWavefront(idx);
519 for (
size_t idx = 0; idx < manager.
leafCount(); ++idx) {
520 computeWavefront(idx);
528 auto subtractTopology = [&](
const size_t idx) {
529 auto& leaf = mManager.leaf(idx);
530 const auto* maskleaf = mask.probeConstLeaf(leaf.origin());
532 leaf.getValueMask() -= maskleaf->getValueMask();
535 if (this->getThreaded()) {
537 [&](
const tbb::blocked_range<size_t>&
r){
538 for (
size_t idx =
r.begin(); idx <
r.end(); ++idx) {
539 subtractTopology(idx);
544 for (
size_t idx = 0; idx < leafCount; ++idx) {
545 subtractTopology(idx);
553 std::vector<MaskType> nodeMasks;
554 this->copyMasks(nodeMasks);
556 if (this->getThreaded()) {
557 const auto range = mManager.getRange();
558 for (
size_t i = 0; i < iter; ++i) {
562 [&](
const tbb::blocked_range<size_t>&
r) {
565 for (
size_t idx = r.begin(); idx < r.end(); ++idx) {
566 const auto& leaf = mManager.leaf(idx);
567 if (leaf.isEmpty())
continue;
570 cache.
erode(leaf, newMask);
576 [&](
const tbb::blocked_range<size_t>& r){
577 for (
size_t idx = r.begin(); idx < r.end(); ++idx)
578 mManager.leaf(idx).setValueMask(nodeMasks[idx]);
585 for (
size_t i = 0; i < iter; ++i) {
588 for (
size_t idx = 0; idx < leafCount; ++idx) {
589 const auto& leaf = mManager.leaf(idx);
590 if (leaf.isEmpty())
continue;
593 cache.
erode(leaf, newMask);
596 for (
size_t idx = 0; idx < leafCount; ++idx) {
597 mManager.leaf(idx).setValueMask(nodeMasks[idx]);
606 zeroVal<typename TreeType::ValueType>(),
607 this->getThreaded());
608 mManager.rebuild(!this->getThreaded());
612 template <
typename TreeType>
616 const bool preserveMaskLeafNodes)
618 if (iter == 0)
return;
620 const bool threaded = this->getThreaded();
626 auto dilate = [iter, nn,
threaded](
auto& manager,
const bool collapse) {
629 using TreeT =
typename LeafManagerT::TreeType;
631 using LeafT =
typename TreeT::LeafNodeType;
637 TreeT& tree = manager.tree();
642 std::vector<MaskType> nodeMasks;
643 std::vector<std::unique_ptr<LeafT>> nodes;
644 const ValueT& bg = tree.background();
645 const bool steal = iter > 1;
647 for (
size_t i = 0; i < iter; ++i) {
648 if (i > 0) manager.rebuild(!threaded);
650 const size_t leafCount = manager.leafCount();
651 if (leafCount == 0)
return;
660 manager.foreach([&](
auto& leaf,
const size_t idx) {
662 const MaskType& oldMask = nodeMasks[idx];
663 const bool dense = oldMask.isOn();
664 cache.
dilate(leaf, oldMask);
669 accessor.
addTile(1, leaf.origin(), bg,
true);
674 tree.template stealNode<LeafT>(leaf.origin(),
675 zeroVal<ValueT>(),
true));
680 if (nodes.empty())
return;
682 for (
auto& node : nodes) {
683 accessor.
addLeaf(node.release());
693 dilate(mManager, isMask && prune);
694 if (!isMask && prune) {
696 zeroVal<typename TreeType::ValueType>(),
708 std::vector<MaskLeafT*> array;
713 topology.topologyUnion(mManager.tree());
714 array.reserve(mManager.leafCount());
715 topology.stealNodes(array);
717 else if (preserveMaskLeafNodes) {
719 array.resize(mManager.leafCount());
721 [&](
const tbb::blocked_range<size_t>&
r){
722 for (
size_t idx =
r.begin(); idx <
r.end(); ++idx) {
723 array[idx] =
new MaskLeafT(mManager.leaf(idx));
728 array.reserve(mManager.leafCount());
729 mask->stealNodes(array);
733 const size_t numThreads = size_t(tbb::this_task_arena::max_concurrency());
734 const size_t subTreeSize =
math::Max(
size_t(1), array.size()/(2*numThreads));
737 tbb::enumerable_thread_specific<std::unique_ptr<MaskTreeT>>
pool;
739 tbb::parallel_for(tbb::blocked_range<MaskLeafT**>(start, start + array.size(), subTreeSize),
740 [&](
const tbb::blocked_range<MaskLeafT**>&
range) {
742 for (
MaskLeafT** it = range.begin(); it != range.end(); ++it) mask->addLeaf(*it);
745 auto& subtree = pool.local();
746 if (!subtree) subtree = std::move(mask);
751 auto piter = pool.begin();
752 MaskTreeT& subtree = mask ? *mask : **piter++;
753 for (; piter != pool.end(); ++piter) subtree.merge(**piter);
756 if (prune)
tools::prune(subtree, zeroVal<typename MaskTreeT::ValueType>(), threaded);
759 if (!mask) mManager.tree().topologyUnion(subtree,
true);
764 mManager.rebuild(!threaded);
768 template <
typename TreeType>
772 for (
int x = 0;
x <
DIM; ++
x) {
775 if (
Word&
w = mask.template getWord<Word>(n)) {
778 (
Word(
w<<1 | (this->
template gather<0,0,-1>(1, n)>>(DIM-1))) &
779 Word(
w>>1 | (this->
template gather<0,0, 1>(2, n)<<(DIM-1)))));
780 w =
Word(
w & this->gatherFacesXY(
x,
y, 0, n, 3));
786 template <
typename TreeType>
788 Morphology<TreeType>::NodeMaskOp::dilate6(
const MaskType& mask)
790 for (
int x = 0;
x < DIM; ++
x ) {
791 for (
int y = 0, n = (
x << LOG2DIM);
794 if (
const Word
w = mask.template getWord<Word>(n)) {
796 this->mWord = Word(
w | (
w>>1) | (
w<<1));
799 if ( (this->mWord = Word(
w<<(DIM-1))) ) {
800 this->
template scatter< 0, 0,-1>(1,
n);
803 if ( (this->mWord = Word(
w>>(DIM-1))) ) {
804 this->
template scatter< 0, 0, 1>(2,
n);
808 this->scatterFacesXY(
x,
y, 0, n, 3);
814 template <
typename TreeType>
816 Morphology<TreeType>::NodeMaskOp::dilate18(
const MaskType& mask)
819 const Coord origin = this->getOrigin();
820 const Coord orig_mz = origin.offsetBy(0, 0, -DIM);
821 const Coord orig_pz = origin.offsetBy(0, 0, DIM);
822 for (
int x = 0;
x < DIM; ++
x ) {
823 for (
int y = 0, n = (
x << LOG2DIM);
y < DIM; ++
y, ++
n) {
824 if (
const Word
w = mask.template getWord<Word>(n)) {
826 this->mWord = Word(
w | (
w>>1) | (
w<<1));
829 this->scatterFacesXY(
x,
y, 0, n, 3);
831 this->scatterEdgesXY(
x,
y, 0, n, 3);
833 if ( (this->mWord = Word(
w<<(DIM-1))) ) {
835 this->
template scatter< 0, 0,-1>(1,
n);
837 this->scatterFacesXY(
x,
y, 1, n, 11);
839 if ( (this->mWord = Word(
w>>(DIM-1))) ) {
841 this->
template scatter< 0, 0, 1>(2,
n);
843 this->scatterFacesXY(
x,
y, 2, n, 15);
851 template <
typename TreeType>
853 Morphology<TreeType>::NodeMaskOp::dilate26(
const MaskType& mask)
856 const Coord origin = this->getOrigin();
857 const Coord orig_mz = origin.offsetBy(0, 0, -DIM);
858 const Coord orig_pz = origin.offsetBy(0, 0, DIM);
859 for (
int x = 0;
x < DIM; ++
x) {
860 for (
int y = 0, n = (
x << LOG2DIM);
y < DIM; ++
y, ++
n) {
861 if (
const Word
w = mask.template getWord<Word>(n)) {
863 this->mWord = Word(
w | (
w>>1) | (
w<<1));
866 this->scatterFacesXY(
x,
y, 0, n, 3);
867 this->scatterEdgesXY(
x,
y, 0, n, 3);
869 if ( (this->mWord = Word(
w<<(DIM-1))) ) {
871 this->
template scatter< 0, 0,-1>(1,
n);
873 this->scatterFacesXY(
x,
y, 1, n, 11);
874 this->scatterEdgesXY(
x,
y, 1, n, 11);
876 if ( (this->mWord = Word(
w>>(DIM-1))) ) {
878 this->
template scatter< 0, 0, 1>(2,
n);
880 this->scatterFacesXY(
x,
y, 2, n, 19);
881 this->scatterEdgesXY(
x,
y, 2, n, 19);
888 template<
typename TreeType>
890 Morphology<TreeType>::NodeMaskOp::scatterFacesXY(
int x,
int y,
int i1,
int n,
int i2)
894 this->scatter(i1, n-DIM);
896 this->
template scatter<-1, 0, 0>(
i2,
n);
900 this->scatter(i1, n+DIM);
902 this->
template scatter< 1, 0, 0>(i2+1,
n);
906 this->scatter(i1, n-1);
908 this->
template scatter< 0,-1, 0>(i2+2,
n);
912 this->scatter(i1, n+1);
914 this->
template scatter< 0, 1, 0>(i2+3,
n);
919 template<
typename TreeType>
921 Morphology<TreeType>::NodeMaskOp::scatterEdgesXY(
int x,
int y,
int i1,
int n,
int i2)
925 this->scatter(i1, n-DIM-1);
927 this->
template scatter< 0,-1, 0>(i2+2, n-DIM);
930 this->scatter(i1, n-DIM+1);
932 this->
template scatter< 0, 1, 0>(i2+3, n-DIM);
936 this->
template scatter<-1, 0, 0>(
i2 , n+1);
938 this->
template scatter<-1, 1, 0>(i2+7,
n );
941 this->
template scatter<-1, 0, 0>(
i2 , n-1);
943 this->
template scatter<-1,-1, 0>(i2+4,
n );
948 this->scatter(i1, n+DIM-1);
950 this->
template scatter< 0,-1, 0>(i2+2, n+DIM);
953 this->scatter(i1, n+DIM+1);
955 this->
template scatter< 0, 1, 0>(i2+3, n+DIM);
959 this->
template scatter< 1, 0, 0>(i2+1, n-1);
961 this->
template scatter< 1,-1, 0>(i2+6,
n );
964 this->
template scatter< 1, 0, 0>(i2+1, n+1);
966 this->
template scatter< 1, 1, 0>(i2+5,
n );
972 template<
typename TreeType>
974 Morphology<TreeType>::NodeMaskOp::gatherFacesXY(
int x,
int y,
int i1,
int n,
int i2)
978 this->gather(i1, n - DIM) :
979 this->template gather<-1,0,0>(i2, n);
982 w = Word(w & (x < DIM - 1 ?
983 this->gather(i1, n + DIM) :
984 this->
template gather<1,0,0>(i2 + 1, n)));
987 w = Word(w & (y > 0 ?
988 this->gather(i1, n - 1) :
989 this->
template gather<0,-1,0>(i2 + 2, n)));
992 w = Word(w & (y < DIM - 1 ?
993 this->gather(i1, n + 1) :
994 this->
template gather<0,1,0>(i2+3, n)));
1000 template<
typename TreeType>
1002 Morphology<TreeType>::NodeMaskOp::gatherEdgesXY(
int x,
int y,
int i1,
int n,
int i2)
1007 w &= y > 0 ? this->gather(i1, n-DIM-1) :
1008 this->template gather< 0,-1, 0>(i2+2, n-DIM);
1009 w &= y < DIM-1 ? this->gather(i1, n-DIM+1) :
1010 this->template gather< 0, 1, 0>(i2+3, n-DIM);
1012 w &= y < DIM-1 ? this->
template gather<-1, 0, 0>(
i2 , n+1):
1013 this->
template gather<-1, 1, 0>(i2+7, n );
1014 w &= y > 0 ? this->
template gather<-1, 0, 0>(
i2 , n-1):
1015 this->
template gather<-1,-1, 0>(i2+4, n );
1018 w &= y > 0 ? this->gather(i1, n+DIM-1) :
1019 this->template gather< 0,-1, 0>(i2+2, n+DIM);
1020 w &= y < DIM-1 ? this->gather(i1, n+DIM+1) :
1021 this->template gather< 0, 1, 0>(i2+3, n+DIM);
1023 w &= y > 0 ? this->
template gather< 1, 0, 0>(i2+1, n-1):
1024 this->
template gather< 1,-1, 0>(i2+6, n );
1025 w &= y < DIM-1 ? this->
template gather< 1, 0, 0>(i2+1, n+1):
1026 this->
template gather< 1, 1, 0>(i2+5, n );
1040 namespace morph_internal {
1041 template <
typename T>
struct Adapter {
1043 static TreeType&
get(
T& tree) {
return tree; }
1044 static void sync(T&) {}
1046 template <
typename T>
1047 struct Adapter<openvdb::tree::LeafManager<T>> {
1049 static TreeType&
get(openvdb::tree::LeafManager<T>& M) {
return M.tree(); }
1050 static void sync(openvdb::tree::LeafManager<T>& M) { M.rebuild(); }
1056 template<
typename TreeOrLeafManagerT>
1058 const int iterations,
1063 using AdapterT = morph_internal::Adapter<TreeOrLeafManagerT>;
1064 using TreeT =
typename AdapterT::TreeType;
1067 if (iterations <= 0)
return;
1073 morph.
dilateVoxels(static_cast<size_t>(iterations), nn,
false);
1088 tree.voxelizeActiveTiles();
1089 AdapterT::sync(treeOrLeafM);
1094 morph.
dilateVoxels(static_cast<size_t>(iterations), nn,
true);
1098 morph.
dilateVoxels(static_cast<size_t>(iterations), nn,
false);
1115 topology.topologyUnion(tree);
1116 topology.voxelizeActiveTiles();
1120 morph.
dilateVoxels(static_cast<size_t>(iterations), nn,
true);
1122 tree.topologyUnion(topology,
true);
1128 tools::prune(tree, zeroVal<typename TreeT::ValueType>(), threaded);
1129 AdapterT::sync(treeOrLeafM);
1133 template<
typename TreeOrLeafManagerT>
1135 const int iterations,
1140 using AdapterT = morph_internal::Adapter<TreeOrLeafManagerT>;
1141 using TreeT =
typename AdapterT::TreeType;
1144 if (iterations <= 0)
return;
1152 topology.topologyUnion(tree);
1153 topology.voxelizeActiveTiles();
1158 morph.
erodeVoxels(static_cast<size_t>(iterations), nn,
false);
1163 tools::prune(topology, zeroVal<typename MaskT::ValueType>(), threaded);
1164 tree.topologyIntersection(topology);
1165 AdapterT::sync(treeOrLeafM);
1173 if (tree.hasActiveTiles()) {
1174 tree.voxelizeActiveTiles();
1175 AdapterT::sync(treeOrLeafM);
1182 morph.
erodeVoxels(static_cast<size_t>(iterations), nn,
false);
1191 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
1193 #ifdef OPENVDB_INSTANTIATE_MORPHOLOGY
1197 #define _FUNCTION(TreeT) \
1198 void dilateActiveValues(TreeT&, \
1199 const int, const NearestNeighbors, const TilePolicy, const bool)
1203 #define _FUNCTION(TreeT) \
1204 void dilateActiveValues(tree::LeafManager<TreeT>&, \
1205 const int, const NearestNeighbors, const TilePolicy, const bool)
1209 #define _FUNCTION(TreeT) \
1210 void erodeActiveValues(TreeT&, \
1211 const int, const NearestNeighbors, const TilePolicy, const bool)
1215 #define _FUNCTION(TreeT) \
1216 void erodeActiveValues(tree::LeafManager<TreeT>&, \
1217 const int, const NearestNeighbors, const TilePolicy, const bool)
1221 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
1228 #endif // OPENVDB_TOOLS_MORPHOLOGY_HAS_BEEN_INCLUDED
bool isValueOn(const Coord &xyz) const
Return the active state of the voxel at the given coordinates.
GLsizei const GLfloat * value
#define OPENVDB_USE_VERSION_NAMESPACE
**But if you need a or simply need to know when the task has note that the like this
The Value Accessor Implementation and API methods. The majoirty of the API matches the API of a compa...
#define OPENVDB_ASSERT(X)
#define OPENVDB_ALL_TREE_INSTANTIATE(Function)
size_t leafCount() const
Return the number of leaf nodes.
GLint GLint GLsizei GLint GLenum GLenum type
ValueAccessors are designed to help accelerate accesses into the OpenVDB Tree structures by storing c...
ImageBuf OIIO_API dilate(const ImageBuf &src, int width=3, int height=-1, ROI roi={}, int nthreads=0)
GT_API const UT_StringHolder topology
Implementation of topological activation/deactivation.
__hostdev__ void setOrigin(const T &ijk)
OIIO_UTIL_API void parallel_for(int32_t begin, int32_t end, function_view< void(int32_t)> task, paropt opt=0)
Defined various multi-threaded utility functions for trees.
auto get(const UT_ARTIterator< T > &it) -> decltype(it.key())
RangeType getRange(size_t grainsize=1) const
Return a tbb::blocked_range of leaf array indices.
LeafNodeT * touchLeaf(const Coord &xyz)
Returns the leaf node that contains voxel (x, y, z) and if it doesn't exist, create it...
LeafNodeT * probeLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains the voxel coordinate xyz. If no LeafNode exists...
GLubyte GLubyte GLubyte GLubyte w
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...
LeafType & leaf(size_t leafIdx) const
Return a pointer to the leaf node at index leafIdx in the array.
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Attribute-owned data structure for points. Point attributes are stored in leaf nodes and ordered by v...
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.
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
void addTile(Index level, const Coord &xyz, const ValueType &value, bool state)
Add a tile at the specified tree level that contains the coordinate xyz, possibly deleting existing n...
#define OPENVDB_THROW(exception, message)
**Note that the tasks the is the thread number *for the pool