OpenVDB  12.0.0
PointRasterizeSDFImpl.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: Apache-2.0
3 //
4 /// @author Nick Avramoussis
5 ///
6 /// @file PointRasterizeSDFImpl.h
7 ///
8 
9 #ifndef OPENVDB_POINTS_RASTERIZE_SDF_IMPL_HAS_BEEN_INCLUDED
10 #define OPENVDB_POINTS_RASTERIZE_SDF_IMPL_HAS_BEEN_INCLUDED
11 
12 namespace openvdb {
14 namespace OPENVDB_VERSION_NAME {
15 namespace points {
16 
17 /// @cond OPENVDB_DOCS_INTERNAL
18 
19 namespace rasterize_sdf_internal
20 {
21 
22 /// @brief Define a fixed index space radius for point rasterization
23 template <typename ValueT>
24 struct FixedRadius
25 {
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; }
30 private:
31  const ValueT mR;
32 };
33 
34 /// @brief Define a fixed narrow band radius for point rasterization
35 /// @details Pass in an index space radius (relative to a PointDataGrid voxel
36 /// size) and the desired half band width of the target surface. The minimum
37 /// radius of influence is clamped to zero.
38 template <typename ValueT>
39 struct FixedBandRadius : public FixedRadius<ValueT>
40 {
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) {}
47 
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; }
54 private:
55  const ValueT mMinSearchIS, mMaxSearchIS;
56  const ValueT mMinSearchSqIS, mMaxSearchSqIS;
57 };
58 
59 /// @brief A varying per point radius with an optional scale
60 template <typename ValueT, typename CodecT = UnknownCodec>
61 struct VaryingRadius
62 {
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) {}
68 
69  inline void reset(const PointDataTree::LeafNodeType& leaf)
70  {
71  mRHandle.reset(new RadiusHandleT(leaf.constAttributeArray(mRIdx)));
72  }
73 
74  /// @brief Compute a fixed radius for a specific point
75  inline const FixedRadius<ValueT> eval(const Index id) const
76  {
77  OPENVDB_ASSERT(mRHandle);
78  return FixedRadius<ValueT>(mRHandle->get(id) * mScale);
79  }
80 
81 private:
82  const size_t mRIdx;
83  typename RadiusHandleT::UniquePtr mRHandle;
84  const ValueT mScale;
85 };
86 
87 /// @brief A varying per point narrow band radius with an optional scale
88 template <typename ValueT, typename CodecT = UnknownCodec>
89 struct VaryingBandRadius : public VaryingRadius<ValueT, CodecT>
90 {
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) {}
95 
96  inline const FixedBandRadius<ValueT> eval(const Index id) const
97  {
98  const auto r = this->BaseT::eval(id).get();
99  return FixedBandRadius<ValueT>(ValueT(r), mHalfBand);
100  }
101 
102 private:
103  const ValueT mHalfBand;
104 };
105 
106 /// @brief Base class for SDF transfers which consolidates member data and
107 /// some required interface methods.
108 /// @note Derives from TransformTransfer for automatic transformation support
109 /// and VolumeTransfer<...T> for automatic buffer setup
110 template <typename SdfT,
111  typename PositionCodecT,
112  typename RadiusType,
113  bool CPG>
114 struct SignedDistanceFieldTransfer :
115  public TransformTransfer,
116  public std::conditional<CPG,
117  VolumeTransfer<typename SdfT::TreeType, Int64Tree>,
118  VolumeTransfer<typename SdfT::TreeType>>::type
119 {
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");
126 
127  using VolumeTransferT = typename std::conditional<CPG,
128  VolumeTransfer<TreeT, Int64Tree>,
129  VolumeTransfer<TreeT>>::type;
130 
131  using PositionHandleT = AttributeHandle<Vec3f, PositionCodecT>;
132 
133  // typically the max radius of all points rounded up
134  inline Int32 range(const Coord&, size_t) const { return Int32(mMaxKernelWidth); }
135 
136  inline bool startPointLeaf(const PointDataTree::LeafNodeType& leaf)
137  {
138  mPosition.reset(new PositionHandleT(leaf.constAttributeArray(mPIdx)));
139  mRadius.reset(leaf);
140  // if CPG, store leaf id in upper 32 bits of mask
141  if (CPG) mPLeafMask = (Index64(mIds->find(&leaf)->second) << 32);
142  return true;
143  }
144 
145 protected:
146  /// @brief Constructor to use when a closet point grid is not in use
147  template <bool EnableT = CPG>
148  SignedDistanceFieldTransfer(const size_t pidx,
149  const size_t width,
150  const RadiusType& rt,
151  const math::Transform& source,
152  SdfT& surface,
153  Int64Tree* cpg,
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)
158  , mPIdx(pidx)
159  , mPosition()
160  , mMaxKernelWidth(width)
161  , mRadius(rt)
162  , mBackground(surface.background())
163  , mDx(surface.voxelSize()[0])
164  , mIds(ids)
165  , mPLeafMask(0) {
166  OPENVDB_ASSERT(cpg && ids);
167  }
168 
169  /// @brief Constructor to use when a closet point grid is in use
170  template <bool EnableT = CPG>
171  SignedDistanceFieldTransfer(const size_t pidx,
172  const size_t width,
173  const RadiusType& rt,
174  const math::Transform& source,
175  SdfT& surface,
176  Int64Tree*,
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())
181  , mPIdx(pidx)
182  , mPosition()
183  , mMaxKernelWidth(width)
184  , mRadius(rt)
185  , mBackground(surface.background())
186  , mDx(surface.voxelSize()[0])
187  , mIds(nullptr)
188  , mPLeafMask(0) {}
189 
190  SignedDistanceFieldTransfer(const SignedDistanceFieldTransfer& other)
191  : TransformTransfer(other)
192  , VolumeTransferT(other)
193  , mPIdx(other.mPIdx)
194  , mPosition()
195  , mMaxKernelWidth(other.mMaxKernelWidth)
196  , mRadius(other.mRadius)
197  , mBackground(other.mBackground)
198  , mDx(other.mDx)
199  , mIds(other.mIds)
200  , mPLeafMask(0) {}
201 
202 protected:
203  const size_t mPIdx;
204  typename PositionHandleT::UniquePtr mPosition;
205  const size_t mMaxKernelWidth;
206  RadiusType mRadius;
207  const ValueT mBackground;
208  const double mDx;
209  const std::unordered_map<const PointDataTree::LeafNodeType*, Index>* mIds;
210  Index64 mPLeafMask;
211 };
212 
213 /// @brief The transfer implementation for spherical stamping of narrow band
214 /// radius values.
215 template <typename SdfT,
216  typename PositionCodecT,
217  typename RadiusType,
218  bool CPG>
219 struct SphericalTransfer :
220  public SignedDistanceFieldTransfer<SdfT, PositionCodecT, RadiusType, CPG>
221 {
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;
227  // The precision of the kernel arithmetic
228  using RealT = double;
229 
230  SphericalTransfer(const size_t pidx,
231  const size_t width,
232  const RadiusType& rt,
233  const math::Transform& source,
234  SdfT& surface,
235  Int64Tree* cpg = nullptr,
236  const std::unordered_map<const PointDataTree::LeafNodeType*, Index>* ids = nullptr)
237  : BaseT(pidx, width, rt, source, surface, cpg, ids) {}
238 
239  /// @brief For each point, stamp a sphere with a given radius by running
240  /// over all intersecting voxels and calculating if this point is closer
241  /// than the currently held distance value. Note that the default value
242  /// of the surface buffer should be the background value of the surface.
243  inline void rasterizePoint(const Coord& ijk,
244  const Index id,
245  const CoordBBox& bounds)
246  {
247  // @note GCC 6.3 warns here in optimised builds
248 #if defined(__GNUC__) && !defined(__clang__)
249 #pragma GCC diagnostic push
250 #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
251 #endif
252  Vec3d P = ijk.asVec3d() + Vec3d(this->mPosition->get(id));
253 #if defined(__GNUC__) && !defined(__clang__)
254 #pragma GCC diagnostic pop
255 #endif
256  P = this->transformSourceToTarget(P);
257 
258  const auto& r = this->mRadius.eval(id);
259  const RealT max = r.max();
260 
261  CoordBBox intersectBox(Coord::floor(P - max), Coord::ceil(P + max));
262  intersectBox.intersect(bounds);
263  if (intersectBox.empty()) return;
264 
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>());
268 
269  // If min2 == 0.0, then the index space radius is equal to or less than
270  // the desired half band. In this case each sphere interior always needs
271  // to be filled with distance values as we won't ever reach the negative
272  // background value. If, however, a point overlaps a voxel coord exactly,
273  // x2y2z2 will be 0.0. Forcing min2 to be less than zero here avoids
274  // incorrectly setting these voxels to inactive -background values as
275  // x2y2z2 will never be < 0.0. We still want the lteq logic in the
276  // (x2y2z2 <= min2) check as this is valid when min2 > 0.0.
277  const RealT min2 = r.minSq() == 0.0 ? -1.0 : r.minSq();
278  const RealT max2 = r.maxSq();
279 
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); // unsigned bit shift mult
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 + /*k*/(c.z() & (DIM-1u));
290  if (!mask.isOn(offset)) continue; // inside existing level set or not in range
291 
292  const RealT x2y2z2 = static_cast<RealT>(x2y2 + math::Pow2(c.z() - P[2]));
293  if (x2y2z2 >= max2) continue; //outside narrow band of particle in positive direction
294  if (x2y2z2 <= min2) { //outside narrow band of the particle in negative direction. can disable this to fill interior
295  data[offset] = -(this->mBackground);
296  mask.setOff(offset);
297  continue;
298  }
299 
300  const ValueT d = ValueT(this->mDx * (math::Sqrt(x2y2z2) - r.get())); // back to world space
301  ValueT& v = data[offset];
302  if (d < v) {
303  v = d;
305  if (CPG) cpg[offset] = Int64(this->mPLeafMask | Index64(id));
307  // transfer attributes - we can't use this here as the exposed
308  // function signatures take vector of attributes (i.e. an unbounded
309  // size). If we instead clamped the attribute transfer to a fixed
310  // amount of attributes we could get rid of the closest point logic
311  // entirely. @todo consider adding another signature for increased
312  // performance
313  // this->foreach([&](auto* buffer, const size_t idx) {
314  // using Type = typename std::remove_pointer<decltype(buffer)>::type;
315  // buffer[offset] = mAttributes.template get<Type>(idx);
316  // })
317  }
318  }
319  }
320  } // outer sdf voxel
321  } // point idx
322 
323  /// @brief Allow early termination if all voxels in the surface have been
324  /// deactivated (all interior)
325  inline bool endPointLeaf(const PointDataTree::LeafNodeType&)
326  {
327  // If the mask is off, terminate rasterization
328  return !(this->template mask<0>()->isOff());
329  }
330 
331  inline bool finalize(const Coord&, const size_t)
332  {
333  // loop over voxels in the outer cube diagonals which won't have been
334  // hit by point rasterizations - these will be on because of the mask
335  // fill technique and need to be turned off.
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());
340  }
341  // apply sdf mask to other grids
342  if (CPG) *(this->template mask<CPG ? 1 : 0>()) = mask;
343  return true;
344  }
345 };
346 
347 /// @brief The transfer implementation for averaging of positions followed by
348 /// spherical stamping.
349 template <typename SdfT,
350  typename PositionCodecT,
351  typename RadiusType,
352  bool CPG>
353 struct AveragePositionTransfer :
354  public SignedDistanceFieldTransfer<SdfT, PositionCodecT, RadiusType, CPG>
355 {
356  using BaseT = SignedDistanceFieldTransfer<SdfT, PositionCodecT, RadiusType, CPG>;
357  using typename BaseT::TreeT;
358  using typename BaseT::ValueT;
359 
360  using VolumeTransferT = typename std::conditional<CPG,
361  VolumeTransfer<typename SdfT::TreeType, Int64Tree>,
362  VolumeTransfer<typename SdfT::TreeType>>::type;
363 
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;
367  // The precision of the kernel arithmetic
368  using RealT = double;
369 
370  struct PosRadPair
371  {
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);
378  RealT R = 0.0;
379  };
380 
381  AveragePositionTransfer(const size_t pidx,
382  const size_t width,
383  const RadiusType& rt,
384  const RealT search,
385  const math::Transform& source,
386  SdfT& surface,
387  Int64Tree* cpg = nullptr,
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)
392  , mWeights()
393  , mDist() {}
394 
395  AveragePositionTransfer(const AveragePositionTransfer& other)
396  : BaseT(other)
397  , mMaxSearchIS(other.mMaxSearchIS)
398  , mMaxSearchSqIS(other.mMaxSearchSqIS)
399  , mWeights()
400  , mDist() {}
401 
402  inline void initialize(const Coord& origin, const size_t idx, const CoordBBox& bounds)
403  {
404  // init buffers
405  this->BaseT::initialize(origin, idx, bounds);
406  mWeights.assign(NUM_VALUES, PosRadPair());
407  if (CPG) mDist.assign(NUM_VALUES, std::numeric_limits<float>::max());
408  // We use the surface buffer to store the intermediate weights as
409  // defined by the sum of k(|x−xj|/R), where k(s) = max(0,(1−s^2)^3)
410  // and R is the maximum search distance. The active buffer currently
411  // holds background values. We could simply subtract the background away
412  // from the final result - however if the background value increases
413  // beyond 1, progressively larger floating point instabilities can be
414  // observed with the weight calculation. Instead, reset all active
415  // values to zero
416  // @todo The surface buffer may not be at RealT precision. Should we
417  // enforce this by storing the weights in another vector?
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);
422  }
423  }
424 
425  inline void rasterizePoint(const Coord& ijk,
426  const Index id,
427  const CoordBBox& bounds)
428  {
429 #if defined(__GNUC__) && !defined(__clang__)
430 #pragma GCC diagnostic push
431 #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
432 #endif
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
436 #endif
437  const Vec3d P = this->targetTransform().worldToIndex(PWS);
438 
439  CoordBBox intersectBox(Coord::floor(P - mMaxSearchIS),
440  Coord::ceil(P + mMaxSearchIS));
441  intersectBox.intersect(bounds);
442  if (intersectBox.empty()) return;
443 
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>());
447 
448  // index space radius
449  const auto& r = this->mRadius.eval(id);
450  const RealT rad = r.get();
451  const RealT invsq = 1.0 / mMaxSearchSqIS;
452 
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); // unsigned bit shift mult
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; //outside search distance
464  const Index offset = ij + /*k*/(c.z() & (DIM-1u));
465  if (!mask.isOn(offset)) continue; // inside existing level set or not in range
466 
467  // @note this algorithm is unable to deactivate voxels within
468  // a computed narrow band during rasterization as all points must
469  // visit their affected voxels.
470 
471  if (CPG) {
472  // CPG still computed directly with each individual point
473  // @note Because voxels can't be discarded, it may be faster to
474  // do this as a post process (and avoid the sqrt per lookup)
475  // @note No need to scale back to world space
476  const float dist = static_cast<float>(math::Sqrt(x2y2z2) - r.get());
477  auto& d = mDist[offset];
478  if (dist < d) {
479  d = dist;
481  cpg[offset] = Int64(this->mPLeafMask | Index64(id));
483  }
484  }
485 
486  x2y2z2 *= invsq; // x2y2z2 = (x - xi) / R
487  // k(s) = max(0,(1−s^2)^3). note that the max is unecessary
488  // as we early terminate above with x2y2z2 >= mMaxSearchSqIS
489  x2y2z2 = math::Pow3(1.0 - x2y2z2);
490  OPENVDB_ASSERT(x2y2z2 >= 0.0);
491  // @todo The surface buffer may not be at RealT precision. Should we
492  // enforce this by storing the weights in another vector?
494  data[offset] += x2y2z2;
496  auto& wt = mWeights[offset];
497  wt.addP(PWS * x2y2z2);
498  wt.addR(rad * x2y2z2);
499  }
500  }
501  } // outer sdf voxel
502  } // point idx
503 
504  inline bool endPointLeaf(const PointDataTree::LeafNodeType&) { return true; }
505 
506  inline bool finalize(const Coord& origin, const size_t)
507  {
508  auto& mask = *(this->template mask<0>());
509  auto* const data = this->template buffer<0>();
510 
511  for (auto iter = mask.beginOn(); iter; ++iter) {
512  const Index idx = iter.pos();
513  auto& w = data[idx];
514  // if background, voxel was out of range. Guaranteed to be outside as
515  // all interior voxels will have at least a single point contribution
516  if (w == 0.0f) {
517  mask.setOff(idx);
518  w = this->mBackground;
519  }
520  else {
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); // wi
524  auto& wt = mWeights[idx];
525  wt.multP(wi); // sum of weighted positions
526  wt.multR(wi * this->mDx); // sum of weighted radii (scale to ws)
527  wt.addP(-ws); // (x - xi) (instead doing (-x + xi))
528  w = static_cast<typename SdfT::ValueType>(wt.length()); // (x - xi) - r
529  // clamp active region and value range to requested narrow band
530  if (std::fabs(w) >= this->mBackground) {
531  w = std::copysign(this->mBackground, w);
532  mask.setOff(idx);
533  }
534  }
535  }
536 
537  // apply sdf mask to other grids
538  if (CPG) *(this->template mask<CPG ? 1 : 0>()) = mask;
539  return true;
540  }
541 
542 private:
543  const RealT mMaxSearchIS, mMaxSearchSqIS;
544  std::vector<PosRadPair> mWeights;
545  std::vector<float> mDist;
546 };
547 
548 /// @brief Base class for surfacing mask initialization
549 template <typename MaskTreeT = MaskTree,
550  typename InterrupterT = util::NullInterrupter>
551 struct SurfaceMaskOp
552 {
553  inline void join(SurfaceMaskOp& other)
554  {
555  if (mMask->leafCount() > other.mMask->leafCount()) {
556  mMask->topologyUnion(*other.mMask);
557  other.mMask.reset();
558  }
559  else {
560  other.mMask->topologyUnion(*mMask);
561  mMask.reset(other.mMask.release());
562  }
563 
564  if (mMaskOff->leafCount() > other.mMaskOff->leafCount()) {
565  mMaskOff->topologyUnion(*other.mMaskOff);
566  other.mMaskOff.reset();
567  }
568  else {
569  other.mMaskOff->topologyUnion(*mMaskOff);
570  mMaskOff.reset(other.mMaskOff.release());
571  }
572  }
573 
574  std::unique_ptr<MaskTreeT> mask() { return std::move(mMask); }
575  std::unique_ptr<MaskTreeT> maskoff() { return std::move(mMaskOff); }
576 
577 protected:
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]))) {}
588 
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) {}
596 
597  // @brief Sparse fill a tree with activated bounding boxes expanded from
598  // each active voxel.
599  // @note This method used to fill from each individual voxel. Whilst more
600  // accurate, this was slower in comparison to using the active node bounds.
601  // As the rasterization is so fast (discarding of voxels out of range)
602  // this overzealous activation results in far superior performance here.
603  template <typename LeafT>
604  inline void fill(const LeafT& leaf, const int dist, const bool active)
605  {
606  auto doFill = [&](CoordBBox b) {
607  b = mSurfaceTransform.worldToIndexCellCentered(
608  mPointsTransform.indexToWorld(b));
609  b.expand(dist);
610  if (active) mMask->sparseFill(b, true, true);
611  else mMaskOff->sparseFill(b, true, true);
612  };
613 
614  const auto& mask = leaf.getValueMask();
615  if (mask.isOn()) {
616  doFill(leaf.getNodeBoundingBox());
617  }
618  else {
619  CoordBBox bounds;
620  for (auto iter = mask.beginOn(); iter; ++iter) {
621  bounds.expand(leaf.offsetToLocalCoord(iter.pos()));
622  }
623  if (!bounds.empty()) {
624  bounds.translate(leaf.origin());
625  doFill(bounds);
626  }
627  }
628  }
629 
630  inline bool interrupted()
631  {
632  if (util::wasInterrupted(mInterrupter)) {
633  thread::cancelGroupExecution();
634  return true;
635  }
636  return false;
637  }
638 
639 protected:
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;
645  // add the size of a single points voxel at the surface resolution
646  // to the distance limits - this is becasue we don't query exact point
647  // positions in SurfaceMaskOp, which means we have to
648  // actviate the "worst case scenario" bounds. The closer the points
649  // voxel size is to the target sdf voxel size the better the topology
650  // estimate
651  const int mVoxelOffset;
652 };
653 
654 /// @brief Initializes a fixed activity mask
655 template <typename MaskTreeT,
656  typename InterrupterT = util::NullInterrupter>
657 struct FixedSurfaceMaskOp
658  : public SurfaceMaskOp<MaskTreeT, InterrupterT>
659 {
660  using BaseT = SurfaceMaskOp<MaskTreeT, InterrupterT>;
661  using LeafManagerT =
662  tree::LeafManager<const points::PointDataTree>;
663 
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)
670  , mMin(), mMax() {
671  // calculate the min interior cube area of activity. this is the side
672  // of the largest possible cube that fits into the radius "min":
673  // d = 2r -> 2r = 3x^2 -> x = 2r / sqrt(3)
674  // finally half it to the new radius
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;
677  mMin = math::Max(0, d);
678  mMax = static_cast<int>(math::RoundUp(maxBandRadius));
679  mMax += this->mVoxelOffset;
680  }
681 
682  FixedSurfaceMaskOp(const FixedSurfaceMaskOp& other, tbb::split)
683  : BaseT(other), mMin(other.mMin), mMax(other.mMax) {}
684 
685  void operator()(const typename LeafManagerT::LeafRange& range)
686  {
687  if (this->interrupted()) return;
688  for (auto leafIter = range.begin(); leafIter; ++leafIter) {
689  this->fill(*leafIter, mMax, true);
690  }
691  if (mMin <= 0) return;
692  for (auto leafIter = range.begin(); leafIter; ++leafIter) {
693  this->fill(*leafIter, mMin, false);
694  }
695  }
696 
697 private:
698  int mMin, mMax;
699 };
700 
701 /// @brief Initializes a variable activity mask
702 template <typename RadiusTreeT,
703  typename MaskTreeT,
704  typename InterrupterT = util::NullInterrupter>
705 struct VariableSurfaceMaskOp
706  : public SurfaceMaskOp<MaskTreeT, InterrupterT>
707 {
708  using BaseT = SurfaceMaskOp<MaskTreeT, InterrupterT>;
709  using LeafManagerT =
710  tree::LeafManager<const points::PointDataTree>;
711 
712  VariableSurfaceMaskOp(const math::Transform& pointsTransform,
713  const math::Transform& surfaceTransform,
714  const RadiusTreeT& min,
715  const RadiusTreeT& max,
716  const Real scale,
717  const Real halfband,
718  InterrupterT* interrupter = nullptr)
719  : BaseT(pointsTransform, surfaceTransform, interrupter)
720  , mMin(min), mMax(max), mScale(scale), mHalfband(halfband) {}
721 
722  VariableSurfaceMaskOp(const VariableSurfaceMaskOp& other, tbb::split)
723  : BaseT(other), mMin(other.mMin), mMax(other.mMax)
724  , mScale(other.mScale), mHalfband(other.mHalfband) {}
725 
726  void operator()(const typename LeafManagerT::LeafRange& range)
727  {
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);
733  }
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);
739  }
740  }
741 
742 private:
743  inline int maxDist(const typename RadiusTreeT::ValueType& max) const
744  {
745  // max radius in index space
746  const Real v = (Real(max) * mScale) + mHalfband;
747  int d = int(math::RoundUp(v)); // next voxel
748  d += this->mVoxelOffset; // final offset
749  return d;
750  }
751 
752  inline int minDist(const typename RadiusTreeT::ValueType& min) const
753  {
754  // min radius in index space
755  Real v = math::Max(0.0, (Real(min) * mScale) - mHalfband);
756  v = math::Max(0.0, Real(v - this->mVoxelOffset)); // final offset
757  // calculate the min interior cube area of activity. this is the side
758  // of the largest possible cube that fits into the radius "v":
759  // d = 2r -> 2r = 3x^2 -> x = 2r / sqrt(3)
760  // finally half it to the new radius
761  v = ((2.0 * v) / std::sqrt(3.0)) / 2.0;
762  return static_cast<int>(v); // round down (v is always >= 0)
763  }
764 
765 private:
766  const RadiusTreeT& mMin;
767  const RadiusTreeT& mMax;
768  const Real mScale, mHalfband;
769 };
770 
771 template <typename SdfT, typename InterrupterT, typename PointDataGridT>
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)
779 {
780  using LeafManagerT = tree::LeafManager<const typename PointDataGridT::TreeType>;
781  using MaskTreeT = typename SdfT::TreeType::template ValueConverter<ValueMask>::Type;
782 
783  typename SdfT::Ptr surface = SdfT::create(bg);
784  surface->setTransform(transform);
785  surface->setGridClass(GRID_LEVEL_SET);
786 
787  FixedSurfaceMaskOp<MaskTreeT, InterrupterT> op(points.transform(),
788  *transform, minBandRadius, maxBandRadius, interrupter);
789 
790  if (interrupter) interrupter->start("Generating uniform surface topology");
791 
792  LeafManagerT leafManager(points.tree());
793  tbb::parallel_reduce(leafManager.leafRange(), op);
794  auto mask = op.mask();
795 
796  if (minBandRadius > 0.0) {
797  auto maskoff = op.maskoff();
798  mask->topologyDifference(*maskoff);
799  // union will copy empty nodes so prune them
800  tools::pruneInactive(*mask);
801  surface->tree().topologyUnion(*mask);
802  mask.reset();
803  // set maskoff values to -background
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; });
808  }
809  };
810  tools::foreach(surface->beginValueOff(), setOffOp,
811  /*thread=*/true, /*shared=*/false);
812  maskoff.reset();
813  }
814  else {
815  surface->tree().topologyUnion(*mask);
816  mask.reset();
817  }
818 
819  surface->tree().voxelizeActiveTiles();
820  if (interrupter) interrupter->end();
821  return surface;
822 }
823 
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,
834  const Real scale,
835  const Real halfband,
836  InterrupterT* interrupter)
837 {
838  using LeafManagerT = tree::LeafManager<const typename PointDataGridT::TreeType>;
839  using MaskTreeT = typename SdfT::TreeType::template ValueConverter<ValueMask>::Type;
840 
841  typename SdfT::Ptr surface = SdfT::create(bg);
842  surface->setTransform(transform);
843  surface->setGridClass(GRID_LEVEL_SET);
844 
845  VariableSurfaceMaskOp<RadiusTreeT, MaskTreeT, InterrupterT>
846  op(points.transform(), *transform, min, max, scale, halfband, interrupter);
847 
848  if (interrupter) interrupter->start("Generating variable surface topology");
849 
850  LeafManagerT leafManager(points.tree());
851  tbb::parallel_reduce(leafManager.leafRange(), op);
852  auto mask = op.mask();
853  auto maskoff = op.maskoff();
854 
855  if (!maskoff->empty()) {
856  mask->topologyDifference(*maskoff);
857  // union will copy empty nodes so prune them
858  tools::pruneInactive(*mask);
859  surface->tree().topologyUnion(*mask);
860  // set maskoff values to -background
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; });
865  }
866  };
867  tools::foreach(surface->beginValueOff(), setOffOp,
868  /*thread=*/true, /*shared=*/false);
869  }
870  else {
871  surface->tree().topologyUnion(*mask);
872  }
873 
874  mask.reset();
875  maskoff.reset();
876  surface->tree().voxelizeActiveTiles();
877  if (interrupter) interrupter->end();
878  return surface;
879 }
880 
881 template <typename PointDataTreeT,
882  typename AttributeTypes>
883 inline GridPtrVec
884 transferAttributes(const tree::LeafManager<const PointDataTreeT>& manager,
885  const std::vector<std::string>& attributes,
886  const Int64Tree& cpg,
887  const math::Transform::Ptr transform)
888 {
889  OPENVDB_ASSERT(manager.leafCount() != 0);
890  // masking uses upper 32 bits for leaf node id
891  // @note we can use a point list impl to support larger counts
892  // if necessary but this is far faster
893  OPENVDB_ASSERT(manager.leafCount() <
895 
896  // linearise cpg to avoid having to probe data
897  const tree::LeafManager<const Int64Tree> cpmanager(cpg);
898 
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>;
902 
903  // init topology
904  tree.topologyUnion(cpg);
905  tree::LeafManager<TreeType> lm(tree);
906 
907  // init values
908  lm.foreach([&manager, &cpmanager, attrIdx]
909  (auto& leaf, const size_t idx) {
910  auto voxel = leaf.beginValueOn();
911  if (!voxel) return;
912 
913  auto* data = leaf.buffer().data();
914  const Int64* ids = cpmanager.leaf(idx).buffer().data();
915  Index prev = Index(ids[voxel.pos()] >> 32);
916  typename HandleT::UniquePtr handle(
917  new HandleT(manager.leaf(prev).constAttributeArray(attrIdx)));
918 
919  for (; voxel; ++voxel) {
920  const Int64 hash = ids[voxel.pos()];
921  const Index lfid = Index(hash >> 32); // upper 32 bits to leaf id
922  const Index ptid = static_cast<Index>(hash); // lower
923  if (lfid != prev) {
924  handle.reset(new HandleT(manager.leaf(lfid).constAttributeArray(attrIdx)));
925  prev = lfid;
926  }
927  data[voxel.pos()] = handle->get(ptid);
928  }
929  });
930  };
931 
932  GridPtrVec grids;
933  grids.reserve(attributes.size());
934  const auto& attrSet = manager.leaf(0).attributeSet();
935  tbb::task_group tasks;
936 
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");
943  }
944 
945  const std::string& type = attrSet.descriptor().valueType(attrIdx);
946  GridBase::Ptr grid = nullptr;
947  AttributeTypes::foreach([&](const auto& v) {
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) {
951  auto typed = Grid<TreeT>::create();
952  grid = typed;
953  typed->setName(name);
954  typed->setTransform(transform);
955  tasks.run([typed, attrIdx, transfer] { transfer(typed->tree(), attrIdx); });
956  grids.emplace_back(grid);
957  }
958  });
959 
960  if (!grid) {
961  OPENVDB_THROW(RuntimeError, "No support for attribute type " + type +
962  " built during closest point surface transfer");
963  }
964  }
965 
966  tasks.wait();
967  return grids;
968 }
969 
970 template <typename SdfT,
971  template <typename, typename, typename, bool> class TransferInterfaceT,
972  typename AttributeTypes,
973  typename InterrupterT,
974  typename PointDataGridT,
975  typename FilterT,
976  typename ...Args>
977 inline GridPtrVec
978 rasterizeSurface(const PointDataGridT& points,
979  const std::vector<std::string>& attributes,
980  const FilterT& filter,
981  SdfT& surface,
982  InterrupterT* interrupter,
983  Args&&... args)
984 {
985  using RadRefT = typename std::tuple_element<1, std::tuple<Args...>>::type;
986  using RadT = typename std::remove_reference<RadRefT>::type;
987 
988  GridPtrVec grids;
989  auto leaf = points.constTree().cbeginLeaf();
990  if (!leaf) return grids;
991 
992  const size_t pidx = leaf->attributeSet().find("P");
993  if (pidx == AttributeSet::INVALID_POS) {
994  OPENVDB_THROW(RuntimeError, "Failed to find position attribute");
995  }
996 
997  // @note Can't split this out into a generic lambda yet as there
998  // are compiler issues with capturing variadic arguments
999  const NamePair& ptype = leaf->attributeSet().descriptor().type(pidx);
1000 
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);
1007  }
1008  else {
1009  using TransferT = TransferInterfaceT<SdfT, UnknownCodec, RadT, false>;
1010  TransferT transfer(pidx, args...);
1011  rasterize<PointDataGridT, TransferT, FilterT, InterrupterT>
1012  (points, transfer, filter, interrupter);
1013  }
1014  }
1015  else {
1016  Int64Tree cpg;
1017  cpg.topologyUnion(surface.tree());
1018  tree::LeafManager<const PointDataTree> manager(points.tree());
1019  // map point leaf nodes to their linear id
1020  // @todo sorted vector of leaf ptr-> index pair then lookup with binary search?
1021  std::unordered_map<const PointDataTree::LeafNodeType*, Index> ids;
1022  manager.foreach([&](auto& leaf, size_t idx) { ids[&leaf] = Index(idx); }, false);
1023 
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);
1029  }
1030  else {
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);
1035  }
1036 
1037  ids.clear();
1038  tools::pruneInactive(cpg);
1039  // Build attribute transfer grids
1040  grids = transferAttributes
1041  <typename PointDataGrid::TreeType, AttributeTypes>
1042  (manager, attributes, cpg, surface.transformPtr());
1043  }
1044 
1045  return grids;
1046 }
1047 
1048 } // namespace rasterize_sdf_internal
1049 
1050 /// @endcond
1051 
1052 ///////////////////////////////////////////////////
1053 ///////////////////////////////////////////////////
1054 
1055 
1056 template <typename PointDataGridT,
1057  typename SdfT,
1058  typename FilterT,
1059  typename InterrupterT>
1060 typename SdfT::Ptr
1061 rasterizeSpheres(const PointDataGridT& points,
1062  const Real radius,
1063  const Real halfband,
1064  math::Transform::Ptr transform,
1065  const FilterT& filter,
1066  InterrupterT* interrupter)
1067 {
1068  auto grids =
1069  rasterizeSpheres<PointDataGridT, TypeList<>, SdfT, FilterT, InterrupterT>
1070  (points, radius, {}, halfband, transform, filter, interrupter);
1071  return StaticPtrCast<SdfT>(grids.front());
1072 }
1073 
1074 template <typename PointDataGridT,
1075  typename RadiusT,
1076  typename SdfT,
1077  typename FilterT,
1078  typename InterrupterT>
1079 typename SdfT::Ptr
1080 rasterizeSpheres(const PointDataGridT& points,
1081  const std::string& radius,
1082  const Real scale,
1083  const Real halfband,
1084  math::Transform::Ptr transform,
1085  const FilterT& filter,
1086  InterrupterT* interrupter)
1087 {
1088  auto grids =
1089  rasterizeSpheres<PointDataGridT, TypeList<>, RadiusT, SdfT, FilterT, InterrupterT>
1090  (points, radius, {}, scale, halfband, transform, filter, interrupter);
1091  return StaticPtrCast<SdfT>(grids.front());
1092 }
1093 
1094 template <typename PointDataGridT,
1095  typename AttributeTypes,
1096  typename SdfT,
1097  typename FilterT,
1098  typename InterrupterT>
1099 GridPtrVec
1100 rasterizeSpheres(const PointDataGridT& points,
1101  const Real radius,
1102  const std::vector<std::string>& attributes,
1103  const Real halfband,
1104  math::Transform::Ptr transform,
1105  const FilterT& filter,
1106  InterrupterT* interrupter)
1107 {
1108  using namespace rasterize_sdf_internal;
1109 
1110  if (!transform) transform = points.transform().copy();
1111  const Real vs = transform->voxelSize()[0];
1112  const float background = static_cast<float>(vs * halfband);
1113 
1114  // search distance at the SDF transform, including its half band
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));
1120 
1121  typename SdfT::Ptr surface =
1122  initFixedSdf<SdfT, InterrupterT>
1123  (points, transform, background, minBandRadius, maxBandRadius, interrupter);
1124 
1125  if (interrupter) interrupter->start("Rasterizing particles to level set using constant Spheres");
1126 
1127  GridPtrVec grids =
1128  rasterizeSurface<SdfT, SphericalTransfer, AttributeTypes, InterrupterT>
1129  (points, attributes, filter, *surface, interrupter,
1130  width, rad, points.transform(), *surface); // args
1131 
1132  tools::pruneLevelSet(surface->tree());
1133  if (interrupter) interrupter->end();
1134  grids.insert(grids.begin(), surface);
1135  return grids;
1136 }
1137 
1138 template <typename PointDataGridT,
1139  typename AttributeTypes,
1140  typename RadiusT,
1141  typename SdfT,
1142  typename FilterT,
1143  typename InterrupterT>
1144 GridPtrVec
1145 rasterizeSpheres(const PointDataGridT& points,
1146  const std::string& radius,
1147  const std::vector<std::string>& attributes,
1148  const Real scale,
1149  const Real halfband,
1150  math::Transform::Ptr transform,
1151  const FilterT& filter,
1152  InterrupterT* interrupter)
1153 {
1154  using namespace rasterize_sdf_internal;
1155  using PointDataTreeT = typename PointDataGridT::TreeType;
1156  using RadTreeT = typename PointDataTreeT::template ValueConverter<RadiusT>::Type;
1157 
1158  if (!transform) transform = points.transform().copy();
1159  const Real vs = transform->voxelSize()[0];
1160  const float background = static_cast<float>(vs * halfband);
1161 
1162  RadiusT min(0), max(0);
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());
1166 
1167  // search distance at the SDF transform
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);
1173  mintree.reset();
1174  maxtree.reset();
1175 
1176  auto leaf = points.constTree().cbeginLeaf();
1177  if (!leaf) return GridPtrVec(1, surface);
1178 
1179  // max possible index space radius
1180  const size_t width = static_cast<size_t>
1181  (math::RoundUp((Real(max) * indexSpaceScale) + Real(halfband)));
1182 
1183  const size_t ridx = leaf->attributeSet().find(radius);
1184  if (ridx == AttributeSet::INVALID_POS) {
1185  OPENVDB_THROW(RuntimeError, "Failed to find radius attribute \"" + radius + "\"");
1186  }
1187 
1188  VaryingBandRadius<RadiusT> rad(ridx, RadiusT(halfband), indexSpaceScale);
1189  if (interrupter) interrupter->start("Rasterizing particles to level set using variable Spheres");
1190 
1191  GridPtrVec grids =
1192  rasterizeSurface<SdfT, SphericalTransfer, AttributeTypes, InterrupterT>
1193  (points, attributes, filter, *surface, interrupter,
1194  width, rad, points.transform(), *surface); // args
1195 
1196  if (interrupter) interrupter->end();
1197  tools::pruneLevelSet(surface->tree());
1198  grids.insert(grids.begin(), surface);
1199  return grids;
1200 }
1201 
1202 ///////////////////////////////////////////////////
1203 
1204 template <typename PointDataGridT,
1205  typename SdfT,
1206  typename FilterT,
1207  typename InterrupterT>
1208 typename SdfT::Ptr
1209 rasterizeSmoothSpheres(const PointDataGridT& points,
1210  const Real radius,
1211  const Real searchRadius,
1212  const Real halfband,
1213  math::Transform::Ptr transform,
1214  const FilterT& filter,
1215  InterrupterT* interrupter)
1216 {
1217  auto grids =
1218  rasterizeSmoothSpheres<PointDataGridT, TypeList<>, SdfT, FilterT, InterrupterT>
1219  (points, radius, searchRadius, {}, halfband, transform, filter, interrupter);
1220  return StaticPtrCast<SdfT>(grids.front());
1221 }
1222 
1223 template <typename PointDataGridT,
1224  typename RadiusT,
1225  typename SdfT,
1226  typename FilterT,
1227  typename InterrupterT>
1228 typename SdfT::Ptr
1229 rasterizeSmoothSpheres(const PointDataGridT& points,
1230  const std::string& radius,
1231  const Real radiusScale,
1232  const Real searchRadius,
1233  const Real halfband,
1234  math::Transform::Ptr transform,
1235  const FilterT& filter,
1236  InterrupterT* interrupter)
1237 {
1238  auto grids =
1239  rasterizeSmoothSpheres<PointDataGridT, TypeList<>, RadiusT, SdfT, FilterT, InterrupterT>
1240  (points, radius, radiusScale, searchRadius, {}, halfband, transform, filter, interrupter);
1241  return StaticPtrCast<SdfT>(grids.front());
1242 }
1243 
1244 template <typename PointDataGridT,
1245  typename AttributeTypes,
1246  typename SdfT,
1247  typename FilterT,
1248  typename InterrupterT>
1249 GridPtrVec
1250 rasterizeSmoothSpheres(const PointDataGridT& points,
1251  const Real radius,
1252  const Real searchRadius,
1253  const std::vector<std::string>& attributes,
1254  const Real halfband,
1255  math::Transform::Ptr transform,
1256  const FilterT& filter,
1257  InterrupterT* interrupter)
1258 {
1259  using namespace rasterize_sdf_internal;
1260 
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();
1267 
1268  typename SdfT::Ptr surface =
1269  initFixedSdf<SdfT, InterrupterT>
1270  (points, transform, background, /*min*/0.0, max, interrupter);
1271 
1272  auto leaf = points.constTree().cbeginLeaf();
1273  if (!leaf) return GridPtrVec(1, surface);
1274 
1275  // max possible index space search radius
1276  const size_t width = static_cast<size_t>(math::RoundUp(indexSpaceSearch));
1277 
1278  const FixedRadius<Real> rad(radius / vs);
1279  if (interrupter) interrupter->start("Rasterizing particles to level set using constant Zhu-Bridson");
1280 
1281  GridPtrVec grids =
1282  rasterizeSurface<SdfT, AveragePositionTransfer, AttributeTypes, InterrupterT>
1283  (points, attributes, filter, *surface, interrupter,
1284  width, rad, indexSpaceSearch, points.transform(), *surface); // args
1285 
1286  tools::pruneInactive(surface->tree());
1287  if (interrupter) interrupter->end();
1288  grids.insert(grids.begin(), surface);
1289  return grids;
1290 }
1291 
1292 template <typename PointDataGridT,
1293  typename AttributeTypes,
1294  typename RadiusT,
1295  typename SdfT,
1296  typename FilterT,
1297  typename InterrupterT>
1298 GridPtrVec
1299 rasterizeSmoothSpheres(const PointDataGridT& points,
1300  const std::string& radius,
1301  const Real radiusScale,
1302  const Real searchRadius,
1303  const std::vector<std::string>& attributes,
1304  const Real halfband,
1305  math::Transform::Ptr transform,
1306  const FilterT& filter,
1307  InterrupterT* interrupter)
1308 {
1309  using namespace rasterize_sdf_internal;
1310  using PointDataTreeT = typename PointDataGridT::TreeType;
1311  using RadTreeT = typename PointDataTreeT::template ValueConverter<RadiusT>::Type;
1312 
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;
1317 
1318  RadiusT min(0), max(0);
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());
1322 
1323  // search distance at the SDF transform
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);
1329  mintree.reset();
1330  maxtree.reset();
1331 
1332  auto leaf = points.constTree().cbeginLeaf();
1333  if (!leaf) return GridPtrVec(1, surface);
1334 
1335  // max possible index space search radius
1336  const size_t width = static_cast<size_t>(math::RoundUp(indexSpaceSearch));
1337 
1338  const size_t ridx = leaf->attributeSet().find(radius);
1339  if (ridx == AttributeSet::INVALID_POS) {
1340  OPENVDB_THROW(RuntimeError, "Failed to find radius attribute");
1341  }
1342 
1343  VaryingRadius<RadiusT> rad(ridx, indexSpaceScale);
1344  if (interrupter) interrupter->start("Rasterizing particles to level set using variable Zhu-Bridson");
1345 
1346  GridPtrVec grids =
1347  rasterizeSurface<SdfT, AveragePositionTransfer, AttributeTypes, InterrupterT>
1348  (points, attributes, filter, *surface, interrupter,
1349  width, rad, indexSpaceSearch, points.transform(), *surface); // args
1350 
1351  tools::pruneInactive(surface->tree());
1352  if (interrupter) interrupter->end();
1353  grids.insert(grids.begin(), surface);
1354  return grids;
1355 }
1356 
1357 } // namespace points
1358 } // namespace OPENVDB_VERSION_NAME
1359 } // namespace openvdb
1360 
1361 #endif //OPENVDB_POINTS_RASTERIZE_SDF_IMPL_HAS_BEEN_INCLUDED
tree::Tree4< int64_t, 5, 4, 3 >::Type Int64Tree
Definition: openvdb.h:57
#define OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN
Bracket code with OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN/_END, to inhibit warnings about type conve...
Definition: Platform.h:221
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
void pruneLevelSet(TreeT &tree, bool threaded=true, size_t grainSize=1)
Reduce the memory footprint of a tree by replacing nodes whose values are all inactive with inactive ...
Definition: Prune.h:390
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
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:110
Vec3< double > Vec3d
Definition: Vec3.h:665
OutGridT XformOp & op
Definition: ValueTransformer.h:139
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:49
#define OPENVDB_NO_TYPE_CONVERSION_WARNING_END
Definition: Platform.h:222
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
Definition: Types.h:455
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
SharedPtr< Transform > Ptr
Definition: Transform.h:42
__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
void pruneInactive(TreeT &tree, bool threaded=true, size_t grainSize=1)
Reduce the memory footprint of a tree by replacing with background tiles any nodes whose values are a...
Definition: Prune.h:355
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:106
#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