14 #ifndef OPENVDB_TOOLS_FILTER_HAS_BEEN_INCLUDED
15 #define OPENVDB_TOOLS_FILTER_HAS_BEEN_INCLUDED
17 #include <tbb/parallel_for.h>
30 #include <type_traits>
31 #include <tbb/concurrent_vector.h>
39 template<
typename GridT,
40 typename MaskT =
typename GridT::template ValueConverter<float>::Type,
41 typename InterruptT = util::NullInterrupter>
48 using LeafType =
typename TreeType::LeafNodeType;
52 using RangeType =
typename LeafManagerType::LeafRange;
54 static_assert(std::is_floating_point<AlphaType>::value,
55 "openvdb::tools::Filter requires a mask grid with floating-point values");
60 Filter(GridT& grid, InterruptT* interrupt =
nullptr)
63 , mInterrupter(interrupt)
77 , mInterrupter(other.mInterrupter)
79 , mGrainSize(other.mGrainSize)
80 , mMinMask(other.mMinMask)
81 , mMaxMask(other.mMaxMask)
82 , mInvertMask(other.mInvertMask)
83 , mTiles(other.mTiles) {}
131 void mean(
int width = 1,
int iterations = 1,
const MaskType* mask =
nullptr);
140 void gaussian(
int width = 1,
int iterations = 1,
const MaskType* mask =
nullptr);
148 void median(
int width = 1,
int iterations = 1,
const MaskType* mask =
nullptr);
153 void offset(ValueType offset,
const MaskType* mask =
nullptr);
161 if (mTask) mTask(
const_cast<Filter*
>(
this), range);
166 using LeafT =
typename TreeType::LeafNodeType;
167 using VoxelIterT =
typename LeafT::ValueOnIter;
168 using VoxelCIterT =
typename LeafT::ValueOnCIter;
170 using LeafIterT =
typename RangeType::Iterator;
173 void cook(LeafManagerType& leafs);
175 template<
size_t Axis>
177 Avg(
const GridT* grid,
Int32 w): acc(grid->tree()), width(w), frac(1.f/float(2*w+1)) {}
178 inline ValueType operator()(
Coord xyz);
179 typename GridT::ConstAccessor acc;
185 template <
typename AvgT>
186 void doBox(
const RangeType& r,
Int32 w);
187 void doBoxX(
const RangeType& r,
Int32 w) { this->doBox<Avg<0> >(r,w); }
188 void doBoxY(
const RangeType& r,
Int32 w) { this->doBox<Avg<1> >(r,w); }
189 void doBoxZ(
const RangeType& r,
Int32 w) { this->doBox<Avg<2> >(r,w); }
190 void doMedian(
const RangeType&,
int);
191 void doOffset(
const RangeType&, ValueType);
196 typename std::function<void (Filter*,
const RangeType&)> mTask;
197 InterruptT* mInterrupter;
198 const MaskType* mMask;
200 AlphaType mMinMask, mMaxMask;
209 namespace filter_internal {
211 template<
typename TreeT>
217 using MaskT =
typename TreeT::template ValueConverter<ValueMask>::Type;
219 Voxelizer(TreeT& tree,
const bool allNeighbours,
const size_t grainSize)
222 , mGrainSize(grainSize)
223 , mOp(tree, mVoxelTopology, allNeighbours ? 26 : 6) {}
233 if (!mOp.tree().hasActiveTiles())
return 0;
236 for (
int i = 0; i < width; i += int(TreeT::LeafNodeType::DIM), ++count) {
237 if (i > 0) mManager->rebuild();
238 mManager->foreachBottomUp(mOp, mGrainSize > 0, mGrainSize);
239 mOp.tree().topologyUnion(mVoxelTopology);
253 mVoxelTopology.topologyUnion(mOp.tree());
254 mManager.reset(
new NodeManagerT(mOp.tree()));
258 struct CreateVoxelMask
260 using LeafT =
typename TreeT::LeafNodeType;
261 using RootT =
typename TreeT::RootNodeType;
263 CreateVoxelMask(TreeT& tree, MaskT& mask,
const size_t NN)
264 : mTree(tree), mVoxelTopology(mask), mNeighbours(NN) {}
266 TreeT& tree() {
return mTree; }
270 void operator()(
const LeafT&)
const { assert(
false); }
272 void operator()(
const RootT& node)
const
274 using ChildT =
typename RootT::ChildNodeType;
275 static constexpr
Int32 CHILDDIM =
Int32(ChildT::DIM);
276 static constexpr
Int32 LEAFDIM =
Int32(LeafT::DIM);
277 const Tester op(mTree, mNeighbours);
280 [&](
const Coord& ijk,
286 Int32& a = offset[axis1];
287 Int32& b = offset[axis2];
288 for (a = 0; a < CHILDDIM; a+=LEAFDIM) {
289 for (b = 0; b < CHILDDIM; b+=LEAFDIM) {
290 const Coord childijk = ijk + offset;
291 if (op.test(childijk, val)) {
292 mVoxelTopology.touchLeaf(childijk);
297 offset.reset(CHILDDIM-1);
298 for (a = 0; a < CHILDDIM; a+=LEAFDIM) {
299 for (b = 0; b < CHILDDIM; b+=LEAFDIM) {
300 const Coord childijk = ijk + offset;
301 if (op.test(childijk, val)) {
302 mVoxelTopology.touchLeaf(childijk);
308 for (
auto iter = node.cbeginValueOn(); iter; ++iter) {
309 const Coord& ijk = iter.getCoord();
312 step(ijk, 0, 1, *iter);
313 step(ijk, 0, 2, *iter);
314 step(ijk, 1, 2, *iter);
318 template<
typename NodeT>
319 void operator()(
const NodeT& node)
const
321 using ChildT =
typename NodeT::ChildNodeType;
322 static constexpr
Int32 CHILDDIM =
Int32(ChildT::DIM);
323 static constexpr
Int32 LEAFDIM =
Int32(LeafT::DIM);
331 std::vector<Coord>& coords)
334 Int32& a = offset[axis1];
335 Int32& b = offset[axis2];
336 for (a = 0; a < CHILDDIM; a+=LEAFDIM) {
337 for (b = 0; b < CHILDDIM; b+=LEAFDIM) {
338 const Coord childijk = ijk + offset;
339 if (op.test(childijk, val)) {
340 coords.emplace_back(childijk);
345 offset.reset(CHILDDIM-1);
346 for (a = 0; a < CHILDDIM; a+=LEAFDIM) {
347 for (b = 0; b < CHILDDIM; b+=LEAFDIM) {
348 const Coord childijk = ijk + offset;
349 if (op.test(childijk, val)) {
350 coords.emplace_back(childijk);
376 if (CHILDDIM == LEAFDIM) {
384 std::vector<char> flags(NodeT::NUM_VALUES,
char(0));
385 tbb::parallel_for(tbb::blocked_range<size_t>(0, NodeT::NUM_VALUES),
386 [&](
const tbb::blocked_range<size_t>& range) {
387 const Tester op(mTree, mNeighbours);
388 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
389 if (node.isValueMaskOn(
Index(n))) {
391 const Coord ijk = node.offsetToGlobalCoord(
Index(n));
392 flags[n] = op.test(ijk, node.getValue(ijk));
399 for (
auto iter = flags.begin(); iter != flags.end(); ++iter, ++idx) {
400 if (*iter) mVoxelTopology.touchLeaf(node.offsetToGlobalCoord(idx));
410 tbb::concurrent_vector<Coord> nodes;
411 tbb::parallel_for(tbb::blocked_range<size_t>(0, NodeT::NUM_VALUES),
412 [&](
const tbb::blocked_range<size_t>& range)
414 const Tester op(mTree, mNeighbours);
415 std::vector<Coord> coords;
417 for (
size_t n = range.begin(), N = range.end(); n < N; ++n) {
418 if (!node.isValueMaskOn(
Index(n)))
continue;
420 const Coord ijk = node.offsetToGlobalCoord(
Index(n));
421 const auto& val = node.getValue(ijk);
424 step(op, ijk, 0, 1, val, coords);
425 step(op, ijk, 0, 2, val, coords);
426 step(op, ijk, 1, 2, val, coords);
429 if (!coords.empty()) {
430 std::copy(coords.begin(), coords.end(),
431 nodes.grow_by(coords.size()));
437 for (
const auto& coord : nodes) {
438 mVoxelTopology.touchLeaf(coord);
446 Tester(
const TreeT& tree,
const size_t NN)
447 : mAcc(tree), mNeighbours(NN) {}
449 inline bool test(
const Coord& ijk,
450 const typename TreeT::ValueType& val)
const
452 static constexpr
Int32 LEAFDIM =
Int32(LeafT::DIM);
454 for (
size_t i = 0; i < mNeighbours; ++i, ++NN) {
455 Coord neighbour(*NN);
456 neighbour.x() *= LEAFDIM;
457 neighbour.y() *= LEAFDIM;
458 neighbour.z() *= LEAFDIM;
461 if (mAcc.getValue(neighbour) != val ||
462 mAcc.probeConstLeaf(neighbour)) {
469 const tree::ValueAccessor<const TreeT> mAcc;
470 const size_t mNeighbours;
475 MaskT& mVoxelTopology;
476 const size_t mNeighbours;
480 MaskT mVoxelTopology;
481 std::unique_ptr<NodeManagerT> mManager;
482 const size_t mGrainSize;
487 template<
typename T>
static inline void accum(T& sum, T addend) { sum += addend; }
489 inline void accum(
bool& sum,
bool addend) { sum = sum || addend; }
496 template<
typename Gr
idT,
typename MaskT,
typename InterruptT>
497 template<
size_t Axis>
498 inline typename GridT::ValueType
499 Filter<GridT, MaskT, InterruptT>::Avg<Axis>::operator()(
Coord xyz)
501 ValueType sum = zeroVal<ValueType>();
505 ValueType value =
static_cast<ValueType
>(sum * frac);
514 template<
typename Gr
idT,
typename MaskT,
typename InterruptT>
518 if (iterations <= 0)
return;
521 const bool serial = mGrainSize == 0;
523 if (mInterrupter) mInterrupter->start(
"Applying mean filter");
525 std::unique_ptr<filter_internal::Voxelizer<TreeType>> voxelizer;
526 if (this->getProcessTiles()) {
529 const bool allNeighbours = iterations > 1;
533 (mGrid->tree(), allNeighbours, mGrainSize));
534 if (!voxelizer->run(w)) voxelizer.reset();
541 for (
int i=0; i<iterations && !this->
wasInterrupted(); ++i, dist+=w) {
542 if (i > 0 && voxelizer) {
545 const int remain = dist - iter * int(TreeType::LeafNodeType::DIM);
547 const int searches = voxelizer->run(remain);
548 if (searches == 0) voxelizer.reset();
549 else leafs.rebuild(serial);
554 mTask = std::bind(&Filter::doBoxX, std::placeholders::_1, std::placeholders::_2, w);
558 mTask = std::bind(&Filter::doBoxZ, std::placeholders::_1, std::placeholders::_2, w);
560 mTask = std::bind(&Filter::doBoxY, std::placeholders::_1, std::placeholders::_2, w);
564 if (mInterrupter) mInterrupter->end();
568 template<
typename Gr
idT,
typename MaskT,
typename InterruptT>
572 if (iterations <= 0)
return;
575 const bool serial = mGrainSize == 0;
577 if (mInterrupter) mInterrupter->start(
"Applying Gaussian filter");
579 std::unique_ptr<filter_internal::Voxelizer<TreeType>> voxelizer;
580 if (this->getProcessTiles()) {
583 const bool allNeighbours = iterations > 1;
588 (mGrid->tree(), allNeighbours, mGrainSize));
589 if (!voxelizer->run(w*4)) voxelizer.reset();
596 for (
int i=0; i<iterations; ++i, dist+=(w*4)) {
597 if (i > 0 && voxelizer) {
600 const int remain = dist - iter * int(TreeType::LeafNodeType::DIM);
602 const int searches = voxelizer->run(remain);
603 if (searches == 0) voxelizer.reset();
604 else leafs.rebuild(serial);
610 mTask = std::bind(&Filter::doBoxX, std::placeholders::_1, std::placeholders::_2, w);
614 mTask = std::bind(&Filter::doBoxZ, std::placeholders::_1, std::placeholders::_2, w);
616 mTask = std::bind(&Filter::doBoxY, std::placeholders::_1, std::placeholders::_2, w);
621 if (mInterrupter) mInterrupter->end();
625 template<
typename Gr
idT,
typename MaskT,
typename InterruptT>
629 if (iterations <= 0)
return;
632 const bool serial = mGrainSize == 0;
634 if (mInterrupter) mInterrupter->start(
"Applying median filter");
636 std::unique_ptr<filter_internal::Voxelizer<TreeType>> voxelizer;
637 if (this->getProcessTiles()) {
641 (mGrid->tree(),
true, mGrainSize));
642 if (!voxelizer->run(w)) voxelizer.reset();
647 mTask = std::bind(&Filter::doMedian, std::placeholders::_1, std::placeholders::_2, w);
651 for (
int i=0; i<iterations && !this->
wasInterrupted(); ++i, dist+=w) {
652 if (i > 0 && voxelizer) {
655 const int remain = dist - iter * int(TreeType::LeafNodeType::DIM);
657 const int searches = voxelizer->run(remain);
658 if (searches == 0) voxelizer.reset();
659 else leafs.rebuild(serial);
667 if (mInterrupter) mInterrupter->end();
671 template<
typename Gr
idT,
typename MaskT,
typename InterruptT>
677 if (mInterrupter) mInterrupter->start(
"Applying offset");
679 if (this->getProcessTiles()) {
683 NodeManagerT manager(mGrid->tree());
686 manager.foreachBottomUp([&](
auto& node) {
688 AlphaMaskT alpha(*mGrid, *mMask, mMinMask, mMaxMask, mInvertMask);
690 for (
auto iter = node.beginValueOn(); iter; ++iter) {
691 if (!alpha(iter.getCoord(), a, b))
continue;
693 iter.modifyValue([&](
ValueType& v) { v += a*value; });
699 manager.foreachBottomUp([&](
auto& node) {
701 for (
auto iter = node.beginValueOn(); iter; ++iter) {
702 iter.modifyValue([&](
ValueType& v) { v += value; });
709 mTask = std::bind(&Filter::doOffset, std::placeholders::_1, std::placeholders::_2, value);
712 if (mInterrupter) mInterrupter->end();
721 template<
typename Gr
idT,
typename MaskT,
typename InterruptT>
726 tbb::parallel_for(leafs.leafRange(mGrainSize), *
this);
728 (*this)(leafs.leafRange());
730 leafs.swapLeafBuffer(1, mGrainSize==0);
735 template<
typename Gr
idT,
typename MaskT,
typename InterruptT>
736 template <
typename AvgT>
738 Filter<GridT, MaskT, InterruptT>::doBox(
const RangeType& range,
Int32 w)
743 typename AlphaMaskT::FloatType a, b;
744 AlphaMaskT alpha(*mGrid, *mMask, mMinMask, mMaxMask, mInvertMask);
745 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
746 BufferT& buffer = leafIter.buffer(1);
747 for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
748 const Coord xyz = iter.getCoord();
749 if (alpha(xyz, a, b)) {
751 const ValueType value(b*(*iter) + a*avg(xyz));
753 buffer.setValue(iter.pos(), value);
758 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
759 BufferT& buffer = leafIter.buffer(1);
760 for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
761 buffer.setValue(iter.pos(), avg(iter.getCoord()));
769 template<
typename Gr
idT,
typename MaskT,
typename InterruptT>
771 Filter<GridT, MaskT, InterruptT>::doMedian(
const RangeType& range,
int width)
774 typename math::DenseStencil<GridType> stencil(*mGrid, width);
776 typename AlphaMaskT::FloatType a, b;
777 AlphaMaskT alpha(*mGrid, *mMask, mMinMask, mMaxMask, mInvertMask);
778 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
779 BufferT& buffer = leafIter.buffer(1);
780 for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
781 if (alpha(iter.getCoord(), a, b)) {
782 stencil.moveTo(iter);
784 ValueType value(b*(*iter) + a*stencil.median());
786 buffer.setValue(iter.pos(), value);
791 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
792 BufferT& buffer = leafIter.buffer(1);
793 for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
794 stencil.moveTo(iter);
795 buffer.setValue(iter.pos(), stencil.median());
803 template<
typename Gr
idT,
typename MaskT,
typename InterruptT>
805 Filter<GridT, MaskT, InterruptT>::doOffset(
const RangeType& range, ValueType offset)
809 typename AlphaMaskT::FloatType a, b;
810 AlphaMaskT alpha(*mGrid, *mMask, mMinMask, mMaxMask, mInvertMask);
811 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
812 for (VoxelIterT iter = leafIter->beginValueOn(); iter; ++iter) {
813 if (alpha(iter.getCoord(), a, b)) {
815 ValueType value(*iter + a*offset);
817 iter.setValue(value);
822 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
823 for (VoxelIterT iter = leafIter->beginValueOn(); iter; ++iter) {
824 iter.setValue(*iter + offset);
831 template<
typename Gr
idT,
typename MaskT,
typename InterruptT>
836 tbb::task::self().cancel_group_execution();
847 #endif // OPENVDB_TOOLS_FILTER_HAS_BEEN_INCLUDED