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