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