150namespace morphology {
153template<
typename TreeType>
160 using MaskTreeT =
typename TreeType::template ValueConverter<ValueMask>::Type;
165 : mManagerPtr(new tree::LeafManager<TreeType>(tree))
166 , mManager(*mManagerPtr)
170 : mManagerPtr(nullptr)
178 inline void setThreaded(
const bool threaded) { mThreaded = threaded; }
190 void erodeVoxels(
const size_t iter,
192 const bool prune =
false);
205 void dilateVoxels(
const size_t iter,
207 const bool prune =
false,
208 const bool preserveMaskLeafNodes =
false);
216 if (masks.size() < mManager.leafCount()) {
217 masks.resize(mManager.leafCount());
220 if (this->getThreaded()) {
222 tbb::parallel_for(mManager.getRange(),
223 [&](
const tbb::blocked_range<size_t>& r){
224 for (size_t idx = r.begin(); idx < r.end(); ++idx)
225 masks[idx] = mManager.leaf(idx).getValueMask();
229 for (
size_t idx = 0; idx < mManager.leafCount(); ++idx) {
230 masks[idx] = mManager.leaf(idx).getValueMask();
244 static const Int32 LOG2DIM =
static_cast<Int32>(LeafType::LOG2DIM);
247 using Word =
typename std::conditional<LOG2DIM == 3, uint8_t,
248 typename std::conditional<LOG2DIM == 4, uint16_t,
249 typename std::conditional<LOG2DIM == 5, uint32_t,
250 typename std::conditional<LOG2DIM == 6, uint64_t,
251 void>::type>::type>::type>::type;
253 static_assert(!std::is_same<Word, void>::value,
254 "Unsupported Node Dimension for node mask dilation/erosion");
260 , mAccessor(&accessor)
276 const MaskType mask = leaf.getValueMask();
277 this->dilate(leaf, mask);
295 mNeighbors[0] = &(leaf.getValueMask());
296 this->setOrigin(leaf.origin());
300 case NN_FACE : { this->dilate6(mask);
return; }
302 assert(
false &&
"Unknown op during dilation.");
return;
320 MaskType mask = leaf.getValueMask();
321 this->erode(leaf, mask);
341 mNeighbors[0] =
const_cast<MaskType*
>(&leaf.getValueMask());
342 this->setOrigin(leaf.origin());
346 case NN_FACE : { this->erode6(mask);
return; }
348 assert(
false &&
"Unknown op during erosion.");
return;
363 void dilate6(
const MaskType& mask);
364 void dilate18(
const MaskType& mask);
365 void dilate26(
const MaskType& mask);
366 void erode6(MaskType& mask);
375 inline void setOrigin(
const Coord& origin) { mOrigin = &origin; }
376 inline const Coord& getOrigin()
const {
return *mOrigin; }
377 inline void clear() { std::fill(mNeighbors.begin(), mNeighbors.end(),
nullptr); }
379 inline void scatter(
size_t n,
int indx)
381 assert(n < mNeighbors.size());
382 assert(mNeighbors[n]);
383 mNeighbors[n]->template getWord<Word>(indx) |= mWord;
386 template<
int DX,
int DY,
int DZ>
387 inline void scatter(
size_t n,
int indx)
389 assert(n < mNeighbors.size());
390 if (!mNeighbors[n]) {
391 mNeighbors[n] = this->getNeighbor<DX,DY,DZ,true>();
393 assert(mNeighbors[n]);
394 this->scatter(n, indx - (DIM - 1)*(DY + DX*DIM));
396 inline Word gather(
size_t n,
int indx)
398 assert(n < mNeighbors.size());
399 return mNeighbors[n]->template getWord<Word>(indx);
401 template<
int DX,
int DY,
int DZ>
402 inline Word gather(
size_t n,
int indx)
404 assert(n < mNeighbors.size());
405 if (!mNeighbors[n]) {
406 mNeighbors[n] = this->getNeighbor<DX,DY,DZ,false>();
408 return this->gather(n, indx - (DIM -1)*(DY + DX*DIM));
411 void scatterFacesXY(
int x,
int y,
int i1,
int n,
int i2);
412 void scatterEdgesXY(
int x,
int y,
int i1,
int n,
int i2);
413 Word gatherFacesXY(
int x,
int y,
int i1,
int n,
int i2);
415 Word gatherEdgesXY(
int x,
int y,
int i1,
int n,
int i2);
417 template<
int DX,
int DY,
int DZ,
bool Create>
418 inline MaskType* getNeighbor()
420 const Coord xyz = mOrigin->offsetBy(DX*DIM, DY*DIM, DZ*DIM);
421 auto* leaf = mAccessor->probeLeaf(xyz);
422 if (leaf)
return &(leaf->getValueMask());
423 if (mAccessor->isValueOn(xyz))
return &mOnTile;
424 if (!Create)
return &mOffTile;
425 leaf = mAccessor->touchLeaf(xyz);
426 return &(leaf->getValueMask());
430 const Coord* mOrigin;
431 std::vector<MaskType*> mNeighbors;
432 AccessorType*
const mAccessor;
434 MaskType mOnTile, mOffTile;
439 std::unique_ptr<tree::LeafManager<TreeType>> mManagerPtr;
440 tree::LeafManager<TreeType>& mManager;
445template <
typename TreeT>
446typename std::enable_if<std::is_same<TreeT, typename TreeT::template ValueConverter<ValueMask>::Type>::value,
447 typename TreeT::template ValueConverter<ValueMask>::Type*>::type
450template <
typename TreeT>
451typename std::enable_if<!std::is_same<TreeT, typename TreeT::template ValueConverter<ValueMask>::Type>::value,
452 typename TreeT::template ValueConverter<ValueMask>::Type*>::type
456template <
typename TreeType>
461 if (iter == 0)
return;
462 const size_t leafCount = mManager.leafCount();
463 if (leafCount == 0)
return;
464 auto& tree = mManager.tree();
494 auto computeWavefront = [&](
const size_t idx) {
495 auto& leaf = manager.
leaf(idx);
496 auto& nodemask = leaf.getValueMask();
497 if (
const auto* original = tree.probeConstLeaf(leaf.origin())) {
498 nodemask ^= original->getValueMask();
505 assert(!nodemask.isOn());
509 if (this->getThreaded()) {
510 tbb::parallel_for(manager.
getRange(),
511 [&](
const tbb::blocked_range<size_t>& r){
512 for (size_t idx = r.begin(); idx < r.end(); ++idx) {
513 computeWavefront(idx);
518 for (
size_t idx = 0; idx < manager.leafCount(); ++idx) {
519 computeWavefront(idx);
524 m.dilateVoxels(iter, nn,
false);
527 auto subtractTopology = [&](
const size_t idx) {
528 auto& leaf = mManager.leaf(idx);
529 const auto* maskleaf = mask.probeConstLeaf(leaf.origin());
531 leaf.getValueMask() -= maskleaf->getValueMask();
534 if (this->getThreaded()) {
535 tbb::parallel_for(mManager.getRange(),
536 [&](
const tbb::blocked_range<size_t>& r){
537 for (size_t idx = r.begin(); idx < r.end(); ++idx) {
538 subtractTopology(idx);
543 for (
size_t idx = 0; idx < leafCount; ++idx) {
544 subtractTopology(idx);
552 std::vector<MaskType> nodeMasks;
553 this->copyMasks(nodeMasks);
555 if (this->getThreaded()) {
556 const auto range = mManager.getRange();
557 for (
size_t i = 0; i < iter; ++i) {
560 tbb::parallel_for(range,
561 [&](
const tbb::blocked_range<size_t>& r) {
562 AccessorType accessor(tree);
563 NodeMaskOp cache(accessor, nn);
564 for (
size_t idx = r.begin(); idx < r.end(); ++idx) {
565 const auto& leaf = mManager.leaf(idx);
566 if (leaf.isEmpty())
continue;
568 MaskType& newMask = nodeMasks[idx];
569 cache.erode(leaf, newMask);
574 tbb::parallel_for(range,
575 [&](
const tbb::blocked_range<size_t>& r){
576 for (
size_t idx = r.begin(); idx < r.end(); ++idx)
577 mManager.leaf(idx).setValueMask(nodeMasks[idx]);
582 AccessorType accessor(tree);
583 NodeMaskOp cache(accessor, nn);
584 for (
size_t i = 0; i < iter; ++i) {
587 for (
size_t idx = 0; idx < leafCount; ++idx) {
588 const auto& leaf = mManager.leaf(idx);
589 if (leaf.isEmpty())
continue;
591 MaskType& newMask = nodeMasks[idx];
592 cache.erode(leaf, newMask);
595 for (
size_t idx = 0; idx < leafCount; ++idx) {
596 mManager.leaf(idx).setValueMask(nodeMasks[idx]);
604 tools::prune(mManager.tree(),
605 zeroVal<typename TreeType::ValueType>(),
606 this->getThreaded());
607 mManager.rebuild(!this->getThreaded());
611template <
typename TreeType>
615 const bool preserveMaskLeafNodes)
617 if (iter == 0)
return;
619 const bool threaded = this->getThreaded();
625 auto dilate = [iter, nn, threaded](
auto& manager,
const bool collapse) {
627 using LeafManagerT =
typename std::decay<
decltype(manager)>::type;
628 using TreeT =
typename LeafManagerT::TreeType;
629 using ValueT =
typename TreeT::ValueType;
630 using LeafT =
typename TreeT::LeafNodeType;
636 TreeT& tree = manager.tree();
641 std::vector<MaskType> nodeMasks;
642 std::vector<std::unique_ptr<LeafT>> nodes;
643 const ValueT& bg = tree.background();
644 const bool steal = iter > 1;
646 for (
size_t i = 0; i < iter; ++i) {
647 if (i > 0) manager.rebuild(!threaded);
649 const size_t leafCount = manager.leafCount();
650 if (leafCount == 0)
return;
659 manager.foreach([&](
auto& leaf,
const size_t idx) {
661 const MaskType& oldMask = nodeMasks[idx];
662 const bool dense = oldMask.isOn();
663 cache.
dilate(leaf, oldMask);
668 accessor.
addTile(1, leaf.origin(), bg,
true);
673 tree.template stealNode<LeafT>(leaf.origin(),
674 zeroVal<ValueT>(),
true));
679 if (nodes.empty())
return;
681 for (
auto& node : nodes) {
682 accessor.
addLeaf(node.release());
691 constexpr bool isMask = std::is_same<TreeType, MaskTreeT>::value;
692 dilate(mManager, isMask &&
prune);
693 if (!isMask &&
prune) {
694 tools::prune(mManager.tree(),
695 zeroVal<typename TreeType::ValueType>(),
707 std::vector<MaskLeafT*> array;
712 topology.topologyUnion(mManager.tree());
713 array.reserve(mManager.leafCount());
714 topology.stealNodes(array);
716 else if (preserveMaskLeafNodes) {
718 array.resize(mManager.leafCount());
719 tbb::parallel_for(mManager.getRange(),
720 [&](
const tbb::blocked_range<size_t>& r){
721 for (size_t idx = r.begin(); idx < r.end(); ++idx) {
722 array[idx] = new MaskLeafT(mManager.leaf(idx));
727 array.reserve(mManager.leafCount());
728 mask->stealNodes(array);
732 const size_t numThreads = size_t(tbb::this_task_arena::max_concurrency());
733 const size_t subTreeSize = math::Max(
size_t(1), array.size()/(2*numThreads));
736 tbb::enumerable_thread_specific<std::unique_ptr<MaskTreeT>> pool;
738 tbb::parallel_for(tbb::blocked_range<MaskLeafT**>(start, start + array.size(), subTreeSize),
739 [&](
const tbb::blocked_range<MaskLeafT**>& range) {
740 std::unique_ptr<MaskTreeT> mask(new MaskTreeT);
741 for (MaskLeafT** it = range.begin(); it != range.end(); ++it) mask->addLeaf(*it);
742 tree::LeafManager<MaskTreeT> manager(*mask, range.begin(), range.end());
743 dilate(manager, prune);
744 auto& subtree = pool.local();
745 if (!subtree) subtree = std::move(mask);
746 else subtree->merge(*mask, MERGE_ACTIVE_STATES);
750 auto piter = pool.begin();
751 MaskTreeT& subtree = mask ? *mask : **piter++;
752 for (; piter != pool.end(); ++piter) subtree.merge(**piter);
755 if (
prune) tools::prune(subtree, zeroVal<typename MaskTreeT::ValueType>(), threaded);
758 if (!mask) mManager.tree().topologyUnion(subtree,
true);
763 mManager.rebuild(!threaded);
767template <
typename TreeType>
769Morphology<TreeType>::NodeMaskOp::erode6(MaskType& mask)
771 for (
int x = 0; x < DIM; ++x) {
772 for (
int y = 0, n = (x << LOG2DIM); y < DIM; ++y, ++n) {
774 if (Word& w = mask.template getWord<Word>(n)) {
777 (Word(w<<1 | (this->
template gather<0,0,-1>(1, n)>>(DIM-1))) &
778 Word(w>>1 | (this->
template gather<0,0, 1>(2, n)<<(DIM-1)))));
779 w = Word(w & this->gatherFacesXY(x, y, 0, n, 3));
785template <
typename TreeType>
787Morphology<TreeType>::NodeMaskOp::dilate6(
const MaskType& mask)
789 for (
int x = 0; x < DIM; ++x ) {
790 for (
int y = 0, n = (x << LOG2DIM);
793 if (
const Word w = mask.template getWord<Word>(n)) {
795 this->mWord = Word(w | (w>>1) | (w<<1));
798 if ( (this->mWord = Word(w<<(DIM-1))) ) {
799 this->
template scatter< 0, 0,-1>(1, n);
802 if ( (this->mWord = Word(w>>(DIM-1))) ) {
803 this->
template scatter< 0, 0, 1>(2, n);
807 this->scatterFacesXY(x, y, 0, n, 3);
813template <
typename TreeType>
815Morphology<TreeType>::NodeMaskOp::dilate18(
const MaskType& mask)
818 const Coord origin = this->getOrigin();
819 const Coord orig_mz = origin.offsetBy(0, 0, -DIM);
820 const Coord orig_pz = origin.offsetBy(0, 0, DIM);
821 for (
int x = 0; x < DIM; ++x ) {
822 for (
int y = 0, n = (x << LOG2DIM); y < DIM; ++y, ++n) {
823 if (
const Word w = mask.template getWord<Word>(n)) {
825 this->mWord = Word(w | (w>>1) | (w<<1));
826 this->setOrigin(origin);
828 this->scatterFacesXY(x, y, 0, n, 3);
830 this->scatterEdgesXY(x, y, 0, n, 3);
832 if ( (this->mWord = Word(w<<(DIM-1))) ) {
833 this->setOrigin(origin);
834 this->
template scatter< 0, 0,-1>(1, n);
835 this->setOrigin(orig_mz);
836 this->scatterFacesXY(x, y, 1, n, 11);
838 if ( (this->mWord = Word(w>>(DIM-1))) ) {
839 this->setOrigin(origin);
840 this->
template scatter< 0, 0, 1>(2, n);
841 this->setOrigin(orig_pz);
842 this->scatterFacesXY(x, y, 2, n, 15);
850template <
typename TreeType>
852Morphology<TreeType>::NodeMaskOp::dilate26(
const MaskType& mask)
855 const Coord origin = this->getOrigin();
856 const Coord orig_mz = origin.offsetBy(0, 0, -DIM);
857 const Coord orig_pz = origin.offsetBy(0, 0, DIM);
858 for (
int x = 0; x < DIM; ++x) {
859 for (
int y = 0, n = (x << LOG2DIM); y < DIM; ++y, ++n) {
860 if (
const Word w = mask.template getWord<Word>(n)) {
862 this->mWord = Word(w | (w>>1) | (w<<1));
863 this->setOrigin(origin);
865 this->scatterFacesXY(x, y, 0, n, 3);
866 this->scatterEdgesXY(x, y, 0, n, 3);
868 if ( (this->mWord = Word(w<<(DIM-1))) ) {
869 this->setOrigin(origin);
870 this->
template scatter< 0, 0,-1>(1, n);
871 this->setOrigin(orig_mz);
872 this->scatterFacesXY(x, y, 1, n, 11);
873 this->scatterEdgesXY(x, y, 1, n, 11);
875 if ( (this->mWord = Word(w>>(DIM-1))) ) {
876 this->setOrigin(origin);
877 this->
template scatter< 0, 0, 1>(2, n);
878 this->setOrigin(orig_pz);
879 this->scatterFacesXY(x, y, 2, n, 19);
880 this->scatterEdgesXY(x, y, 2, n, 19);
887template<
typename TreeType>
889Morphology<TreeType>::NodeMaskOp::scatterFacesXY(
int x,
int y,
int i1,
int n,
int i2)
893 this->scatter(i1, n-DIM);
895 this->
template scatter<-1, 0, 0>(i2, n);
899 this->scatter(i1, n+DIM);
901 this->
template scatter< 1, 0, 0>(i2+1, n);
905 this->scatter(i1, n-1);
907 this->
template scatter< 0,-1, 0>(i2+2, n);
911 this->scatter(i1, n+1);
913 this->
template scatter< 0, 1, 0>(i2+3, n);
918template<
typename TreeType>
920Morphology<TreeType>::NodeMaskOp::scatterEdgesXY(
int x,
int y,
int i1,
int n,
int i2)
924 this->scatter(i1, n-DIM-1);
926 this->
template scatter< 0,-1, 0>(i2+2, n-DIM);
929 this->scatter(i1, n-DIM+1);
931 this->
template scatter< 0, 1, 0>(i2+3, n-DIM);
935 this->
template scatter<-1, 0, 0>(i2 , n+1);
937 this->
template scatter<-1, 1, 0>(i2+7, n );
940 this->
template scatter<-1, 0, 0>(i2 , n-1);
942 this->
template scatter<-1,-1, 0>(i2+4, n );
947 this->scatter(i1, n+DIM-1);
949 this->
template scatter< 0,-1, 0>(i2+2, n+DIM);
952 this->scatter(i1, n+DIM+1);
954 this->
template scatter< 0, 1, 0>(i2+3, n+DIM);
958 this->
template scatter< 1, 0, 0>(i2+1, n-1);
960 this->
template scatter< 1,-1, 0>(i2+6, n );
963 this->
template scatter< 1, 0, 0>(i2+1, n+1);
965 this->
template scatter< 1, 1, 0>(i2+5, n );
971template<
typename TreeType>
972inline typename Morphology<TreeType>::NodeMaskOp::Word
973Morphology<TreeType>::NodeMaskOp::gatherFacesXY(
int x,
int y,
int i1,
int n,
int i2)
977 this->gather(i1, n - DIM) :
978 this->template gather<-1,0,0>(i2, n);
981 w = Word(w & (x < DIM - 1 ?
982 this->gather(i1, n + DIM) :
983 this->template gather<1,0,0>(i2 + 1, n)));
986 w = Word(w & (y > 0 ?
987 this->gather(i1, n - 1) :
988 this->template gather<0,-1,0>(i2 + 2, n)));
991 w = Word(w & (y < DIM - 1 ?
992 this->gather(i1, n + 1) :
993 this->template gather<0,1,0>(i2+3, n)));
999template<
typename TreeType>
1000inline typename Morphology<TreeType>::NodeMaskOp::Word
1001Morphology<TreeType>::NodeMaskOp::gatherEdgesXY(
int x,
int y,
int i1,
int n,
int i2)
1006 w &= y > 0 ? this->gather(i1, n-DIM-1) :
1007 this->template gather< 0,-1, 0>(i2+2, n-DIM);
1008 w &= y < DIM-1 ? this->gather(i1, n-DIM+1) :
1009 this->template gather< 0, 1, 0>(i2+3, n-DIM);
1011 w &= y < DIM-1 ? this->
template gather<-1, 0, 0>(i2 , n+1):
1012 this->template gather<-1, 1, 0>(i2+7, n );
1013 w &= y > 0 ? this->
template gather<-1, 0, 0>(i2 , n-1):
1014 this->template gather<-1,-1, 0>(i2+4, n );
1017 w &= y > 0 ? this->gather(i1, n+DIM-1) :
1018 this->template gather< 0,-1, 0>(i2+2, n+DIM);
1019 w &= y < DIM-1 ? this->gather(i1, n+DIM+1) :
1020 this->template gather< 0, 1, 0>(i2+3, n+DIM);
1022 w &= y > 0 ? this->
template gather< 1, 0, 0>(i2+1, n-1):
1023 this->template gather< 1,-1, 0>(i2+6, n );
1024 w &= y < DIM-1 ? this->
template gather< 1, 0, 0>(i2+1, n+1):
1025 this->template gather< 1, 1, 0>(i2+5, n );