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