OpenVDB  9.0.1
FastSweeping.h
Go to the documentation of this file.
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
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 {
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
69  // Update voxels corresponding to an sdf/fog values that are less than a given 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,
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>
459 {
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;
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  ~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,
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  size_t sweepingVoxelCount() const { return mSweepingVoxelCount; }
659 
660  /// @brief Return the number of voxels that defined the boundary condition.
661  size_t boundaryVoxelCount() const { return mBoundaryVoxelCount; }
662 
663  /// @brief Return true if there are voxels and boundaries to solve for
664  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>
716  : mSdfGrid(nullptr), mExtGrid(nullptr), mSweepingVoxelCount(0), mBoundaryVoxelCount(0), mSweepDirection(FastSweepingDomain::SWEEP_ALL), mIsInputSdf(true)
717 {
718 }
719 
720 template <typename SdfGridT, typename ExtValueT>
722 {
723  mSdfGrid.reset();
724  mExtGrid.reset();
725  mSweepMask.clear();
726  if (mExtGridInput) mExtGridInput.reset();
727  mSweepingVoxelCount = mBoundaryVoxelCount = 0;
728  mSweepDirection = FastSweepingDomain::SWEEP_ALL;
729  mIsInputSdf = true;
730 }
731 
732 template <typename SdfGridT, typename ExtValueT>
734 {
735  // replace any inactive leaf nodes with tiles and voxelize any active tiles
736 
737  pruneInactive(mSweepMask);
738  mSweepMask.voxelizeActiveTiles();
739 
740  using LeafManagerT = tree::LeafManager<SweepMaskTreeT>;
741  using LeafT = typename SweepMaskTreeT::LeafNodeType;
742  LeafManagerT leafManager(mSweepMask);
743 
744  mSweepMaskLeafOrigins.resize(leafManager.leafCount());
745  std::atomic<size_t> sweepingVoxelCount{0};
746  auto kernel = [&](const LeafT& leaf, size_t leafIdx) {
747  mSweepMaskLeafOrigins[leafIdx] = leaf.origin();
748  sweepingVoxelCount += leaf.onVoxelCount();
749  };
750  leafManager.foreach(kernel, /*threaded=*/true, /*grainsize=*/1024);
751 
752  mBoundaryVoxelCount = 0;
753  mSweepingVoxelCount = sweepingVoxelCount;
754  if (mSdfGrid) {
755  const size_t totalCount = mSdfGrid->constTree().activeVoxelCount();
756  assert( totalCount >= mSweepingVoxelCount );
757  mBoundaryVoxelCount = totalCount - mSweepingVoxelCount;
758  }
759 }// FastSweeping::computeSweepMaskLeafOrigins
760 
761 template <typename SdfGridT, typename ExtValueT>
762 bool FastSweeping<SdfGridT, ExtValueT>::initSdf(const SdfGridT &fogGrid, SdfValueT isoValue, bool isInputSdf)
763 {
764  this->clear();
765  mSdfGrid = fogGrid.deepCopy();//very fast
766  mIsInputSdf = isInputSdf;
767  InitSdf kernel(*this);
768  kernel.run(isoValue);
769  return this->isValid();
770 }
771 
772 template <typename SdfGridT, typename ExtValueT>
773 template <typename OpT>
774 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  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  this->clear();
784  mSdfGrid = fogGrid.deepCopy();//very fast
785  mExtGrid = createGrid<ExtGridT>( background );
786  mSweepDirection = mode;
787  mIsInputSdf = isInputSdf;
788  if (mSweepDirection != FastSweepingDomain::SWEEP_ALL) {
789  mExtGridInput = extGrid->deepCopy();
790  }
791  mExtGrid->topologyUnion( *mSdfGrid );//very fast
792  InitExt<OpT> kernel(*this);
793  kernel.run(isoValue, op);
794  return this->isValid();
795 }
796 
797 
798 template <typename SdfGridT, typename ExtValueT>
799 bool FastSweeping<SdfGridT, ExtValueT>::initDilate(const SdfGridT &sdfGrid, int dilate, NearestNeighbors nn, FastSweepingDomain mode)
800 {
801  this->clear();
802  mSdfGrid = sdfGrid.deepCopy();//very fast
803  mSweepDirection = mode;
804  DilateKernel kernel(*this);
805  kernel.run(dilate, nn);
806  return this->isValid();
807 }
808 
809 template <typename SdfGridT, typename ExtValueT>
810 template<typename MaskTreeT>
811 bool FastSweeping<SdfGridT, ExtValueT>::initMask(const SdfGridT &sdfGrid, const Grid<MaskTreeT> &mask, bool ignoreActiveTiles)
812 {
813  this->clear();
814  mSdfGrid = sdfGrid.deepCopy();//very fast
815 
816  if (mSdfGrid->transform() != mask.transform()) {
817  OPENVDB_THROW(RuntimeError, "FastSweeping: Mask not aligned with the grid!");
818  }
819 
820  if (mask.getGridClass() == GRID_LEVEL_SET) {
821  using T = typename MaskTreeT::template ValueConverter<bool>::Type;
822  typename Grid<T>::Ptr tmp = sdfInteriorMask(mask);//might have active tiles
823  tmp->tree().voxelizeActiveTiles();//multi-threaded
824  MaskKernel<T> kernel(*this);
825  kernel.run(tmp->tree());//multi-threaded
826  } else {
827  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  T tmp(mask.tree(), false, TopologyCopy());//multi-threaded
833  tmp.voxelizeActiveTiles(true);//multi-threaded
834  MaskKernel<T> kernel(*this);
835  kernel.run(tmp);//multi-threaded
836  }
837  }
838  return this->isValid();
839 }// FastSweeping::initMask
840 
841 template <typename SdfGridT, typename ExtValueT>
842 void FastSweeping<SdfGridT, ExtValueT>::sweep(int nIter, bool finalize)
843 {
844  if (!mSdfGrid) {
845  OPENVDB_THROW(RuntimeError, "FastSweeping::sweep called before initialization!");
846  }
847  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  if (this->boundaryVoxelCount() == 0) {
852  OPENVDB_THROW(RuntimeError, "FastSweeping: No boundary voxels found!");
853  } 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  std::deque<SweepingKernel> kernels;
859  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  tasks.run([&] { kernels[0].computeVoxelSlices([](const Coord &a){ return a[0]+a[1]+a[2]; });/*+++ & ---*/ });
869  tasks.run([&] { kernels[1].computeVoxelSlices([](const Coord &a){ return a[0]+a[1]-a[2]; });/*++- & --+*/ });
870  tasks.run([&] { kernels[2].computeVoxelSlices([](const Coord &a){ return a[0]-a[1]+a[2]; });/*+-+ & -+-*/ });
871  tasks.run([&] { kernels[3].computeVoxelSlices([](const Coord &a){ return a[0]-a[1]-a[2]; });/*+-- & -++*/ });
872  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  for (int i = 0; i < nIter; ++i) {
881  for (SweepingKernel& kernel : kernels) kernel.sweep();
882  }
883 
884  if (finalize) {
885 #ifdef BENCHMARK_FAST_SWEEPING
886  util::CpuTimer timer("Computing extrema values");
887 #endif
888  MinMaxKernel kernel;
889  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  changeAsymmetricLevelSetBackground(mSdfGrid->tree(), e.max(), e.min());//multi-threaded
896 
897 #ifdef BENCHMARK_FAST_SWEEPING
898  timer.stop();
899 #endif
900  }
901 }// 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 {
910  using LeafRange = typename LeafMgr::LeafRange;
911  MinMaxKernel() : mMin(std::numeric_limits<SdfValueT>::max()), mMax(-mMin) {}
912  MinMaxKernel(MinMaxKernel& other, tbb::split) : mMin(other.mMin), mMax(other.mMax) {}
913 
914  math::MinMax<SdfValueT> run(const SdfGridT &grid)
915  {
916  LeafMgr mgr(grid.tree());// super fast
917  tbb::parallel_reduce(mgr.leafRange(), *this);
918  return math::MinMax<SdfValueT>(mMin, mMax);
919  }
920 
921  void operator()(const LeafRange& r)
922  {
923  for (auto leafIter = r.begin(); leafIter; ++leafIter) {
924  for (auto voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
925  const SdfValueT v = *voxelIter;
926  if (v < mMin) mMin = v;
927  if (v > mMax) mMax = v;
928  }
929  }
930  }
931 
932  void join(const MinMaxKernel& other)
933  {
934  if (other.mMin < mMin) mMin = other.mMin;
935  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 {
949  : mParent(&parent),
950  mBackground(parent.mSdfGrid->background())
951  {
952  mSdfGridInput = mParent->mSdfGrid->deepCopy();
953  }
954  DilateKernel(const DilateKernel &parent) = default;// for tbb::parallel_for
955  DilateKernel& operator=(const DilateKernel&) = delete;
956 
957  void run(int dilation, NearestNeighbors nn)
958  {
959 #ifdef BENCHMARK_FAST_SWEEPING
960  util::CpuTimer timer("Construct LeafManager");
961 #endif
962  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  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  for (int i=0, d = dilation/delta; i<d; ++i) dilateActiveValues(mgr, delta, nn, IGNORE_TILES);
977  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  mParent->mSweepMask.clear();
986  mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
987 
989  using LeafT = typename SdfGridT::TreeType::LeafNodeType;
990 
991  const FastSweepingDomain mode = mParent->mSweepDirection;
992 
993  LeafManagerT leafManager(mParent->mSdfGrid->tree());
994 
995  auto kernel = [&](LeafT& leaf, size_t /*leafIdx*/) {
996  static const SdfValueT Unknown = std::numeric_limits<SdfValueT>::max();
997  const SdfValueT background = mBackground;//local copy
998  auto* maskLeaf = mParent->mSweepMask.probeLeaf(leaf.origin());
999  SdfConstAccT sdfInputAcc(mSdfGridInput->tree());
1000  assert(maskLeaf);
1001  for (auto voxelIter = leaf.beginValueOn(); voxelIter; ++voxelIter) {
1002  const SdfValueT value = *voxelIter;
1003  SdfValueT inputValue;
1004  const Coord ijk = voxelIter.getCoord();
1005 
1006  if (math::Abs(value) < background) {// disable boundary voxels from the mask tree
1007  maskLeaf->setValueOff(voxelIter.pos());
1008  } else {
1009  switch (mode) {
1011  voxelIter.setValue(value > 0 ? Unknown : -Unknown);
1012  break;
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;
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  leafManager.foreach( kernel );
1037 
1038  // cache the leaf node origins for fast lookup in the sweeping kernels
1039  mParent->computeSweepMaskLeafOrigins();
1040 
1041 #ifdef BENCHMARK_FAST_SWEEPING
1042  timer.stop();
1043 #endif
1044  }// FastSweeping::DilateKernel::run
1045 
1046  // Private member data of DilateKernel
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 {
1058  InitSdf(FastSweeping &parent): mParent(&parent),
1059  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  void run(SdfValueT isoValue)
1064  {
1065  mIsoValue = isoValue;
1066  mAboveSign = mParent->mIsInputSdf ? SdfValueT(1) : SdfValueT(-1);
1067  SdfTreeT &tree = mSdfGrid->tree();//sdf
1068  const bool hasActiveTiles = tree.hasActiveTiles();
1069 
1070  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  mParent->mSweepMask.clear();
1078  mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1079 
1080  {// Process all voxels
1081  tree::LeafManager<SdfTreeT> mgr(tree, 1);// we need one auxiliary buffer
1082  tbb::parallel_for(mgr.leafRange(32), *this);//multi-threaded
1083  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  tree::NodeManager<SdfTreeT, SdfTreeT::RootNodeType::LEVEL-1> mgr(tree);
1091  mgr.foreachBottomUp(*this);//multi-threaded
1092  tree.root().setBackground(std::numeric_limits<SdfValueT>::max(), false);
1093  if (hasActiveTiles) tree.voxelizeActiveTiles();//multi-threaded
1094 
1095  // cache the leaf node origins for fast lookup in the sweeping kernels
1096 
1097  mParent->computeSweepMaskLeafOrigins();
1098  }// FastSweeping::InitSdf::run
1099 
1100  void operator()(const LeafRange& r) const
1101  {
1102  SweepMaskAccT sweepMaskAcc(mParent->mSweepMask);
1103  math::GradStencil<SdfGridT, false> stencil(*mSdfGrid);
1104  const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();//local copy
1105  const SdfValueT h = mAboveSign*static_cast<SdfValueT>(mSdfGrid->voxelSize()[0]);//Voxel size
1106  for (auto leafIter = r.begin(); leafIter; ++leafIter) {
1107  SdfValueT* sdf = leafIter.buffer(1).data();
1108  for (auto voxelIter = leafIter->beginValueAll(); voxelIter; ++voxelIter) {
1109  const SdfValueT value = *voxelIter;
1110  const bool isAbove = value > isoValue;
1111  if (!voxelIter.isValueOn()) {// inactive voxels
1112  sdf[voxelIter.pos()] = isAbove ? above : -above;
1113  } else {// active voxels
1114  const Coord ijk = voxelIter.getCoord();
1115  stencil.moveTo(ijk, value);
1116  const auto mask = stencil.intersectionMask( isoValue );
1117  if (mask.none()) {// most common case
1118  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  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  for (int i=0; i<6;) {
1128  SdfValueT d = std::numeric_limits<SdfValueT>::max(), d2;
1129  if (mask.test(i++)) d = math::Abs(delta/(value-stencil.getValue(i)));
1130  if (mask.test(i++)) {
1131  d2 = math::Abs(delta/(value-stencil.getValue(i)));
1132  if (d2 < d) d = d2;
1133  }
1134  if (d < std::numeric_limits<SdfValueT>::max()) sum += 1/(d*d);
1135  }
1136  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  }// FastSweeping::InitSdf::operator(const LeafRange&)
1143 
1144  template<typename RootOrInternalNodeT>
1145  void operator()(const RootOrInternalNodeT& node) const
1146  {
1147  const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();
1148  for (auto it = node.cbeginValueAll(); it; ++it) {
1149  SdfValueT& v = const_cast<SdfValueT&>(*it);
1150  v = v > isoValue ? above : -above;
1151  }//loop over all tiles
1152  }// FastSweeping::InitSdf::operator()(const RootOrInternalNodeT&)
1153 
1154  // Public member data
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 {
1168  using OpPoolT = tbb::enumerable_thread_specific<OpT>;
1169  InitExt(FastSweeping &parent) : mParent(&parent),
1170  mOpPool(nullptr), mSdfGrid(parent.mSdfGrid.get()),
1171  mExtGrid(parent.mExtGrid.get()), mIsoValue(0), mAboveSign(0) {}
1172  InitExt(const InitExt&) = default;// for tbb::parallel_for
1173  InitExt& operator=(const InitExt&) = delete;
1174  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  if (!mExtGrid) {
1178  OPENVDB_THROW(RuntimeError, "FastSweeping::InitExt expected an extension grid!");
1179  }
1180 
1181  mAboveSign = mParent->mIsInputSdf ? SdfValueT(1) : SdfValueT(-1);
1182  mIsoValue = isoValue;
1183  auto &tree1 = mSdfGrid->tree();
1184  auto &tree2 = mExtGrid->tree();
1185  const bool hasActiveTiles = tree1.hasActiveTiles();//very fast
1186 
1187  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  mParent->mSweepMask.clear();
1196  mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1197 
1198  {// Process all voxels
1199  // Define thread-local operators
1200  OpPoolT opPool(opPrototype);
1201  mOpPool = &opPool;
1202 
1203  tree::LeafManager<SdfTreeT> mgr(tree1, 1);// we need one auxiliary buffer
1204  tbb::parallel_for(mgr.leafRange(32), *this);//multi-threaded
1205  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  tree::NodeManager<SdfTreeT, SdfTreeT::RootNodeType::LEVEL-1> mgr(tree1);
1213  mgr.foreachBottomUp(*this);//multi-threaded
1214  tree1.root().setBackground(std::numeric_limits<SdfValueT>::max(), false);
1215  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  mParent->computeSweepMaskLeafOrigins();
1227 
1228 #ifdef BENCHMARK_FAST_SWEEPING
1229  timer.stop();
1230 #endif
1231  }// FastSweeping::InitExt::run
1232 
1233  // int implements an update if minD needs to be updated
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
1239  void sumHelper(ExtT& sum2, ExtT ext, bool /* update */, const SdfT& d2) const { sum2 += static_cast<ExtValueT>(d2 * ext); }// non-int implementation
1240 
1242  ExtT extValHelper(ExtT extSum, const SdfT& /* sdfSum */) const { return extSum; }// int implementation
1243 
1245  ExtT extValHelper(ExtT extSum, const SdfT& sdfSum) const {return ExtT((SdfT(1) / sdfSum) * extSum); }// non-int implementation
1246 
1247  void operator()(const LeafRange& r) const
1248  {
1249  ExtAccT acc(mExtGrid->tree());
1250  SweepMaskAccT sweepMaskAcc(mParent->mSweepMask);
1251  math::GradStencil<SdfGridT, false> stencil(*mSdfGrid);
1252  const math::Transform& xform = mExtGrid->transform();
1253  typename OpPoolT::reference op = mOpPool->local();
1254  const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();//local copy
1255  const SdfValueT h = mAboveSign*static_cast<SdfValueT>(mSdfGrid->voxelSize()[0]);//Voxel size
1256  for (auto leafIter = r.begin(); leafIter; ++leafIter) {
1257  SdfValueT *sdf = leafIter.buffer(1).data();
1258  ExtValueT *ext = acc.probeLeaf(leafIter->origin())->buffer().data();//should be safe!
1259  for (auto voxelIter = leafIter->beginValueAll(); voxelIter; ++voxelIter) {
1260  const SdfValueT value = *voxelIter;
1261  const bool isAbove = value > isoValue;
1262  if (!voxelIter.isValueOn()) {// inactive voxels
1263  sdf[voxelIter.pos()] = isAbove ? above : -above;
1264  } else {// active voxels
1265  const Coord ijk = voxelIter.getCoord();
1266  stencil.moveTo(ijk, value);
1267  const auto mask = stencil.intersectionMask( isoValue );
1268  if (mask.none()) {// no zero-crossing neighbors, most common case
1269  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  const SdfValueT delta = value - isoValue;//offset relative to iso-value
1275  if (math::isApproxZero(delta)) {//voxel is on the iso-surface
1276  sdf[voxelIter.pos()] = 0;
1277  ext[voxelIter.pos()] = ExtValueT(op(xform.indexToWorld(ijk)));
1278  } else {//voxel is neighboring the iso-surface
1279  SdfValueT sum1 = 0;
1280  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  for (int n=0, i=0; i<6;) {
1287  SdfValueT d = std::numeric_limits<SdfValueT>::max(), d2;
1288  if (mask.test(i++)) {
1289  d = math::Abs(delta/(value-stencil.getValue(i)));
1290  n = i - 1;
1291  }
1292  if (mask.test(i++)) {
1293  d2 = math::Abs(delta/(value-stencil.getValue(i)));
1294  if (d2 < d) {
1295  d = d2;
1296  n = i - 1;
1297  }
1298  }
1299  if (d < std::numeric_limits<SdfValueT>::max()) {
1300  d2 = 1/(d*d);
1301  sum1 += d2;
1302  const Vec3R xyz(static_cast<SdfValueT>(ijk[0])+d*static_cast<SdfValueT>(FastSweeping::mOffset[n][0]),
1303  static_cast<SdfValueT>(ijk[1])+d*static_cast<SdfValueT>(FastSweeping::mOffset[n][1]),
1304  static_cast<SdfValueT>(ijk[2])+d*static_cast<SdfValueT>(FastSweeping::mOffset[n][2]));
1305  // If current d is less than minD, update minD
1306  sumHelper(sum2, ExtValueT(op(xform.indexToWorld(xyz))), d < minD, d2);
1307  if (d < minD) minD = d;
1308  }
1309  }//look over six cases
1310  ext[voxelIter.pos()] = extValHelper(sum2, sum1);
1311  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  }// FastSweeping::InitExt::operator(const LeafRange& r)
1318 
1319  template<typename RootOrInternalNodeT>
1320  void operator()(const RootOrInternalNodeT& node) const
1321  {
1322  const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();
1323  for (auto it = node.cbeginValueAll(); it; ++it) {
1324  SdfValueT& v = const_cast<SdfValueT&>(*it);
1325  v = v > isoValue ? above : -above;
1326  }//loop over all tiles
1327  }
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 {
1343  MaskKernel(FastSweeping &parent) : mParent(&parent),
1344  mSdfGrid(parent.mSdfGrid.get()) {}
1345  MaskKernel(const MaskKernel &parent) = default;// for tbb::parallel_for
1346  MaskKernel& operator=(const MaskKernel&) = delete;
1347 
1348  void run(const MaskTreeT &mask)
1349  {
1350 #ifdef BENCHMARK_FAST_SWEEPING
1351  util::CpuTimer timer;
1352 #endif
1353  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  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  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  mParent->mSweepMask.clear();
1375  mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1376 
1377  using LeafManagerT = tree::LeafManager<SweepMaskTreeT>;
1378  using LeafT = typename SweepMaskTreeT::LeafNodeType;
1379  LeafManagerT leafManager(mParent->mSweepMask);
1380 
1381  auto kernel = [&](LeafT& leaf, size_t /*leafIdx*/) {
1382  static const SdfValueT Unknown = std::numeric_limits<SdfValueT>::max();
1383  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  SdfValueT *data = acc.probeLeaf(leaf.origin())->buffer().data();
1387  for (auto voxelIter = leaf.beginValueOn(); voxelIter; ++voxelIter) {// mask voxels
1388  if (math::Abs( data[voxelIter.pos()] ) < Unknown ) {
1389  // disable boundary voxels from the mask tree
1390  voxelIter.setValue(false);
1391  }
1392  }
1393  };
1394  leafManager.foreach( kernel );
1395 
1396  // cache the leaf node origins for fast lookup in the sweeping kernels
1397  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  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>
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  const SweepMaskTreeT& maskTree = mParent->mSweepMask;
1427 
1428  using LeafManagerT = typename tree::LeafManager<const SweepMaskTreeT>;
1429  using LeafT = typename SweepMaskTreeT::LeafNodeType;
1430  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  std::vector<int8_t> leafSliceMasks(leafManager.leafCount()*maskRange);
1442  auto kernel1 = [&](const LeafT& leaf, size_t leafIdx) {
1443  const size_t leafOffset = leafIdx * maskRange;
1444  for (auto voxelIter = leaf.cbeginValueOn(); voxelIter; ++voxelIter) {
1445  const Coord ijk = LeafT::offsetToLocalCoord(voxelIter.pos());
1446  leafSliceMasks[leafOffset + hash(ijk) + maskOffset] = uint8_t(1);
1447  }
1448  };
1449  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  tbb::enumerable_thread_specific<ThreadLocalMap> pool;
1456  auto kernel2 = [&](const LeafT& leaf, size_t leafIdx) {
1457  ThreadLocalMap& map = pool.local();
1458  const Coord& origin = leaf.origin();
1459  const int64_t leafKey = hash(origin);
1460  const size_t leafOffset = leafIdx * maskRange;
1461  for (int sliceIdx = 0; sliceIdx < maskRange; sliceIdx++) {
1462  if (leafSliceMasks[leafOffset + sliceIdx] == uint8_t(1)) {
1463  const int64_t voxelSliceKey = leafKey+sliceIdx-maskOffset;
1464  map[voxelSliceKey].emplace_back(leafIdx);
1465  }
1466  }
1467  };
1468  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  for (auto poolIt = pool.begin(); poolIt != pool.end(); ++poolIt) {
1474  const ThreadLocalMap& map = *poolIt;
1475  for (const auto& it : map) {
1476  for (const size_t leafIdx : it.second) {
1477  mVoxelSliceMap[it.first].emplace_back(leafIdx, NodeMaskPtrT());
1478  }
1479  }
1480  }
1481 
1482  // extract the voxel slice keys for random access into the map
1483  mVoxelSliceKeys.reserve(mVoxelSliceMap.size());
1484  for (const auto& it : mVoxelSliceMap) {
1485  mVoxelSliceKeys.push_back(it.first);
1486  }
1487 
1488  // allocate the node masks in parallel, as the map is populated in serial
1489  auto kernel3 = [&](tbb::blocked_range<size_t>& range) {
1490  for (size_t i = range.begin(); i < range.end(); i++) {
1491  const int64_t key = mVoxelSliceKeys[i];
1492  for (auto& it : mVoxelSliceMap[key]) {
1493  it.second = std::make_unique<NodeMaskT>();
1494  }
1495  }
1496  };
1497  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  auto kernel4 = [&](tbb::blocked_range<size_t>& range) {
1505  for (size_t i = range.begin(); i < range.end(); i++) {
1506  const int64_t voxelSliceKey = mVoxelSliceKeys[i];
1507  LeafSliceArray& leafSliceArray = mVoxelSliceMap[voxelSliceKey];
1508  for (LeafSlice& leafSlice : leafSliceArray) {
1509  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  const int64_t leafKey = hash(origin);
1514  for (auto voxelIter = leaf.cbeginValueOn(); voxelIter; ++voxelIter) {
1515  const Index voxelIdx = voxelIter.pos();
1516  const Coord ijk = LeafT::offsetToLocalCoord(voxelIdx);
1517  const int64_t key = leafKey + hash(ijk);
1518  if (key == voxelSliceKey) {
1519  nodeMask->setOn(voxelIdx);
1520  }
1521  }
1522  }
1523  }
1524  };
1525  tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceKeys.size()), kernel4);
1526  }// 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  inline static Coord ijk(const Coord &p, int i) { return p + FastSweeping::mOffset[i]; }
1533  NN() : v(), n() {}
1534  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  inline bool operator<(const NN &rhs) const { return v < rhs.v; }
1537  inline operator bool() const { return v < SdfValueT(1000); }
1538  };// NN
1539 
1540  /// @note Extending an integer field is based on the nearest-neighbor interpolation
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
1546  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
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
1558  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  return ExtT(w*(d1.v*v1 + d2.v*v2 + d3.v*v3));
1560  }// non-int implementation
1561 
1562  void sweep()
1563  {
1564  typename ExtGridT::TreeType *tree2 = mParent->mExtGrid ? &mParent->mExtGrid->tree() : nullptr;
1565  typename ExtGridT::TreeType *tree3 = mParent->mExtGridInput ? &mParent->mExtGridInput->tree() : nullptr;
1566 
1567  const SdfValueT h = static_cast<SdfValueT>(mParent->mSdfGrid->voxelSize()[0]);
1568  const SdfValueT sqrt2h = math::Sqrt(SdfValueT(2))*h;
1569  const FastSweepingDomain mode = mParent->mSweepDirection;
1570  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  if (tree2 && mode != FastSweepingDomain::SWEEP_ALL) assert(tree3);
1576 
1577  const std::vector<Coord>& leafNodeOrigins = mParent->mSweepMaskLeafOrigins;
1578 
1579  int64_t voxelSliceIndex(0);
1580 
1581  auto kernel = [&](const tbb::blocked_range<size_t>& range) {
1582  using LeafT = typename SdfGridT::TreeType::LeafNodeType;
1583 
1584  SdfAccT acc1(mParent->mSdfGrid->tree());
1585  auto acc2 = std::unique_ptr<ExtAccT>(tree2 ? new ExtAccT(*tree2) : nullptr);
1586  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  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  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  const size_t leafIdx = leafSlice.first;
1601  const NodeMaskPtrT& nodeMask = leafSlice.second;
1602 
1603  const Coord& origin = leafNodeOrigins[leafIdx];
1604 
1605  Coord ijk;
1606  for (auto indexIter = nodeMask->beginOn(); indexIter; ++indexIter) {
1607 
1608  // Get coordinate of center point of the FD stencil
1609  ijk = origin + LeafT::offsetToLocalCoord(indexIter.pos());
1610 
1611  // Find the closes neighbors in the three axial directions
1612  d1 = std::min(NN(acc1, ijk, 0), NN(acc1, ijk, 1));
1613  d2 = std::min(NN(acc1, ijk, 2), NN(acc1, ijk, 3));
1614  d3 = std::min(NN(acc1, ijk, 4), NN(acc1, ijk, 5));
1615 
1616  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  SdfValueT &value = const_cast<SdfValueT&>(acc1.getValue(ijk));
1622 
1623  // Extract the sign
1624  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  if (d2 < d1) std::swap(d1, d2);
1631  if (d3 < d2) std::swap(d2, d3);
1632  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  update = d1.v + h;
1638  if (update <= d2.v) {
1639  if (update < absV) {
1640  value = sign * update;
1641  if (acc2) {
1642  // There is an assert upstream to check if mExtGridInput exists if mode != SWEEP_ALL
1643  ExtValueT updateExt = acc2->getValue(d1(ijk));
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
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  acc2->setValue(ijk, updateExt);
1653  }//update ext?
1654  }//update sdf?
1655  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  if (d2.v <= sqrt2h + d1.v) {
1663  D = SdfValueT(2) * h * h - math::Pow2(d1.v - d2.v);// = 2h^2-(d1-d2)^2
1664  update = SdfValueT(0.5) * (d1.v + d2.v + std::sqrt(D));
1665  if (update > d2.v && update <= d3.v) {
1666  if (update < absV) {
1667  value = sign * update;
1668  if (acc2) {
1669  d1.v -= update;
1670  d2.v -= update;
1671  // affine combination of two neighboring extension values
1672  const SdfValueT w = SdfValueT(1)/(d1.v+d2.v);
1673  const ExtValueT v1 = acc2->getValue(d1(ijk));
1674  const ExtValueT v2 = acc2->getValue(d2(ijk));
1675  const ExtValueT extVal = twoNghbr(d1, d2, w, v1, v2);
1676 
1677  ExtValueT updateExt = extVal;
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
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  acc2->setValue(ijk, updateExt);
1687  }//update ext?
1688  }//update sdf?
1689  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  const SdfValueT d123 = d1.v + d2.v + d3.v;
1698  D = d123*d123 - SdfValueT(3)*(d1.v*d1.v + d2.v*d2.v + d3.v*d3.v - h * h);
1699  if (D >= SdfValueT(0)) {// non-negative discriminant
1700  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  if (update < absV) {
1703  value = sign * update;
1704  if (acc2) {
1705  d1.v -= update;
1706  d2.v -= update;
1707  d3.v -= update;
1708  // affine combination of three neighboring extension values
1709  const SdfValueT w = SdfValueT(1)/(d1.v+d2.v+d3.v);
1710  const ExtValueT v1 = acc2->getValue(d1(ijk));
1711  const ExtValueT v2 = acc2->getValue(d2(ijk));
1712  const ExtValueT v3 = acc2->getValue(d3(ijk));
1713  const ExtValueT extVal = threeNghbr(d1, d2, d3, w, v1, v2, v3);
1714 
1715  ExtValueT updateExt = extVal;
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
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  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  for (size_t i = 0; i < mVoxelSliceKeys.size(); i++) {
1737  voxelSliceIndex = mVoxelSliceKeys[i];
1738  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  for (size_t i = mVoxelSliceKeys.size(); i > 0; i--) {
1745  voxelSliceIndex = mVoxelSliceKeys[i-1];
1746  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  }// 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 {
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 {
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 {
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 {
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 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 {
1836  if (fs.initExt(fogGrid, op, background, isoValue, /*isInputSdf*/false, mode, extGrid))
1837  fs.sweep(nIter, /*finalize*/true);
1838  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 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 {
1852  if (fs.initExt(sdfGrid, op, background, isoValue, /*isInputSdf*/true, mode, extGrid))
1853  fs.sweep(nIter, /*finalize*/true);
1854  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 {
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 maskSdf(const GridT &sdfGrid,
1873  const Grid<MaskTreeT> &mask,
1874  bool ignoreActiveTiles,
1875  int nIter)
1876 {
1878  if (fs.initMask(sdfGrid, mask, ignoreActiveTiles)) fs.sweep(nIter);
1879  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
1892 #endif
1893 
1894 #define _FUNCTION(TreeT) \
1895  Grid<TreeT>::Ptr fogToSdf(const Grid<TreeT>&, TreeT::ValueType, int)
1897 #undef _FUNCTION
1898 
1899 #define _FUNCTION(TreeT) \
1900  Grid<TreeT>::Ptr sdfToSdf(const Grid<TreeT>&, TreeT::ValueType, int)
1902 #undef _FUNCTION
1903 
1904 #define _FUNCTION(TreeT) \
1905  Grid<TreeT>::Ptr dilateSdf(const Grid<TreeT>&, int, NearestNeighbors, int, FastSweepingDomain)
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
ExtT threeNghbr(const NN &d1, const NN &d2, const NN &d3, const SdfT &w, const ExtT &v1, const ExtT &v2, const ExtT &v3) const
Definition: FastSweeping.h:1558
ExtT twoNghbr(const NN &d1, const NN &d2, const SdfT &w, const ExtT &v1, const ExtT &v2) const
Definition: FastSweeping.h:1546
~FastSweeping()
Destructor.
Definition: FastSweeping.h:483
LeafNodeT * probeLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains voxel (x, y, z), or nullptr if no such node exists...
Definition: ValueAccessor.h:379
Vec3d indexToWorld(const Vec3d &xyz) const
Apply this transformation to the given coordinates.
Definition: Transform.h:108
SdfGridT::template ValueConverter< ExtValueT >::Type::Ptr sdfToExt(const SdfGridT &sdfGrid, const OpT &op, const ExtValueT &background, typename SdfGridT::ValueType isoValue, int nIter, FastSweepingDomain mode, const typename SdfGridT::template ValueConverter< ExtValueT >::Type::ConstPtr extGrid)
Definition: FastSweeping.h:1811
SharedPtr< Grid > Ptr
Definition: Grid.h:579
LeafRange leafRange(size_t grainsize=1) const
Return a TBB-compatible LeafRange.
Definition: LeafManager.h:345
bool isInputSdf()
Return whether the fast-sweeping input grid a signed distance function or not (fog).
Definition: FastSweeping.h:675
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
Definition: Morphology.h:80
NearestNeighbors
Voxel topology of nearest neighbors.
Definition: Morphology.h:58
void sweep(int nIter=1, bool finalize=true)
Perform nIter iterations of the fast sweeping algorithm.
Definition: FastSweeping.h:842
typename LeafMgr::LeafRange LeafRange
Definition: FastSweeping.h:910
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
SdfGridT::template ValueConverter< ExtValueT >::Type::Ptr sdfToExt(const SdfGridT &sdfGrid, const ExtOpT &op, const ExtValueT &background, typename SdfGridT::ValueType isoValue=0, int nIter=1, FastSweepingDomain mode=FastSweepingDomain::SWEEP_ALL, const typename SdfGridT::template ValueConverter< ExtValueT >::Type::ConstPtr extGrid=nullptr)
Computes the extension of a field (scalar, vector, or int are supported), defined by the specified fu...
bool initSdf(const SdfGridT &sdfGrid, SdfValueT isoValue, bool isInputSdf)
Initializer for input grids that are either a signed distance field or a scalar fog volume...
Definition: FastSweeping.h:762
Templated class to compute the minimum and maximum values.
Definition: Stats.h:30
std::pair< typename SdfGridT::Ptr, typename SdfGridT::template ValueConverter< ExtValueT >::Type::Ptr > sdfToSdfAndExt(const SdfGridT &sdfGrid, const ExtOpT &op, const ExtValueT &background, typename SdfGridT::ValueType isoValue=0, int nIter=1, FastSweepingDomain mode=FastSweepingDomain::SWEEP_ALL, const typename SdfGridT::template ValueConverter< ExtValueT >::Type::ConstPtr extGrid=nullptr)
Computes the signed distance field and the extension of a field (scalar, vector, or int are supported...
Definition: FastSweeping.h:1843
const SdfValueT mBackground
Definition: FastSweeping.h:1048
typename tree::LeafManager< SdfTreeT >::LeafRange LeafRange
Definition: FastSweeping.h:1057
std::bitset< 6 > intersectionMask(const ValueType &isoValue=zeroVal< ValueType >()) const
Return true a bit-mask where the 6 bits indicates if the center of the stencil intersects the iso-con...
Definition: Stencils.h:188
Private class of FastSweeping to perform concurrent fast sweeping in two directions.
Definition: FastSweeping.h:1411
GridT::Ptr sdfToSdf(const GridT &sdfGrid, typename GridT::ValueType isoValue=0, int nIter=1)
Given an existing approximate SDF it solves the Eikonal equation for all its active voxels...
Definition: FastSweeping.h:1784
void changeAsymmetricLevelSetBackground(TreeOrLeafManagerT &tree, const typename TreeOrLeafManagerT::ValueType &outsideWidth, const typename TreeOrLeafManagerT::ValueType &insideWidth, bool threaded=true, size_t grainSize=32)
Replace the background values in all the nodes of a floating-point tree containing a possibly asymmet...
Definition: ChangeBackground.h:218
size_t MinIndex(const Vec3T &v)
Return the index [0,1,2] of the smallest value in a 3D vector.
Definition: Math.h:938
Functions to efficiently compute histograms, extrema (min/max) and statistics (mean, variance, etc.) of grid values.
bool isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance...
Definition: Math.h:350
bool swapLeafBuffer(size_t bufferIdx, bool serial=false)
Swap each leaf node&#39;s buffer with the nth corresponding auxiliary buffer, where n = bufferIdx...
Definition: LeafManager.h:359
Definition: Coord.h:586
Definition: Transform.h:39
void setValueOff(const Coord &xyz, const ValueType &value)
Set the value of the voxel at the given coordinates and mark the voxel as inactive.
Definition: ValueAccessor.h:266
MinMaxKernel(MinMaxKernel &other, tbb::split)
Definition: FastSweeping.h:912
GridT::Ptr dilateSdf(const GridT &sdfGrid, int dilation, NearestNeighbors nn=NN_FACE, int nIter=1, FastSweepingDomain mode=FastSweepingDomain::SWEEP_ALL)
Dilates the narrow band of an existing signed distance field by a specified number of voxels (like ad...
Definition: FastSweeping.h:1859
ExtT twoNghbr(const NN &d1, const NN &d2, const SdfT &, const ExtT &v1, const ExtT &v2) const
Definition: FastSweeping.h:1542
math::MinMax< SdfValueT > run(const SdfGridT &grid)
Definition: FastSweeping.h:914
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:644
void dilateActiveValues(TreeOrLeafManagerT &tree, const int iterations=1, const NearestNeighbors nn=NN_FACE, const TilePolicy mode=PRESERVE_TILES, const bool threaded=true)
Topologically dilate all active values (i.e. both voxels and tiles) in a tree using one of three near...
Definition: Morphology.h:1055
Update all voxels affected by the sweeping algorithm.
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:107
DilateKernel(FastSweeping &parent)
Definition: FastSweeping.h:948
SdfValueT v
Definition: FastSweeping.h:1530
SdfValueT mIsoValue
Definition: FastSweeping.h:1157
static Coord ijk(const Coord &p, int i)
Definition: FastSweeping.h:1532
Definition: LeafManager.h:101
void clear()
Clears all the grids and counters so initialization can be called again.
Definition: FastSweeping.h:721
Computes signed distance values from an initial iso-surface and optionally performs velocity extensio...
Definition: FastSweeping.h:458
TreeType & tree()
Return a reference to this grid&#39;s tree, which might be shared with other grids.
Definition: Grid.h:917
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:764
std::pair< typename FogGridT::Ptr, typename FogGridT::template ValueConverter< ExtValueT >::Type::Ptr > fogToSdfAndExt(const FogGridT &fogGrid, const ExtOpT &op, const ExtValueT &background, typename FogGridT::ValueType isoValue, int nIter=1, FastSweepingDomain mode=FastSweepingDomain::SWEEP_ALL, const typename FogGridT::template ValueConverter< ExtValueT >::Type::ConstPtr extGrid=nullptr)
Computes the signed distance field and the extension of a field (scalar, vector, or int are supported...
Definition: FastSweeping.h:1827
Coord operator()(const Coord &p) const
Definition: FastSweeping.h:1535
FastSweeping & operator=(const FastSweeping &)=delete
Disallow copy assignment.
GridT::Ptr maskSdf(const GridT &sdfGrid, const Grid< MaskTreeT > &mask, bool ignoreActiveTiles=false, int nIter=1)
Fills mask by extending an existing signed distance field into the active values of this input ree of...
Definition: FastSweeping.h:1872
SdfValueT mAboveSign
Definition: FastSweeping.h:1158
To facilitate threading over the nodes of a tree, cache node pointers in linear arrays, one for each level of the tree.
Definition: NodeManager.h:30
size_t boundaryVoxelCount() const
Return the number of voxels that defined the boundary condition.
Definition: FastSweeping.h:661
Implementation of morphological dilation and erosion.
bool isValid() const
Return true if there are voxels and boundaries to solve for.
Definition: FastSweeping.h:664
Definition: Exceptions.h:13
GridOrTreeType::template ValueConverter< bool >::Type::Ptr sdfInteriorMask(const GridOrTreeType &volume, typename GridOrTreeType::ValueType isovalue=lsutilGridZero< GridOrTreeType >())
Threaded method to construct a boolean mask that represents interior regions in a signed distance fie...
Definition: LevelSetUtil.h:2275
void join(const MinMaxKernel &other)
Definition: FastSweeping.h:932
SweepingKernel(FastSweeping &parent)
Definition: FastSweeping.h:1413
ValueT value
Definition: GridBuilder.h:1287
Definition: Exceptions.h:63
FogGridT::template ValueConverter< ExtValueT >::Type::Ptr fogToExt(const FogGridT &fogGrid, const ExtOpT &op, const ExtValueT &background, typename FogGridT::ValueType isoValue, int nIter=1, FastSweepingDomain mode=FastSweepingDomain::SWEEP_ALL, const typename FogGridT::template ValueConverter< ExtValueT >::Type::ConstPtr extGrid=nullptr)
Computes the extension of a field (scalar, vector, or int are supported), defined by the specified fu...
Definition: FastSweeping.h:1795
typename tree::LeafManager< SdfTreeT >::LeafRange LeafRange
Definition: FastSweeping.h:947
void run(SdfValueT isoValue)
Definition: FastSweeping.h:1063
bool initMask(const SdfGridT &sdfGrid, const Grid< MaskTreeT > &mask, bool ignoreActiveTiles=false)
Initializer used for the extension of an existing signed distance field into the active values of an ...
Definition: FastSweeping.h:811
void pruneInactive(TreeT &tree, bool threaded=true, size_t grainSize=1)
Reduce the memory footprint of a tree by replacing with background tiles any nodes whose values are a...
Definition: Prune.h:355
This class manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional auxil...
Definition: LeafManager.h:84
FastSweepingDomain sweepDirection() const
Return whether the sweep update is in all direction (SWEEP_ALL), greater than isovalue (SWEEP_GREATER...
Definition: FastSweeping.h:672
Index32 Index
Definition: Types.h:54
NN(const SdfAccT &a, const Coord &p, int i)
Definition: FastSweeping.h:1534
__hostdev__ uint32_t hash(uint32_t x)
Definition: common.h:13
void computeVoxelSlices(HashOp hash)
Main method that performs concurrent bi-directional sweeps.
Definition: FastSweeping.h:1419
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:28
SdfGridT::ConstPtr mSdfGridInput
Definition: FastSweeping.h:1049
bool initDilate(const SdfGridT &sdfGrid, int dilation, NearestNeighbors nn=NN_FACE, FastSweepingDomain mode=FastSweepingDomain::SWEEP_ALL)
Initializer used when dilating an existing signed distance field.
Definition: FastSweeping.h:799
bool operator<(const NN &rhs) const
Definition: FastSweeping.h:1536
double restart()
Re-start timer.
Definition: CpuTimer.h:150
GridClass getGridClass() const
Return the class of volumetric data (level set, fog volume, etc.) that is stored in this grid...
bool initExt(const SdfGridT &sdfGrid, const ExtOpT &op, const ExtValueT &background, SdfValueT isoValue, bool isInputSdf, FastSweepingDomain mode=FastSweepingDomain::SWEEP_ALL, const typename ExtGridT::ConstPtr extGrid=nullptr)
Initializer used whenever velocity extension is performed in addition to the computation of signed di...
SdfGridT::Ptr sdfGrid()
Returns a shared pointer to the signed distance field computed by this class.
Definition: FastSweeping.h:497
size_t sweepingVoxelCount() const
Return the number of voxels that will be solved for.
Definition: FastSweeping.h:658
InitSdf(FastSweeping &parent)
Definition: FastSweeping.h:1058
double stop() const
Returns and prints time in milliseconds since construction or start was called.
Definition: CpuTimer.h:128
Definition: Types.h:416
FastSweeping * mParent
Definition: FastSweeping.h:1047
void sweep()
Definition: FastSweeping.h:1562
GridT::Ptr fogToSdf(const GridT &fogGrid, typename GridT::ValueType isoValue, int nIter=1)
Converts a scalar fog volume into a signed distance function. Active input voxels with scalar values ...
Definition: FastSweeping.h:1773
void run(const char *ax, openvdb::GridBase &grid, const AttributeBindings &bindings={})
Run a full AX pipeline (parse, compile and execute) on a single OpenVDB Grid.
void operator()(const RootOrInternalNodeT &node) const
Definition: FastSweeping.h:1145
Miscellaneous utility methods that operate primarily or exclusively on level set grids.
MinMaxKernel()
Definition: FastSweeping.h:911
void run(int dilation, NearestNeighbors nn)
Definition: FastSweeping.h:957
void operator()(const LeafRange &r)
Definition: FastSweeping.h:921
#define OPENVDB_REAL_TREE_INSTANTIATE(Function)
Definition: version.h.in:147
void operator()(const LeafRange &r) const
Definition: FastSweeping.h:1100
ExtGridT::Ptr extGrid()
Returns a shared pointer to the extension field computed by this class.
Definition: FastSweeping.h:505
Definition: Stencils.h:1231
Definition: Morphology.h:58
Coord Abs(const Coord &xyz)
Definition: Coord.h:514
Type Pow2(Type x)
Return x2.
Definition: Math.h:551
SdfValueT mMax
Definition: FastSweeping.h:938
ExtGridT::Ptr extGridInput()
Returns a shared pointer to the extension grid input. This is non-NULL if this class is used to exten...
Definition: FastSweeping.h:513
Simple timer for basic profiling.
Definition: CpuTimer.h:66
A LeafManager manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional au...
MeshToVoxelEdgeData::EdgeData Abs(const MeshToVoxelEdgeData::EdgeData &x)
Definition: MeshToVolume.h:3707
const ValueType & getValue(unsigned int pos=0) const
Return the value from the stencil buffer with linear offset pos.
Definition: Stencils.h:97
Definition: ValueAccessor.h:182
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:116
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:103
Private class of FastSweeping to perform multi-threaded initialization.
Definition: FastSweeping.h:945
ExtT threeNghbr(const NN &d1, const NN &d2, const NN &d3, const SdfT &, const ExtT &v1, const ExtT &v2, const ExtT &v3) const
Definition: FastSweeping.h:1550
void changeLevelSetBackground(TreeOrLeafManagerT &tree, const typename TreeOrLeafManagerT::ValueType &halfWidth, bool threaded=true, size_t grainSize=32)
Replace the background value in all the nodes of a floating-point tree containing a symmetric narrow-...
Definition: ChangeBackground.h:234
SdfGridT * mSdfGrid
Definition: FastSweeping.h:1156
Definition: FastSweeping.h:1055
SdfValueT mMin
Definition: FastSweeping.h:938
math::Transform & transform()
Return a reference to this grid&#39;s transform, which might be shared with other grids.
Definition: Grid.h:415
void moveTo(const Coord &ijk)
Initialize the stencil buffer with the values of voxel (i, j, k) and its neighbors.
Definition: Stencils.h:47
FastSweeping * mParent
Definition: FastSweeping.h:1155
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:202
FastSweepingDomain
Fast Sweeping update mode. This is useful to determine narrow-band extension or field extension in on...
Definition: FastSweeping.h:64