GCC Code Coverage Report

Directory: ./
File: openvdb/openvdb/tools/FastSweeping.h
Date: 2022-07-25 17:40:05
Exec Total Coverage
Lines: 409 461 88.7%
Functions: 63 90 70.0%
Branches: 569 1098 51.8%
Line Branch Exec Source
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 //
4 /// @file FastSweeping.h
5 ///
6 /// @author Ken Museth
7 ///
8 /// @brief Defined the six functions {fog,sdf}To{Sdf,Ext,SdfAndExt} in
9 /// addition to the two functions maskSdf and dilateSdf. Sdf denotes
10 /// a signed-distance field (i.e. negative values are inside), fog
11 /// is a scalar fog volume (i.e. higher values are inside), and Ext is
12 /// a field (of arbitrary type) that is extended off the iso-surface.
13 /// All these functions are implemented with the methods in the class
14 /// named FastSweeping.
15 ///
16 /// @note Solves the (simplified) Eikonal Eq: @f$|\nabla \phi|^2 = 1@f$ and
17 /// performs velocity extension, @f$\nabla f\nabla \phi = 0@f$, both
18 /// by means of the fast sweeping algorithm detailed in:
19 /// "A Fast Sweeping Method For Eikonal Equations"
20 /// by H. Zhao, Mathematics of Computation, Vol 74(230), pp 603-627, 2004
21 ///
22 /// @details The algorithm used below for parallel fast sweeping was first published in:
23 /// "New Algorithm for Sparse and Parallel Fast Sweeping: Efficient
24 /// Computation of Sparse Distance Fields" by K. Museth, ACM SIGGRAPH Talk,
25 /// 2017, http://www.museth.org/Ken/Publications_files/Museth_SIG17.pdf
26
27 #ifndef OPENVDB_TOOLS_FASTSWEEPING_HAS_BEEN_INCLUDED
28 #define OPENVDB_TOOLS_FASTSWEEPING_HAS_BEEN_INCLUDED
29
30 //#define BENCHMARK_FAST_SWEEPING
31
32 #include <openvdb/Platform.h>
33 #include <openvdb/math/Math.h> // for Abs() and isExactlyEqual()
34 #include <openvdb/math/Stencils.h> // for GradStencil
35 #include <openvdb/tree/LeafManager.h>
36 #include "LevelSetUtil.h"
37 #include "Morphology.h"
38 #include <openvdb/openvdb.h>
39
40 #include "Statistics.h"
41 #ifdef BENCHMARK_FAST_SWEEPING
42 #include <openvdb/util/CpuTimer.h>
43 #endif
44
45 #include <tbb/parallel_for.h>
46 #include <tbb/enumerable_thread_specific.h>
47 #include <tbb/task_group.h>
48
49 #include <type_traits>// for static_assert
50 #include <cmath>
51 #include <limits>
52 #include <deque>
53 #include <unordered_map>
54 #include <utility>// for std::make_pair
55
56 namespace openvdb {
57 OPENVDB_USE_VERSION_NAMESPACE
58 namespace OPENVDB_VERSION_NAME {
59 namespace tools {
60
61 /// @brief Fast Sweeping update mode. This is useful to determine
62 /// narrow-band extension or field extension in one side
63 /// of a signed distance field.
64 enum class FastSweepingDomain {
65 /// Update all voxels affected by the sweeping algorithm
66 SWEEP_ALL,
67 // Update voxels corresponding to an sdf/fog values that are greater than a given isovalue
68 SWEEP_GREATER_THAN_ISOVALUE,
69 // Update voxels corresponding to an sdf/fog values that are less than a given isovalue
70 SWEEP_LESS_THAN_ISOVALUE
71 };
72
73 /// @brief Converts a scalar fog volume into a signed distance function. Active input voxels
74 /// with scalar values above the given isoValue will have NEGATIVE distance
75 /// values on output, i.e. they are assumed to be INSIDE the iso-surface.
76 ///
77 /// @return A shared pointer to a signed-distance field defined on the active values
78 /// of the input fog volume.
79 ///
80 /// @param fogGrid Scalar (floating-point) volume from which an
81 /// iso-surface can be defined.
82 ///
83 /// @param isoValue A value which defines a smooth iso-surface that
84 /// intersects active voxels in @a fogGrid.
85 ///
86 /// @param nIter Number of iterations of the fast sweeping algorithm.
87 /// Each iteration performs 2^3 = 8 individual sweeps.
88 ///
89 /// @note Strictly speaking a fog volume is normalized to the range [0,1] but this
90 /// method accepts a scalar volume with an arbritary range, as long as the it
91 /// includes the @a isoValue.
92 ///
93 /// @details Topology of output grid is identical to that of the input grid, except
94 /// active tiles in the input grid will be converted to active voxels
95 /// in the output grid!
96 ///
97 /// @warning If @a isoValue does not intersect any active values in
98 /// @a fogGrid then the returned grid has all its active values set to
99 /// plus or minus infinity, depending on if the input values are larger or
100 /// smaller than @a isoValue.
101 template<typename GridT>
102 typename GridT::Ptr
103 fogToSdf(const GridT &fogGrid,
104 typename GridT::ValueType isoValue,
105 int nIter = 1);
106
107 /// @brief Given an existing approximate SDF it solves the Eikonal equation for all its
108 /// active voxels. Active input voxels with a signed distance value above the
109 /// given isoValue will have POSITIVE distance values on output, i.e. they are
110 /// assumed to be OUTSIDE the iso-surface.
111 ///
112 /// @return A shared pointer to a signed-distance field defined on the active values
113 /// of the input sdf volume.
114 ///
115 /// @param sdfGrid An approximate signed distance field to the specified iso-surface.
116 ///
117 /// @param isoValue A value which defines a smooth iso-surface that
118 /// intersects active voxels in @a sdfGrid.
119 ///
120 /// @param nIter Number of iterations of the fast sweeping algorithm.
121 /// Each iteration performs 2^3 = 8 individual sweeps.
122 ///
123 /// @note The only difference between this method and fogToSdf, defined above, is the
124 /// convention of the sign of the output distance field.
125 ///
126 /// @details Topology of output grid is identical to that of the input grid, except
127 /// active tiles in the input grid will be converted to active voxels
128 /// in the output grid!
129 ///
130 /// @warning If @a isoValue does not intersect any active values in
131 /// @a sdfGrid then the returned grid has all its active values set to
132 /// plus or minus infinity, depending on if the input values are larger or
133 /// smaller than @a isoValue.
134 template<typename GridT>
135 typename GridT::Ptr
136 sdfToSdf(const GridT &sdfGrid,
137 typename GridT::ValueType isoValue = 0,
138 int nIter = 1);
139
140 /// @brief Computes the extension of a field (scalar, vector, or int are supported), defined
141 /// by the specified functor, off an iso-surface from an input FOG volume.
142 ///
143 /// @return A shared pointer to the extension field defined from the active values in
144 /// the input fog volume.
145 ///
146 /// @param fogGrid Scalar (floating-point) volume from which an
147 /// iso-surface can be defined.
148 ///
149 /// @param op Functor with signature [](const Vec3R &xyz)->ExtValueT that
150 /// defines the Dirichlet boundary condition, on the iso-surface,
151 /// of the field to be extended.
152 ///
153 /// @param background Background value of return grid with the extension field.
154 ///
155 /// @param isoValue A value which defines a smooth iso-surface that
156 /// intersects active voxels in @a fogGrid.
157 ///
158 /// @param nIter Number of iterations of the fast sweeping algorithm.
159 /// Each iteration performs 2^3 = 8 individual sweeps.
160 ///
161 /// @param mode Determines the mode of updating the extension field. SWEEP_ALL
162 /// will update all voxels of the extension field affected by the
163 /// fast sweeping algorithm. SWEEP_GREATER_THAN_ISOVALUE will update
164 /// all voxels corresponding to fog values that are greater than a given
165 /// isovalue. SWEEP_LESS_THAN_ISOVALUE will update all voxels corresponding
166 /// to fog values that are less than a given isovalue. If a mode other
167 /// than SWEEP_ALL is chosen, a user needs to supply @a extGrid.
168 ///
169 /// @param extGrid Optional parameter required to supply a default value for the extension
170 /// field when SWEEP_GREATER_THAN_ISOVALUE or SWEEP_LESS_THAN_ISOVALUE
171 /// mode is picked for @a mode. When SWEEP_GREATER_THAN_ISOVALUE is supplied
172 /// as an argument for @a mode, the extension field voxel will default
173 /// to the value of the @a extGrid in that position if it corresponds to a fog
174 /// value that is less than the isovalue. Otherwise, the extension
175 /// field voxel value will be computed by the Fast Sweeping algorithm.
176 /// The opposite convention is implemented when SWEEP_LESS_THAN_ISOVALUE
177 /// is supplied as an argument for @a mode.
178 ///
179 /// @note Strictly speaking a fog volume is normalized to the range [0,1] but this
180 /// method accepts a scalar volume with an arbritary range, as long as the it
181 /// includes the @a isoValue.
182 ///
183 /// @details Topology of output grid is identical to that of the input grid, except
184 /// active tiles in the input grid will be converted to active voxels
185 /// in the output grid!
186 ///
187 /// @warning If @a isoValue does not intersect any active values in
188 /// @a fogGrid then the returned grid has all its active values set to
189 /// @a background.
190 template<typename FogGridT, typename ExtOpT, typename ExtValueT>
191 typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr
192 fogToExt(const FogGridT &fogGrid,
193 const ExtOpT &op,
194 const ExtValueT& background,
195 typename FogGridT::ValueType isoValue,
196 int nIter = 1,
197 FastSweepingDomain mode = FastSweepingDomain::SWEEP_ALL,
198 const typename FogGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid = nullptr);
199
200 /// @brief Computes the extension of a field (scalar, vector, or int are supported), defined
201 /// by the specified functor, off an iso-surface from an input SDF volume.
202 ///
203 /// @return A shared pointer to the extension field defined on the active values in the
204 /// input signed distance field.
205 ///
206 /// @param sdfGrid An approximate signed distance field to the specified iso-surface.
207 ///
208 /// @param op Functor with signature [](const Vec3R &xyz)->ExtValueT that
209 /// defines the Dirichlet boundary condition, on the iso-surface,
210 /// of the field to be extended.
211 ///
212 /// @param background Background value of return grid with the extension field.
213 ///
214 /// @param isoValue A value which defines a smooth iso-surface that
215 /// intersects active voxels in @a sdfGrid.
216 ///
217 /// @param nIter Number of iterations of the fast sweeping algorithm.
218 /// Each iteration performs 2^3 = 8 individual sweeps.
219 ///
220 /// @param mode Determines the mode of updating the extension field. SWEEP_ALL
221 /// will update all voxels of the extension field affected by the
222 /// fast sweeping algorithm. SWEEP_GREATER_THAN_ISOVALUE will update
223 /// all voxels corresponding to level set values that are greater than a given
224 /// isovalue. SWEEP_LESS_THAN_ISOVALUE will update all voxels corresponding
225 /// to level set values that are less than a given isovalue. If a mode other
226 /// than SWEEP_ALL is chosen, a user needs to supply @a extGrid.
227 ///
228 /// @param extGrid Optional parameter required to supply a default value for the extension
229 /// field when SWEEP_GREATER_THAN_ISOVALUE or SWEEP_LESS_THAN_ISOVALUE
230 /// mode is picked for @a mode. When SWEEP_GREATER_THAN_ISOVALUE is supplied
231 /// as an argument for @a mode, the extension field voxel will default
232 /// to the value of the @a extGrid in that position if it corresponds to a level-set
233 /// value that is less than the isovalue. Otherwise, the extension
234 /// field voxel value will be computed by the Fast Sweeping algorithm.
235 /// The opposite convention is implemented when SWEEP_LESS_THAN_ISOVALUE
236 /// is supplied as an argument for @a mode.
237 ///
238 /// @note The only difference between this method and fogToExt, defined above, is the
239 /// convention of the sign of the signed distance field.
240 ///
241 /// @details Topology of output grid is identical to that of the input grid, except
242 /// active tiles in the input grid will be converted to active voxels
243 /// in the output grid!
244 ///
245 /// @warning If @a isoValue does not intersect any active values in
246 /// @a sdfGrid then the returned grid has all its active values set to
247 /// @a background.
248 template<typename SdfGridT, typename ExtOpT, typename ExtValueT>
249 typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr
250 sdfToExt(const SdfGridT &sdfGrid,
251 const ExtOpT &op,
252 const ExtValueT &background,
253 typename SdfGridT::ValueType isoValue = 0,
254 int nIter = 1,
255 FastSweepingDomain mode = FastSweepingDomain::SWEEP_ALL,
256 const typename SdfGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid = nullptr);
257
258 /// @brief Computes the signed distance field and the extension of a field (scalar, vector, or
259 /// int are supported), defined by the specified functor, off an iso-surface from an input
260 /// FOG volume.
261 ///
262 /// @return An pair of two shared pointers to respectively the SDF and extension field
263 ///
264 /// @param fogGrid Scalar (floating-point) volume from which an
265 /// iso-surface can be defined.
266 ///
267 /// @param op Functor with signature [](const Vec3R &xyz)->ExtValueT that
268 /// defines the Dirichlet boundary condition, on the iso-surface,
269 /// of the field to be extended.
270 ///
271 /// @param background Background value of return grid with the extension field.
272 ///
273 /// @param isoValue A value which defines a smooth iso-surface that
274 /// intersects active voxels in @a fogGrid.
275 ///
276 /// @param nIter Number of iterations of the fast sweeping algorithm.
277 /// Each iteration performs 2^3 = 8 individual sweeps.
278 ///
279 /// @param mode Determines the mode of updating the extension field. SWEEP_ALL
280 /// will update all voxels of the extension field affected by the
281 /// fast sweeping algorithm. SWEEP_GREATER_THAN_ISOVALUE will update
282 /// all voxels corresponding to fog values that are greater than a given
283 /// isovalue. SWEEP_LESS_THAN_ISOVALUE will update all voxels corresponding
284 /// to fog values that are less than a given isovalue. If a mode other
285 /// than SWEEP_ALL is chosen, a user needs to supply @a extGrid.
286 ///
287 /// @param extGrid Optional parameter required to supply a default value for the extension
288 /// field when SWEEP_GREATER_THAN_ISOVALUE or SWEEP_LESS_THAN_ISOVALUE
289 /// mode is picked for @a mode. When SWEEP_GREATER_THAN_ISOVALUE is supplied
290 /// as an argument for @a mode, the extension field voxel will default
291 /// to the value of the @a extGrid in that position if it corresponds to a fog
292 /// value that is less than the isovalue. Otherwise, the extension
293 /// field voxel value will be computed by the Fast Sweeping algorithm.
294 /// The opposite convention is implemented when SWEEP_LESS_THAN_ISOVALUE
295 /// is supplied as an argument for @a mode.
296 ///
297 /// @note Strictly speaking a fog volume is normalized to the range [0,1] but this
298 /// method accepts a scalar volume with an arbritary range, as long as the it
299 /// includes the @a isoValue.
300 ///
301 /// @details Topology of output grids are identical to that of the input grid, except
302 /// active tiles in the input grid will be converted to active voxels
303 /// in the output grids!
304 ///
305 /// @warning If @a isoValue does not intersect any active values in
306 /// @a fogGrid then a pair of the following grids is returned: The first
307 /// is a signed distance grid with its active values set to plus or minus
308 /// infinity depending of whether its input values are above or below @a isoValue.
309 /// The second grid, which represents the extension field, has all its active
310 /// values set to @a background.
311 template<typename FogGridT, typename ExtOpT, typename ExtValueT>
312 std::pair<typename FogGridT::Ptr, typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr>
313 fogToSdfAndExt(const FogGridT &fogGrid,
314 const ExtOpT &op,
315 const ExtValueT &background,
316 typename FogGridT::ValueType isoValue,
317 int nIter = 1,
318 FastSweepingDomain mode = FastSweepingDomain::SWEEP_ALL,
319 const typename FogGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid = nullptr);
320
321 /// @brief Computes the signed distance field and the extension of a field (scalar, vector, or
322 /// int are supported), defined by the specified functor, off an iso-surface from an input
323 /// SDF volume.
324 ///
325 /// @return A pair of two shared pointers to respectively the SDF and extension field
326 ///
327 /// @param sdfGrid Scalar (floating-point) volume from which an
328 /// iso-surface can be defined.
329 ///
330 /// @param op Functor with signature [](const Vec3R &xyz)->ExtValueT that
331 /// defines the Dirichlet boundary condition, on the iso-surface,
332 /// of the field to be extended.
333 ///
334 /// @param background Background value of return grid with the extension field.
335 ///
336 /// @param isoValue A value which defines a smooth iso-surface that
337 /// intersects active voxels in @a sdfGrid.
338 ///
339 /// @param nIter Number of iterations of the fast sweeping algorithm.
340 /// Each iteration performs 2^3 = 8 individual sweeps.
341 ///
342 /// @param mode Determines the mode of updating the extension field. SWEEP_ALL
343 /// will update all voxels of the extension field affected by the
344 /// fast sweeping algorithm. SWEEP_GREATER_THAN_ISOVALUE will update
345 /// all voxels corresponding to level set values that are greater than a given
346 /// isovalue. SWEEP_LESS_THAN_ISOVALUE will update all voxels corresponding
347 /// to level set values that are less than a given isovalue. If a mode other
348 /// than SWEEP_ALL is chosen, a user needs to supply @a extGrid.
349 ///
350 /// @param extGrid Optional parameter required to supply a default value for the extension
351 /// field when SWEEP_GREATER_THAN_ISOVALUE or SWEEP_LESS_THAN_ISOVALUE
352 /// mode is picked for @a mode. When SWEEP_GREATER_THAN_ISOVALUE is supplied
353 /// as an argument for @a mode, the extension field voxel will default
354 /// to the value of the @a extGrid in that position if it corresponds to a level-set
355 /// value that is less than the isovalue. Otherwise, the extension
356 /// field voxel value will be computed by the Fast Sweeping algorithm.
357 /// The opposite convention is implemented when SWEEP_LESS_THAN_ISOVALUE
358 /// is supplied as an argument for @a mode.
359 ///
360 /// @note Strictly speaking a fog volume is normalized to the range [0,1] but this
361 /// method accepts a scalar volume with an arbritary range, as long as the it
362 /// includes the @a isoValue.
363 ///
364 /// @details Topology of output grids are identical to that of the input grid, except
365 /// active tiles in the input grid will be converted to active voxels
366 /// in the output grids!
367 ///
368 /// @warning If @a isoValue does not intersect any active values in
369 /// @a sdfGrid then a pair of the following grids is returned: The first
370 /// is a signed distance grid with its active values set to plus or minus
371 /// infinity depending of whether its input values are above or below @a isoValue.
372 /// The second grid, which represents the extension field, has all its active
373 /// values set to @a background.
374 template<typename SdfGridT, typename ExtOpT, typename ExtValueT>
375 std::pair<typename SdfGridT::Ptr, typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr>
376 sdfToSdfAndExt(const SdfGridT &sdfGrid,
377 const ExtOpT &op,
378 const ExtValueT &background,
379 typename SdfGridT::ValueType isoValue = 0,
380 int nIter = 1,
381 FastSweepingDomain mode = FastSweepingDomain::SWEEP_ALL,
382 const typename SdfGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid = nullptr);
383
384 /// @brief Dilates the narrow band of an existing signed distance field by
385 /// a specified number of voxels (like adding "onion-rings").
386 ///
387 /// @note This operation is not to be confused with morphological dilation
388 /// of a level set, which is implemented in LevelSetFilter::offset,
389 /// and involves actual interface tracking of the narrow band.
390 ///
391 /// @return A shared pointer to the dilated signed distance field.
392 ///
393 /// @param sdfGrid Input signed distance field to be dilated.
394 ///
395 /// @param dilation Numer of voxels that the narrow band of the input SDF will be dilated.
396 ///
397 /// @param nn Stencil-pattern used for dilation
398 ///
399 /// @param nIter Number of iterations of the fast sweeping algorithm.
400 /// Each iteration performs 2^3 = 8 individual sweeps.
401 ///
402 /// @param mode Determines the direction of the dilation. SWEEP_ALL
403 /// will dilate in both sides of the signed distance function,
404 /// SWEEP_GREATER_THAN_ISOVALUE will dilate in the positive
405 /// side of the iso-surface, SWEEP_LESS_THAN_ISOVALUE will dilate
406 /// in the negative side of the iso-surface.
407 ///
408 /// @details Topology will change as a result of this dilation. E.g. if
409 /// sdfGrid has a width of 3 and @a dilation = 6 then the grid
410 /// returned by this method is a narrow band signed distance field
411 /// with a total width of 9 units.
412 template<typename GridT>
413 typename GridT::Ptr
414 dilateSdf(const GridT &sdfGrid,
415 int dilation,
416 NearestNeighbors nn = NN_FACE,
417 int nIter = 1,
418 FastSweepingDomain mode = FastSweepingDomain::SWEEP_ALL);
419
420 /// @brief Fills mask by extending an existing signed distance field into
421 /// the active values of this input ree of arbitrary value type.
422 ///
423 /// @return A shared pointer to the masked signed distance field.
424 ///
425 /// @param sdfGrid Input signed distance field to be extended into the mask.
426 ///
427 /// @param mask Mask used to identify the topology of the output SDF.
428 /// Note this mask is assume to overlap with the sdfGrid.
429 ///
430 /// @param ignoreActiveTiles If false, active tiles in the mask are treated
431 /// as active voxels. Else they are ignored.
432 ///
433 /// @param nIter Number of iterations of the fast sweeping algorithm.
434 /// Each iteration performs 2^3 = 8 individual sweeps.
435 ///
436 /// @details Topology of the output SDF is determined by the union of the active
437 /// voxels (or optionally values) in @a sdfGrid and @a mask.
438 template<typename GridT, typename MaskTreeT>
439 typename GridT::Ptr
440 maskSdf(const GridT &sdfGrid,
441 const Grid<MaskTreeT> &mask,
442 bool ignoreActiveTiles = false,
443 int nIter = 1);
444
445 ////////////////////////////////////////////////////////////////////////////////
446 /// @brief Computes signed distance values from an initial iso-surface and
447 /// optionally performs velocity extension at the same time. This is
448 /// done by means of a novel sparse and parallel fast sweeping
449 /// algorithm based on a first order Godunov's scheme.
450 ///
451 /// Solves: @f$|\nabla \phi|^2 = 1 @f$
452 ///
453 /// @warning Note, it is important to call one of the initialization methods before
454 /// called the sweep function. Failure to do so will throw a RuntimeError.
455 /// Consider instead call one of the many higher-level free-standing functions
456 /// defined above!
457 template<typename SdfGridT, typename ExtValueT = typename SdfGridT::ValueType>
458 class FastSweeping
459 {
460 static_assert(std::is_floating_point<typename SdfGridT::ValueType>::value,
461 "FastSweeping requires SdfGridT to have floating-point values");
462 // Defined types related to the signed distance (or fog) grid
463 using SdfValueT = typename SdfGridT::ValueType;
464 using SdfTreeT = typename SdfGridT::TreeType;
465 using SdfAccT = tree::ValueAccessor<SdfTreeT, false>;//don't register accessors
466 using SdfConstAccT = typename tree::ValueAccessor<const SdfTreeT, false>;//don't register accessors
467
468 // define types related to the extension field
469 using ExtGridT = typename SdfGridT::template ValueConverter<ExtValueT>::Type;
470 using ExtTreeT = typename ExtGridT::TreeType;
471 using ExtAccT = tree::ValueAccessor<ExtTreeT, false>;
472
473 // define types related to the tree that masks out the active voxels to be solved for
474 using SweepMaskTreeT = typename SdfTreeT::template ValueConverter<ValueMask>::Type;
475 using SweepMaskAccT = tree::ValueAccessor<SweepMaskTreeT, false>;//don't register accessors
476
477 public:
478
479 /// @brief Constructor
480 FastSweeping();
481
482 /// @brief Destructor.
483 32 ~FastSweeping() { this->clear(); }
484
485 /// @brief Disallow copy construction.
486 FastSweeping(const FastSweeping&) = delete;
487
488 /// @brief Disallow copy assignment.
489 FastSweeping& operator=(const FastSweeping&) = delete;
490
491 /// @brief Returns a shared pointer to the signed distance field computed
492 /// by this class.
493 ///
494 /// @warning This shared pointer might point to NULL if the grid has not been
495 /// initialize (by one of the init methods) or computed (by the sweep
496 /// method).
497 typename SdfGridT::Ptr sdfGrid() { return mSdfGrid; }
498
499 /// @brief Returns a shared pointer to the extension field computed
500 /// by this class.
501 ///
502 /// @warning This shared pointer might point to NULL if the grid has not been
503 /// initialize (by one of the init methods) or computed (by the sweep
504 /// method).
505 typename ExtGridT::Ptr extGrid() { return mExtGrid; }
506
507 /// @brief Returns a shared pointer to the extension grid input. This is non-NULL
508 /// if this class is used to extend a field with a non-default sweep direction.
509 ///
510 /// @warning This shared pointer might point to NULL. This is non-NULL
511 /// if this class is used to extend a field with a non-default sweep direction,
512 /// i.e. SWEEP_LESS_THAN_ISOVALUE or SWEEP_GREATER_THAN_ISOVALUE.
513 typename ExtGridT::Ptr extGridInput() { return mExtGridInput; }
514
515 /// @brief Initializer for input grids that are either a signed distance
516 /// field or a scalar fog volume.
517 ///
518 /// @return True if the initialization succeeded.
519 ///
520 /// @param sdfGrid Input scalar grid that represents an existing signed distance
521 /// field or a fog volume (signified by @a isInputSdf).
522 ///
523 /// @param isoValue Iso-value to be used to define the Dirichlet boundary condition
524 /// of the fast sweeping algorithm (typically 0 for sdfs and a
525 /// positive value for fog volumes).
526 ///
527 /// @param isInputSdf Used to determine if @a sdfGrid is a sigend distance field (true)
528 /// or a scalar fog volume (false).
529 ///
530 /// @details This, or any of ther other initialization methods, should be called
531 /// before any call to sweep(). Failure to do so will throw a RuntimeError.
532 ///
533 /// @warning Note, if this method fails, i.e. returns false, a subsequent call
534 /// to sweep will trow a RuntimeError. Instead call clear and try again.
535 bool initSdf(const SdfGridT &sdfGrid, SdfValueT isoValue, bool isInputSdf);
536
537 /// @brief Initializer used whenever velocity extension is performed in addition
538 /// to the computation of signed distance fields.
539 ///
540 /// @return True if the initialization succeeded.
541 ///
542 ///
543 /// @param sdfGrid Input scalar grid that represents an existing signed distance
544 /// field or a fog volume (signified by @a isInputSdf).
545 ///
546 /// @param op Functor with signature [](const Vec3R &xyz)->ExtValueT that
547 /// defines the Dirichlet boundary condition, on the iso-surface,
548 /// of the field to be extended. Strictly the return type of this functor
549 /// is only required to be convertible to ExtValueT!
550 ///
551 /// @param background Background value of return grid with the extension field.
552 ///
553 /// @param isoValue Iso-value to be used for the boundary condition of the fast
554 /// sweeping algorithm (typically 0 for sdfs and a positive value
555 /// for fog volumes).
556 ///
557 /// @param isInputSdf Used to determine if @a sdfGrid is a sigend distance field (true)
558 /// or a scalar fog volume (false).
559 ///
560 /// @param mode Determines the mode of updating the extension field. SWEEP_ALL
561 /// will update all voxels of the extension field affected by the
562 /// fast sweeping algorithm. SWEEP_GREATER_THAN_ISOVALUE will update
563 /// all voxels corresponding to fog values that are greater than a given
564 /// isovalue. SWEEP_LESS_THAN_ISOVALUE will update all voxels corresponding
565 /// to fog values that are less than a given isovalue. If a mode other
566 /// than SWEEP_ALL is chosen, a user needs to supply @a extGrid.
567 ///
568 /// @param extGrid Optional parameter required to supply a default value for the extension
569 /// field when SWEEP_GREATER_THAN_ISOVALUE or SWEEP_LESS_THAN_ISOVALUE
570 /// mode is picked for @a mode. When SWEEP_GREATER_THAN_ISOVALUE is supplied
571 /// as an argument for @a mode, the extension field voxel will default
572 /// to the value of the @a extGrid in that position if it corresponds to a level-set
573 /// value that is less than the isovalue. Otherwise, the extension
574 /// field voxel value will be computed by the Fast Sweeping algorithm.
575 /// The opposite convention is implemented when SWEEP_LESS_THAN_ISOVALUE
576 /// is supplied as an argument for @a mode.
577 ///
578 /// @details This, or any of ther other initialization methods, should be called
579 /// before any call to sweep(). Failure to do so will throw a RuntimeError.
580 ///
581 /// @warning Note, if this method fails, i.e. returns false, a subsequent call
582 /// to sweep will trow a RuntimeError. Instead call clear and try again.
583 template <typename ExtOpT>
584 bool initExt(const SdfGridT &sdfGrid,
585 const ExtOpT &op,
586 const ExtValueT &background,
587 SdfValueT isoValue,
588 bool isInputSdf,
589 FastSweepingDomain mode = FastSweepingDomain::SWEEP_ALL,
590 const typename ExtGridT::ConstPtr extGrid = nullptr);
591
592 /// @brief Initializer used when dilating an existing signed distance field.
593 ///
594 /// @return True if the initialization succeeded.
595 ///
596 /// @param sdfGrid Input signed distance field to to be dilated.
597 ///
598 /// @param dilation Numer of voxels that the input SDF will be dilated.
599 ///
600 /// @param nn Stencil-pattern used for dilation
601 ///
602 /// @param mode Determines the direction of the dilation. SWEEP_ALL
603 /// will dilate in both sides of the signed distance function,
604 /// SWEEP_GREATER_THAN_ISOVALUE will dilate in the positive
605 /// side of the iso-surface, SWEEP_LESS_THAN_ISOVALUE will dilate
606 /// in the negative side of the iso-surface.
607 ///
608 /// @details This, or any of ther other initialization methods, should be called
609 /// before any call to sweep(). Failure to do so will throw a RuntimeError.
610 ///
611 /// @warning Note, if this method fails, i.e. returns false, a subsequent call
612 /// to sweep will trow a RuntimeError. Instead call clear and try again.
613 bool initDilate(const SdfGridT &sdfGrid,
614 int dilation,
615 NearestNeighbors nn = NN_FACE,
616 FastSweepingDomain mode = FastSweepingDomain::SWEEP_ALL);
617
618 /// @brief Initializer used for the extension of an existing signed distance field
619 /// into the active values of an input mask of arbitrary value type.
620 ///
621 /// @return True if the initialization succeeded.
622 ///
623 /// @param sdfGrid Input signed distance field to be extended into the mask.
624 ///
625 /// @param mask Mask used to identify the topology of the output SDF.
626 /// Note this mask is assume to overlap with the sdfGrid.
627 ///
628 /// @param ignoreActiveTiles If false, active tiles in the mask are treated
629 /// as active voxels. Else they are ignored.
630 ///
631 /// @details This, or any of ther other initialization methods, should be called
632 /// before any call to sweep(). Failure to do so will throw a RuntimeError.
633 ///
634 /// @warning Note, if this method fails, i.e. returns false, a subsequent call
635 /// to sweep will trow a RuntimeError. Instead call clear and try again.
636 template<typename MaskTreeT>
637 bool initMask(const SdfGridT &sdfGrid, const Grid<MaskTreeT> &mask, bool ignoreActiveTiles = false);
638
639 /// @brief Perform @a nIter iterations of the fast sweeping algorithm.
640 ///
641 /// @param nIter Number of iterations of the fast sweeping algorithm.
642 /// Each iteration performs 2^3 = 8 individual sweeps.
643 ///
644 /// @param finalize If true the (possibly asymmetric) inside and outside values of the
645 /// resulting signed distance field are properly set. Unless you're
646 /// an expert this should remain true!
647 ///
648 /// @throw RuntimeError if sweepingVoxelCount() or boundaryVoxelCount() return zero.
649 /// This might happen if none of the initialization methods above were called
650 /// or if that initialization failed.
651 void sweep(int nIter = 1,
652 bool finalize = true);
653
654 /// @brief Clears all the grids and counters so initialization can be called again.
655 void clear();
656
657 /// @brief Return the number of voxels that will be solved for.
658
5/10
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 times.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
 5 size_t sweepingVoxelCount() const { return mSweepingVoxelCount; }
659
660 /// @brief Return the number of voxels that defined the boundary condition.
661
4/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
 4 size_t boundaryVoxelCount() const { return mBoundaryVoxelCount; }
662
663 /// @brief Return true if there are voxels and boundaries to solve for
664
14/28
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 1 times.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 1 times.
✓ Branch 20 taken 1 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✓ Branch 23 taken 1 times.
✓ Branch 24 taken 1 times.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✓ Branch 27 taken 1 times.
 8 bool isValid() const { return mSweepingVoxelCount > 0 && mBoundaryVoxelCount > 0; }
665
666 /// @brief Return whether the sweep update is in all direction (SWEEP_ALL),
667 /// greater than isovalue (SWEEP_GREATER_THAN_ISOVALUE), or less than isovalue
668 /// (SWEEP_LESS_THAN_ISOVALUE).
669 ///
670 /// @note SWEEP_GREATER_THAN_ISOVALUE and SWEEP_LESS_THAN_ISOVALUE modes are used
671 /// in dilating the narrow-band of a levelset or in extending a field.
672 FastSweepingDomain sweepDirection() const { return mSweepDirection; }
673
674 /// @brief Return whether the fast-sweeping input grid a signed distance function or not (fog).
675 bool isInputSdf() { return mIsInputSdf; }
676
677 private:
678
679 /// @brief Private method to prune the sweep mask and cache leaf origins.
680 void computeSweepMaskLeafOrigins();
681
682 // Private utility classes
683 template<typename>
684 struct MaskKernel;// initialization to extend a SDF into a mask
685 template<typename>
686 struct InitExt;
687 struct InitSdf;
688 struct DilateKernel;// initialization to dilate a SDF
689 struct MinMaxKernel;
690 struct SweepingKernel;// performs the actual concurrent sparse fast sweeping
691
692 // Define the topology (i.e. stencil) of the neighboring grid points
693 static const Coord mOffset[6];// = {{-1,0,0},{1,0,0},{0,-1,0},{0,1,0},{0,0,-1},{0,0,1}};
694
695 // Private member data of FastSweeping
696 typename SdfGridT::Ptr mSdfGrid;
697 typename ExtGridT::Ptr mExtGrid;
698 typename ExtGridT::Ptr mExtGridInput; // optional: only used in extending a field in one direction
699 SweepMaskTreeT mSweepMask; // mask tree containing all non-boundary active voxels, in the case of dilation, does not include active voxel
700 std::vector<Coord> mSweepMaskLeafOrigins; // cache of leaf node origins for mask tree
701 size_t mSweepingVoxelCount, mBoundaryVoxelCount;
702 FastSweepingDomain mSweepDirection; // only used in dilate and extending a field
703 bool mIsInputSdf;
704 };// FastSweeping
705
706 ////////////////////////////////////////////////////////////////////////////////
707
708 // Static member data initialization
709 template <typename SdfGridT, typename ExtValueT>
710 const Coord FastSweeping<SdfGridT, ExtValueT>::mOffset[6] = {{-1,0,0},{1,0,0},
711 {0,-1,0},{0,1,0},
712 {0,0,-1},{0,0,1}};
713
714 template <typename SdfGridT, typename ExtValueT>
715 16 FastSweeping<SdfGridT, ExtValueT>::FastSweeping()
716 16 : mSdfGrid(nullptr), mExtGrid(nullptr), mSweepingVoxelCount(0), mBoundaryVoxelCount(0), mSweepDirection(FastSweepingDomain::SWEEP_ALL), mIsInputSdf(true)
717 {
718 16 }
719
720 template <typename SdfGridT, typename ExtValueT>
721
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 9 times.
 34 void FastSweeping<SdfGridT, ExtValueT>::clear()
722 {
723 mSdfGrid.reset();
724 mExtGrid.reset();
725 34 mSweepMask.clear();
726
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 17 times.
 34 if (mExtGridInput) mExtGridInput.reset();
727 34 mSweepingVoxelCount = mBoundaryVoxelCount = 0;
728 34 mSweepDirection = FastSweepingDomain::SWEEP_ALL;
729 34 mIsInputSdf = true;
730 34 }
731
732 template <typename SdfGridT, typename ExtValueT>
733 16 void FastSweeping<SdfGridT, ExtValueT>::computeSweepMaskLeafOrigins()
734 {
735 // replace any inactive leaf nodes with tiles and voxelize any active tiles
736
737 16 pruneInactive(mSweepMask);
738 mSweepMask.voxelizeActiveTiles();
739
740 using LeafManagerT = tree::LeafManager<SweepMaskTreeT>;
741 using LeafT = typename SweepMaskTreeT::LeafNodeType;
742 32 LeafManagerT leafManager(mSweepMask);
743
744
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 16 mSweepMaskLeafOrigins.resize(leafManager.leafCount());
745 16 std::atomic<size_t> sweepingVoxelCount{0};
746 16871 auto kernel = [&](const LeafT& leaf, size_t leafIdx) {
747 16855 mSweepMaskLeafOrigins[leafIdx] = leaf.origin();
748 sweepingVoxelCount += leaf.onVoxelCount();
749 };
750
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 16 leafManager.foreach(kernel, /*threaded=*/true, /*grainsize=*/1024);
751
752
1/2
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
 16 mBoundaryVoxelCount = 0;
753
1/2
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
 16 mSweepingVoxelCount = sweepingVoxelCount;
754
1/2
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
 16 if (mSdfGrid) {
755
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 16 const size_t totalCount = mSdfGrid->constTree().activeVoxelCount();
756
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
 16 assert( totalCount >= mSweepingVoxelCount );
757 16 mBoundaryVoxelCount = totalCount - mSweepingVoxelCount;
758 }
759 16 }// FastSweeping::computeSweepMaskLeafOrigins
760
761 template <typename SdfGridT, typename ExtValueT>
762 1 bool FastSweeping<SdfGridT, ExtValueT>::initSdf(const SdfGridT &fogGrid, SdfValueT isoValue, bool isInputSdf)
763 {
764 1 this->clear();
765 1 mSdfGrid = fogGrid.deepCopy();//very fast
766 1 mIsInputSdf = isInputSdf;
767 InitSdf kernel(*this);
768 1 kernel.run(isoValue);
769 1 return this->isValid();
770 }
771
772 template <typename SdfGridT, typename ExtValueT>
773 template <typename OpT>
774 6 bool FastSweeping<SdfGridT, ExtValueT>::initExt(const SdfGridT &fogGrid, const OpT &op, const ExtValueT &background, SdfValueT isoValue, bool isInputSdf, FastSweepingDomain mode, const typename ExtGridT::ConstPtr extGrid)
775 {
776
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 6 if (mode != FastSweepingDomain::SWEEP_ALL) {
777 ✗ if (!extGrid)
778 ✗ OPENVDB_THROW(RuntimeError, "FastSweeping::initExt Calling initExt with mode != SWEEP_ALL requires an extension grid!");
779 ✗ if (extGrid->transform() != fogGrid.transform())
780 ✗ OPENVDB_THROW(RuntimeError, "FastSweeping::initExt extension grid input should have the same transform as Fog/SDF grid!");
781 }
782
783 6 this->clear();
784 12 mSdfGrid = fogGrid.deepCopy();//very fast
785 6 mExtGrid = createGrid<ExtGridT>( background );
786 6 mSweepDirection = mode;
787 6 mIsInputSdf = isInputSdf;
788
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 6 if (mSweepDirection != FastSweepingDomain::SWEEP_ALL) {
789 ✗ mExtGridInput = extGrid->deepCopy();
790 }
791 mExtGrid->topologyUnion( *mSdfGrid );//very fast
792 InitExt<OpT> kernel(*this);
793 6 kernel.run(isoValue, op);
794 6 return this->isValid();
795 }
796
797
798 template <typename SdfGridT, typename ExtValueT>
799 1 bool FastSweeping<SdfGridT, ExtValueT>::initDilate(const SdfGridT &sdfGrid, int dilate, NearestNeighbors nn, FastSweepingDomain mode)
800 {
801 1 this->clear();
802 1 mSdfGrid = sdfGrid.deepCopy();//very fast
803 1 mSweepDirection = mode;
804 1 DilateKernel kernel(*this);
805
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 kernel.run(dilate, nn);
806 1 return this->isValid();
807 }
808
809 template <typename SdfGridT, typename ExtValueT>
810 template<typename MaskTreeT>
811 6 bool FastSweeping<SdfGridT, ExtValueT>::initMask(const SdfGridT &sdfGrid, const Grid<MaskTreeT> &mask, bool ignoreActiveTiles)
812 {
813 6 this->clear();
814 12 mSdfGrid = sdfGrid.deepCopy();//very fast
815
816
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 6 if (mSdfGrid->transform() != mask.transform()) {
817 ✗ OPENVDB_THROW(RuntimeError, "FastSweeping: Mask not aligned with the grid!");
818 }
819
820
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 1 times.
 6 if (mask.getGridClass() == GRID_LEVEL_SET) {
821 using T = typename MaskTreeT::template ValueConverter<bool>::Type;
822 4 typename Grid<T>::Ptr tmp = sdfInteriorMask(mask);//might have active tiles
823 tmp->tree().voxelizeActiveTiles();//multi-threaded
824 MaskKernel<T> kernel(*this);
825
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 4 kernel.run(tmp->tree());//multi-threaded
826 } else {
827
2/4
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
 4 if (ignoreActiveTiles || !mask.tree().hasActiveTiles()) {
828 MaskKernel<MaskTreeT> kernel(*this);
829 ✗ kernel.run(mask.tree());//multi-threaded
830 } else {
831 using T = typename MaskTreeT::template ValueConverter<ValueMask>::Type;
832
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 4 T tmp(mask.tree(), false, TopologyCopy());//multi-threaded
833 tmp.voxelizeActiveTiles(true);//multi-threaded
834 MaskKernel<T> kernel(*this);
835
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 2 kernel.run(tmp);//multi-threaded
836 }
837 }
838 6 return this->isValid();
839 }// FastSweeping::initMask
840
841 template <typename SdfGridT, typename ExtValueT>
842
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
 16 void FastSweeping<SdfGridT, ExtValueT>::sweep(int nIter, bool finalize)
843 {
844
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
 16 if (!mSdfGrid) {
845 ✗ OPENVDB_THROW(RuntimeError, "FastSweeping::sweep called before initialization!");
846 }
847
3/6
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 16 if (mExtGrid && mSweepDirection != FastSweepingDomain::SWEEP_ALL && !mExtGridInput) {
848 ✗ OPENVDB_THROW(RuntimeError, "FastSweeping: Trying to extend a field in one direction needs"
849 " a non-null reference extension grid input.");
850 }
851
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
 16 if (this->boundaryVoxelCount() == 0) {
852 ✗ OPENVDB_THROW(RuntimeError, "FastSweeping: No boundary voxels found!");
853
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
 16 } else if (this->sweepingVoxelCount() == 0) {
854 ✗ OPENVDB_THROW(RuntimeError, "FastSweeping: No computing voxels found!");
855 }
856
857 // note: Sweeping kernel is non copy-constructible, so use a deque instead of a vector
858 16 std::deque<SweepingKernel> kernels;
859
3/4
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 8 times.
✓ Branch 3 taken 32 times.
✗ Branch 4 not taken.
 80 for (int i = 0; i < 4; i++) kernels.emplace_back(*this);
860
861 { // compute voxel slices
862 #ifdef BENCHMARK_FAST_SWEEPING
863 util::CpuTimer timer("Computing voxel slices");
864 #endif
865
866 // Exploiting nested parallelism - all voxel slice data is precomputed
867 tbb::task_group tasks;
868
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 81050199 tasks.run([&] { kernels[0].computeVoxelSlices([](const Coord &a){ return a[0]+a[1]+a[2]; });/*+++ & ---*/ });
869
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 81113294 tasks.run([&] { kernels[1].computeVoxelSlices([](const Coord &a){ return a[0]+a[1]-a[2]; });/*++- & --+*/ });
870
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 81112321 tasks.run([&] { kernels[2].computeVoxelSlices([](const Coord &a){ return a[0]-a[1]+a[2]; });/*+-+ & -+-*/ });
871
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 81102069 tasks.run([&] { kernels[3].computeVoxelSlices([](const Coord &a){ return a[0]-a[1]-a[2]; });/*+-- & -++*/ });
872
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 16 tasks.wait();
873
874 #ifdef BENCHMARK_FAST_SWEEPING
875 timer.stop();
876 #endif
877 }
878
879 // perform nIter iterations of bi-directional sweeping in all directions
880
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
 32 for (int i = 0; i < nIter; ++i) {
881
3/4
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 8 times.
✓ Branch 3 taken 32 times.
✗ Branch 4 not taken.
 80 for (SweepingKernel& kernel : kernels) kernel.sweep();
882 }
883
884
1/2
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
 16 if (finalize) {
885 #ifdef BENCHMARK_FAST_SWEEPING
886 util::CpuTimer timer("Computing extrema values");
887 #endif
888 MinMaxKernel kernel;
889
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 16 auto e = kernel.run(*mSdfGrid);//multi-threaded
890 //auto e = extrema(mGrid->beginValueOn());// 100x slower!!!!
891 #ifdef BENCHMARK_FAST_SWEEPING
892 std::cerr << "Min = " << e.min() << " Max = " << e.max() << std::endl;
893 timer.restart("Changing asymmetric background value");
894 #endif
895
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 16 changeAsymmetricLevelSetBackground(mSdfGrid->tree(), e.max(), e.min());//multi-threaded
896
897 #ifdef BENCHMARK_FAST_SWEEPING
898 timer.stop();
899 #endif
900 }
901 16 }// FastSweeping::sweep
902
903 /// Private class of FastSweeping to quickly compute the extrema
904 /// values of the active voxels in the leaf nodes. Several orders
905 /// of magnitude faster than tools::extrema!
906 template <typename SdfGridT, typename ExtValueT>
907 struct FastSweeping<SdfGridT, ExtValueT>::MinMaxKernel
908 {
909 using LeafMgr = tree::LeafManager<const SdfTreeT>;
910 using LeafRange = typename LeafMgr::LeafRange;
911
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
 8 MinMaxKernel() : mMin(std::numeric_limits<SdfValueT>::max()), mMax(-mMin) {}
912 31 MinMaxKernel(MinMaxKernel& other, tbb::split) : mMin(other.mMin), mMax(other.mMax) {}
913
914 16 math::MinMax<SdfValueT> run(const SdfGridT &grid)
915 {
916 16 LeafMgr mgr(grid.tree());// super fast
917
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 16 tbb::parallel_reduce(mgr.leafRange(), *this);
918 16 return math::MinMax<SdfValueT>(mMin, mMax);
919 }
920
921 2228 void operator()(const LeafRange& r)
922 {
923
2/2
✓ Branch 1 taken 20124 times.
✓ Branch 2 taken 1114 times.
 42476 for (auto leafIter = r.begin(); leafIter; ++leafIter) {
924
2/2
✓ Branch 0 taken 5505685 times.
✓ Branch 1 taken 20124 times.
 11051618 for (auto voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
925 11011370 const SdfValueT v = *voxelIter;
926
2/2
✓ Branch 0 taken 2838 times.
✓ Branch 1 taken 5502847 times.
 11011370 if (v < mMin) mMin = v;
927
2/2
✓ Branch 0 taken 264 times.
✓ Branch 1 taken 5505421 times.
 11011370 if (v > mMax) mMax = v;
928 }
929 }
930 2228 }
931
932 void join(const MinMaxKernel& other)
933 {
934
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 21 times.
 31 if (other.mMin < mMin) mMin = other.mMin;
935
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 22 times.
 31 if (other.mMax > mMax) mMax = other.mMax;
936 }
937
938 SdfValueT mMin, mMax;
939 };// FastSweeping::MinMaxKernel
940
941 ////////////////////////////////////////////////////////////////////////////////
942
943 /// Private class of FastSweeping to perform multi-threaded initialization
944 template <typename SdfGridT, typename ExtValueT>
945 struct FastSweeping<SdfGridT, ExtValueT>::DilateKernel
946 {
947 using LeafRange = typename tree::LeafManager<SdfTreeT>::LeafRange;
948 1 DilateKernel(FastSweeping &parent)
949 : mParent(&parent),
950
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 mBackground(parent.mSdfGrid->background())
951 {
952
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 mSdfGridInput = mParent->mSdfGrid->deepCopy();
953 1 }
954 DilateKernel(const DilateKernel &parent) = default;// for tbb::parallel_for
955 DilateKernel& operator=(const DilateKernel&) = delete;
956
957 1 void run(int dilation, NearestNeighbors nn)
958 {
959 #ifdef BENCHMARK_FAST_SWEEPING
960 util::CpuTimer timer("Construct LeafManager");
961 #endif
962 2 tree::LeafManager<SdfTreeT> mgr(mParent->mSdfGrid->tree());// super fast
963
964 #ifdef BENCHMARK_FAST_SWEEPING
965 timer.restart("Changing background value");
966 #endif
967 static const SdfValueT Unknown = std::numeric_limits<SdfValueT>::max();
968
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 changeLevelSetBackground(mgr, Unknown);//multi-threaded
969
970 #ifdef BENCHMARK_FAST_SWEEPING
971 timer.restart("Dilating and updating mgr (parallel)");
972 //timer.restart("Dilating and updating mgr (serial)");
973 #endif
974
975 const int delta = 5;
976
3/4
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 2 for (int i=0, d = dilation/delta; i<d; ++i) dilateActiveValues(mgr, delta, nn, IGNORE_TILES);
977
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 dilateActiveValues(mgr, dilation % delta, nn, IGNORE_TILES);
978 //for (int i=0, n=5, d=dilation/n; i<d; ++i) dilateActiveValues(mgr, n, nn, IGNORE_TILES);
979 //dilateVoxels(mgr, dilation, nn);
980
981 #ifdef BENCHMARK_FAST_SWEEPING
982 timer.restart("Initializing grid and sweep mask");
983 #endif
984
985
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 mParent->mSweepMask.clear();
986
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
987
988 using LeafManagerT = tree::LeafManager<typename SdfGridT::TreeType>;
989 using LeafT = typename SdfGridT::TreeType::LeafNodeType;
990
991 1 const FastSweepingDomain mode = mParent->mSweepDirection;
992
993
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 2 LeafManagerT leafManager(mParent->mSdfGrid->tree());
994
995 5636 auto kernel = [&](LeafT& leaf, size_t /*leafIdx*/) {
996 static const SdfValueT Unknown = std::numeric_limits<SdfValueT>::max();
997 2818 const SdfValueT background = mBackground;//local copy
998 2818 auto* maskLeaf = mParent->mSweepMask.probeLeaf(leaf.origin());
999 SdfConstAccT sdfInputAcc(mSdfGridInput->tree());
1000
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2818 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 2818 assert(maskLeaf);
1001
2/4
✓ Branch 0 taken 953186 times.
✓ Branch 1 taken 2818 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 956004 for (auto voxelIter = leaf.beginValueOn(); voxelIter; ++voxelIter) {
1002 953186 const SdfValueT value = *voxelIter;
1003 SdfValueT inputValue;
1004
1/4
✓ Branch 1 taken 953186 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 953186 const Coord ijk = voxelIter.getCoord();
1005
1006
2/4
✓ Branch 0 taken 270638 times.
✓ Branch 1 taken 682548 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 953186 if (math::Abs(value) < background) {// disable boundary voxels from the mask tree
1007 270638 maskLeaf->setValueOff(voxelIter.pos());
1008 } else {
1009
1/8
✓ Branch 0 taken 682548 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 682548 switch (mode) {
1010 682548 case FastSweepingDomain::SWEEP_ALL :
1011
3/8
✓ Branch 0 taken 273402 times.
✓ Branch 1 taken 409146 times.
✓ Branch 3 taken 682548 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 955950 voxelIter.setValue(value > 0 ? Unknown : -Unknown);
1012 682548 break;
1013 ✗ case FastSweepingDomain::SWEEP_GREATER_THAN_ISOVALUE :
1014 ✗ if (value > 0) voxelIter.setValue(Unknown);
1015 else {
1016 ✗ maskLeaf->setValueOff(voxelIter.pos());
1017 ✗ bool isInputOn = sdfInputAcc.probeValue(ijk, inputValue);
1018 ✗ if ( !isInputOn ) voxelIter.setValueOff();
1019 ✗ else voxelIter.setValue(inputValue);
1020 }
1021 break;
1022 ✗ case FastSweepingDomain::SWEEP_LESS_THAN_ISOVALUE :
1023 ✗ if (value < 0) voxelIter.setValue(-Unknown);
1024 else {
1025 ✗ maskLeaf->setValueOff(voxelIter.pos());
1026 ✗ bool isInputOn = sdfInputAcc.probeValue(ijk, inputValue);
1027 ✗ if ( !isInputOn ) voxelIter.setValueOff();
1028 ✗ else voxelIter.setValue(inputValue);
1029 }
1030 break;
1031 }
1032 }
1033 }
1034 };
1035
1036
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 leafManager.foreach( kernel );
1037
1038 // cache the leaf node origins for fast lookup in the sweeping kernels
1039
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 mParent->computeSweepMaskLeafOrigins();
1040
1041 #ifdef BENCHMARK_FAST_SWEEPING
1042 timer.stop();
1043 #endif
1044 1 }// FastSweeping::DilateKernel::run
1045
1046 // Private member data of DilateKernel
1047 FastSweeping *mParent;
1048 const SdfValueT mBackground;
1049 typename SdfGridT::ConstPtr mSdfGridInput;
1050 };// FastSweeping::DilateKernel
1051
1052 ////////////////////////////////////////////////////////////////////////////////
1053
1054 template <typename SdfGridT, typename ExtValueT>
1055 struct FastSweeping<SdfGridT, ExtValueT>::InitSdf
1056 {
1057 using LeafRange = typename tree::LeafManager<SdfTreeT>::LeafRange;
1058 1 InitSdf(FastSweeping &parent): mParent(&parent),
1059 1 mSdfGrid(parent.mSdfGrid.get()), mIsoValue(0), mAboveSign(0) {}
1060 InitSdf(const InitSdf&) = default;// for tbb::parallel_for
1061 InitSdf& operator=(const InitSdf&) = delete;
1062
1063 1 void run(SdfValueT isoValue)
1064 {
1065 1 mIsoValue = isoValue;
1066
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 mAboveSign = mParent->mIsInputSdf ? SdfValueT(1) : SdfValueT(-1);
1067 1 SdfTreeT &tree = mSdfGrid->tree();//sdf
1068 const bool hasActiveTiles = tree.hasActiveTiles();
1069
1070
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 1 if (mParent->mIsInputSdf && hasActiveTiles) {
1071 ✗ OPENVDB_THROW(RuntimeError, "FastSweeping: A SDF should not have active tiles!");
1072 }
1073
1074 #ifdef BENCHMARK_FAST_SWEEPING
1075 util::CpuTimer timer("Initialize voxels");
1076 #endif
1077 1 mParent->mSweepMask.clear();
1078 1 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1079
1080 {// Process all voxels
1081 2 tree::LeafManager<SdfTreeT> mgr(tree, 1);// we need one auxiliary buffer
1082
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 tbb::parallel_for(mgr.leafRange(32), *this);//multi-threaded
1083
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 mgr.swapLeafBuffer(1);//swap voxel values
1084 }
1085
1086 #ifdef BENCHMARK_FAST_SWEEPING
1087 timer.restart("Initialize tiles - new");
1088 #endif
1089 // Process all tiles
1090 2 tree::NodeManager<SdfTreeT, SdfTreeT::RootNodeType::LEVEL-1> mgr(tree);
1091
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 mgr.foreachBottomUp(*this);//multi-threaded
1092
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 tree.root().setBackground(std::numeric_limits<SdfValueT>::max(), false);
1093
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (hasActiveTiles) tree.voxelizeActiveTiles();//multi-threaded
1094
1095 // cache the leaf node origins for fast lookup in the sweeping kernels
1096
1097
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 mParent->computeSweepMaskLeafOrigins();
1098 1 }// FastSweeping::InitSdf::run
1099
1100 32 void operator()(const LeafRange& r) const
1101 {
1102 32 SweepMaskAccT sweepMaskAcc(mParent->mSweepMask);
1103 32 math::GradStencil<SdfGridT, false> stencil(*mSdfGrid);
1104 32 const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();//local copy
1105
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 32 const SdfValueT h = mAboveSign*static_cast<SdfValueT>(mSdfGrid->voxelSize()[0]);//Voxel size
1106
2/2
✓ Branch 1 taken 788 times.
✓ Branch 2 taken 32 times.
 820 for (auto leafIter = r.begin(); leafIter; ++leafIter) {
1107
1/2
✓ Branch 1 taken 788 times.
✗ Branch 2 not taken.
 788 SdfValueT* sdf = leafIter.buffer(1).data();
1108
2/2
✓ Branch 1 taken 403456 times.
✓ Branch 2 taken 788 times.
 405032 for (auto voxelIter = leafIter->beginValueAll(); voxelIter; ++voxelIter) {
1109 403456 const SdfValueT value = *voxelIter;
1110 403456 const bool isAbove = value > isoValue;
1111
3/4
✓ Branch 1 taken 403456 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 184429 times.
✓ Branch 4 taken 219027 times.
 403456 if (!voxelIter.isValueOn()) {// inactive voxels
1112
1/2
✓ Branch 0 taken 184429 times.
✗ Branch 1 not taken.
 184429 sdf[voxelIter.pos()] = isAbove ? above : -above;
1113 } else {// active voxels
1114
1/2
✓ Branch 1 taken 219027 times.
✗ Branch 2 not taken.
 219027 const Coord ijk = voxelIter.getCoord();
1115 stencil.moveTo(ijk, value);
1116 219027 const auto mask = stencil.intersectionMask( isoValue );
1117
2/2
✓ Branch 0 taken 169867 times.
✓ Branch 1 taken 49160 times.
 219027 if (mask.none()) {// most common case
1118
2/2
✓ Branch 0 taken 20316 times.
✓ Branch 1 taken 149551 times.
 169867 sdf[voxelIter.pos()] = isAbove ? above : -above;
1119 } else {// compute distance to iso-surface
1120 // disable boundary voxels from the mask tree
1121 sweepMaskAcc.setValueOff(ijk);
1122
2/2
✓ Branch 0 taken 24918 times.
✓ Branch 1 taken 24242 times.
 49160 const SdfValueT delta = value - isoValue;//offset relative to iso-value
1123 if (math::isApproxZero(delta)) {//voxel is on the iso-surface
1124 ✗ sdf[voxelIter.pos()] = 0;
1125 } else {//voxel is neighboring the iso-surface
1126 SdfValueT sum = 0;
1127
2/2
✓ Branch 0 taken 147480 times.
✓ Branch 1 taken 49160 times.
 196640 for (int i=0; i<6;) {
1128 SdfValueT d = std::numeric_limits<SdfValueT>::max(), d2;
1129
3/4
✓ Branch 1 taken 147480 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 44358 times.
✓ Branch 4 taken 103122 times.
 147480 if (mask.test(i++)) d = math::Abs(delta/(value-stencil.getValue(i)));
1130
3/4
✓ Branch 1 taken 147480 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 44358 times.
✓ Branch 4 taken 103122 times.
 147480 if (mask.test(i++)) {
1131
1/2
✓ Branch 1 taken 44358 times.
✗ Branch 2 not taken.
 44358 d2 = math::Abs(delta/(value-stencil.getValue(i)));
1132
1/2
✓ Branch 0 taken 44358 times.
✗ Branch 1 not taken.
 44358 if (d2 < d) d = d2;
1133 }
1134
2/2
✓ Branch 0 taken 58764 times.
✓ Branch 1 taken 88716 times.
 147480 if (d < std::numeric_limits<SdfValueT>::max()) sum += 1/(d*d);
1135 }
1136
2/2
✓ Branch 0 taken 24242 times.
✓ Branch 1 taken 24918 times.
 49160 sdf[voxelIter.pos()] = isAbove ? h / math::Sqrt(sum) : -h / math::Sqrt(sum);
1137 }// voxel is neighboring the iso-surface
1138 }// intersecting voxels
1139 }// active voxels
1140 }// loop over voxels
1141 }// loop over leaf nodes
1142 32 }// FastSweeping::InitSdf::operator(const LeafRange&)
1143
1144 template<typename RootOrInternalNodeT>
1145 34 void operator()(const RootOrInternalNodeT& node) const
1146 {
1147 34 const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();
1148
3/3
✓ Branch 0 taken 294116 times.
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 1 times.
 588268 for (auto it = node.cbeginValueAll(); it; ++it) {
1149 SdfValueT& v = const_cast<SdfValueT&>(*it);
1150
2/2
✓ Branch 0 taken 293522 times.
✓ Branch 1 taken 594 times.
 588232 v = v > isoValue ? above : -above;
1151 }//loop over all tiles
1152 34 }// FastSweeping::InitSdf::operator()(const RootOrInternalNodeT&)
1153
1154 // Public member data
1155 FastSweeping *mParent;
1156 SdfGridT *mSdfGrid;//raw pointer, i.e. lock free
1157 SdfValueT mIsoValue;
1158 SdfValueT mAboveSign;//sign of distance values above the iso-value
1159 };// FastSweeping::InitSdf
1160
1161
1162 /// Private class of FastSweeping to perform multi-threaded initialization
1163 template <typename SdfGridT, typename ExtValueT>
1164 template <typename OpT>
1165 struct FastSweeping<SdfGridT, ExtValueT>::InitExt
1166 {
1167 using LeafRange = typename tree::LeafManager<SdfTreeT>::LeafRange;
1168 using OpPoolT = tbb::enumerable_thread_specific<OpT>;
1169 3 InitExt(FastSweeping &parent) : mParent(&parent),
1170 mOpPool(nullptr), mSdfGrid(parent.mSdfGrid.get()),
1171 3 mExtGrid(parent.mExtGrid.get()), mIsoValue(0), mAboveSign(0) {}
1172 InitExt(const InitExt&) = default;// for tbb::parallel_for
1173 InitExt& operator=(const InitExt&) = delete;
1174 6 void run(SdfValueT isoValue, const OpT &opPrototype)
1175 {
1176 static_assert(std::is_convertible<decltype(opPrototype(Vec3d(0))),ExtValueT>::value, "Invalid return type of functor");
1177
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 6 if (!mExtGrid) {
1178 ✗ OPENVDB_THROW(RuntimeError, "FastSweeping::InitExt expected an extension grid!");
1179 }
1180
1181
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2 times.
 6 mAboveSign = mParent->mIsInputSdf ? SdfValueT(1) : SdfValueT(-1);
1182 6 mIsoValue = isoValue;
1183 6 auto &tree1 = mSdfGrid->tree();
1184 auto &tree2 = mExtGrid->tree();
1185 const bool hasActiveTiles = tree1.hasActiveTiles();//very fast
1186
1187
3/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
 6 if (mParent->mIsInputSdf && hasActiveTiles) {
1188 ✗ OPENVDB_THROW(RuntimeError, "FastSweeping: A SDF should not have active tiles!");
1189 }
1190
1191 #ifdef BENCHMARK_FAST_SWEEPING
1192 util::CpuTimer timer("Initialize voxels");
1193 #endif
1194
1195 6 mParent->mSweepMask.clear();
1196 6 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1197
1198 {// Process all voxels
1199 // Define thread-local operators
1200 12 OpPoolT opPool(opPrototype);
1201 6 mOpPool = &opPool;
1202
1203
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 12 tree::LeafManager<SdfTreeT> mgr(tree1, 1);// we need one auxiliary buffer
1204
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 6 tbb::parallel_for(mgr.leafRange(32), *this);//multi-threaded
1205
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 6 mgr.swapLeafBuffer(1);//swap out auxiliary buffer
1206 }
1207
1208 #ifdef BENCHMARK_FAST_SWEEPING
1209 timer.restart("Initialize tiles");
1210 #endif
1211 {// Process all tiles
1212 12 tree::NodeManager<SdfTreeT, SdfTreeT::RootNodeType::LEVEL-1> mgr(tree1);
1213
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 6 mgr.foreachBottomUp(*this);//multi-threaded
1214
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 6 tree1.root().setBackground(std::numeric_limits<SdfValueT>::max(), false);
1215
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2 times.
 6 if (hasActiveTiles) {
1216 #ifdef BENCHMARK_FAST_SWEEPING
1217 timer.restart("Voxelizing active tiles");
1218 #endif
1219 tree1.voxelizeActiveTiles();//multi-threaded
1220 tree2.voxelizeActiveTiles();//multi-threaded
1221 }
1222 }
1223
1224 // cache the leaf node origins for fast lookup in the sweeping kernels
1225
1226 6 mParent->computeSweepMaskLeafOrigins();
1227
1228 #ifdef BENCHMARK_FAST_SWEEPING
1229 timer.stop();
1230 #endif
1231 6 }// FastSweeping::InitExt::run
1232
1233 // int implements an update if minD needs to be updated
1234 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<std::is_same<ExtT, int>::value, int>::type = 0>
1235 void sumHelper(ExtT& sum2, ExtT ext, bool update, const SdfT& /* d2 */) const { if (update) sum2 = ext; }// int implementation
1236
1237 // non-int implements a weighted sum
1238 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<!std::is_same<ExtT, int>::value, int>::type = 0>
1239 465090 void sumHelper(ExtT& sum2, ExtT ext, bool /* update */, const SdfT& d2) const { sum2 += static_cast<ExtValueT>(d2 * ext); }// non-int implementation
1240
1241 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<std::is_same<ExtT, int>::value, int>::type = 0>
1242 ExtT extValHelper(ExtT extSum, const SdfT& /* sdfSum */) const { return extSum; }// int implementation
1243
1244 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<!std::is_same<ExtT, int>::value, int>::type = 0>
1245 466572 ExtT extValHelper(ExtT extSum, const SdfT& sdfSum) const {return ExtT((SdfT(1) / sdfSum) * extSum); }// non-int implementation
1246
1247 576 void operator()(const LeafRange& r) const
1248 {
1249 576 ExtAccT acc(mExtGrid->tree());
1250 576 SweepMaskAccT sweepMaskAcc(mParent->mSweepMask);
1251 576 math::GradStencil<SdfGridT, false> stencil(*mSdfGrid);
1252
1/2
✓ Branch 1 taken 288 times.
✗ Branch 2 not taken.
 576 const math::Transform& xform = mExtGrid->transform();
1253
1/2
✓ Branch 1 taken 288 times.
✗ Branch 2 not taken.
 576 typename OpPoolT::reference op = mOpPool->local();
1254 576 const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();//local copy
1255
1/2
✓ Branch 1 taken 288 times.
✗ Branch 2 not taken.
 576 const SdfValueT h = mAboveSign*static_cast<SdfValueT>(mSdfGrid->voxelSize()[0]);//Voxel size
1256
2/2
✓ Branch 1 taken 8838 times.
✓ Branch 2 taken 288 times.
 18252 for (auto leafIter = r.begin(); leafIter; ++leafIter) {
1257
1/2
✓ Branch 1 taken 8838 times.
✗ Branch 2 not taken.
 17676 SdfValueT *sdf = leafIter.buffer(1).data();
1258
1/2
✓ Branch 1 taken 8838 times.
✗ Branch 2 not taken.
 17676 ExtValueT *ext = acc.probeLeaf(leafIter->origin())->buffer().data();//should be safe!
1259
2/2
✓ Branch 1 taken 4525056 times.
✓ Branch 2 taken 8838 times.
 9085464 for (auto voxelIter = leafIter->beginValueAll(); voxelIter; ++voxelIter) {
1260 9050112 const SdfValueT value = *voxelIter;
1261 9050112 const bool isAbove = value > isoValue;
1262
3/4
✓ Branch 1 taken 4525056 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2798049 times.
✓ Branch 4 taken 1727007 times.
 9050112 if (!voxelIter.isValueOn()) {// inactive voxels
1263
2/2
✓ Branch 0 taken 1371819 times.
✓ Branch 1 taken 1426230 times.
 5596098 sdf[voxelIter.pos()] = isAbove ? above : -above;
1264 } else {// active voxels
1265
1/2
✓ Branch 1 taken 1727007 times.
✗ Branch 2 not taken.
 3454014 const Coord ijk = voxelIter.getCoord();
1266 stencil.moveTo(ijk, value);
1267 3454014 const auto mask = stencil.intersectionMask( isoValue );
1268
2/2
✓ Branch 0 taken 1260135 times.
✓ Branch 1 taken 466872 times.
 3454014 if (mask.none()) {// no zero-crossing neighbors, most common case
1269
2/2
✓ Branch 0 taken 544552 times.
✓ Branch 1 taken 715583 times.
 2520270 sdf[voxelIter.pos()] = isAbove ? above : -above;
1270 // the ext grid already has its active values set to the background value
1271 } else {// compute distance to iso-surface
1272 // disable boundary voxels from the mask tree
1273 sweepMaskAcc.setValueOff(ijk);
1274
2/2
✓ Branch 0 taken 232686 times.
✓ Branch 1 taken 234186 times.
 933744 const SdfValueT delta = value - isoValue;//offset relative to iso-value
1275 if (math::isApproxZero(delta)) {//voxel is on the iso-surface
1276
1/2
✓ Branch 1 taken 300 times.
✗ Branch 2 not taken.
 600 sdf[voxelIter.pos()] = 0;
1277 600 ext[voxelIter.pos()] = ExtValueT(op(xform.indexToWorld(ijk)));
1278 } else {//voxel is neighboring the iso-surface
1279 SdfValueT sum1 = 0;
1280 417412 ExtValueT sum2 = zeroVal<ExtValueT>();
1281 // minD is used to update sum2 in the integer case,
1282 // where we choose the value of the extension grid corresponding to
1283 // the smallest value of d in the 6 direction neighboring the current
1284 // voxel.
1285 SdfValueT minD = std::numeric_limits<SdfValueT>::max();
1286
2/2
✓ Branch 0 taken 1399716 times.
✓ Branch 1 taken 466572 times.
 3732576 for (int n=0, i=0; i<6;) {
1287 SdfValueT d = std::numeric_limits<SdfValueT>::max(), d2;
1288
3/4
✓ Branch 1 taken 1399716 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 420732 times.
✓ Branch 4 taken 978984 times.
 2799432 if (mask.test(i++)) {
1289 841464 d = math::Abs(delta/(value-stencil.getValue(i)));
1290 n = i - 1;
1291 }
1292
3/4
✓ Branch 1 taken 1399716 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 420732 times.
✓ Branch 4 taken 978984 times.
 2799432 if (mask.test(i++)) {
1293
1/2
✓ Branch 1 taken 420732 times.
✗ Branch 2 not taken.
 841464 d2 = math::Abs(delta/(value-stencil.getValue(i)));
1294
1/2
✓ Branch 0 taken 420732 times.
✗ Branch 1 not taken.
 841464 if (d2 < d) {
1295 d = d2;
1296 n = i - 1;
1297 }
1298 }
1299
2/2
✓ Branch 0 taken 558252 times.
✓ Branch 1 taken 841464 times.
 2799432 if (d < std::numeric_limits<SdfValueT>::max()) {
1300 1682928 d2 = 1/(d*d);
1301
1/2
✓ Branch 1 taken 841464 times.
✗ Branch 2 not taken.
 1682928 sum1 += d2;
1302 1682928 const Vec3R xyz(static_cast<SdfValueT>(ijk[0])+d*static_cast<SdfValueT>(FastSweeping::mOffset[n][0]),
1303 1682928 static_cast<SdfValueT>(ijk[1])+d*static_cast<SdfValueT>(FastSweeping::mOffset[n][1]),
1304
1/2
✓ Branch 1 taken 841464 times.
✗ Branch 2 not taken.
 1682928 static_cast<SdfValueT>(ijk[2])+d*static_cast<SdfValueT>(FastSweeping::mOffset[n][2]));
1305 // If current d is less than minD, update minD
1306 1682928 sumHelper(sum2, ExtValueT(op(xform.indexToWorld(xyz))), d < minD, d2);
1307
2/2
✓ Branch 0 taken 632212 times.
✓ Branch 1 taken 209252 times.
 1682928 if (d < minD) minD = d;
1308 }
1309 }//look over six cases
1310
2/2
✓ Branch 0 taken 104972 times.
✓ Branch 1 taken 103734 times.
 933144 ext[voxelIter.pos()] = extValHelper(sum2, sum1);
1311
2/2
✓ Branch 0 taken 234186 times.
✓ Branch 1 taken 232386 times.
 933144 sdf[voxelIter.pos()] = isAbove ? h / math::Sqrt(sum1) : -h / math::Sqrt(sum1);
1312 }// voxel is neighboring the iso-surface
1313 }// intersecting voxels
1314 }// active voxels
1315 }// loop over voxels
1316 }// loop over leaf nodes
1317 576 }// FastSweeping::InitExt::operator(const LeafRange& r)
1318
1319 template<typename RootOrInternalNodeT>
1320 102 void operator()(const RootOrInternalNodeT& node) const
1321 {
1322 102 const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();
1323
3/3
✓ Branch 0 taken 875874 times.
✓ Branch 1 taken 48 times.
✓ Branch 2 taken 3 times.
 1751856 for (auto it = node.cbeginValueAll(); it; ++it) {
1324 SdfValueT& v = const_cast<SdfValueT&>(*it);
1325
2/2
✓ Branch 0 taken 306132 times.
✓ Branch 1 taken 569742 times.
 1751748 v = v > isoValue ? above : -above;
1326 }//loop over all tiles
1327 102 }
1328 // Public member data
1329 FastSweeping *mParent;
1330 OpPoolT *mOpPool;
1331 SdfGridT *mSdfGrid;
1332 ExtGridT *mExtGrid;
1333 SdfValueT mIsoValue;
1334 SdfValueT mAboveSign;//sign of distance values above the iso-value
1335 };// FastSweeping::InitExt
1336
1337 /// Private class of FastSweeping to perform multi-threaded initialization
1338 template <typename SdfGridT, typename ExtValueT>
1339 template <typename MaskTreeT>
1340 struct FastSweeping<SdfGridT, ExtValueT>::MaskKernel
1341 {
1342 using LeafRange = typename tree::LeafManager<const MaskTreeT>::LeafRange;
1343 3 MaskKernel(FastSweeping &parent) : mParent(&parent),
1344
2/8
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
 3 mSdfGrid(parent.mSdfGrid.get()) {}
1345 MaskKernel(const MaskKernel &parent) = default;// for tbb::parallel_for
1346 MaskKernel& operator=(const MaskKernel&) = delete;
1347
1348 6 void run(const MaskTreeT &mask)
1349 {
1350 #ifdef BENCHMARK_FAST_SWEEPING
1351 util::CpuTimer timer;
1352 #endif
1353 6 auto &lsTree = mSdfGrid->tree();
1354
1355 static const SdfValueT Unknown = std::numeric_limits<SdfValueT>::max();
1356
1357 #ifdef BENCHMARK_FAST_SWEEPING
1358 timer.restart("Changing background value");
1359 #endif
1360 6 changeLevelSetBackground(lsTree, Unknown);//multi-threaded
1361
1362 #ifdef BENCHMARK_FAST_SWEEPING
1363 timer.restart("Union with mask");//multi-threaded
1364 #endif
1365 lsTree.topologyUnion(mask);//multi-threaded
1366
1367 // ignore active tiles since the input grid is assumed to be a level set
1368 12 tree::LeafManager<const MaskTreeT> mgr(mask);// super fast
1369
1370 #ifdef BENCHMARK_FAST_SWEEPING
1371 timer.restart("Initializing grid and sweep mask");
1372 #endif
1373
1374
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 6 mParent->mSweepMask.clear();
1375
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 6 mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1376
1377 using LeafManagerT = tree::LeafManager<SweepMaskTreeT>;
1378 using LeafT = typename SweepMaskTreeT::LeafNodeType;
1379
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 12 LeafManagerT leafManager(mParent->mSweepMask);
1380
1381 12990 auto kernel = [&](LeafT& leaf, size_t /*leafIdx*/) {
1382 static const SdfValueT Unknown = std::numeric_limits<SdfValueT>::max();
1383 6492 SdfAccT acc(mSdfGrid->tree());
1384 // The following hack is safe due to the topology union in
1385 // init and the fact that SdfValueT is known to be a floating point!
1386 6492 SdfValueT *data = acc.probeLeaf(leaf.origin())->buffer().data();
1387
4/6
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 374964 times.
✓ Branch 3 taken 1781 times.
✓ Branch 4 taken 1623245 times.
✓ Branch 5 taken 4711 times.
 2004701 for (auto voxelIter = leaf.beginValueOn(); voxelIter; ++voxelIter) {// mask voxels
1388
4/6
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 270638 times.
✓ Branch 3 taken 104326 times.
✓ Branch 4 taken 541276 times.
✓ Branch 5 taken 1081969 times.
 1998209 if (math::Abs( data[voxelIter.pos()] ) < Unknown ) {
1389 // disable boundary voxels from the mask tree
1390 811914 voxelIter.setValue(false);
1391 }
1392 }
1393 };
1394
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 6 leafManager.foreach( kernel );
1395
1396 // cache the leaf node origins for fast lookup in the sweeping kernels
1397
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 6 mParent->computeSweepMaskLeafOrigins();
1398
1399 #ifdef BENCHMARK_FAST_SWEEPING
1400 timer.stop();
1401 #endif
1402 }// FastSweeping::MaskKernel::run
1403
1404 // Private member data of MaskKernel
1405 FastSweeping *mParent;
1406 SdfGridT *mSdfGrid;//raw pointer, i.e. lock free
1407 };// FastSweeping::MaskKernel
1408
1409 /// @brief Private class of FastSweeping to perform concurrent fast sweeping in two directions
1410 template <typename SdfGridT, typename ExtValueT>
1411 struct FastSweeping<SdfGridT, ExtValueT>::SweepingKernel
1412 {
1413 32 SweepingKernel(FastSweeping &parent) : mParent(&parent) {}
1414 SweepingKernel(const SweepingKernel&) = delete;
1415 SweepingKernel& operator=(const SweepingKernel&) = delete;
1416
1417 /// Main method that performs concurrent bi-directional sweeps
1418 template<typename HashOp>
1419 64 void computeVoxelSlices(HashOp hash)
1420 {
1421 #ifdef BENCHMARK_FAST_SWEEPING
1422 util::CpuTimer timer;
1423 #endif
1424
1425 // mask of the active voxels to be solved for, i.e. excluding boundary voxels
1426 64 const SweepMaskTreeT& maskTree = mParent->mSweepMask;
1427
1428 using LeafManagerT = typename tree::LeafManager<const SweepMaskTreeT>;
1429 using LeafT = typename SweepMaskTreeT::LeafNodeType;
1430 128 LeafManagerT leafManager(maskTree);
1431
1432 // compute the leaf node slices that have active voxels in them
1433 // the sliding window of the has keys is -14 to 21 (based on an 8x8x8 leaf node
1434 // and the extrema hash values i-j-k and i+j+k), but we use a larger mask window here to
1435 // easily accommodate any leaf dimension. The mask offset is used to be able to
1436 // store this in a fixed-size byte array
1437 constexpr int maskOffset = LeafT::DIM * 3;
1438 constexpr int maskRange = maskOffset * 2;
1439
1440 // mark each possible slice in each leaf node that has one or more active voxels in it
1441
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 64 std::vector<int8_t> leafSliceMasks(leafManager.leafCount()*maskRange);
1442 15763308 auto kernel1 = [&](const LeafT& leaf, size_t leafIdx) {
1443 67420 const size_t leafOffset = leafIdx * maskRange;
1444
16/16
✓ Branch 0 taken 545134 times.
✓ Branch 1 taken 4025 times.
✓ Branch 2 taken 545134 times.
✓ Branch 3 taken 4025 times.
✓ Branch 4 taken 545134 times.
✓ Branch 5 taken 4025 times.
✓ Branch 6 taken 545134 times.
✓ Branch 7 taken 4025 times.
✓ Branch 8 taken 3361967 times.
✓ Branch 9 taken 12830 times.
✓ Branch 10 taken 3361967 times.
✓ Branch 11 taken 12830 times.
✓ Branch 12 taken 3361967 times.
✓ Branch 13 taken 12830 times.
✓ Branch 14 taken 3361967 times.
✓ Branch 15 taken 12830 times.
 15695824 for (auto voxelIter = leaf.cbeginValueOn(); voxelIter; ++voxelIter) {
1445 15628404 const Coord ijk = LeafT::offsetToLocalCoord(voxelIter.pos());
1446 15628404 leafSliceMasks[leafOffset + hash(ijk) + maskOffset] = uint8_t(1);
1447 }
1448 };
1449
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 64 leafManager.foreach( kernel1 );
1450
1451 // compute the voxel slice map using a thread-local-storage hash map
1452 // the key of the hash map is the slice index of the voxel coord (ijk.x() + ijk.y() + ijk.z())
1453 // the values are an array of indices for every leaf that has active voxels with this slice index
1454 using ThreadLocalMap = std::unordered_map</*voxelSliceKey=*/int64_t, /*leafIdx=*/std::deque<size_t>>;
1455
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 64 tbb::enumerable_thread_specific<ThreadLocalMap> pool;
1456 67484 auto kernel2 = [&](const LeafT& leaf, size_t leafIdx) {
1457 ThreadLocalMap& map = pool.local();
1458 const Coord& origin = leaf.origin();
1459 67420 const int64_t leafKey = hash(origin);
1460 67420 const size_t leafOffset = leafIdx * maskRange;
1461
16/16
✓ Branch 0 taken 193200 times.
✓ Branch 1 taken 4025 times.
✓ Branch 2 taken 193200 times.
✓ Branch 3 taken 4025 times.
✓ Branch 4 taken 193200 times.
✓ Branch 5 taken 4025 times.
✓ Branch 6 taken 193200 times.
✓ Branch 7 taken 4025 times.
✓ Branch 8 taken 615840 times.
✓ Branch 9 taken 12830 times.
✓ Branch 10 taken 615840 times.
✓ Branch 11 taken 12830 times.
✓ Branch 12 taken 615840 times.
✓ Branch 13 taken 12830 times.
✓ Branch 14 taken 615840 times.
✓ Branch 15 taken 12830 times.
 3303580 for (int sliceIdx = 0; sliceIdx < maskRange; sliceIdx++) {
1462
16/16
✓ Branch 0 taken 56585 times.
✓ Branch 1 taken 136615 times.
✓ Branch 2 taken 56585 times.
✓ Branch 3 taken 136615 times.
✓ Branch 4 taken 56585 times.
✓ Branch 5 taken 136615 times.
✓ Branch 6 taken 56549 times.
✓ Branch 7 taken 136651 times.
✓ Branch 8 taken 213899 times.
✓ Branch 9 taken 401941 times.
✓ Branch 10 taken 213832 times.
✓ Branch 11 taken 402008 times.
✓ Branch 12 taken 213828 times.
✓ Branch 13 taken 402012 times.
✓ Branch 14 taken 213792 times.
✓ Branch 15 taken 402048 times.
 3236160 if (leafSliceMasks[leafOffset + sliceIdx] == uint8_t(1)) {
1463 1081655 const int64_t voxelSliceKey = leafKey+sliceIdx-maskOffset;
1464 1081655 map[voxelSliceKey].emplace_back(leafIdx);
1465 }
1466 }
1467 };
1468
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 64 leafManager.foreach( kernel2 );
1469
1470 // combine into a single ordered map keyed by the voxel slice key
1471 // note that this is now stored in a map ordered by voxel slice key,
1472 // so sweep slices can be processed in order
1473
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 32 times.
 128 for (auto poolIt = pool.begin(); poolIt != pool.end(); ++poolIt) {
1474 const ThreadLocalMap& map = *poolIt;
1475
2/2
✓ Branch 0 taken 7172 times.
✓ Branch 1 taken 32 times.
 14408 for (const auto& it : map) {
1476
2/2
✓ Branch 0 taken 1081655 times.
✓ Branch 1 taken 7172 times.
 2177654 for (const size_t leafIdx : it.second) {
1477
4/6
✓ Branch 1 taken 1081655 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1081655 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1068377 times.
✓ Branch 7 taken 13278 times.
 4326620 mVoxelSliceMap[it.first].emplace_back(leafIdx, NodeMaskPtrT());
1478 }
1479 }
1480 }
1481
1482 // extract the voxel slice keys for random access into the map
1483
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 64 mVoxelSliceKeys.reserve(mVoxelSliceMap.size());
1484
2/2
✓ Branch 0 taken 7172 times.
✓ Branch 1 taken 32 times.
 14408 for (const auto& it : mVoxelSliceMap) {
1485
1/2
✓ Branch 1 taken 7172 times.
✗ Branch 2 not taken.
 14344 mVoxelSliceKeys.push_back(it.first);
1486 }
1487
1488 // allocate the node masks in parallel, as the map is populated in serial
1489
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 15428 auto kernel3 = [&](tbb::blocked_range<size_t>& range) {
1490
16/16
✓ Branch 0 taken 357 times.
✓ Branch 1 taken 128 times.
✓ Branch 2 taken 357 times.
✓ Branch 3 taken 128 times.
✓ Branch 4 taken 357 times.
✓ Branch 5 taken 128 times.
✓ Branch 6 taken 357 times.
✓ Branch 7 taken 128 times.
✓ Branch 8 taken 1433 times.
✓ Branch 9 taken 896 times.
✓ Branch 10 taken 1439 times.
✓ Branch 11 taken 896 times.
✓ Branch 12 taken 1439 times.
✓ Branch 13 taken 896 times.
✓ Branch 14 taken 1433 times.
✓ Branch 15 taken 896 times.
 11268 for (size_t i = range.begin(); i < range.end(); i++) {
1491 7172 const int64_t key = mVoxelSliceKeys[i];
1492
24/32
✓ Branch 1 taken 357 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 56585 times.
✓ Branch 4 taken 357 times.
✓ Branch 6 taken 357 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 56585 times.
✓ Branch 9 taken 357 times.
✓ Branch 11 taken 357 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 56585 times.
✓ Branch 14 taken 357 times.
✓ Branch 16 taken 357 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 56549 times.
✓ Branch 19 taken 357 times.
✓ Branch 21 taken 1433 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 213899 times.
✓ Branch 24 taken 1433 times.
✓ Branch 26 taken 1439 times.
✗ Branch 27 not taken.
✓ Branch 28 taken 213832 times.
✓ Branch 29 taken 1439 times.
✓ Branch 31 taken 1439 times.
✗ Branch 32 not taken.
✓ Branch 33 taken 213828 times.
✓ Branch 34 taken 1439 times.
✓ Branch 36 taken 1433 times.
✗ Branch 37 not taken.
✓ Branch 38 taken 213792 times.
✓ Branch 39 taken 1433 times.
 1095999 for (auto& it : mVoxelSliceMap[key]) {
1493 it.second = std::make_unique<NodeMaskT>();
1494 }
1495 }
1496 };
1497
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 64 tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceKeys.size()), kernel3);
1498
1499 // each voxel slice contains a leafIdx-nodeMask pair,
1500 // this routine populates these node masks to select only the active voxels
1501 // from the mask tree that have the same voxel slice key
1502 // TODO: a small optimization here would be to union this leaf node mask with
1503 // a pre-computed one for this particular slice pattern
1504
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 11389 auto kernel4 = [&](tbb::blocked_range<size_t>& range) {
1505
16/16
✓ Branch 0 taken 357 times.
✓ Branch 1 taken 128 times.
✓ Branch 2 taken 357 times.
✓ Branch 3 taken 132 times.
✓ Branch 4 taken 357 times.
✓ Branch 5 taken 128 times.
✓ Branch 6 taken 357 times.
✓ Branch 7 taken 128 times.
✓ Branch 8 taken 1433 times.
✓ Branch 9 taken 896 times.
✓ Branch 10 taken 1439 times.
✓ Branch 11 taken 896 times.
✓ Branch 12 taken 1439 times.
✓ Branch 13 taken 949 times.
✓ Branch 14 taken 1433 times.
✓ Branch 15 taken 896 times.
 11325 for (size_t i = range.begin(); i < range.end(); i++) {
1506 7172 const int64_t voxelSliceKey = mVoxelSliceKeys[i];
1507 7172 LeafSliceArray& leafSliceArray = mVoxelSliceMap[voxelSliceKey];
1508
16/16
✓ Branch 0 taken 56585 times.
✓ Branch 1 taken 357 times.
✓ Branch 2 taken 56585 times.
✓ Branch 3 taken 357 times.
✓ Branch 4 taken 56585 times.
✓ Branch 5 taken 357 times.
✓ Branch 6 taken 56549 times.
✓ Branch 7 taken 357 times.
✓ Branch 8 taken 213899 times.
✓ Branch 9 taken 1433 times.
✓ Branch 10 taken 213832 times.
✓ Branch 11 taken 1439 times.
✓ Branch 12 taken 213828 times.
✓ Branch 13 taken 1439 times.
✓ Branch 14 taken 213792 times.
✓ Branch 15 taken 1433 times.
 1088827 for (LeafSlice& leafSlice : leafSliceArray) {
1509
8/16
✗ Branch 0 not taken.
✓ Branch 1 taken 56585 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 56585 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 56585 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 56549 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 213899 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 213832 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 213828 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 213792 times.
 1081655 const size_t leafIdx = leafSlice.first;
1510 NodeMaskPtrT& nodeMask = leafSlice.second;
1511 const LeafT& leaf = leafManager.leaf(leafIdx);
1512 const Coord& origin = leaf.origin();
1513 1081655 const int64_t leafKey = hash(origin);
1514
16/16
✓ Branch 0 taken 9429389 times.
✓ Branch 1 taken 56585 times.
✓ Branch 2 taken 9429389 times.
✓ Branch 3 taken 56585 times.
✓ Branch 4 taken 9429389 times.
✓ Branch 5 taken 56585 times.
✓ Branch 6 taken 9422104 times.
✓ Branch 7 taken 56549 times.
✓ Branch 8 taken 67478216 times.
✓ Branch 9 taken 213899 times.
✓ Branch 10 taken 67488535 times.
✓ Branch 11 taken 213832 times.
✓ Branch 12 taken 67489512 times.
✓ Branch 13 taken 213828 times.
✓ Branch 14 taken 67433774 times.
✓ Branch 15 taken 213792 times.
 308681963 for (auto voxelIter = leaf.cbeginValueOn(); voxelIter; ++voxelIter) {
1515 const Index voxelIdx = voxelIter.pos();
1516 307600308 const Coord ijk = LeafT::offsetToLocalCoord(voxelIdx);
1517 307600308 const int64_t key = leafKey + hash(ijk);
1518
16/16
✓ Branch 0 taken 545134 times.
✓ Branch 1 taken 8884255 times.
✓ Branch 2 taken 545134 times.
✓ Branch 3 taken 8884255 times.
✓ Branch 4 taken 545134 times.
✓ Branch 5 taken 8884255 times.
✓ Branch 6 taken 545134 times.
✓ Branch 7 taken 8876970 times.
✓ Branch 8 taken 3361967 times.
✓ Branch 9 taken 64116249 times.
✓ Branch 10 taken 3361967 times.
✓ Branch 11 taken 64126568 times.
✓ Branch 12 taken 3361967 times.
✓ Branch 13 taken 64127545 times.
✓ Branch 14 taken 3361967 times.
✓ Branch 15 taken 64071807 times.
 307600308 if (key == voxelSliceKey) {
1519 15628404 nodeMask->setOn(voxelIdx);
1520 }
1521 }
1522 }
1523 }
1524 };
1525
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 64 tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceKeys.size()), kernel4);
1526 64 }// FastSweeping::SweepingKernel::computeVoxelSlices
1527
1528 // Private struct for nearest neighbor grid points (very memory light!)
1529 struct NN {
1530 SdfValueT v;
1531 int n;
1532 199635633 inline static Coord ijk(const Coord &p, int i) { return p + FastSweeping::mOffset[i]; }
1533 1598891 NN() : v(), n() {}
1534 375081696 NN(const SdfAccT &a, const Coord &p, int i) : v(math::Abs(a.getValue(ijk(p,i)))), n(i) {}
1535 ✗ inline Coord operator()(const Coord &p) const { return ijk(p, n); }
1536
12/12
✓ Branch 0 taken 1983754 times.
✓ Branch 1 taken 2377318 times.
✓ Branch 2 taken 1983991 times.
✓ Branch 3 taken 2377081 times.
✓ Branch 4 taken 1986202 times.
✓ Branch 5 taken 2374870 times.
✓ Branch 6 taken 11515067 times.
✓ Branch 7 taken 15380669 times.
✓ Branch 8 taken 12218056 times.
✓ Branch 9 taken 14677680 times.
✓ Branch 10 taken 12212823 times.
✓ Branch 11 taken 14682913 times.
 93770424 inline bool operator<(const NN &rhs) const { return v < rhs.v; }
1537 31957969 inline operator bool() const { return v < SdfValueT(1000); }
1538 };// NN
1539
1540 /// @note Extending an integer field is based on the nearest-neighbor interpolation
1541 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<std::is_same<ExtT, int>::value, int>::type = 0>
1542 ExtT twoNghbr(const NN& d1, const NN& d2, const SdfT& /* w */, const ExtT& v1, const ExtT& v2) const { return d1.v < d2.v ? v1 : v2; }// int implementation
1543
1544 /// @note Extending a non-integer field is based on a weighted interpolation
1545 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<!std::is_same<ExtT, int>::value, int>::type = 0>
1546 1012251 ExtT twoNghbr(const NN& d1, const NN& d2, const SdfT& w, const ExtT& v1, const ExtT& v2) const { return ExtT(w*(d1.v*v1 + d2.v*v2)); }// non-int implementation
1547
1548 /// @note Extending an integer field is based on the nearest-neighbor interpolation
1549 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<std::is_same<ExtT, int>::value, int>::type = 0>
1550 ExtT threeNghbr(const NN& d1, const NN& d2, const NN& d3, const SdfT& /* w */, const ExtT& v1, const ExtT& v2, const ExtT& v3) const {
1551 math::Vec3<SdfT> d(d1.v, d2.v, d3.v);
1552 math::Vec3<ExtT> v(v1, v2, v3);
1553 return v[static_cast<int>(math::MinIndex(d))];
1554 }// int implementation
1555
1556 /// @note Extending a non-integer field is based on a weighted interpolation
1557 template<typename ExtT = ExtValueT, typename SdfT = SdfValueT, typename std::enable_if<!std::is_same<ExtT, int>::value, int>::type = 0>
1558 942297 ExtT threeNghbr(const NN& d1, const NN& d2, const NN& d3, const SdfT& w, const ExtT& v1, const ExtT& v2, const ExtT& v3) const {
1559 3172041 return ExtT(w*(d1.v*v1 + d2.v*v2 + d3.v*v3));
1560 }// non-int implementation
1561
1562 64 void sweep()
1563 {
1564
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 20 times.
 64 typename ExtGridT::TreeType *tree2 = mParent->mExtGrid ? &mParent->mExtGrid->tree() : nullptr;
1565
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
 64 typename ExtGridT::TreeType *tree3 = mParent->mExtGridInput ? &mParent->mExtGridInput->tree() : nullptr;
1566
1567 64 const SdfValueT h = static_cast<SdfValueT>(mParent->mSdfGrid->voxelSize()[0]);
1568 64 const SdfValueT sqrt2h = math::Sqrt(SdfValueT(2))*h;
1569 64 const FastSweepingDomain mode = mParent->mSweepDirection;
1570 64 const bool isInputSdf = mParent->mIsInputSdf;
1571
1572 // If we are using an extension in one direction, we need a reference grid
1573 // for the default value of the extension for the voxels that are not
1574 // intended to be updated by the sweeping algorithm.
1575
3/6
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 12 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 64 if (tree2 && mode != FastSweepingDomain::SWEEP_ALL) assert(tree3);
1576
1577 64 const std::vector<Coord>& leafNodeOrigins = mParent->mSweepMaskLeafOrigins;
1578
1579 64 int64_t voxelSliceIndex(0);
1580
1581 1598955 auto kernel = [&](const tbb::blocked_range<size_t>& range) {
1582 using LeafT = typename SdfGridT::TreeType::LeafNodeType;
1583
1584 1598891 SdfAccT acc1(mParent->mSdfGrid->tree());
1585
3/4
✓ Branch 0 taken 354976 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 505864 times.
✓ Branch 5 taken 738051 times.
 1598891 auto acc2 = std::unique_ptr<ExtAccT>(tree2 ? new ExtAccT(*tree2) : nullptr);
1586
2/8
✗ Branch 0 not taken.
✓ Branch 1 taken 354976 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1243915 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 1598891 auto acc3 = std::unique_ptr<ExtAccT>(tree3 ? new ExtAccT(*tree3) : nullptr);
1587 SdfValueT absV, sign, update, D;
1588 NN d1, d2, d3;//distance values and coordinates of closest neighbor points
1589
1590
2/4
✓ Branch 1 taken 354976 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1243915 times.
✗ Branch 5 not taken.
 1598891 const LeafSliceArray& leafSliceArray = mVoxelSliceMap[voxelSliceIndex];
1591
1592 // Solves Godunov's scheme: [x-d1]^2 + [x-d2]^2 + [x-d3]^2 = h^2
1593 // where [X] = (X>0?X:0) and ai=min(di+1,di-1)
1594
4/4
✓ Branch 0 taken 452608 times.
✓ Branch 1 taken 354976 times.
✓ Branch 2 taken 1710702 times.
✓ Branch 3 taken 1243915 times.
 3762201 for (size_t i = range.begin(); i < range.end(); ++i) {
1595
1596 // iterate over all leafs in the slice and extract the leaf
1597 // and node mask for each slice pattern
1598
1599 const LeafSlice& leafSlice = leafSliceArray[i];
1600 2163310 const size_t leafIdx = leafSlice.first;
1601 const NodeMaskPtrT& nodeMask = leafSlice.second;
1602
1603 const Coord& origin = leafNodeOrigins[leafIdx];
1604
1605 Coord ijk;
1606
4/4
✓ Branch 1 taken 4361072 times.
✓ Branch 2 taken 452608 times.
✓ Branch 4 taken 26895736 times.
✓ Branch 5 taken 1710702 times.
 35583428 for (auto indexIter = nodeMask->beginOn(); indexIter; ++indexIter) {
1607
1608 // Get coordinate of center point of the FD stencil
1609 31256808 ijk = origin + LeafT::offsetToLocalCoord(indexIter.pos());
1610
1611 // Find the closes neighbors in the three axial directions
1612
4/8
✓ Branch 1 taken 4361072 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4361072 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 26895736 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 26895736 times.
✗ Branch 11 not taken.
 31256808 d1 = std::min(NN(acc1, ijk, 0), NN(acc1, ijk, 1));
1613
4/8
✓ Branch 1 taken 4361072 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4361072 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 26895736 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 26895736 times.
✗ Branch 11 not taken.
 31256808 d2 = std::min(NN(acc1, ijk, 2), NN(acc1, ijk, 3));
1614
4/8
✓ Branch 1 taken 4361072 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4361072 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 26895736 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 26895736 times.
✗ Branch 11 not taken.
 31256808 d3 = std::min(NN(acc1, ijk, 4), NN(acc1, ijk, 5));
1615
1616
12/12
✓ Branch 0 taken 94933 times.
✓ Branch 1 taken 4266139 times.
✓ Branch 2 taken 54095 times.
✓ Branch 3 taken 40838 times.
✓ Branch 4 taken 43169 times.
✓ Branch 5 taken 10926 times.
✓ Branch 6 taken 309397 times.
✓ Branch 7 taken 26586339 times.
✓ Branch 8 taken 242736 times.
✓ Branch 9 taken 66661 times.
✓ Branch 10 taken 226729 times.
✓ Branch 11 taken 16007 times.
 31256808 if (!(d1 || d2 || d3)) continue;//no valid neighbors
1617
1618 // Get the center point of the FD stencil (assumed to be an active voxel)
1619 // Note this const_cast is normally unsafe but by design we know the tree
1620 // to be static, of floating-point type and containing active voxels only!
1621
2/4
✓ Branch 1 taken 4317903 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26669007 times.
✗ Branch 5 not taken.
 30986910 SdfValueT &value = const_cast<SdfValueT&>(acc1.getValue(ijk));
1622
1623 // Extract the sign
1624
4/4
✓ Branch 0 taken 2082522 times.
✓ Branch 1 taken 2235381 times.
✓ Branch 2 taken 18533751 times.
✓ Branch 3 taken 8135256 times.
 30986910 sign = value >= SdfValueT(0) ? SdfValueT(1) : SdfValueT(-1);
1625
1626 // Absolute value
1627 absV = math::Abs(value);
1628
1629 // sort values so d1 <= d2 <= d3
1630
4/4
✓ Branch 0 taken 2140268 times.
✓ Branch 1 taken 2177635 times.
✓ Branch 2 taken 11722692 times.
✓ Branch 3 taken 14946315 times.
 30986910 if (d2 < d1) std::swap(d1, d2);
1631
4/4
✓ Branch 0 taken 2852247 times.
✓ Branch 1 taken 1465656 times.
✓ Branch 2 taken 17335589 times.
✓ Branch 3 taken 9333418 times.
 30986910 if (d3 < d2) std::swap(d2, d3);
1632
4/4
✓ Branch 0 taken 1428016 times.
✓ Branch 1 taken 2889887 times.
✓ Branch 2 taken 7804584 times.
✓ Branch 3 taken 18864423 times.
 30986910 if (d2 < d1) std::swap(d1, d2);
1633
1634 // Test if there is a solution depending on ONE of the neighboring voxels
1635 // if d2 - d1 >= h => d2 >= d1 + h then:
1636 // (x-d1)^2=h^2 => x = d1 + h
1637 30986910 update = d1.v + h;
1638
4/4
✓ Branch 0 taken 159183 times.
✓ Branch 1 taken 4158720 times.
✓ Branch 2 taken 1457454 times.
✓ Branch 3 taken 25211553 times.
 30986910 if (update <= d2.v) {
1639
4/4
✓ Branch 0 taken 159102 times.
✓ Branch 1 taken 81 times.
✓ Branch 2 taken 1449225 times.
✓ Branch 3 taken 8229 times.
 1616637 if (update < absV) {
1640 1608327 value = sign * update;
1641
3/4
✓ Branch 0 taken 159102 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 395058 times.
✓ Branch 3 taken 1054167 times.
 1608327 if (acc2) {
1642 // There is an assert upstream to check if mExtGridInput exists if mode != SWEEP_ALL
1643
2/4
✓ Branch 1 taken 159102 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 395058 times.
✗ Branch 5 not taken.
 554160 ExtValueT updateExt = acc2->getValue(d1(ijk));
1644
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 159102 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 395058 times.
 554160 if (mode == FastSweepingDomain::SWEEP_GREATER_THAN_ISOVALUE) {
1645 ✗ if (isInputSdf) updateExt = (value >= SdfValueT(0)) ? acc2->getValue(d1(ijk)) : acc3->getValue(ijk);
1646 ✗ else updateExt = (value <= SdfValueT(0)) ? acc2->getValue(d1(ijk)) : acc3->getValue(ijk);
1647 } // SWEEP_GREATER_THAN_ISOVALUE
1648
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 159102 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 395058 times.
 554160 else if (mode == FastSweepingDomain::SWEEP_LESS_THAN_ISOVALUE) {
1649 ✗ if (isInputSdf) updateExt = (value <= SdfValueT(0)) ? acc2->getValue(d1(ijk)) : acc3->getValue(ijk);
1650 ✗ else updateExt = (value >= SdfValueT(0)) ? acc2->getValue(d1(ijk)) : acc3->getValue(ijk);
1651 } // SWEEP_LESS_THAN_ISOVALUE
1652
2/4
✓ Branch 1 taken 159102 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 395058 times.
✗ Branch 5 not taken.
 554160 acc2->setValue(ijk, updateExt);
1653 }//update ext?
1654 }//update sdf?
1655 1616637 continue;
1656 }// one neighbor case
1657
1658 // Test if there is a solution depending on TWO of the neighboring voxels
1659 // (x-d1)^2 + (x-d2)^2 = h^2
1660 //D = SdfValueT(2) * h * h - math::Pow2(d1.v - d2.v);// = 2h^2-(d1-d2)^2
1661 //if (D >= SdfValueT(0)) {// non-negative discriminant
1662
2/4
✓ Branch 0 taken 4158720 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 25211553 times.
✗ Branch 3 not taken.
 29370273 if (d2.v <= sqrt2h + d1.v) {
1663
2/4
✓ Branch 0 taken 4158720 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 25211553 times.
✗ Branch 3 not taken.
 29370273 D = SdfValueT(2) * h * h - math::Pow2(d1.v - d2.v);// = 2h^2-(d1-d2)^2
1664
2/4
✓ Branch 0 taken 4158720 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 25211553 times.
✗ Branch 3 not taken.
 29370273 update = SdfValueT(0.5) * (d1.v + d2.v + std::sqrt(D));
1665
6/8
✓ Branch 0 taken 4158720 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 366951 times.
✓ Branch 3 taken 3791769 times.
✓ Branch 4 taken 25211553 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2962176 times.
✓ Branch 7 taken 22249377 times.
 29370273 if (update > d2.v && update <= d3.v) {
1666
4/4
✓ Branch 0 taken 302835 times.
✓ Branch 1 taken 64116 times.
✓ Branch 2 taken 2743484 times.
✓ Branch 3 taken 218692 times.
 3329127 if (update < absV) {
1667 3046319 value = sign * update;
1668
3/4
✓ Branch 0 taken 302835 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 709416 times.
✓ Branch 3 taken 2034068 times.
 3046319 if (acc2) {
1669 1012251 d1.v -= update;
1670 1012251 d2.v -= update;
1671 // affine combination of two neighboring extension values
1672 1012251 const SdfValueT w = SdfValueT(1)/(d1.v+d2.v);
1673
4/8
✓ Branch 1 taken 302835 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 302835 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 709416 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 709416 times.
✗ Branch 11 not taken.
 1012251 const ExtValueT v1 = acc2->getValue(d1(ijk));
1674
2/4
✓ Branch 1 taken 302835 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 709416 times.
✗ Branch 5 not taken.
 1012251 const ExtValueT v2 = acc2->getValue(d2(ijk));
1675 302835 const ExtValueT extVal = twoNghbr(d1, d2, w, v1, v2);
1676
1677 1012251 ExtValueT updateExt = extVal;
1678
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 302835 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 709416 times.
 1012251 if (mode == FastSweepingDomain::SWEEP_GREATER_THAN_ISOVALUE) {
1679 ✗ if (isInputSdf) updateExt = (value >= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1680 ✗ else updateExt = (value <= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1681 } // SWEEP_GREATER_THAN_ISOVALUE
1682
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 302835 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 709416 times.
 1012251 else if (mode == FastSweepingDomain::SWEEP_LESS_THAN_ISOVALUE) {
1683 ✗ if (isInputSdf) updateExt = (value <= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1684 ✗ else updateExt = (value >= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1685 } // SWEEP_LESS_THAN_ISOVALUE
1686
2/4
✓ Branch 1 taken 302835 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 709416 times.
✗ Branch 5 not taken.
 1012251 acc2->setValue(ijk, updateExt);
1687 }//update ext?
1688 }//update sdf?
1689 3329127 continue;
1690 }//test for two neighbor case
1691 }//test for non-negative determinant
1692
1693 // Test if there is a solution depending on THREE of the neighboring voxels
1694 // (x-d1)^2 + (x-d2)^2 + (x-d3)^2 = h^2
1695 // 3x^2 - 2(d1 + d2 + d3)x + d1^2 + d2^2 + d3^2 = h^2
1696 // ax^2 + bx + c=0, a=3, b=-2(d1+d2+d3), c=d1^2 + d2^2 + d3^2 - h^2
1697 26041146 const SdfValueT d123 = d1.v + d2.v + d3.v;
1698 26041146 D = d123*d123 - SdfValueT(3)*(d1.v*d1.v + d2.v*d2.v + d3.v*d3.v - h * h);
1699
2/4
✓ Branch 0 taken 3791769 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 22249377 times.
✗ Branch 3 not taken.
 26041146 if (D >= SdfValueT(0)) {// non-negative discriminant
1700 26041146 update = SdfValueT(1.0/3.0) * (d123 + std::sqrt(D));//always passes test
1701 //if (update > d3.v) {//disabled due to round-off errors
1702
4/4
✓ Branch 0 taken 942297 times.
✓ Branch 1 taken 2849472 times.
✓ Branch 2 taken 8533983 times.
✓ Branch 3 taken 13715394 times.
 26041146 if (update < absV) {
1703 9476280 value = sign * update;
1704
3/4
✓ Branch 0 taken 942297 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2229744 times.
✓ Branch 3 taken 6304239 times.
 9476280 if (acc2) {
1705 3172041 d1.v -= update;
1706 3172041 d2.v -= update;
1707 3172041 d3.v -= update;
1708 // affine combination of three neighboring extension values
1709 3172041 const SdfValueT w = SdfValueT(1)/(d1.v+d2.v+d3.v);
1710
4/8
✓ Branch 1 taken 942297 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 942297 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2229744 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2229744 times.
✗ Branch 11 not taken.
 3172041 const ExtValueT v1 = acc2->getValue(d1(ijk));
1711
4/8
✓ Branch 1 taken 942297 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 942297 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2229744 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2229744 times.
✗ Branch 11 not taken.
 3172041 const ExtValueT v2 = acc2->getValue(d2(ijk));
1712
2/4
✓ Branch 1 taken 942297 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2229744 times.
✗ Branch 5 not taken.
 3172041 const ExtValueT v3 = acc2->getValue(d3(ijk));
1713 942297 const ExtValueT extVal = threeNghbr(d1, d2, d3, w, v1, v2, v3);
1714
1715 3172041 ExtValueT updateExt = extVal;
1716
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 942297 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2229744 times.
 3172041 if (mode == FastSweepingDomain::SWEEP_GREATER_THAN_ISOVALUE) {
1717 ✗ if (isInputSdf) updateExt = (value >= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1718 ✗ else updateExt = (value <= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1719 } // SWEEP_GREATER_THAN_ISOVALUE
1720
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 942297 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2229744 times.
 3172041 else if (mode == FastSweepingDomain::SWEEP_LESS_THAN_ISOVALUE) {
1721 ✗ if (isInputSdf) updateExt = (value <= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1722 ✗ else updateExt = (value >= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1723 } // SWEEP_LESS_THAN_ISOVALUE
1724
2/4
✓ Branch 1 taken 942297 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2229744 times.
✗ Branch 5 not taken.
 3172041 acc2->setValue(ijk, updateExt);
1725 }//update ext?
1726 }//update sdf?
1727 }//test for non-negative determinant
1728 }//loop over coordinates
1729 }
1730 };
1731
1732 #ifdef BENCHMARK_FAST_SWEEPING
1733 util::CpuTimer timer("Forward sweep");
1734 #endif
1735
1736
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 7172 times.
 14408 for (size_t i = 0; i < mVoxelSliceKeys.size(); i++) {
1737 14344 voxelSliceIndex = mVoxelSliceKeys[i];
1738 14344 tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceMap[voxelSliceIndex].size()), kernel);
1739 }
1740
1741 #ifdef BENCHMARK_FAST_SWEEPING
1742 timer.restart("Backward sweeps");
1743 #endif
1744
2/2
✓ Branch 0 taken 7172 times.
✓ Branch 1 taken 32 times.
 14408 for (size_t i = mVoxelSliceKeys.size(); i > 0; i--) {
1745 14344 voxelSliceIndex = mVoxelSliceKeys[i-1];
1746 14344 tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceMap[voxelSliceIndex].size()), kernel);
1747 }
1748
1749 #ifdef BENCHMARK_FAST_SWEEPING
1750 timer.stop();
1751 #endif
1752 64 }// FastSweeping::SweepingKernel::sweep
1753
1754 private:
1755 using NodeMaskT = typename SweepMaskTreeT::LeafNodeType::NodeMaskType;
1756 using NodeMaskPtrT = std::unique_ptr<NodeMaskT>;
1757 // using a unique ptr for the NodeMask allows for parallel allocation,
1758 // but makes this class not copy-constructible
1759 using LeafSlice = std::pair</*leafIdx=*/size_t, /*leafMask=*/NodeMaskPtrT>;
1760 using LeafSliceArray = std::deque<LeafSlice>;
1761 using VoxelSliceMap = std::map</*voxelSliceKey=*/int64_t, LeafSliceArray>;
1762
1763 // Private member data of SweepingKernel
1764 FastSweeping *mParent;
1765 VoxelSliceMap mVoxelSliceMap;
1766 std::vector<int64_t> mVoxelSliceKeys;
1767 };// FastSweeping::SweepingKernel
1768
1769 ////////////////////////////////////////////////////////////////////////////////
1770
1771 template<typename GridT>
1772 typename GridT::Ptr
1773 ✗ fogToSdf(const GridT &fogGrid,
1774 typename GridT::ValueType isoValue,
1775 int nIter)
1776 {
1777 ✗ FastSweeping<GridT> fs;
1778 ✗ if (fs.initSdf(fogGrid, isoValue, /*isInputSdf*/false)) fs.sweep(nIter);
1779 ✗ return fs.sdfGrid();
1780 }
1781
1782 template<typename GridT>
1783 typename GridT::Ptr
1784 ✗ sdfToSdf(const GridT &sdfGrid,
1785 typename GridT::ValueType isoValue,
1786 int nIter)
1787 {
1788 ✗ FastSweeping<GridT> fs;
1789 ✗ if (fs.initSdf(sdfGrid, isoValue, /*isInputSdf*/true)) fs.sweep(nIter);
1790 ✗ return fs.sdfGrid();
1791 }
1792
1793 template<typename FogGridT, typename ExtOpT, typename ExtValueT>
1794 typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr
1795 fogToExt(const FogGridT &fogGrid,
1796 const ExtOpT &op,
1797 const ExtValueT& background,
1798 typename FogGridT::ValueType isoValue,
1799 int nIter,
1800 FastSweepingDomain mode,
1801 const typename FogGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid)
1802 {
1803 FastSweeping<FogGridT, ExtValueT> fs;
1804 if (fs.initExt(fogGrid, op, background, isoValue, /*isInputSdf*/false, mode, extGrid))
1805 fs.sweep(nIter, /*finalize*/true);
1806 return fs.extGrid();
1807 }
1808
1809 template<typename SdfGridT, typename OpT, typename ExtValueT>
1810 typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr
1811 sdfToExt(const SdfGridT &sdfGrid,
1812 const OpT &op,
1813 const ExtValueT &background,
1814 typename SdfGridT::ValueType isoValue,
1815 int nIter,
1816 FastSweepingDomain mode,
1817 const typename SdfGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid)
1818 {
1819 FastSweeping<SdfGridT, ExtValueT> fs;
1820 if (fs.initExt(sdfGrid, op, background, isoValue, /*isInputSdf*/true, mode, extGrid))
1821 fs.sweep(nIter, /*finalize*/true);
1822 return fs.extGrid();
1823 }
1824
1825 template<typename FogGridT, typename ExtOpT, typename ExtValueT>
1826 std::pair<typename FogGridT::Ptr, typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr>
1827 1 fogToSdfAndExt(const FogGridT &fogGrid,
1828 const ExtOpT &op,
1829 const ExtValueT &background,
1830 typename FogGridT::ValueType isoValue,
1831 int nIter,
1832 FastSweepingDomain mode,
1833 const typename FogGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid)
1834 {
1835 1 FastSweeping<FogGridT, ExtValueT> fs;
1836
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 2 if (fs.initExt(fogGrid, op, background, isoValue, /*isInputSdf*/false, mode, extGrid))
1837
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fs.sweep(nIter, /*finalize*/true);
1838
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 return std::make_pair(fs.sdfGrid(), fs.extGrid());
1839 }
1840
1841 template<typename SdfGridT, typename ExtOpT, typename ExtValueT>
1842 std::pair<typename SdfGridT::Ptr, typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr>
1843 4 sdfToSdfAndExt(const SdfGridT &sdfGrid,
1844 const ExtOpT &op,
1845 const ExtValueT &background,
1846 typename SdfGridT::ValueType isoValue,
1847 int nIter,
1848 FastSweepingDomain mode,
1849 const typename SdfGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid)
1850 {
1851 4 FastSweeping<SdfGridT, ExtValueT> fs;
1852
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 8 if (fs.initExt(sdfGrid, op, background, isoValue, /*isInputSdf*/true, mode, extGrid))
1853
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 4 fs.sweep(nIter, /*finalize*/true);
1854
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 return std::make_pair(fs.sdfGrid(), fs.extGrid());
1855 }
1856
1857 template<typename GridT>
1858 typename GridT::Ptr
1859 ✗ dilateSdf(const GridT &sdfGrid,
1860 int dilation,
1861 NearestNeighbors nn,
1862 int nIter,
1863 FastSweepingDomain mode)
1864 {
1865 ✗ FastSweeping<GridT> fs;
1866 ✗ if (fs.initDilate(sdfGrid, dilation, nn, /*sweep direction*/ mode)) fs.sweep(nIter);
1867 ✗ return fs.sdfGrid();
1868 }
1869
1870 template<typename GridT, typename MaskTreeT>
1871 typename GridT::Ptr
1872 3 maskSdf(const GridT &sdfGrid,
1873 const Grid<MaskTreeT> &mask,
1874 bool ignoreActiveTiles,
1875 int nIter)
1876 {
1877 6 FastSweeping<GridT> fs;
1878 3 if (fs.initMask(sdfGrid, mask, ignoreActiveTiles)) fs.sweep(nIter);
1879 3 return fs.sdfGrid();
1880 }
1881
1882
1883 ////////////////////////////////////////
1884
1885
1886 // Explicit Template Instantiation
1887
1888 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
1889
1890 #ifdef OPENVDB_INSTANTIATE_FASTSWEEPING
1891 #include <openvdb/util/ExplicitInstantiation.h>
1892 #endif
1893
1894 #define _FUNCTION(TreeT) \
1895 Grid<TreeT>::Ptr fogToSdf(const Grid<TreeT>&, TreeT::ValueType, int)
1896 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
1897 #undef _FUNCTION
1898
1899 #define _FUNCTION(TreeT) \
1900 Grid<TreeT>::Ptr sdfToSdf(const Grid<TreeT>&, TreeT::ValueType, int)
1901 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
1902 #undef _FUNCTION
1903
1904 #define _FUNCTION(TreeT) \
1905 Grid<TreeT>::Ptr dilateSdf(const Grid<TreeT>&, int, NearestNeighbors, int, FastSweepingDomain)
1906 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
1907 #undef _FUNCTION
1908
1909 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
1910
1911
1912 } // namespace tools
1913 } // namespace OPENVDB_VERSION_NAME
1914 } // namespace openvdb
1915
1916 #endif // OPENVDB_TOOLS_FASTSWEEPING_HAS_BEEN_INCLUDED
1917