GCC Code Coverage Report

Directory: ./
File: openvdb/openvdb/points/impl/PointRasterizeSDFImpl.h
Date: 2022-07-25 17:40:05
Exec Total Coverage
Lines: 378 391 96.7%
Functions: 100 144 69.4%
Branches: 298 630 47.3%
Line Branch Exec Source
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-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 {
13 OPENVDB_USE_VERSION_NAMESPACE
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 4600 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 160008 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 1370 FixedBandRadius(const ValueT ris, const ValueT hb)
42 : FixedRadius<ValueT>(ris)
43 2753 , mMinSearchIS(math::Max(ValueT(0.0), ris - hb))
44 1383 , mMaxSearchIS(ris + hb)
45 1370 , mMinSearchSqIS(mMinSearchIS*mMinSearchIS)
46 21 , 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 2836 inline ValueT minSq() const { return mMinSearchSqIS; }
52 2328 inline ValueT max() const { return mMaxSearchIS; }
53 2836 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 8 struct VaryingRadius
62 {
63 using RadiusHandleT = AttributeHandle<ValueT, CodecT>;
64 19 VaryingRadius(const size_t ridx, const ValueT scale = 1.0)
65
4/8
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 5 times.
 8 : mRIdx(ridx), mRHandle(), mScale(scale) {}
66 175 VaryingRadius(const VaryingRadius& other)
67
2/8
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 10 taken 5 times.
✗ Branch 11 not taken.
 8 : mRIdx(other.mRIdx), mRHandle(), mScale(other.mScale) {}
68
69 3070 inline void reset(const PointDataTree::LeafNodeType& leaf)
70 {
71
1/2
✓ Branch 3 taken 3070 times.
✗ Branch 4 not taken.
 3070 mRHandle.reset(new RadiusHandleT(leaf.constAttributeArray(mRIdx)));
72 3070 }
73
74 /// @brief Compute a fixed radius for a specific point
75
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3219 times.
 3219 inline const FixedRadius<ValueT> eval(const Index id) const
76 {
77
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3219 times.
 3219 assert(mRHandle);
78 3219 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
2/8
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 10 taken 8 times.
✗ Branch 11 not taken.
 22 struct VaryingBandRadius : public VaryingRadius<ValueT, CodecT>
90 {
91 using BaseT = VaryingRadius<ValueT, CodecT>;
92 11 VaryingBandRadius(const size_t ridx, const ValueT halfband,
93 const ValueT scale = 1.0)
94 11 : BaseT(ridx, scale), mHalfBand(halfband) {}
95
96 1349 inline const FixedBandRadius<ValueT> eval(const Index id) const
97 {
98 1349 const auto r = this->BaseT::eval(id).get();
99
2/2
✓ Branch 0 taken 753 times.
✓ Branch 1 taken 596 times.
 1349 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 2395 inline Int32 range(const Coord&, size_t) const { return Int32(mMaxKernelWidth); }
135
136 16398 inline bool startPointLeaf(const PointDataTree::LeafNodeType& leaf)
137 {
138
1/2
✓ Branch 3 taken 8199 times.
✗ Branch 4 not taken.
 16398 mPosition.reset(new PositionHandleT(leaf.constAttributeArray(mPIdx)));
139 6140 mRadius.reset(leaf);
140 // if CPG, store leaf id in upper 32 bits of mask
141 3138 if (CPG) mPLeafMask = (Index64(mIds->find(&leaf)->second) << 32);
142 16398 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 24 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 24 , mDx(surface.voxelSize()[0])
164 , mIds(ids)
165
1/2
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
 24 , mPLeafMask(0) {
166
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
 24 assert(cpg && ids);
167 }
168
169 /// @brief Constructor to use when a closet point grid is in use
170 template <bool EnableT = CPG>
171
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
 74 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 74 , mDx(surface.voxelSize()[0])
187 , mIds(nullptr)
188
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
 148 , mPLeafMask(0) {}
189
190 615 SignedDistanceFieldTransfer(const SignedDistanceFieldTransfer& other)
191 : TransformTransfer(other)
192 , VolumeTransferT(other)
193 615 , mPIdx(other.mPIdx)
194 , mPosition()
195 615 , mMaxKernelWidth(other.mMaxKernelWidth)
196 , mRadius(other.mRadius)
197 615 , mBackground(other.mBackground)
198 615 , mDx(other.mDx)
199 615 , mIds(other.mIds)
200 615 , 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 17 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 30 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 30 : 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 7354 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 7354 Vec3d P = ijk.asVec3d() + Vec3d(this->mPosition->get(id));
253 #if defined(__GNUC__) && !defined(__clang__)
254 #pragma GCC diagnostic pop
255 #endif
256 7354 P = this->transformSourceToTarget(P);
257
258 2698 const auto& r = this->mRadius.eval(id);
259 2698 const RealT max = r.max();
260
261 7354 CoordBBox intersectBox(Coord::floor(P - max), Coord::ceil(P + max));
262 7354 intersectBox.intersect(bounds);
263 1682 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
2/2
✓ Branch 0 taken 1484 times.
✓ Branch 1 taken 1352 times.
 5672 const RealT min2 = r.minSq() == 0.0 ? -1.0 : r.minSq();
278 1948 const RealT max2 = r.maxSq();
279
280 const Coord& a(intersectBox.min());
281 const Coord& b(intersectBox.max());
282
2/2
✓ Branch 0 taken 15437 times.
✓ Branch 1 taken 2836 times.
 36546 for (Coord c = a; c.x() <= b.x(); ++c.x()) {
283 30874 const RealT x2 = static_cast<RealT>(math::Pow2(c.x() - P[0]));
284 30874 const Index i = ((c.x() & (DIM-1u)) << 2*LOG2DIM); // unsigned bit shift mult
285
2/2
✓ Branch 0 taken 86267 times.
✓ Branch 1 taken 15437 times.
 203408 for (c.y() = a.y(); c.y() <= b.y(); ++c.y()) {
286 172534 const RealT x2y2 = static_cast<RealT>(x2 + math::Pow2(c.y() - P[1]));
287 172534 const Index ij = i + ((c.y() & (DIM-1u)) << LOG2DIM);
288
2/2
✓ Branch 0 taken 487347 times.
✓ Branch 1 taken 86267 times.
 1147228 for (c.z() = a.z(); c.z() <= b.z(); ++c.z()) {
289 974694 const Index offset = ij + /*k*/(c.z() & (DIM-1u));
290
2/2
✓ Branch 1 taken 13932 times.
✓ Branch 2 taken 473415 times.
 974694 if (!mask.isOn(offset)) continue; // inside existing level set or not in range
291
292
2/2
✓ Branch 0 taken 319444 times.
✓ Branch 1 taken 153971 times.
 946830 const RealT x2y2z2 = static_cast<RealT>(x2y2 + math::Pow2(c.z() - P[2]));
293
2/2
✓ Branch 0 taken 319444 times.
✓ Branch 1 taken 153971 times.
 946830 if (x2y2z2 >= max2) continue; //outside narrow band of particle in positive direction
294
2/2
✓ Branch 0 taken 4219 times.
✓ Branch 1 taken 149752 times.
 307942 if (x2y2z2 <= min2) { //outside narrow band of the particle in negative direction. can disable this to fill interior
295 8438 data[offset] = -(this->mBackground);
296 8438 mask.setOff(offset);
297 8438 continue;
298 }
299
300
2/2
✓ Branch 0 taken 94168 times.
✓ Branch 1 taken 55584 times.
 299504 const ValueT d = ValueT(this->mDx * (math::Sqrt(x2y2z2) - r.get())); // back to world space
301 299504 ValueT& v = data[offset];
302
2/2
✓ Branch 0 taken 94168 times.
✓ Branch 1 taken 55584 times.
 299504 if (d < v) {
303 188336 v = d;
304 OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN
305 23628 if (CPG) cpg[offset] = Int64(this->mPLeafMask | Index64(id));
306 OPENVDB_NO_TYPE_CONVERSION_WARNING_END
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 2666 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
2/2
✓ Branch 1 taken 291426 times.
✓ Branch 2 taken 1333 times.
 588184 for (auto iter = mask.beginOn(); iter; ++iter) {
339
2/2
✓ Branch 0 taken 222447 times.
✓ Branch 1 taken 68979 times.
 582852 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 2666 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 222407 template <typename S> inline void addR(const S r) { R += r; }
374 101607 template <typename S> inline void multR(const S w) { R *= w; }
375 template <typename S> inline void multP(const S w) { P *= w; }
376 203214 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 19 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 19 , mMaxSearchSqIS(search*search)
392 , mWeights()
393
10/64
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✓ Branch 31 taken 1 times.
✗ Branch 32 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✓ Branch 43 taken 1 times.
✗ Branch 44 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✓ Branch 49 taken 4 times.
✗ Branch 50 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✓ Branch 55 taken 1 times.
✗ Branch 56 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
✗ Branch 64 not taken.
✗ Branch 65 not taken.
✓ Branch 67 taken 1 times.
✗ Branch 68 not taken.
✗ Branch 70 not taken.
✗ Branch 71 not taken.
✗ Branch 73 not taken.
✗ Branch 74 not taken.
✗ Branch 76 not taken.
✗ Branch 77 not taken.
✓ Branch 79 taken 1 times.
✗ Branch 80 not taken.
✗ Branch 82 not taken.
✗ Branch 83 not taken.
✓ Branch 85 taken 7 times.
✗ Branch 86 not taken.
✗ Branch 88 not taken.
✗ Branch 89 not taken.
✓ Branch 91 taken 1 times.
✗ Branch 92 not taken.
✗ Branch 94 not taken.
✗ Branch 95 not taken.
 19 , mDist() {}
394
395 19 AveragePositionTransfer(const AveragePositionTransfer& other)
396 : BaseT(other)
397 19 , mMaxSearchIS(other.mMaxSearchIS)
398 19 , mMaxSearchSqIS(other.mMaxSearchSqIS)
399 , mWeights()
400 6 , mDist() {}
401
402 2124 inline void initialize(const Coord& origin, const size_t idx, const CoordBBox& bounds)
403 {
404 // init buffers
405 2124 this->BaseT::initialize(origin, idx, bounds);
406
2/2
✓ Branch 0 taken 2124 times.
✓ Branch 1 taken 1062 times.
 6372 mWeights.assign(NUM_VALUES, PosRadPair());
407 768 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
2/2
✓ Branch 1 taken 216409 times.
✓ Branch 2 taken 1062 times.
 437066 for (auto iter = mask.beginOn(); iter; ++iter) {
421 432818 data[iter.pos()] = ValueT(0);
422 }
423 }
424
425 5478 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 10956 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 5478 CoordBBox intersectBox(Coord::floor(P - mMaxSearchIS),
440 5478 Coord::ceil(P + mMaxSearchIS));
441 5478 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 3740 const auto& r = this->mRadius.eval(id);
450 3740 const RealT rad = r.get();
451 5478 const RealT invsq = 1.0 / mMaxSearchSqIS;
452
453 const Coord& a(intersectBox.min());
454 const Coord& b(intersectBox.max());
455
2/2
✓ Branch 0 taken 14109 times.
✓ Branch 1 taken 2739 times.
 33696 for (Coord c = a; c.x() <= b.x(); ++c.x()) {
456 28218 const RealT x2 = static_cast<RealT>(math::Pow2(c.x() - P[0]));
457 28218 const Index i = ((c.x() & (DIM-1u)) << 2*LOG2DIM); // unsigned bit shift mult
458
2/2
✓ Branch 0 taken 74889 times.
✓ Branch 1 taken 14109 times.
 177996 for (c.y() = a.y(); c.y() <= b.y(); ++c.y()) {
459 149778 const RealT x2y2 = static_cast<RealT>(x2 + math::Pow2(c.y() - P[1]));
460 149778 const Index ij = i + ((c.y() & (DIM-1u)) << LOG2DIM);
461
2/2
✓ Branch 0 taken 410914 times.
✓ Branch 1 taken 74889 times.
 971606 for (c.z() = a.z(); c.z() <= b.z(); ++c.z()) {
462
2/2
✓ Branch 0 taken 182113 times.
✓ Branch 1 taken 228801 times.
 821828 RealT x2y2z2 = static_cast<RealT>(x2y2 + math::Pow2(c.z() - P[2]));
463
2/2
✓ Branch 0 taken 182113 times.
✓ Branch 1 taken 228801 times.
 821828 if (x2y2z2 >= mMaxSearchSqIS) continue; //outside search distance
464 457602 const Index offset = ij + /*k*/(c.z() & (DIM-1u));
465
2/2
✓ Branch 1 taken 6394 times.
✓ Branch 2 taken 222407 times.
 457602 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 201216 const float dist = static_cast<float>(math::Sqrt(x2y2z2) - r.get());
477
2/2
✓ Branch 0 taken 51999 times.
✓ Branch 1 taken 48609 times.
 201216 auto& d = mDist[offset];
478
2/2
✓ Branch 0 taken 51999 times.
✓ Branch 1 taken 48609 times.
 201216 if (dist < d) {
479 103998 d = dist;
480 OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN
481 103998 cpg[offset] = Int64(this->mPLeafMask | Index64(id));
482 OPENVDB_NO_TYPE_CONVERSION_WARNING_END
483 }
484 }
485
486 444814 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
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 222407 times.
 444814 x2y2z2 = math::Pow3(1.0 - x2y2z2);
490
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 222407 times.
 444814 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?
493 OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN
494 444814 data[offset] += x2y2z2;
495 OPENVDB_NO_TYPE_CONVERSION_WARNING_END
496 auto& wt = mWeights[offset];
497 444814 wt.addP(PWS * x2y2z2);
498 444814 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 2124 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
2/2
✓ Branch 1 taken 216409 times.
✓ Branch 2 taken 1062 times.
 437066 for (auto iter = mask.beginOn(); iter; ++iter) {
512 const Index idx = iter.pos();
513 432818 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
2/2
✓ Branch 0 taken 114802 times.
✓ Branch 1 taken 101607 times.
 432818 if (w == 0.0f) {
517 229604 mask.setOff(idx);
518 229604 w = this->mBackground;
519 }
520 else {
521 203214 const Coord ijk = origin + TreeT::LeafNodeType::offsetToLocalCoord(idx);
522 const Vec3d ws = this->targetTransform().indexToWorld(ijk);
523 203214 const RealT wi = RealT(1.0) / RealT(w); // wi
524 auto& wt = mWeights[idx];
525 wt.multP(wi); // sum of weighted positions
526 203214 wt.multR(wi * this->mDx); // sum of weighted radii (scale to ws)
527 wt.addP(-ws); // (x - xi) (instead doing (-x + xi))
528
2/2
✓ Branch 1 taken 63909 times.
✓ Branch 2 taken 37698 times.
 203214 w = static_cast<typename SdfT::ValueType>(wt.length()); // (x - xi) - r
529 // clamp active region and value range to requested narrow band
530
2/2
✓ Branch 0 taken 63909 times.
✓ Branch 1 taken 37698 times.
 203214 if (std::fabs(w) >= this->mBackground) {
531 127818 w = std::copysign(this->mBackground, w);
532 127818 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 2124 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 56 inline void join(SurfaceMaskOp& other)
554 {
555
2/2
✓ Branch 2 taken 30 times.
✓ Branch 3 taken 26 times.
 56 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
2/2
✓ Branch 2 taken 9 times.
✓ Branch 3 taken 47 times.
 56 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 56 }
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 57 SurfaceMaskOp(const math::Transform& points,
579 const math::Transform& surface,
580 InterrupterT* interrupter = nullptr)
581 57 : mMask(new MaskTreeT)
582
1/2
✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
 57 , mMaskOff(new MaskTreeT)
583 , mPointsTransform(points)
584 , mSurfaceTransform(surface)
585 , mInterrupter(interrupter)
586
1/4
✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 57 , mVoxelOffset(static_cast<int>(math::RoundUp(points.voxelSize()[0] /
587
1/2
✓ Branch 1 taken 57 times.
✗ Branch 2 not taken.
 114 surface.voxelSize()[0]))) {}
588
589 56 SurfaceMaskOp(const SurfaceMaskOp& other)
590 56 : mMask(new MaskTreeT)
591
1/2
✓ Branch 1 taken 56 times.
✗ Branch 2 not taken.
 56 , mMaskOff(new MaskTreeT)
592 56 , mPointsTransform(other.mPointsTransform)
593 56 , mSurfaceTransform(other.mSurfaceTransform)
594 56 , mInterrupter(other.mInterrupter)
595 56 , 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 221 inline void fill(const LeafT& leaf, const int dist, const bool active)
605 {
606 442 auto doFill = [&](CoordBBox b) {
607 221 b = mSurfaceTransform.worldToIndexCellCentered(
608 221 mPointsTransform.indexToWorld(b));
609
2/2
✓ Branch 0 taken 189 times.
✓ Branch 1 taken 32 times.
 221 b.expand(dist);
610
2/2
✓ Branch 0 taken 189 times.
✓ Branch 1 taken 32 times.
 410 if (active) mMask->sparseFill(b, true, true);
611 64 else mMaskOff->sparseFill(b, true, true);
612 };
613
614 const auto& mask = leaf.getValueMask();
615
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 221 times.
 221 if (mask.isOn()) {
616 ✗ doFill(leaf.getNodeBoundingBox());
617 }
618 else {
619 221 CoordBBox bounds;
620
2/2
✓ Branch 1 taken 221 times.
✓ Branch 2 taken 221 times.
 663 for (auto iter = mask.beginOn(); iter; ++iter) {
621 221 bounds.expand(leaf.offsetToLocalCoord(iter.pos()));
622 }
623 if (!bounds.empty()) {
624 bounds.translate(leaf.origin());
625 221 doFill(bounds);
626 }
627 }
628 221 }
629
630 189 inline bool interrupted()
631 {
632
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 189 times.
 189 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 73 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 34 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 34 , 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
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 32 times.
 34 const Real min = ((2.0 * minBandRadius) / std::sqrt(3.0)) / 2.0;
676 34 const int d = static_cast<int>(math::RoundDown(min)) - this->mVoxelOffset;
677
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 32 times.
 34 mMin = math::Max(0, d);
678 34 mMax = static_cast<int>(math::RoundUp(maxBandRadius));
679 34 mMax += this->mVoxelOffset;
680 34 }
681
682 39 FixedSurfaceMaskOp(const FixedSurfaceMaskOp& other, tbb::split)
683 39 : BaseT(other), mMin(other.mMin), mMax(other.mMax) {}
684
685 128 void operator()(const typename LeafManagerT::LeafRange& range)
686 {
687
1/2
✓ Branch 1 taken 128 times.
✗ Branch 2 not taken.
 128 if (this->interrupted()) return;
688
2/2
✓ Branch 1 taken 128 times.
✓ Branch 2 taken 128 times.
 256 for (auto leafIter = range.begin(); leafIter; ++leafIter) {
689 128 this->fill(*leafIter, mMax, true);
690 }
691
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 112 times.
 128 if (mMin <= 0) return;
692
2/2
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 16 times.
 32 for (auto leafIter = range.begin(); leafIter; ++leafIter) {
693 16 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 40 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 23 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
1/2
✓ Branch 1 taken 23 times.
✗ Branch 2 not taken.
 23 , mMin(min), mMax(max), mScale(scale), mHalfband(halfband) {}
721
722 17 VariableSurfaceMaskOp(const VariableSurfaceMaskOp& other, tbb::split)
723 17 : BaseT(other), mMin(other.mMin), mMax(other.mMax)
724 17 , mScale(other.mScale), mHalfband(other.mHalfband) {}
725
726 61 void operator()(const typename LeafManagerT::LeafRange& range)
727 {
728
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 61 times.
 61 if (this->interrupted()) return;
729 61 const tree::ValueAccessor<const RadiusTreeT> maxacc(mMax);
730
2/2
✓ Branch 1 taken 61 times.
✓ Branch 2 taken 61 times.
 122 for (auto leafIter = range.begin(); leafIter; ++leafIter) {
731
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 const int max = this->maxDist(maxacc.getValue(leafIter->origin()));
732
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 this->fill(*leafIter, max, true);
733 }
734
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 const tree::ValueAccessor<const RadiusTreeT> minacc(mMin);
735
2/2
✓ Branch 1 taken 61 times.
✓ Branch 2 taken 61 times.
 122 for (auto leafIter = range.begin(); leafIter; ++leafIter) {
736
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 const int min = this->minDist(minacc.getValue(leafIter->origin()));
737
2/2
✓ Branch 0 taken 45 times.
✓ Branch 1 taken 16 times.
 61 if (min <= 0) continue;
738
1/2
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
 16 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 61 const Real v = (Real(max) * mScale) + mHalfband;
747 61 int d = int(math::RoundUp(v)); // next voxel
748 61 d += this->mVoxelOffset; // final offset
749 return d;
750 }
751
752 61 inline int minDist(const typename RadiusTreeT::ValueType& min) const
753 {
754 // min radius in index space
755
2/2
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 37 times.
 61 Real v = math::Max(0.0, (Real(min) * mScale) - mHalfband);
756
2/2
✓ Branch 0 taken 23 times.
✓ Branch 1 taken 38 times.
 61 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 61 v = ((2.0 * v) / std::sqrt(3.0)) / 2.0;
762 61 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 34 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 34 typename SdfT::Ptr surface = SdfT::create(bg);
784
2/6
✓ Branch 2 taken 34 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 34 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 68 surface->setTransform(transform);
785
1/2
✓ Branch 1 taken 34 times.
✗ Branch 2 not taken.
 34 surface->setGridClass(GRID_LEVEL_SET);
786
787
1/2
✓ Branch 1 taken 34 times.
✗ Branch 2 not taken.
 34 FixedSurfaceMaskOp<MaskTreeT, InterrupterT> op(points.transform(),
788 *transform, minBandRadius, maxBandRadius, interrupter);
789
790
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 34 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 34 if (interrupter) interrupter->start("Generating uniform surface topology");
791
792
1/2
✓ Branch 1 taken 34 times.
✗ Branch 2 not taken.
 68 LeafManagerT leafManager(points.tree());
793
3/4
✓ Branch 1 taken 34 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 7 times.
✓ Branch 4 taken 27 times.
 34 tbb::parallel_reduce(leafManager.leafRange(), op);
794 auto mask = op.mask();
795
796
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 27 times.
 34 if (minBandRadius > 0.0) {
797 auto maskoff = op.maskoff();
798 mask->topologyDifference(*maskoff);
799 // union will copy empty nodes so prune them
800
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 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
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
 7 auto setOffOp = [acc](auto& iter) {
806
2/2
✓ Branch 2 taken 432 times.
✓ Branch 3 taken 1909813 times.
 1910245 if (acc.isValueOn(iter.getCoord())) {
807 864 iter.modifyValue([](auto& v) { v = -v; });
808 }
809 };
810
2/6
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 7 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 7 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
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 34 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 34 if (interrupter) interrupter->end();
821 34 return surface;
822 }
823
824 template <typename SdfT,
825 typename InterrupterT,
826 typename PointDataGridT,
827 typename RadiusTreeT>
828 inline typename SdfT::Ptr
829 23 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 23 typename SdfT::Ptr surface = SdfT::create(bg);
842
2/6
✓ Branch 2 taken 23 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 23 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 46 surface->setTransform(transform);
843
1/2
✓ Branch 1 taken 23 times.
✗ Branch 2 not taken.
 23 surface->setGridClass(GRID_LEVEL_SET);
844
845 VariableSurfaceMaskOp<RadiusTreeT, MaskTreeT, InterrupterT>
846 op(points.transform(), *transform, min, max, scale, halfband, interrupter);
847
848
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 23 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 23 if (interrupter) interrupter->start("Generating variable surface topology");
849
850
1/2
✓ Branch 1 taken 23 times.
✗ Branch 2 not taken.
 46 LeafManagerT leafManager(points.tree());
851
1/2
✓ Branch 1 taken 23 times.
✗ Branch 2 not taken.
 23 tbb::parallel_reduce(leafManager.leafRange(), op);
852 auto mask = op.mask();
853 auto maskoff = op.maskoff();
854
855
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 19 times.
 23 if (!maskoff->empty()) {
856 mask->topologyDifference(*maskoff);
857 // union will copy empty nodes so prune them
858
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 tools::pruneInactive(*mask);
859 surface->tree().topologyUnion(*mask);
860 // set maskoff values to -background
861 tree::ValueAccessor<const MaskTreeT> acc(*maskoff);
862
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 4 auto setOffOp = [acc](auto& iter) {
863
2/2
✓ Branch 2 taken 432 times.
✓ Branch 3 taken 1292454 times.
 1292886 if (acc.isValueOn(iter.getCoord())) {
864 864 iter.modifyValue([](auto& v) { v = -v; });
865 }
866 };
867
2/6
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 4 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
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 23 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 23 if (interrupter) interrupter->end();
878 23 return surface;
879 }
880
881 template <typename PointDataTreeT,
882 typename AttributeTypes>
883 inline GridPtrVec
884
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
 24 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
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
 24 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
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
 24 assert(manager.leafCount() <
894 size_t(std::numeric_limits<Index>::max()));
895
896 // linearise cpg to avoid having to probe data
897 48 const tree::LeafManager<const Int64Tree> cpmanager(cpg);
898
899
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 40 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 8 tree::LeafManager<TreeType> lm(tree);
906
907 // init values
908
2/6
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 344 lm.foreach([&manager, &cpmanager, attrIdx]
909 (auto& leaf, const size_t idx) {
910 auto voxel = leaf.beginValueOn();
911
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 56 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 56 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 112 if (!voxel) return;
912
913 112 auto* data = leaf.buffer().data();
914 112 const Int64* ids = cpmanager.leaf(idx).buffer().data();
915
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 56 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 56 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 112 Index prev = Index(ids[voxel.pos()] >> 32);
916
2/6
✓ Branch 2 taken 56 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 56 times.
✗ Branch 8 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 112 typename HandleT::UniquePtr handle(
917 new HandleT(manager.leaf(prev).constAttributeArray(attrIdx)));
918
919
4/6
✓ Branch 0 taken 4307 times.
✓ Branch 1 taken 56 times.
✓ Branch 2 taken 4307 times.
✓ Branch 3 taken 56 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 8726 for (; voxel; ++voxel) {
920 8614 const Int64 hash = ids[voxel.pos()];
921 8614 const Index lfid = Index(hash >> 32); // upper 32 bits to leaf id
922 8614 const Index ptid = static_cast<Index>(hash); // lower
923
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 4307 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4307 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 8614 if (lfid != prev) {
924 ✗ handle.reset(new HandleT(manager.leaf(lfid).constAttributeArray(attrIdx)));
925 prev = lfid;
926 }
927
2/6
✓ Branch 1 taken 4307 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 4307 times.
✗ Branch 6 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 8614 data[voxel.pos()] = handle->get(ptid);
928 }
929 });
930 };
931
932 GridPtrVec grids;
933
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 24 grids.reserve(attributes.size());
934 const auto& attrSet = manager.leaf(0).attributeSet();
935 tbb::task_group tasks;
936
937
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 4 times.
 40 for (const auto& name : attributes) {
938
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 32 const size_t attrIdx = attrSet.find(name);
939
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 16 times.
 32 if (attrIdx == points::AttributeSet::INVALID_POS) continue;
940
3/6
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 16 times.
 32 if (attrSet.get(attrIdx)->stride() != 1) {
941 ✗ OPENVDB_THROW(RuntimeError, "Transfer of attribute " + name +
942 " not supported since it is strided");
943 }
944
945
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 32 const std::string& type = attrSet.descriptor().valueType(attrIdx);
946 GridBase::Ptr grid = nullptr;
947
3/4
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 8 times.
 68 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
8/12
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 4 times.
✓ Branch 7 taken 4 times.
✓ Branch 8 taken 4 times.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 4 times.
 32 if (!grid && typeNameAsString<ValueType>() == type) {
951 auto typed = Grid<TreeT>::create();
952 grid = typed;
953
2/6
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 8 typed->setName(name);
954
2/6
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 16 typed->setTransform(transform);
955
2/12
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 9 taken 4 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
 16 tasks.run([typed, attrIdx, transfer] { transfer(typed->tree(), attrIdx); });
956
2/6
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 8 grids.emplace_back(grid);
957 }
958 });
959
960
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 8 times.
 24 if (!grid) {
961
4/10
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 8 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 8 times.
✗ Branch 12 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
 80 OPENVDB_THROW(RuntimeError, "No support for attribute type " + type +
962 " built during closest point surface transfer");
963 }
964 }
965
966
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 tasks.wait();
967 8 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
1/2
✓ Branch 1 taken 51 times.
✗ Branch 2 not taken.
 102 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
2/2
✓ Branch 0 taken 49 times.
✓ Branch 1 taken 2 times.
 102 if (!leaf) return grids;
991
992
2/4
✓ Branch 1 taken 49 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 49 times.
✗ Branch 5 not taken.
 98 const size_t pidx = leaf->attributeSet().find("P");
993
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 49 times.
 98 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
1/2
✓ Branch 1 taken 49 times.
✗ Branch 2 not taken.
 98 const NamePair& ptype = leaf->attributeSet().descriptor().type(pidx);
1000
1001
2/2
✓ Branch 0 taken 37 times.
✓ Branch 1 taken 12 times.
 98 if (attributes.empty()) {
1002
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
 74 if (ptype.second == NullCodec::name()) {
1003 using TransferT = TransferInterfaceT<SdfT, NullCodec, RadT, false>;
1004
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
 100 TransferT transfer(pidx, args...);
1005 rasterize<PointDataGridT, TransferT, FilterT, InterrupterT>
1006
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
 74 (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 48 Int64Tree cpg;
1017 cpg.topologyUnion(surface.tree());
1018
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 48 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
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 78 manager.foreach([&](auto& leaf, size_t idx) { ids[&leaf] = Index(idx); }, false);
1023
1024
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 24 if (ptype.second == NullCodec::name()) {
1025 using TransferT = TransferInterfaceT<SdfT, NullCodec, RadT, true>;
1026
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 36 TransferT transfer(pidx, args..., &cpg, &ids);
1027 rasterize<PointDataGridT, TransferT, FilterT, InterrupterT>
1028
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 24 (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
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 24 tools::pruneInactive(cpg);
1039 // Build attribute transfer grids
1040
3/4
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 8 times.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 56 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 34 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
2/6
✓ Branch 1 taken 17 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 17 times.
✗ Branch 4 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 102 auto grids =
1069 rasterizeSpheres<PointDataGridT, TypeList<>, SdfT, FilterT, InterrupterT>
1070 (points, radius, {}, halfband, transform, filter, interrupter);
1071 34 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 18 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
2/6
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 9 times.
✗ Branch 4 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 54 auto grids =
1089 rasterizeSpheres<PointDataGridT, TypeList<>, RadiusT, SdfT, FilterT, InterrupterT>
1090 (points, radius, {}, scale, halfband, transform, filter, interrupter);
1091 18 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
2/2
✓ Branch 0 taken 17 times.
✓ Branch 1 taken 4 times.
 42 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
2/2
✓ Branch 0 taken 17 times.
✓ Branch 1 taken 4 times.
 76 if (!transform) transform = points.transform().copy();
1111 42 const Real vs = transform->voxelSize()[0];
1112 42 const float background = static_cast<float>(vs * halfband);
1113
1114 // search distance at the SDF transform, including its half band
1115
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 14 times.
 42 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 42 const size_t width = static_cast<size_t>(math::RoundUp(maxBandRadius));
1120
1121
1/2
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
 42 typename SdfT::Ptr surface =
1122 initFixedSdf<SdfT, InterrupterT>
1123 (points, transform, background, minBandRadius, maxBandRadius, interrupter);
1124
1125
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 21 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 42 if (interrupter) interrupter->start("Rasterizing particles to level set using constant Spheres");
1126
1127
2/2
✓ Branch 1 taken 19 times.
✓ Branch 2 taken 2 times.
 42 GridPtrVec grids =
1128 rasterizeSurface<SdfT, SphericalTransfer, AttributeTypes, InterrupterT>
1129 (points, attributes, filter, *surface, interrupter,
1130 width, rad, points.transform(), *surface); // args
1131
1132
1/2
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
 38 tools::pruneLevelSet(surface->tree());
1133
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 19 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 38 if (interrupter) interrupter->end();
1134
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 19 times.
✓ Branch 2 taken 19 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 38 grids.insert(grids.begin(), surface);
1135 38 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
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 5 times.
 26 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
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 5 times.
 42 if (!transform) transform = points.transform().copy();
1159 26 const Real vs = transform->voxelSize()[0];
1160 26 const float background = static_cast<float>(vs * halfband);
1161
1162 26 RadiusT min(0), max(0);
1163
1/2
✓ Branch 3 taken 13 times.
✗ Branch 4 not taken.
 52 typename RadTreeT::Ptr mintree(new RadTreeT), maxtree(new RadTreeT);
1164 points::evalMinMax<RadiusT, UnknownCodec>
1165
1/2
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
 26 (points.tree(), radius, min, max, filter, mintree.get(), maxtree.get());
1166
1167 // search distance at the SDF transform
1168 26 const RadiusT indexSpaceScale = RadiusT(scale / vs);
1169
3/6
✓ Branch 2 taken 13 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
 56 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
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 11 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
 38 if (!leaf) return GridPtrVec(1, surface);
1178
1179 // max possible index space radius
1180
1/2
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
 22 const size_t width = static_cast<size_t>
1181
1/2
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
 22 (math::RoundUp((Real(max) * indexSpaceScale) + Real(halfband)));
1182
1183
1/2
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
 22 const size_t ridx = leaf->attributeSet().find(radius);
1184
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
 22 if (ridx == AttributeSet::INVALID_POS) {
1185 ✗ OPENVDB_THROW(RuntimeError, "Failed to find radius attribute \"" + radius + "\"");
1186 }
1187
1188
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
 22 VaryingBandRadius<RadiusT> rad(ridx, RadiusT(halfband), indexSpaceScale);
1189
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 22 if (interrupter) interrupter->start("Rasterizing particles to level set using variable Spheres");
1190
1191
2/2
✓ Branch 1 taken 9 times.
✓ Branch 2 taken 2 times.
 40 GridPtrVec grids =
1192 rasterizeSurface<SdfT, SphericalTransfer, AttributeTypes, InterrupterT>
1193 (points, attributes, filter, *surface, interrupter,
1194 width, rad, points.transform(), *surface); // args
1195
1196
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 18 if (interrupter) interrupter->end();
1197
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 18 tools::pruneLevelSet(surface->tree());
1198
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 18 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 18 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
2/6
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 9 times.
✗ Branch 4 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 54 auto grids =
1218 rasterizeSmoothSpheres<PointDataGridT, TypeList<>, SdfT, FilterT, InterrupterT>
1219 (points, radius, searchRadius, {}, halfband, transform, filter, interrupter);
1220 18 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 12 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
2/6
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 36 auto grids =
1239 rasterizeSmoothSpheres<PointDataGridT, TypeList<>, RadiusT, SdfT, FilterT, InterrupterT>
1240 (points, radius, radiusScale, searchRadius, {}, halfband, transform, filter, interrupter);
1241 12 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
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 2 times.
 26 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
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 2 times.
 48 if (!transform) transform = points.transform().copy();
1262 26 const Real vs = transform->voxelSize()[0];
1263 26 const float background = static_cast<float>(vs * halfband);
1264 26 const Real indexSpaceSearch = searchRadius / vs;
1265 26 const FixedBandRadius<Real> bands(radius / vs, halfband);
1266 const Real max = bands.max();
1267
1268
2/4
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
 52 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
4/6
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 11 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
 38 if (!leaf) return GridPtrVec(1, surface);
1274
1275 // max possible index space search radius
1276
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
 22 const size_t width = static_cast<size_t>(math::RoundUp(indexSpaceSearch));
1277
1278 const FixedRadius<Real> rad(radius / vs);
1279
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 22 if (interrupter) interrupter->start("Rasterizing particles to level set using constant Zhu-Bridson");
1280
1281
2/2
✓ Branch 1 taken 9 times.
✓ Branch 2 taken 2 times.
 40 GridPtrVec grids =
1282 rasterizeSurface<SdfT, AveragePositionTransfer, AttributeTypes, InterrupterT>
1283 (points, attributes, filter, *surface, interrupter,
1284 width, rad, indexSpaceSearch, points.transform(), *surface); // args
1285
1286
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 18 tools::pruneInactive(surface->tree());
1287
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 18 if (interrupter) interrupter->end();
1288
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 18 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 10 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 8 if (!transform) transform = points.transform().copy();
1314 10 const Real vs = transform->voxelSize()[0];
1315 10 const float background = static_cast<float>(vs * halfband);
1316 10 const Real indexSpaceSearch = searchRadius / vs;
1317
1318 10 RadiusT min(0), max(0);
1319 20 typename RadTreeT::Ptr mintree(new RadTreeT), maxtree(new RadTreeT);
1320 points::evalMinMax<RadiusT, UnknownCodec>
1321 10 (points.tree(), radius, min, max, filter, mintree.get(), maxtree.get());
1322
1323 // search distance at the SDF transform
1324 10 const RadiusT indexSpaceScale = RadiusT(radiusScale / vs);
1325 22 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 16 if (!leaf) return GridPtrVec(1, surface);
1334
1335 // max possible index space search radius
1336 8 const size_t width = static_cast<size_t>(math::RoundUp(indexSpaceSearch));
1337
1338 8 const size_t ridx = leaf->attributeSet().find(radius);
1339 8 if (ridx == AttributeSet::INVALID_POS) {
1340 ✗ OPENVDB_THROW(RuntimeError, "Failed to find radius attribute");
1341 }
1342
1343 VaryingRadius<RadiusT> rad(ridx, indexSpaceScale);
1344 8 if (interrupter) interrupter->start("Rasterizing particles to level set using variable Zhu-Bridson");
1345
1346 14 GridPtrVec grids =
1347 rasterizeSurface<SdfT, AveragePositionTransfer, AttributeTypes, InterrupterT>
1348 (points, attributes, filter, *surface, interrupter,
1349 width, rad, indexSpaceSearch, points.transform(), *surface); // args
1350
1351 6 tools::pruneInactive(surface->tree());
1352 6 if (interrupter) interrupter->end();
1353 6 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
1362