9 #ifndef OPENVDB_POINTS_RASTERIZE_SDF_IMPL_HAS_BEEN_INCLUDED 10 #define OPENVDB_POINTS_RASTERIZE_SDF_IMPL_HAS_BEEN_INCLUDED 19 namespace rasterize_sdf_internal
23 template <
typename ValueT>
26 FixedRadius(
const ValueT ris) : mR(ris) {}
27 inline void reset(
const PointDataTree::LeafNodeType&)
const {}
28 inline const FixedRadius& eval(
const Index)
const {
return *
this; }
29 inline ValueT
get()
const {
return mR; }
38 template <
typename ValueT>
39 struct FixedBandRadius :
public FixedRadius<ValueT>
41 FixedBandRadius(
const ValueT ris,
const ValueT hb)
42 : FixedRadius<ValueT>(ris)
43 , mMinSearchIS(math::
Max(ValueT(0.0), ris - hb))
44 , mMaxSearchIS(ris + hb)
45 , mMinSearchSqIS(mMinSearchIS*mMinSearchIS)
46 , mMaxSearchSqIS(mMaxSearchIS*mMaxSearchIS) {}
48 inline void reset(
const PointDataTree::LeafNodeType&)
const {}
49 inline const FixedBandRadius& eval(
const Index)
const {
return *
this; }
50 inline ValueT
min()
const {
return mMinSearchIS; }
51 inline ValueT minSq()
const {
return mMinSearchSqIS; }
52 inline ValueT
max()
const {
return mMaxSearchIS; }
53 inline ValueT maxSq()
const {
return mMaxSearchSqIS; }
55 const ValueT mMinSearchIS, mMaxSearchIS;
56 const ValueT mMinSearchSqIS, mMaxSearchSqIS;
60 template <
typename ValueT,
typename CodecT = UnknownCodec>
63 using RadiusHandleT = AttributeHandle<ValueT, CodecT>;
64 VaryingRadius(
const size_t ridx,
const ValueT
scale = 1.0)
65 : mRIdx(ridx), mRHandle(), mScale(
scale) {}
66 VaryingRadius(
const VaryingRadius& other)
67 : mRIdx(other.mRIdx), mRHandle(), mScale(other.mScale) {}
69 inline void reset(
const PointDataTree::LeafNodeType& leaf)
71 mRHandle.reset(
new RadiusHandleT(leaf.constAttributeArray(mRIdx)));
75 inline const FixedRadius<ValueT> eval(
const Index
id)
const 78 return FixedRadius<ValueT>(mRHandle->get(
id) * mScale);
83 typename RadiusHandleT::UniquePtr mRHandle;
88 template <
typename ValueT,
typename CodecT = UnknownCodec>
89 struct VaryingBandRadius :
public VaryingRadius<ValueT, CodecT>
91 using BaseT = VaryingRadius<ValueT, CodecT>;
92 VaryingBandRadius(
const size_t ridx,
const ValueT halfband,
93 const ValueT
scale = 1.0)
94 : BaseT(ridx,
scale), mHalfBand(halfband) {}
96 inline const FixedBandRadius<ValueT> eval(
const Index
id)
const 98 const auto r = this->BaseT::eval(
id).get();
99 return FixedBandRadius<ValueT>(ValueT(r), mHalfBand);
103 const ValueT mHalfBand;
110 template <
typename SdfT,
111 typename PositionCodecT,
114 struct SignedDistanceFieldTransfer :
115 public TransformTransfer,
116 public std::conditional<CPG,
117 VolumeTransfer<typename SdfT::TreeType, Int64Tree>,
118 VolumeTransfer<typename SdfT::TreeType>>::type
120 using TreeT =
typename SdfT::TreeType;
121 using ValueT =
typename TreeT::ValueType;
122 static_assert(std::is_floating_point<ValueT>::value,
123 "Spherical transfer only supports floating point values.");
124 static_assert(!std::is_reference<RadiusType>::value && !std::is_pointer<RadiusType>::value,
125 "Templated radius type must not be a reference or pointer");
127 using VolumeTransferT =
typename std::conditional<CPG,
128 VolumeTransfer<TreeT, Int64Tree>,
129 VolumeTransfer<TreeT>>::type;
131 using PositionHandleT = AttributeHandle<Vec3f, PositionCodecT>;
134 inline Int32 range(
const Coord&,
size_t)
const {
return Int32(mMaxKernelWidth); }
136 inline bool startPointLeaf(
const PointDataTree::LeafNodeType& leaf)
138 mPosition.reset(
new PositionHandleT(leaf.constAttributeArray(mPIdx)));
141 if (CPG) mPLeafMask = (
Index64(mIds->find(&leaf)->second) << 32);
147 template <
bool EnableT = CPG>
148 SignedDistanceFieldTransfer(
const size_t pidx,
150 const RadiusType& rt,
151 const math::Transform& source,
154 const std::unordered_map<const PointDataTree::LeafNodeType*, Index>* ids,
155 typename std::enable_if<EnableT>::type* = 0)
156 : TransformTransfer(source, surface.transform())
157 , VolumeTransferT(&(surface.tree()), cpg)
160 , mMaxKernelWidth(width)
162 , mBackground(surface.background())
163 , mDx(surface.voxelSize()[0])
170 template <
bool EnableT = CPG>
171 SignedDistanceFieldTransfer(
const size_t pidx,
173 const RadiusType& rt,
174 const math::Transform& source,
177 const std::unordered_map<const PointDataTree::LeafNodeType*, Index>*,
178 typename std::enable_if<!EnableT>::type* = 0)
179 : TransformTransfer(source, surface.transform())
180 , VolumeTransferT(surface.tree())
183 , mMaxKernelWidth(width)
185 , mBackground(surface.background())
186 , mDx(surface.voxelSize()[0])
190 SignedDistanceFieldTransfer(
const SignedDistanceFieldTransfer& other)
191 : TransformTransfer(other)
192 , VolumeTransferT(other)
195 , mMaxKernelWidth(other.mMaxKernelWidth)
196 , mRadius(other.mRadius)
197 , mBackground(other.mBackground)
204 typename PositionHandleT::UniquePtr mPosition;
205 const size_t mMaxKernelWidth;
207 const ValueT mBackground;
209 const std::unordered_map<const PointDataTree::LeafNodeType*, Index>* mIds;
215 template <
typename SdfT,
216 typename PositionCodecT,
219 struct SphericalTransfer :
220 public SignedDistanceFieldTransfer<SdfT, PositionCodecT, RadiusType, CPG>
222 using BaseT = SignedDistanceFieldTransfer<SdfT, PositionCodecT, RadiusType, CPG>;
223 using typename BaseT::TreeT;
224 using typename BaseT::ValueT;
225 static const Index DIM = TreeT::LeafNodeType::DIM;
226 static const Index LOG2DIM = TreeT::LeafNodeType::LOG2DIM;
228 using RealT = double;
230 SphericalTransfer(
const size_t pidx,
232 const RadiusType& rt,
233 const math::Transform& source,
236 const std::unordered_map<const PointDataTree::LeafNodeType*, Index>* ids =
nullptr)
237 : BaseT(pidx, width, rt, source, surface, cpg, ids) {}
243 inline void rasterizePoint(
const Coord& ijk,
245 const CoordBBox& bounds)
248 #if defined(__GNUC__) && !defined(__clang__) 249 #pragma GCC diagnostic push 250 #pragma GCC diagnostic ignored "-Wmaybe-uninitialized" 252 Vec3d P = ijk.asVec3d() +
Vec3d(this->mPosition->get(
id));
253 #if defined(__GNUC__) && !defined(__clang__) 254 #pragma GCC diagnostic pop 256 P = this->transformSourceToTarget(P);
258 const auto& r = this->mRadius.eval(
id);
259 const RealT
max = r.max();
261 CoordBBox intersectBox(Coord::floor(P - max), Coord::ceil(P + max));
262 intersectBox.intersect(bounds);
263 if (intersectBox.empty())
return;
265 auto*
const data = this->
template buffer<0>();
266 auto*
const cpg = CPG ? this->
template buffer<CPG ? 1 : 0>() :
nullptr;
267 auto& mask = *(this->
template mask<0>());
277 const RealT min2 = r.minSq() == 0.0 ? -1.0 : r.minSq();
278 const RealT max2 = r.maxSq();
280 const Coord& a(intersectBox.min());
281 const Coord& b(intersectBox.max());
282 for (Coord c = a; c.x() <= b.x(); ++c.x()) {
283 const RealT x2 =
static_cast<RealT
>(
math::Pow2(c.x() - P[0]));
284 const Index i = ((c.x() & (DIM-1u)) << 2*LOG2DIM);
285 for (c.y() = a.y(); c.y() <= b.y(); ++c.y()) {
286 const RealT x2y2 =
static_cast<RealT
>(x2 +
math::Pow2(c.y() - P[1]));
287 const Index ij = i + ((c.y() & (DIM-1u)) << LOG2DIM);
288 for (c.z() = a.z(); c.z() <= b.z(); ++c.z()) {
289 const Index offset = ij + (c.z() & (DIM-1u));
290 if (!mask.isOn(offset))
continue;
292 const RealT x2y2z2 =
static_cast<RealT
>(x2y2 +
math::Pow2(c.z() - P[2]));
293 if (x2y2z2 >= max2)
continue;
294 if (x2y2z2 <= min2) {
295 data[offset] = -(this->mBackground);
300 const ValueT d = ValueT(this->mDx * (
math::Sqrt(x2y2z2) - r.get()));
301 ValueT& v = data[offset];
305 if (CPG) cpg[offset] =
Int64(this->mPLeafMask |
Index64(
id));
325 inline bool endPointLeaf(
const PointDataTree::LeafNodeType&)
328 return !(this->
template mask<0>()->isOff());
331 inline bool finalize(
const Coord&,
const size_t)
336 auto& mask = *(this->
template mask<0>());
337 auto*
const data = this->
template buffer<0>();
338 for (
auto iter = mask.beginOn(); iter; ++iter) {
339 if (data[iter.pos()] == this->mBackground) mask.setOff(iter.pos());
342 if (CPG) *(this->
template mask<CPG ? 1 : 0>()) = mask;
349 template <
typename SdfT,
350 typename PositionCodecT,
353 struct AveragePositionTransfer :
354 public SignedDistanceFieldTransfer<SdfT, PositionCodecT, RadiusType, CPG>
356 using BaseT = SignedDistanceFieldTransfer<SdfT, PositionCodecT, RadiusType, CPG>;
357 using typename BaseT::TreeT;
358 using typename BaseT::ValueT;
360 using VolumeTransferT =
typename std::conditional<CPG,
361 VolumeTransfer<typename SdfT::TreeType, Int64Tree>,
362 VolumeTransfer<typename SdfT::TreeType>>::type;
364 static const Index DIM = TreeT::LeafNodeType::DIM;
365 static const Index LOG2DIM = TreeT::LeafNodeType::LOG2DIM;
366 static const Index NUM_VALUES = TreeT::LeafNodeType::NUM_VALUES;
368 using RealT = double;
372 template <
typename S>
inline void addP(
const math::Vec3<S>& v) { P += v; }
373 template <
typename S>
inline void addR(
const S r) { R += r; }
374 template <
typename S>
inline void multR(
const S w) { R *= w; }
375 template <
typename S>
inline void multP(
const S w) { P *= w; }
376 inline RealT length()
const {
return P.length() - R; }
377 math::Vec3<RealT> P = math::Vec3<RealT>(0.0);
381 AveragePositionTransfer(
const size_t pidx,
383 const RadiusType& rt,
385 const math::Transform& source,
388 const std::unordered_map<const PointDataTree::LeafNodeType*, Index>* ids =
nullptr)
389 : BaseT(pidx, width, rt, source, surface, cpg, ids)
390 , mMaxSearchIS(search)
391 , mMaxSearchSqIS(search*search)
395 AveragePositionTransfer(
const AveragePositionTransfer& other)
397 , mMaxSearchIS(other.mMaxSearchIS)
398 , mMaxSearchSqIS(other.mMaxSearchSqIS)
402 inline void initialize(
const Coord& origin,
const size_t idx,
const CoordBBox& bounds)
406 mWeights.assign(NUM_VALUES, PosRadPair());
418 auto*
const data = this->
template buffer<0>();
419 const auto& mask = *(this->
template mask<0>());
420 for (
auto iter = mask.beginOn(); iter; ++iter) {
421 data[iter.pos()] = ValueT(0);
425 inline void rasterizePoint(
const Coord& ijk,
427 const CoordBBox& bounds)
429 #if defined(__GNUC__) && !defined(__clang__) 430 #pragma GCC diagnostic push 431 #pragma GCC diagnostic ignored "-Wmaybe-uninitialized" 433 const Vec3d PWS = this->sourceTransform().indexToWorld(ijk.asVec3d() +
Vec3d(this->mPosition->get(
id)));
434 #if defined(__GNUC__) && !defined(__clang__) 435 #pragma GCC diagnostic pop 437 const Vec3d P = this->targetTransform().worldToIndex(PWS);
439 CoordBBox intersectBox(Coord::floor(P - mMaxSearchIS),
440 Coord::ceil(P + mMaxSearchIS));
441 intersectBox.intersect(bounds);
442 if (intersectBox.empty())
return;
444 auto*
const data = this->
template buffer<0>();
445 auto*
const cpg = CPG ? this->
template buffer<CPG ? 1 : 0>() :
nullptr;
446 const auto& mask = *(this->
template mask<0>());
449 const auto& r = this->mRadius.eval(
id);
450 const RealT rad = r.get();
451 const RealT invsq = 1.0 / mMaxSearchSqIS;
453 const Coord& a(intersectBox.min());
454 const Coord& b(intersectBox.max());
455 for (Coord c = a; c.x() <= b.x(); ++c.x()) {
456 const RealT x2 =
static_cast<RealT
>(
math::Pow2(c.x() - P[0]));
457 const Index i = ((c.x() & (DIM-1u)) << 2*LOG2DIM);
458 for (c.y() = a.y(); c.y() <= b.y(); ++c.y()) {
459 const RealT x2y2 =
static_cast<RealT
>(x2 +
math::Pow2(c.y() - P[1]));
460 const Index ij = i + ((c.y() & (DIM-1u)) << LOG2DIM);
461 for (c.z() = a.z(); c.z() <= b.z(); ++c.z()) {
462 RealT x2y2z2 =
static_cast<RealT
>(x2y2 +
math::Pow2(c.z() - P[2]));
463 if (x2y2z2 >= mMaxSearchSqIS)
continue;
464 const Index offset = ij + (c.z() & (DIM-1u));
465 if (!mask.isOn(offset))
continue;
476 const float dist =
static_cast<float>(
math::Sqrt(x2y2z2) - r.get());
477 auto& d = mDist[offset];
494 data[offset] += x2y2z2;
496 auto& wt = mWeights[offset];
497 wt.addP(PWS * x2y2z2);
498 wt.addR(rad * x2y2z2);
504 inline bool endPointLeaf(
const PointDataTree::LeafNodeType&) {
return true; }
506 inline bool finalize(
const Coord& origin,
const size_t)
508 auto& mask = *(this->
template mask<0>());
509 auto*
const data = this->
template buffer<0>();
511 for (
auto iter = mask.beginOn(); iter; ++iter) {
512 const Index idx = iter.pos();
518 w = this->mBackground;
521 const Coord ijk = origin + TreeT::LeafNodeType::offsetToLocalCoord(idx);
522 const Vec3d ws = this->targetTransform().indexToWorld(ijk);
523 const RealT wi = RealT(1.0) / RealT(w);
524 auto& wt = mWeights[idx];
526 wt.multR(wi * this->mDx);
528 w =
static_cast<typename SdfT::ValueType
>(wt.length());
530 if (std::fabs(w) >= this->mBackground) {
531 w = std::copysign(this->mBackground, w);
538 if (CPG) *(this->
template mask<CPG ? 1 : 0>()) = mask;
543 const RealT mMaxSearchIS, mMaxSearchSqIS;
544 std::vector<PosRadPair> mWeights;
545 std::vector<float> mDist;
549 template <
typename MaskTreeT =
MaskTree,
550 typename InterrupterT = util::NullInterrupter>
553 inline void join(SurfaceMaskOp& other)
555 if (mMask->leafCount() > other.mMask->leafCount()) {
556 mMask->topologyUnion(*other.mMask);
560 other.mMask->topologyUnion(*mMask);
561 mMask.reset(other.mMask.release());
564 if (mMaskOff->leafCount() > other.mMaskOff->leafCount()) {
565 mMaskOff->topologyUnion(*other.mMaskOff);
566 other.mMaskOff.reset();
569 other.mMaskOff->topologyUnion(*mMaskOff);
570 mMaskOff.reset(other.mMaskOff.release());
574 std::unique_ptr<MaskTreeT> mask() {
return std::move(mMask); }
575 std::unique_ptr<MaskTreeT> maskoff() {
return std::move(mMaskOff); }
578 SurfaceMaskOp(
const math::Transform& points,
579 const math::Transform& surface,
580 InterrupterT* interrupter =
nullptr)
581 : mMask(new MaskTreeT)
582 , mMaskOff(new MaskTreeT)
583 , mPointsTransform(points)
584 , mSurfaceTransform(surface)
585 , mInterrupter(interrupter)
586 , mVoxelOffset(static_cast<int>(math::
RoundUp(points.voxelSize()[0] /
587 surface.voxelSize()[0]))) {}
589 SurfaceMaskOp(
const SurfaceMaskOp& other)
590 : mMask(new MaskTreeT)
591 , mMaskOff(new MaskTreeT)
592 , mPointsTransform(other.mPointsTransform)
593 , mSurfaceTransform(other.mSurfaceTransform)
594 , mInterrupter(other.mInterrupter)
595 , mVoxelOffset(other.mVoxelOffset) {}
603 template <
typename LeafT>
604 inline void fill(
const LeafT& leaf,
const int dist,
const bool active)
606 auto doFill = [&](CoordBBox b) {
607 b = mSurfaceTransform.worldToIndexCellCentered(
608 mPointsTransform.indexToWorld(b));
610 if (active) mMask->sparseFill(b,
true,
true);
611 else mMaskOff->sparseFill(b,
true,
true);
614 const auto& mask = leaf.getValueMask();
616 doFill(leaf.getNodeBoundingBox());
620 for (
auto iter = mask.beginOn(); iter; ++iter) {
621 bounds.expand(leaf.offsetToLocalCoord(iter.pos()));
623 if (!bounds.empty()) {
624 bounds.translate(leaf.origin());
630 inline bool interrupted()
633 thread::cancelGroupExecution();
640 std::unique_ptr<MaskTreeT> mMask;
641 std::unique_ptr<MaskTreeT> mMaskOff;
642 const math::Transform& mPointsTransform;
643 const math::Transform& mSurfaceTransform;
644 InterrupterT* mInterrupter;
651 const int mVoxelOffset;
655 template <
typename MaskTreeT,
656 typename InterrupterT = util::NullInterrupter>
657 struct FixedSurfaceMaskOp
658 :
public SurfaceMaskOp<MaskTreeT, InterrupterT>
660 using BaseT = SurfaceMaskOp<MaskTreeT, InterrupterT>;
662 tree::LeafManager<const points::PointDataTree>;
664 FixedSurfaceMaskOp(
const math::Transform& points,
665 const math::Transform& surface,
666 const double minBandRadius,
667 const double maxBandRadius,
668 InterrupterT* interrupter =
nullptr)
669 : BaseT(points, surface, interrupter)
675 const Real min = ((2.0 * minBandRadius) / std::sqrt(3.0)) / 2.0;
676 const int d =
static_cast<int>(
math::RoundDown(min)) - this->mVoxelOffset;
679 mMax += this->mVoxelOffset;
682 FixedSurfaceMaskOp(
const FixedSurfaceMaskOp& other,
tbb::split)
683 : BaseT(other), mMin(other.mMin), mMax(other.mMax) {}
685 void operator()(
const typename LeafManagerT::LeafRange& range)
687 if (this->interrupted())
return;
688 for (
auto leafIter = range.begin(); leafIter; ++leafIter) {
689 this->fill(*leafIter, mMax,
true);
691 if (mMin <= 0)
return;
692 for (
auto leafIter = range.begin(); leafIter; ++leafIter) {
693 this->fill(*leafIter, mMin,
false);
702 template <
typename RadiusTreeT,
704 typename InterrupterT = util::NullInterrupter>
705 struct VariableSurfaceMaskOp
706 :
public SurfaceMaskOp<MaskTreeT, InterrupterT>
708 using BaseT = SurfaceMaskOp<MaskTreeT, InterrupterT>;
710 tree::LeafManager<const points::PointDataTree>;
712 VariableSurfaceMaskOp(
const math::Transform& pointsTransform,
713 const math::Transform& surfaceTransform,
714 const RadiusTreeT& min,
715 const RadiusTreeT& max,
718 InterrupterT* interrupter =
nullptr)
719 : BaseT(pointsTransform, surfaceTransform, interrupter)
720 , mMin(min), mMax(max), mScale(scale), mHalfband(halfband) {}
722 VariableSurfaceMaskOp(
const VariableSurfaceMaskOp& other,
tbb::split)
723 : BaseT(other), mMin(other.mMin), mMax(other.mMax)
724 , mScale(other.mScale), mHalfband(other.mHalfband) {}
726 void operator()(
const typename LeafManagerT::LeafRange& range)
728 if (this->interrupted())
return;
729 const tree::ValueAccessor<const RadiusTreeT> maxacc(mMax);
730 for (
auto leafIter = range.begin(); leafIter; ++leafIter) {
731 const int max = this->maxDist(maxacc.getValue(leafIter->origin()));
732 this->fill(*leafIter, max,
true);
734 const tree::ValueAccessor<const RadiusTreeT> minacc(mMin);
735 for (
auto leafIter = range.begin(); leafIter; ++leafIter) {
736 const int min = this->minDist(minacc.getValue(leafIter->origin()));
737 if (min <= 0)
continue;
738 this->fill(*leafIter, min,
false);
743 inline int maxDist(
const typename RadiusTreeT::ValueType& max)
const 746 const Real v = (
Real(max) * mScale) + mHalfband;
748 d += this->mVoxelOffset;
752 inline int minDist(
const typename RadiusTreeT::ValueType& min)
const 761 v = ((2.0 * v) / std::sqrt(3.0)) / 2.0;
762 return static_cast<int>(v);
766 const RadiusTreeT& mMin;
767 const RadiusTreeT& mMax;
768 const Real mScale, mHalfband;
771 template <
typename SdfT,
typename InterrupterT,
typename Po
intDataGr
idT>
772 inline typename SdfT::Ptr
773 initFixedSdf(
const PointDataGridT& points,
774 math::Transform::Ptr transform,
775 const typename SdfT::ValueType bg,
776 const double minBandRadius,
777 const double maxBandRadius,
778 InterrupterT* interrupter)
780 using LeafManagerT = tree::LeafManager<const typename PointDataGridT::TreeType>;
781 using MaskTreeT =
typename SdfT::TreeType::template ValueConverter<ValueMask>::Type;
783 typename SdfT::Ptr surface = SdfT::create(bg);
784 surface->setTransform(transform);
787 FixedSurfaceMaskOp<MaskTreeT, InterrupterT>
op(points.transform(),
788 *transform, minBandRadius, maxBandRadius, interrupter);
790 if (interrupter) interrupter->start(
"Generating uniform surface topology");
792 LeafManagerT leafManager(points.tree());
793 tbb::parallel_reduce(leafManager.leafRange(),
op);
794 auto mask =
op.mask();
796 if (minBandRadius > 0.0) {
797 auto maskoff =
op.maskoff();
798 mask->topologyDifference(*maskoff);
801 surface->tree().topologyUnion(*mask);
804 tree::ValueAccessor<const MaskTreeT> acc(*maskoff);
805 auto setOffOp = [acc](
auto& iter) {
806 if (acc.isValueOn(iter.getCoord())) {
807 iter.modifyValue([](
auto& v) { v = -v; });
815 surface->tree().topologyUnion(*mask);
819 surface->tree().voxelizeActiveTiles();
820 if (interrupter) interrupter->end();
824 template <
typename SdfT,
825 typename InterrupterT,
826 typename PointDataGridT,
827 typename RadiusTreeT>
828 inline typename SdfT::Ptr
829 initVariableSdf(
const PointDataGridT& points,
830 math::Transform::Ptr transform,
831 const typename SdfT::ValueType bg,
832 const RadiusTreeT& min,
833 const RadiusTreeT& max,
836 InterrupterT* interrupter)
838 using LeafManagerT = tree::LeafManager<const typename PointDataGridT::TreeType>;
839 using MaskTreeT =
typename SdfT::TreeType::template ValueConverter<ValueMask>::Type;
841 typename SdfT::Ptr surface = SdfT::create(bg);
842 surface->setTransform(transform);
845 VariableSurfaceMaskOp<RadiusTreeT, MaskTreeT, InterrupterT>
846 op(points.transform(), *transform,
min,
max,
scale, halfband, interrupter);
848 if (interrupter) interrupter->start(
"Generating variable surface topology");
850 LeafManagerT leafManager(points.tree());
851 tbb::parallel_reduce(leafManager.leafRange(),
op);
852 auto mask =
op.mask();
853 auto maskoff =
op.maskoff();
855 if (!maskoff->empty()) {
856 mask->topologyDifference(*maskoff);
859 surface->tree().topologyUnion(*mask);
861 tree::ValueAccessor<const MaskTreeT> acc(*maskoff);
862 auto setOffOp = [acc](
auto& iter) {
863 if (acc.isValueOn(iter.getCoord())) {
864 iter.modifyValue([](
auto& v) { v = -v; });
871 surface->tree().topologyUnion(*mask);
876 surface->tree().voxelizeActiveTiles();
877 if (interrupter) interrupter->end();
881 template <
typename PointDataTreeT,
884 transferAttributes(
const tree::LeafManager<const PointDataTreeT>& manager,
885 const std::vector<std::string>& attributes,
887 const math::Transform::Ptr transform)
897 const tree::LeafManager<const Int64Tree> cpmanager(cpg);
899 auto transfer = [&](
auto& tree,
const size_t attrIdx) {
900 using TreeType =
typename std::decay<decltype(tree)>::type;
901 using HandleT = AttributeHandle<typename TreeType::ValueType>;
904 tree.topologyUnion(cpg);
905 tree::LeafManager<TreeType> lm(tree);
908 lm.foreach([&manager, &cpmanager, attrIdx]
909 (
auto& leaf,
const size_t idx) {
910 auto voxel = leaf.beginValueOn();
913 auto* data = leaf.buffer().data();
914 const Int64* ids = cpmanager.leaf(idx).buffer().data();
916 typename HandleT::UniquePtr handle(
917 new HandleT(manager.leaf(prev).constAttributeArray(attrIdx)));
919 for (; voxel; ++voxel) {
924 handle.reset(
new HandleT(manager.leaf(lfid).constAttributeArray(attrIdx)));
927 data[voxel.pos()] = handle->get(ptid);
933 grids.reserve(attributes.size());
934 const auto& attrSet = manager.leaf(0).attributeSet();
935 tbb::task_group tasks;
937 for (
const auto& name : attributes) {
938 const size_t attrIdx = attrSet.find(name);
939 if (attrIdx == points::AttributeSet::INVALID_POS)
continue;
940 if (attrSet.get(attrIdx)->stride() != 1) {
941 OPENVDB_THROW(RuntimeError,
"Transfer of attribute " + name +
942 " not supported since it is strided");
945 const std::string& type = attrSet.descriptor().valueType(attrIdx);
948 using ValueType =
typename std::remove_const<typename std::decay<decltype(v)>::type>::type;
949 using TreeT =
typename PointDataTreeT::template ValueConverter<ValueType>::Type;
950 if (!grid && typeNameAsString<ValueType>() == type) {
953 typed->setName(name);
954 typed->setTransform(transform);
955 tasks.run([typed, attrIdx, transfer] { transfer(typed->tree(), attrIdx); });
956 grids.emplace_back(grid);
961 OPENVDB_THROW(RuntimeError,
"No support for attribute type " + type +
962 " built during closest point surface transfer");
970 template <
typename SdfT,
971 template <
typename,
typename,
typename,
bool>
class TransferInterfaceT,
973 typename InterrupterT,
974 typename PointDataGridT,
978 rasterizeSurface(
const PointDataGridT& points,
979 const std::vector<std::string>& attributes,
980 const FilterT& filter,
982 InterrupterT* interrupter,
985 using RadRefT =
typename std::tuple_element<1, std::tuple<Args...>>::type;
986 using RadT =
typename std::remove_reference<RadRefT>::type;
989 auto leaf = points.constTree().cbeginLeaf();
990 if (!leaf)
return grids;
992 const size_t pidx = leaf->attributeSet().find(
"P");
993 if (pidx == AttributeSet::INVALID_POS) {
994 OPENVDB_THROW(RuntimeError,
"Failed to find position attribute");
999 const NamePair& ptype = leaf->attributeSet().descriptor().type(pidx);
1001 if (attributes.empty()) {
1002 if (ptype.second == NullCodec::name()) {
1003 using TransferT = TransferInterfaceT<SdfT, NullCodec, RadT, false>;
1004 TransferT transfer(pidx, args...);
1005 rasterize<PointDataGridT, TransferT, FilterT, InterrupterT>
1006 (points, transfer, filter, interrupter);
1009 using TransferT = TransferInterfaceT<SdfT, UnknownCodec, RadT, false>;
1010 TransferT transfer(pidx, args...);
1011 rasterize<PointDataGridT, TransferT, FilterT, InterrupterT>
1012 (points, transfer, filter, interrupter);
1017 cpg.topologyUnion(surface.tree());
1018 tree::LeafManager<const PointDataTree> manager(points.tree());
1021 std::unordered_map<const PointDataTree::LeafNodeType*, Index> ids;
1022 manager.foreach([&](
auto& leaf,
size_t idx) { ids[&leaf] =
Index(idx); },
false);
1024 if (ptype.second == NullCodec::name()) {
1025 using TransferT = TransferInterfaceT<SdfT, NullCodec, RadT, true>;
1026 TransferT transfer(pidx, args..., &cpg, &ids);
1027 rasterize<PointDataGridT, TransferT, FilterT, InterrupterT>
1028 (points, transfer, filter, interrupter);
1031 using TransferT = TransferInterfaceT<SdfT, UnknownCodec, RadT, true>;
1032 TransferT transfer(pidx, args..., &cpg, &ids);
1033 rasterize<PointDataGridT, TransferT, FilterT, InterrupterT>
1034 (points, transfer, filter, interrupter);
1040 grids = transferAttributes
1041 <
typename PointDataGrid::TreeType, AttributeTypes>
1042 (manager, attributes, cpg, surface.transformPtr());
1056 template <
typename PointDataGridT,
1059 typename InterrupterT>
1063 const Real halfband,
1065 const FilterT& filter,
1066 InterrupterT* interrupter)
1069 rasterizeSpheres<PointDataGridT, TypeList<>, SdfT, FilterT, InterrupterT>
1070 (points, radius, {}, halfband, transform, filter, interrupter);
1071 return StaticPtrCast<SdfT>(grids.front());
1074 template <
typename PointDataGridT,
1078 typename InterrupterT>
1081 const std::string& radius,
1083 const Real halfband,
1085 const FilterT& filter,
1086 InterrupterT* interrupter)
1089 rasterizeSpheres<PointDataGridT, TypeList<>, RadiusT, SdfT, FilterT, InterrupterT>
1090 (points, radius, {},
scale, halfband, transform, filter, interrupter);
1091 return StaticPtrCast<SdfT>(grids.front());
1094 template <
typename PointDataGridT,
1098 typename InterrupterT>
1102 const std::vector<std::string>& attributes,
1103 const Real halfband,
1105 const FilterT& filter,
1106 InterrupterT* interrupter)
1108 using namespace rasterize_sdf_internal;
1110 if (!transform) transform = points.transform().copy();
1111 const Real vs = transform->voxelSize()[0];
1112 const float background =
static_cast<float>(vs * halfband);
1115 const Real radiusIndexSpace = radius / vs;
1116 const FixedBandRadius<Real> rad(radiusIndexSpace, halfband);
1117 const Real minBandRadius = rad.min();
1118 const Real maxBandRadius = rad.max();
1119 const size_t width =
static_cast<size_t>(
math::RoundUp(maxBandRadius));
1121 typename SdfT::Ptr surface =
1122 initFixedSdf<SdfT, InterrupterT>
1123 (points, transform, background, minBandRadius, maxBandRadius, interrupter);
1125 if (interrupter) interrupter->start(
"Rasterizing particles to level set using constant Spheres");
1128 rasterizeSurface<SdfT, SphericalTransfer, AttributeTypes, InterrupterT>
1129 (points, attributes, filter, *surface, interrupter,
1130 width, rad, points.transform(), *surface);
1133 if (interrupter) interrupter->end();
1134 grids.insert(grids.begin(), surface);
1138 template <
typename PointDataGridT,
1143 typename InterrupterT>
1146 const std::string& radius,
1147 const std::vector<std::string>& attributes,
1149 const Real halfband,
1151 const FilterT& filter,
1152 InterrupterT* interrupter)
1154 using namespace rasterize_sdf_internal;
1155 using PointDataTreeT =
typename PointDataGridT::TreeType;
1156 using RadTreeT =
typename PointDataTreeT::template ValueConverter<RadiusT>::Type;
1158 if (!transform) transform = points.transform().copy();
1159 const Real vs = transform->voxelSize()[0];
1160 const float background =
static_cast<float>(vs * halfband);
1163 typename RadTreeT::Ptr mintree(
new RadTreeT), maxtree(
new RadTreeT);
1164 points::evalMinMax<RadiusT, UnknownCodec>
1165 (points.tree(), radius,
min,
max, filter, mintree.get(), maxtree.get());
1168 const RadiusT indexSpaceScale = RadiusT(scale / vs);
1169 typename SdfT::Ptr surface =
1170 initVariableSdf<SdfT, InterrupterT>
1171 (points, transform, background, *mintree, *maxtree,
1172 indexSpaceScale, halfband, interrupter);
1176 auto leaf = points.constTree().cbeginLeaf();
1180 const size_t width =
static_cast<size_t> 1183 const size_t ridx = leaf->attributeSet().find(radius);
1184 if (ridx == AttributeSet::INVALID_POS) {
1188 VaryingBandRadius<RadiusT> rad(ridx, RadiusT(halfband), indexSpaceScale);
1189 if (interrupter) interrupter->start(
"Rasterizing particles to level set using variable Spheres");
1192 rasterizeSurface<SdfT, SphericalTransfer, AttributeTypes, InterrupterT>
1193 (points, attributes, filter, *surface, interrupter,
1194 width, rad, points.transform(), *surface);
1196 if (interrupter) interrupter->end();
1198 grids.insert(grids.begin(), surface);
1204 template <
typename PointDataGridT,
1207 typename InterrupterT>
1211 const Real searchRadius,
1212 const Real halfband,
1214 const FilterT& filter,
1215 InterrupterT* interrupter)
1218 rasterizeSmoothSpheres<PointDataGridT, TypeList<>, SdfT, FilterT, InterrupterT>
1219 (points, radius, searchRadius, {}, halfband, transform, filter, interrupter);
1220 return StaticPtrCast<SdfT>(grids.front());
1223 template <
typename PointDataGridT,
1227 typename InterrupterT>
1230 const std::string& radius,
1231 const Real radiusScale,
1232 const Real searchRadius,
1233 const Real halfband,
1235 const FilterT& filter,
1236 InterrupterT* interrupter)
1239 rasterizeSmoothSpheres<PointDataGridT, TypeList<>, RadiusT, SdfT, FilterT, InterrupterT>
1240 (points, radius, radiusScale, searchRadius, {}, halfband, transform, filter, interrupter);
1241 return StaticPtrCast<SdfT>(grids.front());
1244 template <
typename PointDataGridT,
1248 typename InterrupterT>
1252 const Real searchRadius,
1253 const std::vector<std::string>& attributes,
1254 const Real halfband,
1256 const FilterT& filter,
1257 InterrupterT* interrupter)
1259 using namespace rasterize_sdf_internal;
1261 if (!transform) transform = points.transform().copy();
1262 const Real vs = transform->voxelSize()[0];
1263 const float background =
static_cast<float>(vs * halfband);
1264 const Real indexSpaceSearch = searchRadius / vs;
1265 const FixedBandRadius<Real> bands(radius / vs, halfband);
1266 const Real max = bands.max();
1268 typename SdfT::Ptr surface =
1269 initFixedSdf<SdfT, InterrupterT>
1270 (points, transform, background, 0.0,
max, interrupter);
1272 auto leaf = points.constTree().cbeginLeaf();
1276 const size_t width =
static_cast<size_t>(
math::RoundUp(indexSpaceSearch));
1278 const FixedRadius<Real> rad(radius / vs);
1279 if (interrupter) interrupter->start(
"Rasterizing particles to level set using constant Zhu-Bridson");
1282 rasterizeSurface<SdfT, AveragePositionTransfer, AttributeTypes, InterrupterT>
1283 (points, attributes, filter, *surface, interrupter,
1284 width, rad, indexSpaceSearch, points.transform(), *surface);
1287 if (interrupter) interrupter->end();
1288 grids.insert(grids.begin(), surface);
1292 template <
typename PointDataGridT,
1297 typename InterrupterT>
1300 const std::string& radius,
1301 const Real radiusScale,
1302 const Real searchRadius,
1303 const std::vector<std::string>& attributes,
1304 const Real halfband,
1306 const FilterT& filter,
1307 InterrupterT* interrupter)
1309 using namespace rasterize_sdf_internal;
1310 using PointDataTreeT =
typename PointDataGridT::TreeType;
1311 using RadTreeT =
typename PointDataTreeT::template ValueConverter<RadiusT>::Type;
1313 if (!transform) transform = points.transform().copy();
1314 const Real vs = transform->voxelSize()[0];
1315 const float background =
static_cast<float>(vs * halfband);
1316 const Real indexSpaceSearch = searchRadius / vs;
1319 typename RadTreeT::Ptr mintree(
new RadTreeT), maxtree(
new RadTreeT);
1320 points::evalMinMax<RadiusT, UnknownCodec>
1321 (points.tree(), radius,
min,
max, filter, mintree.get(), maxtree.get());
1324 const RadiusT indexSpaceScale = RadiusT(radiusScale / vs);
1325 typename SdfT::Ptr surface =
1326 initVariableSdf<SdfT, InterrupterT>
1327 (points, transform, background, *mintree, *maxtree,
1328 indexSpaceScale, halfband, interrupter);
1332 auto leaf = points.constTree().cbeginLeaf();
1336 const size_t width =
static_cast<size_t>(
math::RoundUp(indexSpaceSearch));
1338 const size_t ridx = leaf->attributeSet().find(radius);
1339 if (ridx == AttributeSet::INVALID_POS) {
1343 VaryingRadius<RadiusT> rad(ridx, indexSpaceScale);
1344 if (interrupter) interrupter->start(
"Rasterizing particles to level set using variable Zhu-Bridson");
1347 rasterizeSurface<SdfT, AveragePositionTransfer, AttributeTypes, InterrupterT>
1348 (points, attributes, filter, *surface, interrupter,
1349 width, rad, indexSpaceSearch, points.transform(), *surface);
1352 if (interrupter) interrupter->end();
1353 grids.insert(grids.begin(), surface);
1361 #endif //OPENVDB_POINTS_RASTERIZE_SDF_IMPL_HAS_BEEN_INCLUDED tree::Tree4< int64_t, 5, 4, 3 >::Type Int64Tree
Definition: openvdb.h:57
Type Pow2(Type x)
Return x2.
Definition: Math.h:548
SdfT::Ptr rasterizeSpheres(const PointDataGridT &points, const Real radius, const Real halfband=LEVEL_SET_HALF_WIDTH, math::Transform::Ptr transform=nullptr, const FilterT &filter=NullFilter(), InterrupterT *interrupter=nullptr)
Narrow band sphere stamping with a uniform radius.
Definition: PointRasterizeSDFImpl.h:1061
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
float RoundDown(float x)
Return x rounded down to the nearest integer.
Definition: Math.h:803
uint64_t Index64
Definition: Types.h:53
static Ptr create()
Return a new grid with background value zero.
Definition: Grid.h:1343
SharedPtr< GridBase > Ptr
Definition: Grid.h:80
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:595
int32_t Int32
Definition: Types.h:56
Index32 Index
Definition: Types.h:54
Vec3< double > Vec3d
Definition: Vec3.h:665
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:49
std::vector< GridBase::Ptr > GridPtrVec
Definition: Grid.h:508
Type Pow3(Type x)
Return x3.
Definition: Math.h:552
void split(ContainerT &out, const std::string &in, const char delim)
Definition: Name.h:43
OPENVDB_IMPORT void initialize()
Global registration of native Grid, Transform, Metadata and Point attribute types. Also initializes blosc (if enabled).
#define OPENVDB_ASSERT(X)
Definition: Assert.h:41
NumericAttributeTypes::Append< Vec3AttributeTypes >::Append< Mat3AttributeTypes >::Append< Mat4AttributeTypes >::Append< QuatAttributeTypes >::Append< points::GroupAttributeArray >::Append< points::StringAttributeArray >::Append< points::TypedAttributeArray< bool >> AttributeTypes
The attribute array types which OpenVDB will register by default.
Definition: openvdb.h:186
SdfT::Ptr rasterizeSmoothSpheres(const PointDataGridT &points, const Real radius, const Real searchRadius, const Real halfband=LEVEL_SET_HALF_WIDTH, math::Transform::Ptr transform=nullptr, const FilterT &filter=NullFilter(), InterrupterT *interrupter=nullptr)
Smoothed point distribution based sphere stamping with a uniform radius.
Definition: PointRasterizeSDFImpl.h:1209
MatType scale(const Vec3< typename MatType::value_type > &s)
Return a matrix that scales by s.
Definition: Mat.h:615
void foreach(T &&t, const F &func, std::integer_sequence< size_t, Is... >)
Definition: PointTransfer.h:335
Definition: Exceptions.h:13
int64_t Int64
Definition: Types.h:57
std::pair< Name, Name > NamePair
Definition: AttributeArray.h:40
double Real
Definition: Types.h:60
__hostdev__ uint32_t hash(uint32_t x)
Definition: common.h:14
tree::Tree4< ValueMask, 5, 4, 3 >::Type MaskTree
Definition: openvdb.h:58
float RoundUp(float x)
Return x rounded up to the nearest integer.
Definition: Math.h:787
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:121
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:761
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:218
Definition: Exceptions.h:63