FastSweeping.h

Go to the documentation of this file.

// Copyright Contributors to the OpenVDB Project
// SPDX-License-Identifier: MPL-2.0
//
/// @file    FastSweeping.h
///
/// @author  Ken Museth
///
/// @brief   Defined the six functions {fog,sdf}To{Sdf,Ext,SdfAndExt} in
///          addition to the two functions maskSdf and dilateSdf. Sdf denotes
///          a signed-distance field (i.e. negative values are inside), fog
///          is a scalar fog volume (i.e. higher values are inside), and Ext is
///          a field (of arbitrary type) that is extended off the iso-surface.
///          All these functions are implemented with the methods in the class
///          named FastSweeping.
///
/// @note    Solves the (simplified) Eikonal Eq: @f$|\nabla \phi|^2 = 1@f$ and
///          performs velocity extension,  @f$\nabla f\nabla \phi = 0@f$, both
///          by means of the fast sweeping algorithm detailed in:
///          "A Fast Sweeping Method For Eikonal Equations"
///          by H. Zhao, Mathematics of Computation, Vol 74(230), pp 603-627, 2004
///
/// @details The algorithm used below for parallel fast sweeping was first published in:
///          "New Algorithm for Sparse and Parallel Fast Sweeping: Efficient
///          Computation of Sparse Distance Fields" by K. Museth, ACM SIGGRAPH Talk,
///          2017, http://www.museth.org/Ken/Publications_files/Museth_SIG17.pdf
 
#ifndef OPENVDB_TOOLS_FASTSWEEPING_HAS_BEEN_INCLUDED
#define OPENVDB_TOOLS_FASTSWEEPING_HAS_BEEN_INCLUDED
 
//#define BENCHMARK_FAST_SWEEPING
 
#include <openvdb/openvdb.h>
#include <openvdb/Platform.h>
#include <openvdb/math/Math.h> // for Abs() and isExactlyEqual()
#include <openvdb/math/Stencils.h> // for GradStencil
#include <openvdb/tree/LeafManager.h>
#include <openvdb/tree/NodeManager.h> // for PruneMinMaxFltKernel
 
#include "LevelSetUtil.h"
#include "Morphology.h"
 
#include "Statistics.h"
#ifdef BENCHMARK_FAST_SWEEPING
#include <openvdb/util/CpuTimer.h>
#endif
 
#include <tbb/parallel_for.h>
#include <tbb/enumerable_thread_specific.h>
#include <tbb/task_group.h>
 
#include <type_traits>// for static_assert
#include <cmath>
#include <limits>
#include <deque>
#include <unordered_map>
#include <utility>// for std::make_pair
 
namespace openvdb {
OPENVDB_USE_VERSION_NAMESPACE
namespace OPENVDB_VERSION_NAME {
namespace tools {
 
/// @brief Fast Sweeping update mode. This is useful to determine
///        narrow-band extension or field extension in one side
///        of a signed distance field.
enum class FastSweepingDomain {
    /// Update all voxels affected by the sweeping algorithm
    SWEEP_ALL,
    // Update voxels corresponding to an sdf/fog values that are greater than a given isovalue
    SWEEP_GREATER_THAN_ISOVALUE,
    // Update voxels corresponding to an sdf/fog values that are less than a given isovalue
    SWEEP_LESS_THAN_ISOVALUE
};
 
/// @brief Converts a scalar fog volume into a signed distance function. Active input voxels
///        with scalar values above the given isoValue will have NEGATIVE distance
///        values on output, i.e. they are assumed to be INSIDE the iso-surface.
///
/// @return A shared pointer to a signed-distance field defined on the active values
///         of the input fog volume.
///
/// @param fogGrid  Scalar (floating-point) volume from which an
///                 iso-surface can be defined.
///
/// @param isoValue A value which defines a smooth iso-surface that
///                 intersects active voxels in @a fogGrid.
///
/// @param nIter    Number of iterations of the fast sweeping algorithm.
///                 Each iteration performs 2^3 = 8 individual sweeps.
///
/// @note Strictly speaking a fog volume is normalized to the range [0,1] but this
///       method accepts a scalar volume with an arbritary range, as long as the it
///       includes the @a isoValue.
///
/// @details Topology of output grid is identical to that of the input grid, except
///          active tiles in the input grid will be converted to active voxels
///          in the output grid!
///
/// @warning If @a isoValue does not intersect any active values in
///          @a fogGrid then the returned grid has all its active values set to
///          plus or minus infinity, depending on if the input values are larger or
///          smaller than @a isoValue.
template<typename GridT>
typename GridT::Ptr
fogToSdf(const GridT &fogGrid,
         typename GridT::ValueType isoValue,
         int nIter = 1);
 
/// @brief Given an existing approximate SDF it solves the Eikonal equation for all its
///        active voxels. Active input voxels with a signed distance value above the
///        given isoValue will have POSITIVE distance values on output, i.e. they are
///        assumed to be OUTSIDE the iso-surface.
///
/// @return A shared pointer to a signed-distance field defined on the active values
///         of the input sdf volume.
///
/// @param sdfGrid  An approximate signed distance field to the specified iso-surface.
///
/// @param isoValue A value which defines a smooth iso-surface that
///                 intersects active voxels in @a sdfGrid.
///
/// @param nIter    Number of iterations of the fast sweeping algorithm.
///                 Each iteration performs 2^3 = 8 individual sweeps.
///
/// @note The only difference between this method and fogToSdf, defined above, is the
///       convention of the sign of the output distance field.
///
/// @details Topology of output grid is identical to that of the input grid, except
///          active tiles in the input grid will be converted to active voxels
///          in the output grid!
///
/// @warning If @a isoValue does not intersect any active values in
///          @a sdfGrid then the returned grid has all its active values set to
///          plus or minus infinity, depending on if the input values are larger or
///          smaller than @a isoValue.
template<typename GridT>
typename GridT::Ptr
sdfToSdf(const GridT &sdfGrid,
         typename GridT::ValueType isoValue = 0,
         int nIter = 1);
 
/// @brief Computes the extension of a field (scalar, vector, or int are supported), defined
///        by the specified functor, off an iso-surface from an input FOG volume.
///
/// @return A shared pointer to the extension field defined from the active values in
///         the input fog volume.
///
/// @param fogGrid  Scalar (floating-point) volume from which an
///                 iso-surface can be defined.
///
/// @param op       Functor with signature [](const Vec3R &xyz)->ExtValueT that
///                 defines the Dirichlet boundary condition, on the iso-surface,
///                 of the field to be extended.
///
/// @param background Background value of return grid with the extension field.
///
/// @param isoValue A value which defines a smooth iso-surface that
///                 intersects active voxels in @a fogGrid.
///
/// @param nIter    Number of iterations of the fast sweeping algorithm.
///                 Each iteration performs 2^3 = 8 individual sweeps.
///
/// @param mode     Determines the mode of updating the extension field. SWEEP_ALL
///                 will update all voxels of the extension field affected by the
///                 fast sweeping algorithm. SWEEP_GREATER_THAN_ISOVALUE will update
///                 all voxels corresponding to fog values that are greater than a given
///                 isovalue. SWEEP_LESS_THAN_ISOVALUE will update all voxels corresponding
///                 to fog values that are less than a given isovalue. If a mode other
///                 than SWEEP_ALL is chosen, a user needs to supply @a extGrid.
///
/// @param extGrid  Optional parameter required to supply a default value for the extension
///                 field when SWEEP_GREATER_THAN_ISOVALUE or SWEEP_LESS_THAN_ISOVALUE
///                 mode is picked for @a mode. When SWEEP_GREATER_THAN_ISOVALUE is supplied
///                 as an argument for @a mode, the extension field voxel will default
///                 to the value of the @a extGrid in that position if it corresponds to a fog
///                 value that is less than the isovalue. Otherwise, the extension
///                 field voxel value will be computed by the Fast Sweeping algorithm.
///                 The opposite convention is implemented when SWEEP_LESS_THAN_ISOVALUE
///                 is supplied as an argument for @a mode.
///
/// @note Strictly speaking a fog volume is normalized to the range [0,1] but this
///       method accepts a scalar volume with an arbritary range, as long as the it
///       includes the @a isoValue.
///
/// @details Topology of output grid is identical to that of the input grid, except
///          active tiles in the input grid will be converted to active voxels
///          in the output grid!
///
/// @warning If @a isoValue does not intersect any active values in
///          @a fogGrid then the returned grid has all its active values set to
///          @a background.
template<typename FogGridT, typename ExtOpT, typename ExtValueT>
typename 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);
 
/// @brief Computes the extension of a field (scalar, vector, or int are supported), defined
///        by the specified functor, off an iso-surface from an input SDF volume.
///
/// @return A shared pointer to the extension field defined on the active values in the
///         input signed distance field.
///
/// @param sdfGrid  An approximate signed distance field to the specified iso-surface.
///
/// @param op       Functor with signature [](const Vec3R &xyz)->ExtValueT that
///                 defines the Dirichlet boundary condition, on the iso-surface,
///                 of the field to be extended.
///
/// @param background Background value of return grid with the extension field.
///
/// @param isoValue A value which defines a smooth iso-surface that
///                 intersects active voxels in @a sdfGrid.
///
/// @param nIter    Number of iterations of the fast sweeping algorithm.
///                 Each iteration performs 2^3 = 8 individual sweeps.
///
/// @param mode     Determines the mode of updating the extension field. SWEEP_ALL
///                 will update all voxels of the extension field affected by the
///                 fast sweeping algorithm. SWEEP_GREATER_THAN_ISOVALUE will update
///                 all voxels corresponding to level set values that are greater than a given
///                 isovalue. SWEEP_LESS_THAN_ISOVALUE will update all voxels corresponding
///                 to level set values that are less than a given isovalue. If a mode other
///                 than SWEEP_ALL is chosen, a user needs to supply @a extGrid.
///
/// @param extGrid  Optional parameter required to supply a default value for the extension
///                 field when SWEEP_GREATER_THAN_ISOVALUE or SWEEP_LESS_THAN_ISOVALUE
///                 mode is picked for @a mode. When SWEEP_GREATER_THAN_ISOVALUE is supplied
///                 as an argument for @a mode, the extension field voxel will default
///                 to the value of the @a extGrid in that position if it corresponds to a level-set
///                 value that is less than the isovalue. Otherwise, the extension
///                 field voxel value will be computed by the Fast Sweeping algorithm.
///                 The opposite convention is implemented when SWEEP_LESS_THAN_ISOVALUE
///                 is supplied as an argument for @a mode.
///
/// @note The only difference between this method and fogToExt, defined above, is the
///       convention of the sign of the signed distance field.
///
/// @details Topology of output grid is identical to that of the input grid, except
///          active tiles in the input grid will be converted to active voxels
///          in the output grid!
///
/// @warning If @a isoValue does not intersect any active values in
///          @a sdfGrid then the returned grid has all its active values set to
///          @a background.
template<typename SdfGridT, typename ExtOpT, typename ExtValueT>
typename 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);
 
/// @brief Computes the signed distance field and the extension of a field (scalar, vector, or
///        int are supported), defined by the specified functor, off an iso-surface from an input
///        FOG volume.
///
/// @return An pair of two shared pointers to respectively the SDF and extension field
///
/// @param fogGrid  Scalar (floating-point) volume from which an
///                 iso-surface can be defined.
///
/// @param op       Functor with signature [](const Vec3R &xyz)->ExtValueT that
///                 defines the Dirichlet boundary condition, on the iso-surface,
///                 of the field to be extended.
///
/// @param background Background value of return grid with the extension field.
///
/// @param isoValue A value which defines a smooth iso-surface that
///                 intersects active voxels in @a fogGrid.
///
/// @param nIter    Number of iterations of the fast sweeping algorithm.
///                 Each iteration performs 2^3 = 8 individual sweeps.
///
/// @param mode     Determines the mode of updating the extension field. SWEEP_ALL
///                 will update all voxels of the extension field affected by the
///                 fast sweeping algorithm. SWEEP_GREATER_THAN_ISOVALUE will update
///                 all voxels corresponding to fog values that are greater than a given
///                 isovalue. SWEEP_LESS_THAN_ISOVALUE will update all voxels corresponding
///                 to fog values that are less than a given isovalue. If a mode other
///                 than SWEEP_ALL is chosen, a user needs to supply @a extGrid.
///
/// @param extGrid  Optional parameter required to supply a default value for the extension
///                 field when SWEEP_GREATER_THAN_ISOVALUE or SWEEP_LESS_THAN_ISOVALUE
///                 mode is picked for @a mode. When SWEEP_GREATER_THAN_ISOVALUE is supplied
///                 as an argument for @a mode, the extension field voxel will default
///                 to the value of the @a extGrid in that position if it corresponds to a fog
///                 value that is less than the isovalue. Otherwise, the extension
///                 field voxel value will be computed by the Fast Sweeping algorithm.
///                 The opposite convention is implemented when SWEEP_LESS_THAN_ISOVALUE
///                 is supplied as an argument for @a mode.
///
/// @note Strictly speaking a fog volume is normalized to the range [0,1] but this
///       method accepts a scalar volume with an arbritary range, as long as the it
///       includes the @a isoValue.
///
/// @details Topology of output grids are identical to that of the input grid, except
///          active tiles in the input grid will be converted to active voxels
///          in the output grids!
///
/// @warning If @a isoValue does not intersect any active values in
///          @a fogGrid then a pair of the following grids is returned: The first
///          is a signed distance grid with its active values set to plus or minus
///          infinity depending of whether its input values are above or below @a isoValue.
///          The second grid, which represents the extension field, has all its active
///          values set to @a background.
template<typename FogGridT, typename ExtOpT, typename ExtValueT>
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);
 
/// @brief Computes the signed distance field and the extension of a field (scalar, vector, or
///        int are supported), defined by the specified functor, off an iso-surface from an input
///        SDF volume.
///
/// @return A pair of two shared pointers to respectively the SDF and extension field
///
/// @param sdfGrid  Scalar (floating-point) volume from which an
///                 iso-surface can be defined.
///
/// @param op       Functor with signature [](const Vec3R &xyz)->ExtValueT that
///                 defines the Dirichlet boundary condition, on the iso-surface,
///                 of the field to be extended.
///
/// @param background Background value of return grid with the extension field.
///
/// @param isoValue A value which defines a smooth iso-surface that
///                 intersects active voxels in @a sdfGrid.
///
/// @param nIter    Number of iterations of the fast sweeping algorithm.
///                 Each iteration performs 2^3 = 8 individual sweeps.
///
/// @param mode     Determines the mode of updating the extension field. SWEEP_ALL
///                 will update all voxels of the extension field affected by the
///                 fast sweeping algorithm. SWEEP_GREATER_THAN_ISOVALUE will update
///                 all voxels corresponding to level set values that are greater than a given
///                 isovalue. SWEEP_LESS_THAN_ISOVALUE will update all voxels corresponding
///                 to level set values that are less than a given isovalue. If a mode other
///                 than SWEEP_ALL is chosen, a user needs to supply @a extGrid.
///
/// @param extGrid  Optional parameter required to supply a default value for the extension
///                 field when SWEEP_GREATER_THAN_ISOVALUE or SWEEP_LESS_THAN_ISOVALUE
///                 mode is picked for @a mode. When SWEEP_GREATER_THAN_ISOVALUE is supplied
///                 as an argument for @a mode, the extension field voxel will default
///                 to the value of the @a extGrid in that position if it corresponds to a level-set
///                 value that is less than the isovalue. Otherwise, the extension
///                 field voxel value will be computed by the Fast Sweeping algorithm.
///                 The opposite convention is implemented when SWEEP_LESS_THAN_ISOVALUE
///                 is supplied as an argument for @a mode.
///
/// @note Strictly speaking a fog volume is normalized to the range [0,1] but this
///       method accepts a scalar volume with an arbritary range, as long as the it
///       includes the @a isoValue.
///
/// @details Topology of output grids are identical to that of the input grid, except
///          active tiles in the input grid will be converted to active voxels
///          in the output grids!
///
/// @warning If @a isoValue does not intersect any active values in
///          @a sdfGrid then a pair of the following grids is returned: The first
///          is a signed distance grid with its active values set to plus or minus
///          infinity depending of whether its input values are above or below @a isoValue.
///          The second grid, which represents the extension field, has all its active
///          values set to @a background.
template<typename SdfGridT, typename ExtOpT, typename ExtValueT>
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);
 
/// @brief Dilates the narrow band of an existing signed distance field by
///        a specified number of voxels (like adding "onion-rings").
///
/// @note This operation is not to be confused with morphological dilation
///       of a level set, which is implemented in LevelSetFilter::offset,
///       and involves actual interface tracking of the narrow band.
///
/// @return A shared pointer to the dilated signed distance field.
///
/// @param sdfGrid  Input signed distance field to be dilated.
///
/// @param dilation Numer of voxels that the narrow band of the input SDF will be dilated.
///
/// @param nn       Stencil-pattern used for dilation
///
/// @param nIter    Number of iterations of the fast sweeping algorithm.
///                 Each iteration performs 2^3 = 8 individual sweeps.
///
/// @param mode     Determines the direction of the dilation. SWEEP_ALL
///                 will dilate in both sides of the signed distance function,
///                 SWEEP_GREATER_THAN_ISOVALUE will dilate in the positive
///                 side of the iso-surface, SWEEP_LESS_THAN_ISOVALUE will dilate
///                 in the negative side of the iso-surface.
///
/// @details Topology will change as a result of this dilation. E.g. if
///          sdfGrid has a width of 3 and @a dilation = 6 then the grid
///          returned by this method is a narrow band signed distance field
///          with a total width of 9 units.
template<typename GridT>
typename GridT::Ptr
dilateSdf(const GridT &sdfGrid,
          int dilation,
          NearestNeighbors nn = NN_FACE,
          int nIter = 1,
          FastSweepingDomain mode = FastSweepingDomain::SWEEP_ALL);
 
/// @brief Fills mask by extending an existing signed distance field into
///        the active values of this input ree of arbitrary value type.
///
/// @return A shared pointer to the masked signed distance field.
///
/// @param sdfGrid  Input signed distance field to be extended into the mask.
///
/// @param mask     Mask used to identify the topology of the output SDF.
///                 Note this mask is assume to overlap with the sdfGrid.
///
/// @param ignoreActiveTiles If false, active tiles in the mask are treated
///                 as active voxels. Else they are ignored.
///
/// @param nIter    Number of iterations of the fast sweeping algorithm.
///                 Each iteration performs 2^3 = 8 individual sweeps.
///
/// @details Topology of the output SDF is determined by the union of the active
///          voxels (or optionally values) in @a sdfGrid and @a mask.
template<typename GridT, typename MaskTreeT>
typename GridT::Ptr
maskSdf(const GridT &sdfGrid,
        const Grid<MaskTreeT> &mask,
        bool ignoreActiveTiles = false,
        int nIter = 1);
 
////////////////////////////////////////////////////////////////////////////////
/// @brief Computes signed distance values from an initial iso-surface and
///        optionally performs velocity extension at the same time. This is
///        done by means of a novel sparse and parallel fast sweeping
///        algorithm based on a first order Godunov's scheme.
///
///        Solves: @f$|\nabla \phi|^2 = 1 @f$
///
/// @warning Note, it is important to call one of the initialization methods before
///          called the sweep function. Failure to do so will throw a RuntimeError.