OpenVDB  9.0.0
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 an existing signed distance field by a specified number of voxels
385 ///
386 /// @return A shared pointer to the dilated signed distance field.
387 ///
388 /// @param sdfGrid Input signed distance field to be dilated.
389 ///
390 /// @param dilation Numer of voxels that the input SDF will be dilated.
391 ///
392 /// @param nn Stencil-pattern used for dilation
393 ///
394 /// @param nIter Number of iterations of the fast sweeping algorithm.
395 /// Each iteration performs 2^3 = 8 individual sweeps.
396 ///
397 /// @param mode Determines the direction of the dilation. SWEEP_ALL
398 /// will dilate in both sides of the signed distance function,
399 /// SWEEP_GREATER_THAN_ISOVALUE will dilate in the positive
400 /// side of the iso-surface, SWEEP_LESS_THAN_ISOVALUE will dilate
401 /// in the negative side of the iso-surface.
402 ///
403 /// @details Topology will change as a result of this dilation. E.g. if
404 /// sdfGrid has a width of 3 and @a dilation = 6 then the grid
405 /// returned by this method is a narrow band signed distance field
406 /// with a total width of 9 units.
407 template<typename GridT>
408 typename GridT::Ptr
409 dilateSdf(const GridT &sdfGrid,
410  int dilation,
412  int nIter = 1,
413  FastSweepingDomain mode = FastSweepingDomain::SWEEP_ALL);
414 
415 /// @brief Fills mask by extending an existing signed distance field into
416 /// the active values of this input ree of arbitrary value type.
417 ///
418 /// @return A shared pointer to the masked signed distance field.
419 ///
420 /// @param sdfGrid Input signed distance field to be extended into the mask.
421 ///
422 /// @param mask Mask used to identify the topology of the output SDF.
423 /// Note this mask is assume to overlap with the sdfGrid.
424 ///
425 /// @param ignoreActiveTiles If false, active tiles in the mask are treated
426 /// as active voxels. Else they are ignored.
427 ///
428 /// @param nIter Number of iterations of the fast sweeping algorithm.
429 /// Each iteration performs 2^3 = 8 individual sweeps.
430 ///
431 /// @details Topology of the output SDF is determined by the union of the active
432 /// voxels (or optionally values) in @a sdfGrid and @a mask.
433 template<typename GridT, typename MaskTreeT>
434 typename GridT::Ptr
435 maskSdf(const GridT &sdfGrid,
436  const Grid<MaskTreeT> &mask,
437  bool ignoreActiveTiles = false,
438  int nIter = 1);
439 
440 ////////////////////////////////////////////////////////////////////////////////
441 /// @brief Computes signed distance values from an initial iso-surface and
442 /// optionally performs velocity extension at the same time. This is
443 /// done by means of a novel sparse and parallel fast sweeping
444 /// algorithm based on a first order Godunov's scheme.
445 ///
446 /// Solves: @f$|\nabla \phi|^2 = 1 @f$
447 ///
448 /// @warning Note, it is important to call one of the initialization methods before
449 /// called the sweep function. Failure to do so will throw a RuntimeError.
450 /// Consider instead call one of the many higher-level free-standing functions
451 /// defined above!
452 template<typename SdfGridT, typename ExtValueT = typename SdfGridT::ValueType>
454 {
456  "FastSweeping requires SdfGridT to have floating-point values");
457  // Defined types related to the signed distance (or fog) grid
458  using SdfValueT = typename SdfGridT::ValueType;
459  using SdfTreeT = typename SdfGridT::TreeType;
460  using SdfAccT = tree::ValueAccessor<SdfTreeT, false>;//don't register accessors
461  using SdfConstAccT = typename tree::ValueAccessor<const SdfTreeT, false>;//don't register accessors
462 
463  // define types related to the extension field
464  using ExtGridT = typename SdfGridT::template ValueConverter<ExtValueT>::Type;
465  using ExtTreeT = typename ExtGridT::TreeType;
467 
468  // define types related to the tree that masks out the active voxels to be solved for
469  using SweepMaskTreeT = typename SdfTreeT::template ValueConverter<ValueMask>::Type;
470  using SweepMaskAccT = tree::ValueAccessor<SweepMaskTreeT, false>;//don't register accessors
471 
472 public:
473 
474  /// @brief Constructor
475  FastSweeping();
476 
477  /// @brief Destructor.
478  ~FastSweeping() { this->clear(); }
479 
480  /// @brief Disallow copy construction.
481  FastSweeping(const FastSweeping&) = delete;
482 
483  /// @brief Disallow copy assignment.
485 
486  /// @brief Returns a shared pointer to the signed distance field computed
487  /// by this class.
488  ///
489  /// @warning This shared pointer might point to NULL if the grid has not been
490  /// initialize (by one of the init methods) or computed (by the sweep
491  /// method).
492  typename SdfGridT::Ptr sdfGrid() { return mSdfGrid; }
493 
494  /// @brief Returns a shared pointer to the extension field computed
495  /// by this class.
496  ///
497  /// @warning This shared pointer might point to NULL if the grid has not been
498  /// initialize (by one of the init methods) or computed (by the sweep
499  /// method).
500  typename ExtGridT::Ptr extGrid() { return mExtGrid; }
501 
502  /// @brief Returns a shared pointer to the extension grid input. This is non-NULL
503  /// if this class is used to extend a field with a non-default sweep direction.
504  ///
505  /// @warning This shared pointer might point to NULL. This is non-NULL
506  /// if this class is used to extend a field with a non-default sweep direction,
507  /// i.e. SWEEP_LESS_THAN_ISOVALUE or SWEEP_GREATER_THAN_ISOVALUE.
508  typename ExtGridT::Ptr extGridInput() { return mExtGridInput; }
509 
510  /// @brief Initializer for input grids that are either a signed distance
511  /// field or a scalar fog volume.
512  ///
513  /// @return True if the initialization succeeded.
514  ///
515  /// @param sdfGrid Input scalar grid that represents an existing signed distance
516  /// field or a fog volume (signified by @a isInputSdf).
517  ///
518  /// @param isoValue Iso-value to be used to define the Dirichlet boundary condition
519  /// of the fast sweeping algorithm (typically 0 for sdfs and a
520  /// positive value for fog volumes).
521  ///
522  /// @param isInputSdf Used to determine if @a sdfGrid is a sigend distance field (true)
523  /// or a scalar fog volume (false).
524  ///
525  /// @details This, or any of ther other initialization methods, should be called
526  /// before any call to sweep(). Failure to do so will throw a RuntimeError.
527  ///
528  /// @warning Note, if this method fails, i.e. returns false, a subsequent call
529  /// to sweep will trow a RuntimeError. Instead call clear and try again.
530  bool initSdf(const SdfGridT &sdfGrid, SdfValueT isoValue, bool isInputSdf);
531 
532  /// @brief Initializer used whenever velocity extension is performed in addition
533  /// to the computation of signed distance fields.
534  ///
535  /// @return True if the initialization succeeded.
536  ///
537  ///
538  /// @param sdfGrid Input scalar grid that represents an existing signed distance
539  /// field or a fog volume (signified by @a isInputSdf).
540  ///
541  /// @param op Functor with signature [](const Vec3R &xyz)->ExtValueT that
542  /// defines the Dirichlet boundary condition, on the iso-surface,
543  /// of the field to be extended. Strictly the return type of this functor
544  /// is only required to be convertible to ExtValueT!
545  ///
546  /// @param background Background value of return grid with the extension field.
547  ///
548  /// @param isoValue Iso-value to be used for the boundary condition of the fast
549  /// sweeping algorithm (typically 0 for sdfs and a positive value
550  /// for fog volumes).
551  ///
552  /// @param isInputSdf Used to determine if @a sdfGrid is a sigend distance field (true)
553  /// or a scalar fog volume (false).
554  ///
555  /// @param mode Determines the mode of updating the extension field. SWEEP_ALL
556  /// will update all voxels of the extension field affected by the
557  /// fast sweeping algorithm. SWEEP_GREATER_THAN_ISOVALUE will update
558  /// all voxels corresponding to fog values that are greater than a given
559  /// isovalue. SWEEP_LESS_THAN_ISOVALUE will update all voxels corresponding
560  /// to fog values that are less than a given isovalue. If a mode other
561  /// than SWEEP_ALL is chosen, a user needs to supply @a extGrid.
562  ///
563  /// @param extGrid Optional parameter required to supply a default value for the extension
564  /// field when SWEEP_GREATER_THAN_ISOVALUE or SWEEP_LESS_THAN_ISOVALUE
565  /// mode is picked for @a mode. When SWEEP_GREATER_THAN_ISOVALUE is supplied
566  /// as an argument for @a mode, the extension field voxel will default
567  /// to the value of the @a extGrid in that position if it corresponds to a level-set
568  /// value that is less than the isovalue. Otherwise, the extension
569  /// field voxel value will be computed by the Fast Sweeping algorithm.
570  /// The opposite convention is implemented when SWEEP_LESS_THAN_ISOVALUE
571  /// is supplied as an argument for @a mode.
572  ///
573  /// @details This, or any of ther other initialization methods, should be called
574  /// before any call to sweep(). Failure to do so will throw a RuntimeError.
575  ///
576  /// @warning Note, if this method fails, i.e. returns false, a subsequent call
577  /// to sweep will trow a RuntimeError. Instead call clear and try again.
578  template <typename ExtOpT>
579  bool initExt(const SdfGridT &sdfGrid,
580  const ExtOpT &op,
581  const ExtValueT &background,
582  SdfValueT isoValue,
583  bool isInputSdf,
584  FastSweepingDomain mode = FastSweepingDomain::SWEEP_ALL,
585  const typename ExtGridT::ConstPtr extGrid = nullptr);
586 
587  /// @brief Initializer used when dilating an existing signed distance field.
588  ///
589  /// @return True if the initialization succeeded.
590  ///
591  /// @param sdfGrid Input signed distance field to to be dilated.
592  ///
593  /// @param dilation Numer of voxels that the input SDF will be dilated.
594  ///
595  /// @param nn Stencil-pattern used for dilation
596  ///
597  /// @param mode Determines the direction of the dilation. SWEEP_ALL
598  /// will dilate in both sides of the signed distance function,
599  /// SWEEP_GREATER_THAN_ISOVALUE will dilate in the positive
600  /// side of the iso-surface, SWEEP_LESS_THAN_ISOVALUE will dilate
601  /// in the negative side of the iso-surface.
602  ///
603  /// @details This, or any of ther other initialization methods, should be called
604  /// before any call to sweep(). Failure to do so will throw a RuntimeError.
605  ///
606  /// @warning Note, if this method fails, i.e. returns false, a subsequent call
607  /// to sweep will trow a RuntimeError. Instead call clear and try again.
608  bool initDilate(const SdfGridT &sdfGrid,
609  int dilation,
611  FastSweepingDomain mode = FastSweepingDomain::SWEEP_ALL);
612 
613  /// @brief Initializer used for the extension of an existing signed distance field
614  /// into the active values of an input mask of arbitrary value type.
615  ///
616  /// @return True if the initialization succeeded.
617  ///
618  /// @param sdfGrid Input signed distance field to be extended into the mask.
619  ///
620  /// @param mask Mask used to identify the topology of the output SDF.
621  /// Note this mask is assume to overlap with the sdfGrid.
622  ///
623  /// @param ignoreActiveTiles If false, active tiles in the mask are treated
624  /// as active voxels. Else they are ignored.
625  ///
626  /// @details This, or any of ther other initialization methods, should be called
627  /// before any call to sweep(). Failure to do so will throw a RuntimeError.
628  ///
629  /// @warning Note, if this method fails, i.e. returns false, a subsequent call
630  /// to sweep will trow a RuntimeError. Instead call clear and try again.
631  template<typename MaskTreeT>
632  bool initMask(const SdfGridT &sdfGrid, const Grid<MaskTreeT> &mask, bool ignoreActiveTiles = false);
633 
634  /// @brief Perform @a nIter iterations of the fast sweeping algorithm.
635  ///
636  /// @param nIter Number of iterations of the fast sweeping algorithm.
637  /// Each iteration performs 2^3 = 8 individual sweeps.
638  ///
639  /// @param finalize If true the (possibly asymmetric) inside and outside values of the
640  /// resulting signed distance field are properly set. Unless you're
641  /// an expert this should remain true!
642  ///
643  /// @throw RuntimeError if sweepingVoxelCount() or boundaryVoxelCount() return zero.
644  /// This might happen if none of the initialization methods above were called
645  /// or if that initialization failed.
646  void sweep(int nIter = 1,
647  bool finalize = true);
648 
649  /// @brief Clears all the grids and counters so initialization can be called again.
650  void clear();
651 
652  /// @brief Return the number of voxels that will be solved for.
653  size_t sweepingVoxelCount() const { return mSweepingVoxelCount; }
654 
655  /// @brief Return the number of voxels that defined the boundary condition.
656  size_t boundaryVoxelCount() const { return mBoundaryVoxelCount; }
657 
658  /// @brief Return true if there are voxels and boundaries to solve for
659  bool isValid() const { return mSweepingVoxelCount > 0 && mBoundaryVoxelCount > 0; }
660 
661  /// @brief Return whether the sweep update is in all direction (SWEEP_ALL),
662  /// greater than isovalue (SWEEP_GREATER_THAN_ISOVALUE), or less than isovalue
663  /// (SWEEP_LESS_THAN_ISOVALUE).
664  ///
665  /// @note SWEEP_GREATER_THAN_ISOVALUE and SWEEP_LESS_THAN_ISOVALUE modes are used
666  /// in dilating the narrow-band of a levelset or in extending a field.
667  FastSweepingDomain sweepDirection() const { return mSweepDirection; }
668 
669  /// @brief Return whether the fast-sweeping input grid a signed distance function or not (fog).
670  bool isInputSdf() { return mIsInputSdf; }
671 
672 private:
673 
674  /// @brief Private method to prune the sweep mask and cache leaf origins.
675  void computeSweepMaskLeafOrigins();
676 
677  // Private utility classes
678  template<typename>
679  struct MaskKernel;// initialization to extend a SDF into a mask
680  template<typename>
681  struct InitExt;
682  struct InitSdf;
683  struct DilateKernel;// initialization to dilate a SDF
684  struct MinMaxKernel;
685  struct SweepingKernel;// performs the actual concurrent sparse fast sweeping
686 
687  // Define the topology (i.e. stencil) of the neighboring grid points
688  static const Coord mOffset[6];// = {{-1,0,0},{1,0,0},{0,-1,0},{0,1,0},{0,0,-1},{0,0,1}};
689 
690  // Private member data of FastSweeping
691  typename SdfGridT::Ptr mSdfGrid;
692  typename ExtGridT::Ptr mExtGrid;
693  typename ExtGridT::Ptr mExtGridInput; // optional: only used in extending a field in one direction
694  SweepMaskTreeT mSweepMask; // mask tree containing all non-boundary active voxels, in the case of dilation, does not include active voxel
695  std::vector<Coord> mSweepMaskLeafOrigins; // cache of leaf node origins for mask tree
696  size_t mSweepingVoxelCount, mBoundaryVoxelCount;
697  FastSweepingDomain mSweepDirection; // only used in dilate and extending a field
698  bool mIsInputSdf;
699 };// FastSweeping
700 
701 ////////////////////////////////////////////////////////////////////////////////
702 
703 // Static member data initialization
704 template <typename SdfGridT, typename ExtValueT>
705 const Coord FastSweeping<SdfGridT, ExtValueT>::mOffset[6] = {{-1,0,0},{1,0,0},
706  {0,-1,0},{0,1,0},
707  {0,0,-1},{0,0,1}};
708 
709 template <typename SdfGridT, typename ExtValueT>
711  : mSdfGrid(nullptr), mExtGrid(nullptr), mSweepingVoxelCount(0), mBoundaryVoxelCount(0), mSweepDirection(FastSweepingDomain::SWEEP_ALL), mIsInputSdf(true)
712 {
713 }
714 
715 template <typename SdfGridT, typename ExtValueT>
717 {
718  mSdfGrid.reset();
719  mExtGrid.reset();
720  mSweepMask.clear();
721  if (mExtGridInput) mExtGridInput.reset();
722  mSweepingVoxelCount = mBoundaryVoxelCount = 0;
723  mSweepDirection = FastSweepingDomain::SWEEP_ALL;
724  mIsInputSdf = true;
725 }
726 
727 template <typename SdfGridT, typename ExtValueT>
729 {
730  // replace any inactive leaf nodes with tiles and voxelize any active tiles
731 
732  pruneInactive(mSweepMask);
733  mSweepMask.voxelizeActiveTiles();
734 
735  using LeafManagerT = tree::LeafManager<SweepMaskTreeT>;
736  using LeafT = typename SweepMaskTreeT::LeafNodeType;
737  LeafManagerT leafManager(mSweepMask);
738 
739  mSweepMaskLeafOrigins.resize(leafManager.leafCount());
740  std::atomic<size_t> sweepingVoxelCount{0};
741  auto kernel = [&](const LeafT& leaf, size_t leafIdx) {
742  mSweepMaskLeafOrigins[leafIdx] = leaf.origin();
743  sweepingVoxelCount += leaf.onVoxelCount();
744  };
745  leafManager.foreach(kernel, /*threaded=*/true, /*grainsize=*/1024);
746 
747  mBoundaryVoxelCount = 0;
748  mSweepingVoxelCount = sweepingVoxelCount;
749  if (mSdfGrid) {
750  const size_t totalCount = mSdfGrid->constTree().activeVoxelCount();
751  assert( totalCount >= mSweepingVoxelCount );
752  mBoundaryVoxelCount = totalCount - mSweepingVoxelCount;
753  }
754 }// FastSweeping::computeSweepMaskLeafOrigins
755 
756 template <typename SdfGridT, typename ExtValueT>
757 bool FastSweeping<SdfGridT, ExtValueT>::initSdf(const SdfGridT &fogGrid, SdfValueT isoValue, bool isInputSdf)
758 {
759  this->clear();
760  mSdfGrid = fogGrid.deepCopy();//very fast
761  mIsInputSdf = isInputSdf;
762  InitSdf kernel(*this);
763  kernel.run(isoValue);
764  return this->isValid();
765 }
766 
767 template <typename SdfGridT, typename ExtValueT>
768 template <typename OpT>
769 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)
770 {
771  if (mode != FastSweepingDomain::SWEEP_ALL) {
772  if (!extGrid)
773  OPENVDB_THROW(RuntimeError, "FastSweeping::initExt Calling initExt with mode != SWEEP_ALL requires an extension grid!");
774  if (extGrid->transform() != fogGrid.transform())
775  OPENVDB_THROW(RuntimeError, "FastSweeping::initExt extension grid input should have the same transform as Fog/SDF grid!");
776  }
777 
778  this->clear();
779  mSdfGrid = fogGrid.deepCopy();//very fast
780  mExtGrid = createGrid<ExtGridT>( background );
781  mSweepDirection = mode;
782  mIsInputSdf = isInputSdf;
783  if (mSweepDirection != FastSweepingDomain::SWEEP_ALL) {
784  mExtGridInput = extGrid->deepCopy();
785  }
786  mExtGrid->topologyUnion( *mSdfGrid );//very fast
787  InitExt<OpT> kernel(*this);
788  kernel.run(isoValue, op);
789  return this->isValid();
790 }
791 
792 
793 template <typename SdfGridT, typename ExtValueT>
794 bool FastSweeping<SdfGridT, ExtValueT>::initDilate(const SdfGridT &sdfGrid, int dilate, NearestNeighbors nn, FastSweepingDomain mode)
795 {
796  this->clear();
797  mSdfGrid = sdfGrid.deepCopy();//very fast
798  mSweepDirection = mode;
799  DilateKernel kernel(*this);
800  kernel.run(dilate, nn);
801  return this->isValid();
802 }
803 
804 template <typename SdfGridT, typename ExtValueT>
805 template<typename MaskTreeT>
806 bool FastSweeping<SdfGridT, ExtValueT>::initMask(const SdfGridT &sdfGrid, const Grid<MaskTreeT> &mask, bool ignoreActiveTiles)
807 {
808  this->clear();
809  mSdfGrid = sdfGrid.deepCopy();//very fast
810 
811  if (mSdfGrid->transform() != mask.transform()) {
812  OPENVDB_THROW(RuntimeError, "FastSweeping: Mask not aligned with the grid!");
813  }
814 
815  if (mask.getGridClass() == GRID_LEVEL_SET) {
816  using T = typename MaskTreeT::template ValueConverter<bool>::Type;
817  typename Grid<T>::Ptr tmp = sdfInteriorMask(mask);//might have active tiles
818  tmp->tree().voxelizeActiveTiles();//multi-threaded
819  MaskKernel<T> kernel(*this);
820  kernel.run(tmp->tree());//multi-threaded
821  } else {
822  if (ignoreActiveTiles || !mask.tree().hasActiveTiles()) {
823  MaskKernel<MaskTreeT> kernel(*this);
824  kernel.run(mask.tree());//multi-threaded
825  } else {
826  using T = typename MaskTreeT::template ValueConverter<ValueMask>::Type;
827  T tmp(mask.tree(), false, TopologyCopy());//multi-threaded
828  tmp.voxelizeActiveTiles(true);//multi-threaded
829  MaskKernel<T> kernel(*this);
830  kernel.run(tmp);//multi-threaded
831  }
832  }
833  return this->isValid();
834 }// FastSweeping::initMask
835 
836 template <typename SdfGridT, typename ExtValueT>
837 void FastSweeping<SdfGridT, ExtValueT>::sweep(int nIter, bool finalize)
838 {
839  if (!mSdfGrid) {
840  OPENVDB_THROW(RuntimeError, "FastSweeping::sweep called before initialization!");
841  }
842  if (mExtGrid && mSweepDirection != FastSweepingDomain::SWEEP_ALL && !mExtGridInput) {
843  OPENVDB_THROW(RuntimeError, "FastSweeping: Trying to extend a field in one direction needs"
844  " a non-null reference extension grid input.");
845  }
846  if (this->boundaryVoxelCount() == 0) {
847  OPENVDB_THROW(RuntimeError, "FastSweeping: No boundary voxels found!");
848  } else if (this->sweepingVoxelCount() == 0) {
849  OPENVDB_THROW(RuntimeError, "FastSweeping: No computing voxels found!");
850  }
851 
852  // note: Sweeping kernel is non copy-constructible, so use a deque instead of a vector
853  std::deque<SweepingKernel> kernels;
854  for (int i = 0; i < 4; i++) kernels.emplace_back(*this);
855 
856  { // compute voxel slices
857 #ifdef BENCHMARK_FAST_SWEEPING
858  util::CpuTimer timer("Computing voxel slices");
859 #endif
860 
861  // Exploiting nested parallelism - all voxel slice data is precomputed
862  tbb::task_group tasks;
863  tasks.run([&] { kernels[0].computeVoxelSlices([](const Coord &a){ return a[0]+a[1]+a[2]; });/*+++ & ---*/ });
864  tasks.run([&] { kernels[1].computeVoxelSlices([](const Coord &a){ return a[0]+a[1]-a[2]; });/*++- & --+*/ });
865  tasks.run([&] { kernels[2].computeVoxelSlices([](const Coord &a){ return a[0]-a[1]+a[2]; });/*+-+ & -+-*/ });
866  tasks.run([&] { kernels[3].computeVoxelSlices([](const Coord &a){ return a[0]-a[1]-a[2]; });/*+-- & -++*/ });
867  tasks.wait();
868 
869 #ifdef BENCHMARK_FAST_SWEEPING
870  timer.stop();
871 #endif
872  }
873 
874  // perform nIter iterations of bi-directional sweeping in all directions
875  for (int i = 0; i < nIter; ++i) {
876  for (SweepingKernel& kernel : kernels) kernel.sweep();
877  }
878 
879  if (finalize) {
880 #ifdef BENCHMARK_FAST_SWEEPING
881  util::CpuTimer timer("Computing extrema values");
882 #endif
883  MinMaxKernel kernel;
884  auto e = kernel.run(*mSdfGrid);//multi-threaded
885  //auto e = extrema(mGrid->beginValueOn());// 100x slower!!!!
886 #ifdef BENCHMARK_FAST_SWEEPING
887  std::cerr << "Min = " << e.min() << " Max = " << e.max() << std::endl;
888  timer.restart("Changing asymmetric background value");
889 #endif
890  changeAsymmetricLevelSetBackground(mSdfGrid->tree(), e.max(), e.min());//multi-threaded
891 
892 #ifdef BENCHMARK_FAST_SWEEPING
893  timer.stop();
894 #endif
895  }
896 }// FastSweeping::sweep
897 
898 /// Private class of FastSweeping to quickly compute the extrema
899 /// values of the active voxels in the leaf nodes. Several orders
900 /// of magnitude faster than tools::extrema!
901 template <typename SdfGridT, typename ExtValueT>
902 struct FastSweeping<SdfGridT, ExtValueT>::MinMaxKernel
903 {
905  using LeafRange = typename LeafMgr::LeafRange;
906  MinMaxKernel() : mMin(std::numeric_limits<SdfValueT>::max()), mMax(-mMin) {}
907  MinMaxKernel(MinMaxKernel& other, tbb::split) : mMin(other.mMin), mMax(other.mMax) {}
908 
909  math::MinMax<SdfValueT> run(const SdfGridT &grid)
910  {
911  LeafMgr mgr(grid.tree());// super fast
912  tbb::parallel_reduce(mgr.leafRange(), *this);
913  return math::MinMax<SdfValueT>(mMin, mMax);
914  }
915 
916  void operator()(const LeafRange& r)
917  {
918  for (auto leafIter = r.begin(); leafIter; ++leafIter) {
919  for (auto voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
920  const SdfValueT v = *voxelIter;
921  if (v < mMin) mMin = v;
922  if (v > mMax) mMax = v;
923  }
924  }
925  }
926 
927  void join(const MinMaxKernel& other)
928  {
929  if (other.mMin < mMin) mMin = other.mMin;
930  if (other.mMax > mMax) mMax = other.mMax;
931  }
932 
933  SdfValueT mMin, mMax;
934 };// FastSweeping::MinMaxKernel
935 
936 ////////////////////////////////////////////////////////////////////////////////
937 
938 /// Private class of FastSweeping to perform multi-threaded initialization
939 template <typename SdfGridT, typename ExtValueT>
940 struct FastSweeping<SdfGridT, ExtValueT>::DilateKernel
941 {
944  : mParent(&parent),
945  mBackground(parent.mSdfGrid->background())
946  {
947  mSdfGridInput = mParent->mSdfGrid->deepCopy();
948  }
949  DilateKernel(const DilateKernel &parent) = default;// for tbb::parallel_for
951 
952  void run(int dilation, NearestNeighbors nn)
953  {
954 #ifdef BENCHMARK_FAST_SWEEPING
955  util::CpuTimer timer("Construct LeafManager");
956 #endif
957  tree::LeafManager<SdfTreeT> mgr(mParent->mSdfGrid->tree());// super fast
958 
959 #ifdef BENCHMARK_FAST_SWEEPING
960  timer.restart("Changing background value");
961 #endif
962  static const SdfValueT Unknown = std::numeric_limits<SdfValueT>::max();
963  changeLevelSetBackground(mgr, Unknown);//multi-threaded
964 
965  #ifdef BENCHMARK_FAST_SWEEPING
966  timer.restart("Dilating and updating mgr (parallel)");
967  //timer.restart("Dilating and updating mgr (serial)");
968 #endif
969 
970  const int delta = 5;
971  for (int i=0, d = dilation/delta; i<d; ++i) dilateActiveValues(mgr, delta, nn, IGNORE_TILES);
972  dilateActiveValues(mgr, dilation % delta, nn, IGNORE_TILES);
973  //for (int i=0, n=5, d=dilation/n; i<d; ++i) dilateActiveValues(mgr, n, nn, IGNORE_TILES);
974  //dilateVoxels(mgr, dilation, nn);
975 
976 #ifdef BENCHMARK_FAST_SWEEPING
977  timer.restart("Initializing grid and sweep mask");
978 #endif
979 
980  mParent->mSweepMask.clear();
981  mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
982 
984  using LeafT = typename SdfGridT::TreeType::LeafNodeType;
985 
986  const FastSweepingDomain mode = mParent->mSweepDirection;
987 
988  LeafManagerT leafManager(mParent->mSdfGrid->tree());
989 
990  auto kernel = [&](LeafT& leaf, size_t /*leafIdx*/) {
991  static const SdfValueT Unknown = std::numeric_limits<SdfValueT>::max();
992  const SdfValueT background = mBackground;//local copy
993  auto* maskLeaf = mParent->mSweepMask.probeLeaf(leaf.origin());
994  SdfConstAccT sdfInputAcc(mSdfGridInput->tree());
995  assert(maskLeaf);
996  for (auto voxelIter = leaf.beginValueOn(); voxelIter; ++voxelIter) {
997  const SdfValueT value = *voxelIter;
998  SdfValueT inputValue;
999  const Coord ijk = voxelIter.getCoord();
1000 
1001  if (math::Abs(value) < background) {// disable boundary voxels from the mask tree
1002  maskLeaf->setValueOff(voxelIter.pos());
1003  } else {
1004  switch (mode) {
1006  voxelIter.setValue(value > 0 ? Unknown : -Unknown);
1007  break;
1009  if (value > 0) voxelIter.setValue(Unknown);
1010  else {
1011  maskLeaf->setValueOff(voxelIter.pos());
1012  bool isInputOn = sdfInputAcc.probeValue(ijk, inputValue);
1013  if ( !isInputOn ) voxelIter.setValueOff();
1014  else voxelIter.setValue(inputValue);
1015  }
1016  break;
1018  if (value < 0) voxelIter.setValue(-Unknown);
1019  else {
1020  maskLeaf->setValueOff(voxelIter.pos());
1021  bool isInputOn = sdfInputAcc.probeValue(ijk, inputValue);
1022  if ( !isInputOn ) voxelIter.setValueOff();
1023  else voxelIter.setValue(inputValue);
1024  }
1025  break;
1026  }
1027  }
1028  }
1029  };
1030 
1031  leafManager.foreach( kernel );
1032 
1033  // cache the leaf node origins for fast lookup in the sweeping kernels
1034  mParent->computeSweepMaskLeafOrigins();
1035 
1036 #ifdef BENCHMARK_FAST_SWEEPING
1037  timer.stop();
1038 #endif
1039  }// FastSweeping::DilateKernel::run
1040 
1041  // Private member data of DilateKernel
1043  const SdfValueT mBackground;
1044  typename SdfGridT::ConstPtr mSdfGridInput;
1045 };// FastSweeping::DilateKernel
1046 
1047 ////////////////////////////////////////////////////////////////////////////////
1048 
1049 template <typename SdfGridT, typename ExtValueT>
1050 struct FastSweeping<SdfGridT, ExtValueT>::InitSdf
1051 {
1053  InitSdf(FastSweeping &parent): mParent(&parent),
1054  mSdfGrid(parent.mSdfGrid.get()), mIsoValue(0), mAboveSign(0) {}
1055  InitSdf(const InitSdf&) = default;// for tbb::parallel_for
1056  InitSdf& operator=(const InitSdf&) = delete;
1057 
1058  void run(SdfValueT isoValue)
1059  {
1060  mIsoValue = isoValue;
1061  mAboveSign = mParent->mIsInputSdf ? SdfValueT(1) : SdfValueT(-1);
1062  SdfTreeT &tree = mSdfGrid->tree();//sdf
1063  const bool hasActiveTiles = tree.hasActiveTiles();
1064 
1065  if (mParent->mIsInputSdf && hasActiveTiles) {
1066  OPENVDB_THROW(RuntimeError, "FastSweeping: A SDF should not have active tiles!");
1067  }
1068 
1069 #ifdef BENCHMARK_FAST_SWEEPING
1070  util::CpuTimer timer("Initialize voxels");
1071 #endif
1072  mParent->mSweepMask.clear();
1073  mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1074 
1075  {// Process all voxels
1076  tree::LeafManager<SdfTreeT> mgr(tree, 1);// we need one auxiliary buffer
1077  tbb::parallel_for(mgr.leafRange(32), *this);//multi-threaded
1078  mgr.swapLeafBuffer(1);//swap voxel values
1079  }
1080 
1081 #ifdef BENCHMARK_FAST_SWEEPING
1082  timer.restart("Initialize tiles - new");
1083 #endif
1084  // Process all tiles
1085  tree::NodeManager<SdfTreeT, SdfTreeT::RootNodeType::LEVEL-1> mgr(tree);
1086  mgr.foreachBottomUp(*this);//multi-threaded
1087  tree.root().setBackground(std::numeric_limits<SdfValueT>::max(), false);
1088  if (hasActiveTiles) tree.voxelizeActiveTiles();//multi-threaded
1089 
1090  // cache the leaf node origins for fast lookup in the sweeping kernels
1091 
1092  mParent->computeSweepMaskLeafOrigins();
1093  }// FastSweeping::InitSdf::run
1094 
1095  void operator()(const LeafRange& r) const
1096  {
1097  SweepMaskAccT sweepMaskAcc(mParent->mSweepMask);
1098  math::GradStencil<SdfGridT, false> stencil(*mSdfGrid);
1099  const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();//local copy
1100  const SdfValueT h = mAboveSign*static_cast<SdfValueT>(mSdfGrid->voxelSize()[0]);//Voxel size
1101  for (auto leafIter = r.begin(); leafIter; ++leafIter) {
1102  SdfValueT* sdf = leafIter.buffer(1).data();
1103  for (auto voxelIter = leafIter->beginValueAll(); voxelIter; ++voxelIter) {
1104  const SdfValueT value = *voxelIter;
1105  const bool isAbove = value > isoValue;
1106  if (!voxelIter.isValueOn()) {// inactive voxels
1107  sdf[voxelIter.pos()] = isAbove ? above : -above;
1108  } else {// active voxels
1109  const Coord ijk = voxelIter.getCoord();
1110  stencil.moveTo(ijk, value);
1111  const auto mask = stencil.intersectionMask( isoValue );
1112  if (mask.none()) {// most common case
1113  sdf[voxelIter.pos()] = isAbove ? above : -above;
1114  } else {// compute distance to iso-surface
1115  // disable boundary voxels from the mask tree
1116  sweepMaskAcc.setValueOff(ijk);
1117  const SdfValueT delta = value - isoValue;//offset relative to iso-value
1118  if (math::isApproxZero(delta)) {//voxel is on the iso-surface
1119  sdf[voxelIter.pos()] = 0;
1120  } else {//voxel is neighboring the iso-surface
1121  SdfValueT sum = 0;
1122  for (int i=0; i<6;) {
1123  SdfValueT d = std::numeric_limits<SdfValueT>::max(), d2;
1124  if (mask.test(i++)) d = math::Abs(delta/(value-stencil.getValue(i)));
1125  if (mask.test(i++)) {
1126  d2 = math::Abs(delta/(value-stencil.getValue(i)));
1127  if (d2 < d) d = d2;
1128  }
1129  if (d < std::numeric_limits<SdfValueT>::max()) sum += 1/(d*d);
1130  }
1131  sdf[voxelIter.pos()] = isAbove ? h / math::Sqrt(sum) : -h / math::Sqrt(sum);
1132  }// voxel is neighboring the iso-surface
1133  }// intersecting voxels
1134  }// active voxels
1135  }// loop over voxels
1136  }// loop over leaf nodes
1137  }// FastSweeping::InitSdf::operator(const LeafRange&)
1138 
1139  template<typename RootOrInternalNodeT>
1140  void operator()(const RootOrInternalNodeT& node) const
1141  {
1142  const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();
1143  for (auto it = node.cbeginValueAll(); it; ++it) {
1144  SdfValueT& v = const_cast<SdfValueT&>(*it);
1145  v = v > isoValue ? above : -above;
1146  }//loop over all tiles
1147  }// FastSweeping::InitSdf::operator()(const RootOrInternalNodeT&)
1148 
1149  // Public member data
1151  SdfGridT *mSdfGrid;//raw pointer, i.e. lock free
1152  SdfValueT mIsoValue;
1153  SdfValueT mAboveSign;//sign of distance values above the iso-value
1154 };// FastSweeping::InitSdf
1155 
1156 
1157 /// Private class of FastSweeping to perform multi-threaded initialization
1158 template <typename SdfGridT, typename ExtValueT>
1159 template <typename OpT>
1160 struct FastSweeping<SdfGridT, ExtValueT>::InitExt
1161 {
1162  using LeafRange = typename tree::LeafManager<SdfTreeT>::LeafRange;
1163  using OpPoolT = tbb::enumerable_thread_specific<OpT>;
1164  InitExt(FastSweeping &parent) : mParent(&parent),
1165  mOpPool(nullptr), mSdfGrid(parent.mSdfGrid.get()),
1166  mExtGrid(parent.mExtGrid.get()), mIsoValue(0), mAboveSign(0) {}
1167  InitExt(const InitExt&) = default;// for tbb::parallel_for
1168  InitExt& operator=(const InitExt&) = delete;
1169  void run(SdfValueT isoValue, const OpT &opPrototype)
1170  {
1171  static_assert(std::is_convertible<decltype(opPrototype(Vec3d(0))),ExtValueT>::value, "Invalid return type of functor");
1172  if (!mExtGrid) {
1173  OPENVDB_THROW(RuntimeError, "FastSweeping::InitExt expected an extension grid!");
1174  }
1175 
1176  mAboveSign = mParent->mIsInputSdf ? SdfValueT(1) : SdfValueT(-1);
1177  mIsoValue = isoValue;
1178  auto &tree1 = mSdfGrid->tree();
1179  auto &tree2 = mExtGrid->tree();
1180  const bool hasActiveTiles = tree1.hasActiveTiles();//very fast
1181 
1182  if (mParent->mIsInputSdf && hasActiveTiles) {
1183  OPENVDB_THROW(RuntimeError, "FastSweeping: A SDF should not have active tiles!");
1184  }
1185 
1186 #ifdef BENCHMARK_FAST_SWEEPING
1187  util::CpuTimer timer("Initialize voxels");
1188 #endif
1189 
1190  mParent->mSweepMask.clear();
1191  mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1192 
1193  {// Process all voxels
1194  // Define thread-local operators
1195  OpPoolT opPool(opPrototype);
1196  mOpPool = &opPool;
1197 
1198  tree::LeafManager<SdfTreeT> mgr(tree1, 1);// we need one auxiliary buffer
1199  tbb::parallel_for(mgr.leafRange(32), *this);//multi-threaded
1200  mgr.swapLeafBuffer(1);//swap out auxiliary buffer
1201  }
1202 
1203 #ifdef BENCHMARK_FAST_SWEEPING
1204  timer.restart("Initialize tiles");
1205 #endif
1206  {// Process all tiles
1207  tree::NodeManager<SdfTreeT, SdfTreeT::RootNodeType::LEVEL-1> mgr(tree1);
1208  mgr.foreachBottomUp(*this);//multi-threaded
1209  tree1.root().setBackground(std::numeric_limits<SdfValueT>::max(), false);
1210  if (hasActiveTiles) {
1211 #ifdef BENCHMARK_FAST_SWEEPING
1212  timer.restart("Voxelizing active tiles");
1213 #endif
1214  tree1.voxelizeActiveTiles();//multi-threaded
1215  tree2.voxelizeActiveTiles();//multi-threaded
1216  }
1217  }
1218 
1219  // cache the leaf node origins for fast lookup in the sweeping kernels
1220 
1221  mParent->computeSweepMaskLeafOrigins();
1222 
1223 #ifdef BENCHMARK_FAST_SWEEPING
1224  timer.stop();
1225 #endif
1226  }// FastSweeping::InitExt::run
1227 
1228  // int implements an update if minD needs to be updated
1230  void sumHelper(ExtT& sum2, ExtT ext, bool update, const SdfT& /* d2 */) const { if (update) sum2 = ext; }// int implementation
1231 
1232  // non-int implements a weighted sum
1234  void sumHelper(ExtT& sum2, ExtT ext, bool /* update */, const SdfT& d2) const { sum2 += static_cast<ExtValueT>(d2 * ext); }// non-int implementation
1235 
1237  ExtT extValHelper(ExtT extSum, const SdfT& /* sdfSum */) const { return extSum; }// int implementation
1238 
1240  ExtT extValHelper(ExtT extSum, const SdfT& sdfSum) const {return ExtT((SdfT(1) / sdfSum) * extSum); }// non-int implementation
1241 
1242  void operator()(const LeafRange& r) const
1243  {
1244  ExtAccT acc(mExtGrid->tree());
1245  SweepMaskAccT sweepMaskAcc(mParent->mSweepMask);
1246  math::GradStencil<SdfGridT, false> stencil(*mSdfGrid);
1247  const math::Transform& xform = mExtGrid->transform();
1248  typename OpPoolT::reference op = mOpPool->local();
1249  const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();//local copy
1250  const SdfValueT h = mAboveSign*static_cast<SdfValueT>(mSdfGrid->voxelSize()[0]);//Voxel size
1251  for (auto leafIter = r.begin(); leafIter; ++leafIter) {
1252  SdfValueT *sdf = leafIter.buffer(1).data();
1253  ExtValueT *ext = acc.probeLeaf(leafIter->origin())->buffer().data();//should be safe!
1254  for (auto voxelIter = leafIter->beginValueAll(); voxelIter; ++voxelIter) {
1255  const SdfValueT value = *voxelIter;
1256  const bool isAbove = value > isoValue;
1257  if (!voxelIter.isValueOn()) {// inactive voxels
1258  sdf[voxelIter.pos()] = isAbove ? above : -above;
1259  } else {// active voxels
1260  const Coord ijk = voxelIter.getCoord();
1261  stencil.moveTo(ijk, value);
1262  const auto mask = stencil.intersectionMask( isoValue );
1263  if (mask.none()) {// no zero-crossing neighbors, most common case
1264  sdf[voxelIter.pos()] = isAbove ? above : -above;
1265  // the ext grid already has its active values set to the background value
1266  } else {// compute distance to iso-surface
1267  // disable boundary voxels from the mask tree
1268  sweepMaskAcc.setValueOff(ijk);
1269  const SdfValueT delta = value - isoValue;//offset relative to iso-value
1270  if (math::isApproxZero(delta)) {//voxel is on the iso-surface
1271  sdf[voxelIter.pos()] = 0;
1272  ext[voxelIter.pos()] = ExtValueT(op(xform.indexToWorld(ijk)));
1273  } else {//voxel is neighboring the iso-surface
1274  SdfValueT sum1 = 0;
1275  ExtValueT sum2 = zeroVal<ExtValueT>();
1276  // minD is used to update sum2 in the integer case,
1277  // where we choose the value of the extension grid corresponding to
1278  // the smallest value of d in the 6 direction neighboring the current
1279  // voxel.
1280  SdfValueT minD = std::numeric_limits<SdfValueT>::max();
1281  for (int n=0, i=0; i<6;) {
1282  SdfValueT d = std::numeric_limits<SdfValueT>::max(), d2;
1283  if (mask.test(i++)) {
1284  d = math::Abs(delta/(value-stencil.getValue(i)));
1285  n = i - 1;
1286  }
1287  if (mask.test(i++)) {
1288  d2 = math::Abs(delta/(value-stencil.getValue(i)));
1289  if (d2 < d) {
1290  d = d2;
1291  n = i - 1;
1292  }
1293  }
1295  d2 = 1/(d*d);
1296  sum1 += d2;
1297  const Vec3R xyz(static_cast<SdfValueT>(ijk[0])+d*static_cast<SdfValueT>(FastSweeping::mOffset[n][0]),
1298  static_cast<SdfValueT>(ijk[1])+d*static_cast<SdfValueT>(FastSweeping::mOffset[n][1]),
1299  static_cast<SdfValueT>(ijk[2])+d*static_cast<SdfValueT>(FastSweeping::mOffset[n][2]));
1300  // If current d is less than minD, update minD
1301  sumHelper(sum2, ExtValueT(op(xform.indexToWorld(xyz))), d < minD, d2);
1302  if (d < minD) minD = d;
1303  }
1304  }//look over six cases
1305  ext[voxelIter.pos()] = extValHelper(sum2, sum1);
1306  sdf[voxelIter.pos()] = isAbove ? h / math::Sqrt(sum1) : -h / math::Sqrt(sum1);
1307  }// voxel is neighboring the iso-surface
1308  }// intersecting voxels
1309  }// active voxels
1310  }// loop over voxels
1311  }// loop over leaf nodes
1312  }// FastSweeping::InitExt::operator(const LeafRange& r)
1313 
1314  template<typename RootOrInternalNodeT>
1315  void operator()(const RootOrInternalNodeT& node) const
1316  {
1317  const SdfValueT isoValue = mIsoValue, above = mAboveSign*std::numeric_limits<SdfValueT>::max();
1318  for (auto it = node.cbeginValueAll(); it; ++it) {
1319  SdfValueT& v = const_cast<SdfValueT&>(*it);
1320  v = v > isoValue ? above : -above;
1321  }//loop over all tiles
1322  }
1323  // Public member data
1324  FastSweeping *mParent;
1325  OpPoolT *mOpPool;
1326  SdfGridT *mSdfGrid;
1327  ExtGridT *mExtGrid;
1328  SdfValueT mIsoValue;
1329  SdfValueT mAboveSign;//sign of distance values above the iso-value
1330 };// FastSweeping::InitExt
1331 
1332 /// Private class of FastSweeping to perform multi-threaded initialization
1333 template <typename SdfGridT, typename ExtValueT>
1334 template <typename MaskTreeT>
1335 struct FastSweeping<SdfGridT, ExtValueT>::MaskKernel
1336 {
1337  using LeafRange = typename tree::LeafManager<const MaskTreeT>::LeafRange;
1338  MaskKernel(FastSweeping &parent) : mParent(&parent),
1339  mSdfGrid(parent.mSdfGrid.get()) {}
1340  MaskKernel(const MaskKernel &parent) = default;// for tbb::parallel_for
1341  MaskKernel& operator=(const MaskKernel&) = delete;
1342 
1343  void run(const MaskTreeT &mask)
1344  {
1345 #ifdef BENCHMARK_FAST_SWEEPING
1346  util::CpuTimer timer;
1347 #endif
1348  auto &lsTree = mSdfGrid->tree();
1349 
1350  static const SdfValueT Unknown = std::numeric_limits<SdfValueT>::max();
1351 
1352 #ifdef BENCHMARK_FAST_SWEEPING
1353  timer.restart("Changing background value");
1354 #endif
1355  changeLevelSetBackground(lsTree, Unknown);//multi-threaded
1356 
1357 #ifdef BENCHMARK_FAST_SWEEPING
1358  timer.restart("Union with mask");//multi-threaded
1359 #endif
1360  lsTree.topologyUnion(mask);//multi-threaded
1361 
1362  // ignore active tiles since the input grid is assumed to be a level set
1363  tree::LeafManager<const MaskTreeT> mgr(mask);// super fast
1364 
1365 #ifdef BENCHMARK_FAST_SWEEPING
1366  timer.restart("Initializing grid and sweep mask");
1367 #endif
1368 
1369  mParent->mSweepMask.clear();
1370  mParent->mSweepMask.topologyUnion(mParent->mSdfGrid->constTree());
1371 
1372  using LeafManagerT = tree::LeafManager<SweepMaskTreeT>;
1373  using LeafT = typename SweepMaskTreeT::LeafNodeType;
1374  LeafManagerT leafManager(mParent->mSweepMask);
1375 
1376  auto kernel = [&](LeafT& leaf, size_t /*leafIdx*/) {
1377  static const SdfValueT Unknown = std::numeric_limits<SdfValueT>::max();
1378  SdfAccT acc(mSdfGrid->tree());
1379  // The following hack is safe due to the topology union in
1380  // init and the fact that SdfValueT is known to be a floating point!
1381  SdfValueT *data = acc.probeLeaf(leaf.origin())->buffer().data();
1382  for (auto voxelIter = leaf.beginValueOn(); voxelIter; ++voxelIter) {// mask voxels
1383  if (math::Abs( data[voxelIter.pos()] ) < Unknown ) {
1384  // disable boundary voxels from the mask tree
1385  voxelIter.setValue(false);
1386  }
1387  }
1388  };
1389  leafManager.foreach( kernel );
1390 
1391  // cache the leaf node origins for fast lookup in the sweeping kernels
1392  mParent->computeSweepMaskLeafOrigins();
1393 
1394 #ifdef BENCHMARK_FAST_SWEEPING
1395  timer.stop();
1396 #endif
1397  }// FastSweeping::MaskKernel::run
1398 
1399  // Private member data of MaskKernel
1400  FastSweeping *mParent;
1401  SdfGridT *mSdfGrid;//raw pointer, i.e. lock free
1402 };// FastSweeping::MaskKernel
1403 
1404 /// @brief Private class of FastSweeping to perform concurrent fast sweeping in two directions
1405 template <typename SdfGridT, typename ExtValueT>
1406 struct FastSweeping<SdfGridT, ExtValueT>::SweepingKernel
1407 {
1408  SweepingKernel(FastSweeping &parent) : mParent(&parent) {}
1409  SweepingKernel(const SweepingKernel&) = delete;
1411 
1412  /// Main method that performs concurrent bi-directional sweeps
1413  template<typename HashOp>
1415  {
1416 #ifdef BENCHMARK_FAST_SWEEPING
1417  util::CpuTimer timer;
1418 #endif
1419 
1420  // mask of the active voxels to be solved for, i.e. excluding boundary voxels
1421  const SweepMaskTreeT& maskTree = mParent->mSweepMask;
1422 
1423  using LeafManagerT = typename tree::LeafManager<const SweepMaskTreeT>;
1424  using LeafT = typename SweepMaskTreeT::LeafNodeType;
1425  LeafManagerT leafManager(maskTree);
1426 
1427  // compute the leaf node slices that have active voxels in them
1428  // the sliding window of the has keys is -14 to 21 (based on an 8x8x8 leaf node
1429  // and the extrema hash values i-j-k and i+j+k), but we use a larger mask window here to
1430  // easily accommodate any leaf dimension. The mask offset is used to be able to
1431  // store this in a fixed-size byte array
1432  constexpr int maskOffset = LeafT::DIM * 3;
1433  constexpr int maskRange = maskOffset * 2;
1434 
1435  // mark each possible slice in each leaf node that has one or more active voxels in it
1436  std::vector<int8_t> leafSliceMasks(leafManager.leafCount()*maskRange);
1437  auto kernel1 = [&](const LeafT& leaf, size_t leafIdx) {
1438  const size_t leafOffset = leafIdx * maskRange;
1439  for (auto voxelIter = leaf.cbeginValueOn(); voxelIter; ++voxelIter) {
1440  const Coord ijk = LeafT::offsetToLocalCoord(voxelIter.pos());
1441  leafSliceMasks[leafOffset + hash(ijk) + maskOffset] = uint8_t(1);
1442  }
1443  };
1444  leafManager.foreach( kernel1 );
1445 
1446  // compute the voxel slice map using a thread-local-storage hash map
1447  // the key of the hash map is the slice index of the voxel coord (ijk.x() + ijk.y() + ijk.z())
1448  // the values are an array of indices for every leaf that has active voxels with this slice index
1449  using ThreadLocalMap = std::unordered_map</*voxelSliceKey=*/int64_t, /*leafIdx=*/std::deque<size_t>>;
1450  tbb::enumerable_thread_specific<ThreadLocalMap> pool;
1451  auto kernel2 = [&](const LeafT& leaf, size_t leafIdx) {
1452  ThreadLocalMap& map = pool.local();
1453  const Coord& origin = leaf.origin();
1454  const int64_t leafKey = hash(origin);
1455  const size_t leafOffset = leafIdx * maskRange;
1456  for (int sliceIdx = 0; sliceIdx < maskRange; sliceIdx++) {
1457  if (leafSliceMasks[leafOffset + sliceIdx] == uint8_t(1)) {
1458  const int64_t voxelSliceKey = leafKey+sliceIdx-maskOffset;
1459  map[voxelSliceKey].emplace_back(leafIdx);
1460  }
1461  }
1462  };
1463  leafManager.foreach( kernel2 );
1464 
1465  // combine into a single ordered map keyed by the voxel slice key
1466  // note that this is now stored in a map ordered by voxel slice key,
1467  // so sweep slices can be processed in order
1468  for (auto poolIt = pool.begin(); poolIt != pool.end(); ++poolIt) {
1469  const ThreadLocalMap& map = *poolIt;
1470  for (const auto& it : map) {
1471  for (const size_t leafIdx : it.second) {
1472  mVoxelSliceMap[it.first].emplace_back(leafIdx, NodeMaskPtrT());
1473  }
1474  }
1475  }
1476 
1477  // extract the voxel slice keys for random access into the map
1478  mVoxelSliceKeys.reserve(mVoxelSliceMap.size());
1479  for (const auto& it : mVoxelSliceMap) {
1480  mVoxelSliceKeys.push_back(it.first);
1481  }
1482 
1483  // allocate the node masks in parallel, as the map is populated in serial
1484  auto kernel3 = [&](tbb::blocked_range<size_t>& range) {
1485  for (size_t i = range.begin(); i < range.end(); i++) {
1486  const int64_t key = mVoxelSliceKeys[i];
1487  for (auto& it : mVoxelSliceMap[key]) {
1488  it.second = std::make_unique<NodeMaskT>();
1489  }
1490  }
1491  };
1492  tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceKeys.size()), kernel3);
1493 
1494  // each voxel slice contains a leafIdx-nodeMask pair,
1495  // this routine populates these node masks to select only the active voxels
1496  // from the mask tree that have the same voxel slice key
1497  // TODO: a small optimization here would be to union this leaf node mask with
1498  // a pre-computed one for this particular slice pattern
1499  auto kernel4 = [&](tbb::blocked_range<size_t>& range) {
1500  for (size_t i = range.begin(); i < range.end(); i++) {
1501  const int64_t voxelSliceKey = mVoxelSliceKeys[i];
1502  LeafSliceArray& leafSliceArray = mVoxelSliceMap[voxelSliceKey];
1503  for (LeafSlice& leafSlice : leafSliceArray) {
1504  const size_t leafIdx = leafSlice.first;
1505  NodeMaskPtrT& nodeMask = leafSlice.second;
1506  const LeafT& leaf = leafManager.leaf(leafIdx);
1507  const Coord& origin = leaf.origin();
1508  const int64_t leafKey = hash(origin);
1509  for (auto voxelIter = leaf.cbeginValueOn(); voxelIter; ++voxelIter) {
1510  const Index voxelIdx = voxelIter.pos();
1511  const Coord ijk = LeafT::offsetToLocalCoord(voxelIdx);
1512  const int64_t key = leafKey + hash(ijk);
1513  if (key == voxelSliceKey) {
1514  nodeMask->setOn(voxelIdx);
1515  }
1516  }
1517  }
1518  }
1519  };
1520  tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceKeys.size()), kernel4);
1521  }// FastSweeping::SweepingKernel::computeVoxelSlices
1522 
1523  // Private struct for nearest neighbor grid points (very memory light!)
1524  struct NN {
1525  SdfValueT v;
1526  int n;
1527  inline static Coord ijk(const Coord &p, int i) { return p + FastSweeping::mOffset[i]; }
1528  NN() : v(), n() {}
1529  NN(const SdfAccT &a, const Coord &p, int i) : v(math::Abs(a.getValue(ijk(p,i)))), n(i) {}
1530  inline Coord operator()(const Coord &p) const { return ijk(p, n); }
1531  inline bool operator<(const NN &rhs) const { return v < rhs.v; }
1532  inline operator bool() const { return v < SdfValueT(1000); }
1533  };// NN
1534 
1535  /// @note Extending an integer field is based on the nearest-neighbor interpolation
1537  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
1538 
1539  /// @note Extending a non-integer field is based on a weighted interpolation
1541  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
1542 
1543  /// @note Extending an integer field is based on the nearest-neighbor interpolation
1545  ExtT threeNghbr(const NN& d1, const NN& d2, const NN& d3, const SdfT& /* w */, const ExtT& v1, const ExtT& v2, const ExtT& v3) const {
1546  math::Vec3<SdfT> d(d1.v, d2.v, d3.v);
1547  math::Vec3<ExtT> v(v1, v2, v3);
1548  return v[static_cast<int>(math::MinIndex(d))];
1549  }// int implementation
1550 
1551  /// @note Extending a non-integer field is based on a weighted interpolation
1553  ExtT threeNghbr(const NN& d1, const NN& d2, const NN& d3, const SdfT& w, const ExtT& v1, const ExtT& v2, const ExtT& v3) const {
1554  return ExtT(w*(d1.v*v1 + d2.v*v2 + d3.v*v3));
1555  }// non-int implementation
1556 
1557  void sweep()
1558  {
1559  typename ExtGridT::TreeType *tree2 = mParent->mExtGrid ? &mParent->mExtGrid->tree() : nullptr;
1560  typename ExtGridT::TreeType *tree3 = mParent->mExtGridInput ? &mParent->mExtGridInput->tree() : nullptr;
1561 
1562  const SdfValueT h = static_cast<SdfValueT>(mParent->mSdfGrid->voxelSize()[0]);
1563  const SdfValueT sqrt2h = math::Sqrt(SdfValueT(2))*h;
1564  const FastSweepingDomain mode = mParent->mSweepDirection;
1565  const bool isInputSdf = mParent->mIsInputSdf;
1566 
1567  // If we are using an extension in one direction, we need a reference grid
1568  // for the default value of the extension for the voxels that are not
1569  // intended to be updated by the sweeping algorithm.
1570  if (tree2 && mode != FastSweepingDomain::SWEEP_ALL) assert(tree3);
1571 
1572  const std::vector<Coord>& leafNodeOrigins = mParent->mSweepMaskLeafOrigins;
1573 
1574  int64_t voxelSliceIndex(0);
1575 
1576  auto kernel = [&](const tbb::blocked_range<size_t>& range) {
1577  using LeafT = typename SdfGridT::TreeType::LeafNodeType;
1578 
1579  SdfAccT acc1(mParent->mSdfGrid->tree());
1580  auto acc2 = std::unique_ptr<ExtAccT>(tree2 ? new ExtAccT(*tree2) : nullptr);
1581  auto acc3 = std::unique_ptr<ExtAccT>(tree3 ? new ExtAccT(*tree3) : nullptr);
1582  SdfValueT absV, sign, update, D;
1583  NN d1, d2, d3;//distance values and coordinates of closest neighbor points
1584 
1585  const LeafSliceArray& leafSliceArray = mVoxelSliceMap[voxelSliceIndex];
1586 
1587  // Solves Godunov's scheme: [x-d1]^2 + [x-d2]^2 + [x-d3]^2 = h^2
1588  // where [X] = (X>0?X:0) and ai=min(di+1,di-1)
1589  for (size_t i = range.begin(); i < range.end(); ++i) {
1590 
1591  // iterate over all leafs in the slice and extract the leaf
1592  // and node mask for each slice pattern
1593 
1594  const LeafSlice& leafSlice = leafSliceArray[i];
1595  const size_t leafIdx = leafSlice.first;
1596  const NodeMaskPtrT& nodeMask = leafSlice.second;
1597 
1598  const Coord& origin = leafNodeOrigins[leafIdx];
1599 
1600  Coord ijk;
1601  for (auto indexIter = nodeMask->beginOn(); indexIter; ++indexIter) {
1602 
1603  // Get coordinate of center point of the FD stencil
1604  ijk = origin + LeafT::offsetToLocalCoord(indexIter.pos());
1605 
1606  // Find the closes neighbors in the three axial directions
1607  d1 = std::min(NN(acc1, ijk, 0), NN(acc1, ijk, 1));
1608  d2 = std::min(NN(acc1, ijk, 2), NN(acc1, ijk, 3));
1609  d3 = std::min(NN(acc1, ijk, 4), NN(acc1, ijk, 5));
1610 
1611  if (!(d1 || d2 || d3)) continue;//no valid neighbors
1612 
1613  // Get the center point of the FD stencil (assumed to be an active voxel)
1614  // Note this const_cast is normally unsafe but by design we know the tree
1615  // to be static, of floating-point type and containing active voxels only!
1616  SdfValueT &value = const_cast<SdfValueT&>(acc1.getValue(ijk));
1617 
1618  // Extract the sign
1619  sign = value >= SdfValueT(0) ? SdfValueT(1) : SdfValueT(-1);
1620 
1621  // Absolute value
1622  absV = math::Abs(value);
1623 
1624  // sort values so d1 <= d2 <= d3
1625  if (d2 < d1) std::swap(d1, d2);
1626  if (d3 < d2) std::swap(d2, d3);
1627  if (d2 < d1) std::swap(d1, d2);
1628 
1629  // Test if there is a solution depending on ONE of the neighboring voxels
1630  // if d2 - d1 >= h => d2 >= d1 + h then:
1631  // (x-d1)^2=h^2 => x = d1 + h
1632  update = d1.v + h;
1633  if (update <= d2.v) {
1634  if (update < absV) {
1635  value = sign * update;
1636  if (acc2) {
1637  // There is an assert upstream to check if mExtGridInput exists if mode != SWEEP_ALL
1638  ExtValueT updateExt = acc2->getValue(d1(ijk));
1640  if (isInputSdf) updateExt = (value >= SdfValueT(0)) ? acc2->getValue(d1(ijk)) : acc3->getValue(ijk);
1641  else updateExt = (value <= SdfValueT(0)) ? acc2->getValue(d1(ijk)) : acc3->getValue(ijk);
1642  } // SWEEP_GREATER_THAN_ISOVALUE
1644  if (isInputSdf) updateExt = (value <= SdfValueT(0)) ? acc2->getValue(d1(ijk)) : acc3->getValue(ijk);
1645  else updateExt = (value >= SdfValueT(0)) ? acc2->getValue(d1(ijk)) : acc3->getValue(ijk);
1646  } // SWEEP_LESS_THAN_ISOVALUE
1647  acc2->setValue(ijk, updateExt);
1648  }//update ext?
1649  }//update sdf?
1650  continue;
1651  }// one neighbor case
1652 
1653  // Test if there is a solution depending on TWO of the neighboring voxels
1654  // (x-d1)^2 + (x-d2)^2 = h^2
1655  //D = SdfValueT(2) * h * h - math::Pow2(d1.v - d2.v);// = 2h^2-(d1-d2)^2
1656  //if (D >= SdfValueT(0)) {// non-negative discriminant
1657  if (d2.v <= sqrt2h + d1.v) {
1658  D = SdfValueT(2) * h * h - math::Pow2(d1.v - d2.v);// = 2h^2-(d1-d2)^2
1659  update = SdfValueT(0.5) * (d1.v + d2.v + std::sqrt(D));
1660  if (update > d2.v && update <= d3.v) {
1661  if (update < absV) {
1662  value = sign * update;
1663  if (acc2) {
1664  d1.v -= update;
1665  d2.v -= update;
1666  // affine combination of two neighboring extension values
1667  const SdfValueT w = SdfValueT(1)/(d1.v+d2.v);
1668  const ExtValueT v1 = acc2->getValue(d1(ijk));
1669  const ExtValueT v2 = acc2->getValue(d2(ijk));
1670  const ExtValueT extVal = twoNghbr(d1, d2, w, v1, v2);
1671 
1672  ExtValueT updateExt = extVal;
1674  if (isInputSdf) updateExt = (value >= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1675  else updateExt = (value <= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1676  } // SWEEP_GREATER_THAN_ISOVALUE
1678  if (isInputSdf) updateExt = (value <= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1679  else updateExt = (value >= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1680  } // SWEEP_LESS_THAN_ISOVALUE
1681  acc2->setValue(ijk, updateExt);
1682  }//update ext?
1683  }//update sdf?
1684  continue;
1685  }//test for two neighbor case
1686  }//test for non-negative determinant
1687 
1688  // Test if there is a solution depending on THREE of the neighboring voxels
1689  // (x-d1)^2 + (x-d2)^2 + (x-d3)^2 = h^2
1690  // 3x^2 - 2(d1 + d2 + d3)x + d1^2 + d2^2 + d3^2 = h^2
1691  // ax^2 + bx + c=0, a=3, b=-2(d1+d2+d3), c=d1^2 + d2^2 + d3^2 - h^2
1692  const SdfValueT d123 = d1.v + d2.v + d3.v;
1693  D = d123*d123 - SdfValueT(3)*(d1.v*d1.v + d2.v*d2.v + d3.v*d3.v - h * h);
1694  if (D >= SdfValueT(0)) {// non-negative discriminant
1695  update = SdfValueT(1.0/3.0) * (d123 + std::sqrt(D));//always passes test
1696  //if (update > d3.v) {//disabled due to round-off errors
1697  if (update < absV) {
1698  value = sign * update;
1699  if (acc2) {
1700  d1.v -= update;
1701  d2.v -= update;
1702  d3.v -= update;
1703  // affine combination of three neighboring extension values
1704  const SdfValueT w = SdfValueT(1)/(d1.v+d2.v+d3.v);
1705  const ExtValueT v1 = acc2->getValue(d1(ijk));
1706  const ExtValueT v2 = acc2->getValue(d2(ijk));
1707  const ExtValueT v3 = acc2->getValue(d3(ijk));
1708  const ExtValueT extVal = threeNghbr(d1, d2, d3, w, v1, v2, v3);
1709 
1710  ExtValueT updateExt = extVal;
1712  if (isInputSdf) updateExt = (value >= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1713  else updateExt = (value <= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1714  } // SWEEP_GREATER_THAN_ISOVALUE
1716  if (isInputSdf) updateExt = (value <= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1717  else updateExt = (value >= SdfValueT(0)) ? extVal : acc3->getValue(ijk);
1718  } // SWEEP_LESS_THAN_ISOVALUE
1719  acc2->setValue(ijk, updateExt);
1720  }//update ext?
1721  }//update sdf?
1722  }//test for non-negative determinant
1723  }//loop over coordinates
1724  }
1725  };
1726 
1727 #ifdef BENCHMARK_FAST_SWEEPING
1728  util::CpuTimer timer("Forward sweep");
1729 #endif
1730 
1731  for (size_t i = 0; i < mVoxelSliceKeys.size(); i++) {
1732  voxelSliceIndex = mVoxelSliceKeys[i];
1733  tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceMap[voxelSliceIndex].size()), kernel);
1734  }
1735 
1736 #ifdef BENCHMARK_FAST_SWEEPING
1737  timer.restart("Backward sweeps");
1738 #endif
1739  for (size_t i = mVoxelSliceKeys.size(); i > 0; i--) {
1740  voxelSliceIndex = mVoxelSliceKeys[i-1];
1741  tbb::parallel_for(tbb::blocked_range<size_t>(0, mVoxelSliceMap[voxelSliceIndex].size()), kernel);
1742  }
1743 
1744 #ifdef BENCHMARK_FAST_SWEEPING
1745  timer.stop();
1746 #endif
1747  }// FastSweeping::SweepingKernel::sweep
1748 
1749 private:
1750  using NodeMaskT = typename SweepMaskTreeT::LeafNodeType::NodeMaskType;
1751  using NodeMaskPtrT = std::unique_ptr<NodeMaskT>;
1752  // using a unique ptr for the NodeMask allows for parallel allocation,
1753  // but makes this class not copy-constructible
1754  using LeafSlice = std::pair</*leafIdx=*/size_t, /*leafMask=*/NodeMaskPtrT>;
1755  using LeafSliceArray = std::deque<LeafSlice>;
1756  using VoxelSliceMap = std::map</*voxelSliceKey=*/int64_t, LeafSliceArray>;
1757 
1758  // Private member data of SweepingKernel
1759  FastSweeping *mParent;
1760  VoxelSliceMap mVoxelSliceMap;
1761  std::vector<int64_t> mVoxelSliceKeys;
1762 };// FastSweeping::SweepingKernel
1763 
1764 ////////////////////////////////////////////////////////////////////////////////
1765 
1766 template<typename GridT>
1767 typename GridT::Ptr
1768 fogToSdf(const GridT &fogGrid,
1769  typename GridT::ValueType isoValue,
1770  int nIter)
1771 {
1773  if (fs.initSdf(fogGrid, isoValue, /*isInputSdf*/false)) fs.sweep(nIter);
1774  return fs.sdfGrid();
1775 }
1776 
1777 template<typename GridT>
1778 typename GridT::Ptr
1779 sdfToSdf(const GridT &sdfGrid,
1780  typename GridT::ValueType isoValue,
1781  int nIter)
1782 {
1784  if (fs.initSdf(sdfGrid, isoValue, /*isInputSdf*/true)) fs.sweep(nIter);
1785  return fs.sdfGrid();
1786 }
1787 
1788 template<typename FogGridT, typename ExtOpT, typename ExtValueT>
1789 typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr
1790 fogToExt(const FogGridT &fogGrid,
1791  const ExtOpT &op,
1792  const ExtValueT& background,
1793  typename FogGridT::ValueType isoValue,
1794  int nIter,
1795  FastSweepingDomain mode,
1796  const typename FogGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid)
1797 {
1799  if (fs.initExt(fogGrid, op, background, isoValue, /*isInputSdf*/false, mode, extGrid))
1800  fs.sweep(nIter, /*finalize*/true);
1801  return fs.extGrid();
1802 }
1803 
1804 template<typename SdfGridT, typename OpT, typename ExtValueT>
1805 typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr
1806 sdfToExt(const SdfGridT &sdfGrid,
1807  const OpT &op,
1808  const ExtValueT &background,
1809  typename SdfGridT::ValueType isoValue,
1810  int nIter,
1811  FastSweepingDomain mode,
1812  const typename SdfGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid)
1813 {
1815  if (fs.initExt(sdfGrid, op, background, isoValue, /*isInputSdf*/true, mode, extGrid))
1816  fs.sweep(nIter, /*finalize*/true);
1817  return fs.extGrid();
1818 }
1819 
1820 template<typename FogGridT, typename ExtOpT, typename ExtValueT>
1821 std::pair<typename FogGridT::Ptr, typename FogGridT::template ValueConverter<ExtValueT>::Type::Ptr>
1822 fogToSdfAndExt(const FogGridT &fogGrid,
1823  const ExtOpT &op,
1824  const ExtValueT &background,
1825  typename FogGridT::ValueType isoValue,
1826  int nIter,
1827  FastSweepingDomain mode,
1828  const typename FogGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid)
1829 {
1831  if (fs.initExt(fogGrid, op, background, isoValue, /*isInputSdf*/false, mode, extGrid))
1832  fs.sweep(nIter, /*finalize*/true);
1833  return std::make_pair(fs.sdfGrid(), fs.extGrid());
1834 }
1835 
1836 template<typename SdfGridT, typename ExtOpT, typename ExtValueT>
1837 std::pair<typename SdfGridT::Ptr, typename SdfGridT::template ValueConverter<ExtValueT>::Type::Ptr>
1838 sdfToSdfAndExt(const SdfGridT &sdfGrid,
1839  const ExtOpT &op,
1840  const ExtValueT &background,
1841  typename SdfGridT::ValueType isoValue,
1842  int nIter,
1843  FastSweepingDomain mode,
1844  const typename SdfGridT::template ValueConverter<ExtValueT>::Type::ConstPtr extGrid)
1845 {
1847  if (fs.initExt(sdfGrid, op, background, isoValue, /*isInputSdf*/true, mode, extGrid))
1848  fs.sweep(nIter, /*finalize*/true);
1849  return std::make_pair(fs.sdfGrid(), fs.extGrid());
1850 }
1851 
1852 template<typename GridT>
1853 typename GridT::Ptr
1854 dilateSdf(const GridT &sdfGrid,
1855  int dilation,
1856  NearestNeighbors nn,
1857  int nIter,
1858  FastSweepingDomain mode)
1859 {
1861  if (fs.initDilate(sdfGrid, dilation, nn, /*sweep direction*/ mode)) fs.sweep(nIter);
1862  return fs.sdfGrid();
1863 }
1864 
1865 template<typename GridT, typename MaskTreeT>
1866 typename GridT::Ptr
1867 maskSdf(const GridT &sdfGrid,
1868  const Grid<MaskTreeT> &mask,
1869  bool ignoreActiveTiles,
1870  int nIter)
1871 {
1873  if (fs.initMask(sdfGrid, mask, ignoreActiveTiles)) fs.sweep(nIter);
1874  return fs.sdfGrid();
1875 }
1876 
1877 
1878 ////////////////////////////////////////
1879 
1880 
1881 // Explicit Template Instantiation
1882 
1883 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
1884 
1885 #ifdef OPENVDB_INSTANTIATE_FASTSWEEPING
1887 #endif
1888 
1889 #define _FUNCTION(TreeT) \
1890  Grid<TreeT>::Ptr fogToSdf(const Grid<TreeT>&, TreeT::ValueType, int)
1892 #undef _FUNCTION
1893 
1894 #define _FUNCTION(TreeT) \
1895  Grid<TreeT>::Ptr sdfToSdf(const Grid<TreeT>&, TreeT::ValueType, int)
1897 #undef _FUNCTION
1898 
1899 #define _FUNCTION(TreeT) \
1900  Grid<TreeT>::Ptr dilateSdf(const Grid<TreeT>&, int, NearestNeighbors, int, FastSweepingDomain)
1902 #undef _FUNCTION
1903 
1904 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
1905 
1906 
1907 } // namespace tools
1908 } // namespace OPENVDB_VERSION_NAME
1909 } // namespace openvdb
1910 
1911 #endif // OPENVDB_TOOLS_FASTSWEEPING_HAS_BEEN_INCLUDED
ValueT value
Definition: GridBuilder.h:1287
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...
Miscellaneous utility methods that operate primarily or exclusively on level set grids.
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
Implementation of morphological dilation and erosion.
Functions to efficiently compute histograms, extrema (min/max) and statistics (mean,...
math::Transform & transform()
Return a reference to this grid's transform, which might be shared with other grids.
Definition: Grid.h:415
GridClass getGridClass() const
Return the class of volumetric data (level set, fog volume, etc.) that is stored in this grid.
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:577
TreeType & tree()
Return a reference to this grid's tree, which might be shared with other grids.
Definition: Grid.h:917
SharedPtr< Grid > Ptr
Definition: Grid.h:579
Definition: Exceptions.h:63
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:564
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
const ValueType & getValue(unsigned int pos=0) const
Return the value from the stencil buffer with linear offset pos.
Definition: Stencils.h:97
void moveTo(const Coord &ijk)
Initialize the stencil buffer with the values of voxel (i, j, k) and its neighbors.
Definition: Stencils.h:47
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:25
Definition: Stencils.h:1232
Templated class to compute the minimum and maximum values.
Definition: Stats.h:31
Definition: Vec3.h:24
Computes signed distance values from an initial iso-surface and optionally performs velocity extensio...
Definition: FastSweeping.h:454
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:806
size_t boundaryVoxelCount() const
Return the number of voxels that defined the boundary condition.
Definition: FastSweeping.h:656
FastSweepingDomain sweepDirection() const
Return whether the sweep update is in all direction (SWEEP_ALL), greater than isovalue (SWEEP_GREATER...
Definition: FastSweeping.h:667
bool isInputSdf()
Return whether the fast-sweeping input grid a signed distance function or not (fog).
Definition: FastSweeping.h:670
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:508
FastSweeping()
Constructor.
Definition: FastSweeping.h:710
bool isValid() const
Return true if there are voxels and boundaries to solve for.
Definition: FastSweeping.h:659
size_t sweepingVoxelCount() const
Return the number of voxels that will be solved for.
Definition: FastSweeping.h:653
~FastSweeping()
Destructor.
Definition: FastSweeping.h:478
FastSweeping & operator=(const FastSweeping &)=delete
Disallow copy assignment.
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:757
SdfGridT::Ptr sdfGrid()
Returns a shared pointer to the signed distance field computed by this class.
Definition: FastSweeping.h:492
void sweep(int nIter=1, bool finalize=true)
Perform nIter iterations of the fast sweeping algorithm.
Definition: FastSweeping.h:837
ExtGridT::Ptr extGrid()
Returns a shared pointer to the extension field computed by this class.
Definition: FastSweeping.h:500
void clear()
Clears all the grids and counters so initialization can be called again.
Definition: FastSweeping.h:716
FastSweeping(const FastSweeping &)=delete
Disallow copy construction.
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:794
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...
Definition: LeafManager.h:102
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:85
bool swapLeafBuffer(size_t bufferIdx, bool serial=false)
Swap each leaf node's buffer with the nth corresponding auxiliary buffer, where n = bufferIdx.
Definition: LeafManager.h:359
LeafRange leafRange(size_t grainsize=1) const
Return a TBB-compatible LeafRange.
Definition: LeafManager.h:345
To facilitate threading over the nodes of a tree, cache node pointers in linear arrays,...
Definition: NodeManager.h:531
Definition: ValueAccessor.h:183
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
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition: ValueAccessor.h:219
Simple timer for basic profiling.
Definition: CpuTimer.h:67
double stop() const
Returns and prints time in milliseconds since construction or start was called.
Definition: CpuTimer.h:128
double restart()
Re-start timer.
Definition: CpuTimer.h:150
__hostdev__ uint32_t hash(uint32_t x)
Definition: common.h:13
bool isValid(GridType gridType, GridClass gridClass)
return true if the combination of GridType and GridClass is valid.
Definition: NanoVDB.h:520
void run(const char *ax, openvdb::GridBase &grid)
Run a full AX pipeline (parse, compile and execute) on a single OpenVDB Grid.
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:764
Coord Abs(const Coord &xyz)
Definition: Coord.h:514
size_t MinIndex(const Vec3T &v)
Return the index [0,1,2] of the smallest value in a 3D vector.
Definition: Math.h:938
Type Pow2(Type x)
Return x2.
Definition: Math.h:551
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
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:107
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:103
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:1768
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...
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:1838
NearestNeighbors
Voxel topology of nearest neighbors.
Definition: Morphology.h:58
@ NN_FACE
Definition: Morphology.h:58
GridT::Ptr dilateSdf(const GridT &sdfGrid, int dilation, NearestNeighbors nn=NN_FACE, int nIter=1, FastSweepingDomain mode=FastSweepingDomain::SWEEP_ALL)
Dilates an existing signed distance field by a specified number of voxels.
Definition: FastSweeping.h:1854
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
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:1779
FastSweepingDomain
Fast Sweeping update mode. This is useful to determine narrow-band extension or field extension in on...
Definition: FastSweeping.h:64
@ SWEEP_ALL
Update all voxels affected by the sweeping algorithm.
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
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
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:1822
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:1867
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:1790
MeshToVoxelEdgeData::EdgeData Abs(const MeshToVoxelEdgeData::EdgeData &x)
Definition: MeshToVolume.h:3707
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:1806
@ IGNORE_TILES
Definition: Morphology.h:80
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
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
Index32 Index
Definition: Types.h:54
@ GRID_LEVEL_SET
Definition: Types.h:337
math::Vec3< Real > Vec3R
Definition: Types.h:72
Definition: Exceptions.h:13
Definition: Coord.h:586
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
Defines various finite difference stencils by means of the "curiously recurring template pattern" on ...
Private class of FastSweeping to perform multi-threaded initialization.
Definition: FastSweeping.h:941
DilateKernel & operator=(const DilateKernel &)=delete
const SdfValueT mBackground
Definition: FastSweeping.h:1043
typename tree::LeafManager< SdfTreeT >::LeafRange LeafRange
Definition: FastSweeping.h:942
SdfGridT::ConstPtr mSdfGridInput
Definition: FastSweeping.h:1044
void run(int dilation, NearestNeighbors nn)
Definition: FastSweeping.h:952
DilateKernel(const DilateKernel &parent)=default
DilateKernel(FastSweeping &parent)
Definition: FastSweeping.h:943
FastSweeping * mParent
Definition: FastSweeping.h:1042
Definition: FastSweeping.h:1051
typename tree::LeafManager< SdfTreeT >::LeafRange LeafRange
Definition: FastSweeping.h:1052
void operator()(const RootOrInternalNodeT &node) const
Definition: FastSweeping.h:1140
void operator()(const LeafRange &r) const
Definition: FastSweeping.h:1095
void run(SdfValueT isoValue)
Definition: FastSweeping.h:1058
SdfGridT * mSdfGrid
Definition: FastSweeping.h:1151
SdfValueT mAboveSign
Definition: FastSweeping.h:1153
SdfValueT mIsoValue
Definition: FastSweeping.h:1152
InitSdf & operator=(const InitSdf &)=delete
InitSdf(FastSweeping &parent)
Definition: FastSweeping.h:1053
FastSweeping * mParent
Definition: FastSweeping.h:1150
SdfValueT mMin
Definition: FastSweeping.h:933
MinMaxKernel(MinMaxKernel &other, tbb::split)
Definition: FastSweeping.h:907
void operator()(const LeafRange &r)
Definition: FastSweeping.h:916
void join(const MinMaxKernel &other)
Definition: FastSweeping.h:927
typename LeafMgr::LeafRange LeafRange
Definition: FastSweeping.h:905
math::MinMax< SdfValueT > run(const SdfGridT &grid)
Definition: FastSweeping.h:909
SdfValueT mMax
Definition: FastSweeping.h:933
MinMaxKernel()
Definition: FastSweeping.h:906
SdfValueT v
Definition: FastSweeping.h:1525
Coord operator()(const Coord &p) const
Definition: FastSweeping.h:1530
bool operator<(const NN &rhs) const
Definition: FastSweeping.h:1531
NN(const SdfAccT &a, const Coord &p, int i)
Definition: FastSweeping.h:1529
static Coord ijk(const Coord &p, int i)
Definition: FastSweeping.h:1527
Private class of FastSweeping to perform concurrent fast sweeping in two directions.
Definition: FastSweeping.h:1407
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:1545
ExtT twoNghbr(const NN &d1, const NN &d2, const SdfT &w, const ExtT &v1, const ExtT &v2) const
Definition: FastSweeping.h:1541
void sweep()
Definition: FastSweeping.h:1557
SweepingKernel(const SweepingKernel &)=delete
SweepingKernel(FastSweeping &parent)
Definition: FastSweeping.h:1408
ExtT twoNghbr(const NN &d1, const NN &d2, const SdfT &, const ExtT &v1, const ExtT &v2) const
Definition: FastSweeping.h:1537
SweepingKernel & operator=(const SweepingKernel &)=delete
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:1553
void computeVoxelSlices(HashOp hash)
Main method that performs concurrent bi-directional sweeps.
Definition: FastSweeping.h:1414
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:116
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:202
#define OPENVDB_REAL_TREE_INSTANTIATE(Function)
Definition: version.h.in:147