Initial implementation of anisotropic surfacing, including re-work of…

… PointRasterize and other fixes

Signed-off-by: Nick Avramoussis <4256455+Idclip@users.noreply.github.com>

openvdb/openvdb/points/PointTransfer.h

@@ -113,7 +113,7 @@ namespace points {
 ///         , mHandle(nullptr) {}
 ///
 ///     /// @brief Range in index space of the source points
-///     Int32 range(const Coord&, size_t) const { return Int32(1); }
+///     Vec3i range(const Coord&, size_t) const { return Vec3i(1); }
 ///
 ///     /// @brief Every time we start a new point leaf, init the position array.
 ///     ///   Always return true as we don't skip any leaf nodes.
@@ -416,10 +416,12 @@ struct RasterizePoints
 
     RasterizePoints(const points::PointDataTree& tree,
                     const TransferT& transfer,
+                    const CoordBBox& pointBounds,
                     const PointFilterT& filter = PointFilterT(),
                     InterrupterT* interrupter = nullptr)
         : mPointAccessor(tree)
         , mTransfer(transfer)
+        , mPointBounds(pointBounds)
         , mFilter(filter)
         , mInterrupter(interrupter) {}
 
@@ -449,8 +451,12 @@ struct RasterizePoints
 
         mTransfer.initialize(origin, idx, bounds);
 
-        CoordBBox search = bounds.expandBy(mTransfer.range(origin, idx));
+        CoordBBox search = bounds;
+        const Vec3i range(mTransfer.range(origin, idx));
+        search.min() -= Coord(range);
+        search.max() += Coord(range);
         this->transform<>(search);
+        search.intersect(mPointBounds);
 
         // start the iteration from a leaf origin
         const Coord min = (search.min() & ~(DIM-1));
@@ -539,6 +545,7 @@ struct RasterizePoints
 private:
     const PointDataGrid::ConstAccessor mPointAccessor;
     mutable TransferT mTransfer;
+    const CoordBBox& mPointBounds;
     const PointFilterT& mFilter;
     InterrupterT* mInterrupter;
 };
@@ -566,9 +573,14 @@ rasterize(const PointDataTreeOrGridT& points,
 
     auto& topology = transfer.topology();
     using TreeT = typename std::decay<decltype(topology)>::type;
+
+    // Compute max search bounds
+    CoordBBox bounds;
+    tree.evalLeafBoundingBox(bounds);
+
     tree::LeafManager<TreeT> manager(topology);
     transfer_internal::RasterizePoints<TransferT, TreeT, FilterT, InterrupterT>
-        raster(tree, transfer, filter, interrupter);
+        raster(tree, transfer, bounds, filter, interrupter);
     manager.foreach(raster);
 }
 

openvdb/openvdb/points/PrincipleComponentAnalysis.h

@@ -0,0 +1,139 @@
+// Copyright Contributors to the OpenVDB Project
+// SPDX-License-Identifier: MPL-2.0
+
+/// @author Richard Jones, Nick Avramoussis
+///
+/// @file PrincipleComponentAnalysis.h
+///
+/// @brief  Provides methods to perform principle component analysis (PCA) over
+///   a point set to compute rotational and affine transformations for each
+///   point that represent a their neighborhoods anisotropy. The techniques
+///   and algorithms used here are described in:
+///       [Reconstructing Surfaces of Particle-Based Fluids Using Anisotropic
+///        Kernel - Yu Turk 2010].
+///   The parameters and results of these methods can be combines with the
+///   ellipsoidal surfacing technique found in PointRasterizeSDF.h.
+
+#ifndef OPENVDB_POINTS_POINT_PRINCIPLE_COMPONENT_ANALYSIS_HAS_BEEN_INCLUDED
+#define OPENVDB_POINTS_POINT_PRINCIPLE_COMPONENT_ANALYSIS_HAS_BEEN_INCLUDED
+
+#include <openvdb/openvdb.h>
+#include <openvdb/Grid.h>
+#include <openvdb/Types.h>
+#include <openvdb/math/Coord.h>
+#include <openvdb/thread/Threading.h>
+#include <openvdb/util/NullInterrupter.h>
+#include <openvdb/points/PointAttribute.h>
+#include <openvdb/points/PointGroup.h>
+#include <openvdb/points/PointTransfer.h>
+#include <openvdb/points/PointDataGrid.h>
+
+#include <string>
+#include <vector>
+#include <memory>
+#include <cmath> // std::cbrt
+#include <algorithm> // std::accumulate
+
+namespace openvdb {
+OPENVDB_USE_VERSION_NAMESPACE
+namespace OPENVDB_VERSION_NAME {
+namespace points {
+
+struct PcaSettings;
+struct PcaAttributes;
+
+/// @brief  Calculate  Calculate ellipsoid transformations from the local point
+///   distributions as described in Yu and Turk's 'Reconstructing Fluid Surfaces
+///   with Anisotropic Kernels'.
+/// @param points   The point data grid to analyses
+/// @param settings The PCA settings for controlling the behavior of the
+///   neighborhood searches and the resulting ellipsoidal values
+/// @param attrs    The PCA attributes to create
+/// @param interrupt An optional interrupter.
+template <typename PointDataGridT,
+    typename FilterT = NullFilter,
+    typename InterrupterT = util::NullInterrupter>
+inline void
+pca(PointDataGridT& points,
+    const PcaSettings& settings,
+    const PcaAttributes& attrs,
+    InterrupterT* interrupt = nullptr);
+
+
+/// @brief  Various settings for the neighborhood analysis of point
+///   distributions.
+struct PcaSettings
+{
+    /// @param searchRadius  the world space search radius of the neighborhood
+    ///   around each point. Increasing this value will result in points
+    ///   including more of their neighbors into their ellipsoidal calculations.
+    ///   This may or may not be desirable depending on the point set's
+    ///   distribution and can be tweaked as necessary. Note however that large
+    ///   values will cause the PCA calculation to become exponentially more
+    ///   expensive.
+    float searchRadius = 1.0f;
+    /// @param allowedAnisotropyRatio  the maximum allowed ratio between the
+    ///   components in each ellipse' stretch coefficients such that:
+    /// @code
+    ///     const auto s = stretch.sorted();
+    ///     assert(s[0]/s[2] >= allowedAnisotropyRatio);
+    /// @endcode
+    ///   This parameter effectively clamps the allowed anisotropy, with a
+    ///   value of 1.0f resulting in uniform stretch values (representing a
+    ///   a sphere). Values tending towards zero will allow for greater
+    ///   anisotropy i.e. much more exaggerated stretches along the computed
+    ///   principle axis and corresponding squashes along the others to
+    ///   compensate.
+    /// @note Very small values may cause very thing ellipses to be produced,
+    ///   so a reasonable minimum should be set. Values equal to or less than
+    ///   0.0, or values greater than 1.0 have undefined results.
+    float allowedAnisotropyRatio = 0.25f;
+    /// @param neighbourThreshold  the number of points required to have a full
+    ///   neighbourhood. Should never be 0 (as this would simply result in
+    ///   an identity transformations i.e. a canonical sphere).
+    size_t neighbourThreshold = 20;
+    /// @param averagePositions  the amount (between 0 and 1) to average out
+    ///   the ellipsoid positions
+    float averagePositions = 1.0f;
+};
+
+/// @brief  The persistent attributes created by the PCA methods.
+/// @note   These can be passed to points::rasterizeSdf with the
+///   EllipsoidSettings to perform ellipsoidal surface construction.
+struct PcaAttributes
+{
+    /// @brief  Settings for the "stretch" attribute, a floating point vector
+    ///   attribute which represents the scaling components of each points
+    ///   ellipse.
+    using StretchT = math::Vec3<float>;
+    std::string stretch = "_stretch";
+
+    /// @brief  Settings for the "rotation" attribute, a floating point matrix
+    ///   attribute which represents the rotation of each points ellipse.
+    using RotationT = math::Mat3<float>;
+    std::string rotation = "_rotation";
+
+    /// @brief  Settings for the world space position of every point. This may
+    ///   end up being different to their actual position if the
+    ///   PcaSettings::averagePositions value is not exactly 0. This attribute
+    ///   is used as the points position when calling points::rasterizeSdf.
+    /// @note  This should always be at least at double precision
+    using PosWsT = math::Vec3<double>;
+    std::string positionWS = "_pws";
+
+    /// @brief  inclusion  A point group to create that represents points which
+    ///   have valid ellipsoidal neighborhood. Points outside of this group
+    ///   will have their stretch and rotation attributes set to describe a
+    ///   canonical sphere. Note however that all points, regardless of this
+    ///   groups membership flag, may still have their world space position
+    ///   deformed in relation to their neighboring points.
+    std::string inclusion = "_ellipsoids";
+};
+
+}
+}
+}
+
+#include "impl/PrincipleComponentAnalysisImpl.h"
+
+#endif // OPENVDB_POINTS_POINT_PRINCIPLE_COMPONENT_ANALYSIS_HAS_BEEN_INCLUDED