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