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