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