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