PrincipalComponentAnalysis.h

Go to the documentation of this file.

// Copyright Contributors to the OpenVDB Project
// SPDX-License-Identifier: Apache-2.0
 
/// @author Richard Jones, Nick Avramoussis
///
/// @file PrincipalComponentAnalysis.h
///
/// @brief  Provides methods to perform principal 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_PRINCIPAL_COMPONENT_ANALYSIS_HAS_BEEN_INCLUDED
#define OPENVDB_POINTS_POINT_PRINCIPAL_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/util/Assert.h>
#include <openvdb/simd/Simd.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 <limits>
#include <cmath> // std::cbrt
 
namespace openvdb {
OPENVDB_USE_VERSION_NAMESPACE
namespace OPENVDB_VERSION_NAME {
namespace points {
 
struct PcaSettings;
struct PcaAttributes;
 
/// @brief  Calculate ellipsoid transformations from the local point
///   distributions as described in Yu and Turk's 'Reconstructing Fluid Surfaces
///   with Anisotropic Kernels'. The results are stored on the attributes
///   pointed to by the PcaAttributes. See the PcaSettings and PcaAttributes
///   structs for more details.
/// @warning  This method will throw if the 'strech', 'xform' or 'pws'
///   attributes already exist on this input point set.
/// @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
template <typename PointDataGridT>
inline void
pca(PointDataGridT& points,
    const PcaSettings& settings,
    const PcaAttributes& attrs);
 
 
/// @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 and should be used in conjunction with the max point per
    ///   voxel settings below.
    /// @note  This value is used to normalize, so you should not use maximum
    ///   limit values to try and include all points. If, for example, you
    ///   want include every point against eachother, consider computing the
    ///   extents of the grid and using that as the searchRadius.
    /// @warning  Valid range is [0, inf). Behaviour is undefined when outside
    ///   this range.
    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
    ///   sphere). Values tending towards zero will allow for greater
    ///   anisotropy i.e. much more exaggerated stretches along the computed
    ///   principal 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. Valid range is (0, 1].
    ///   Behaviour is undefined when outside this range.
    float allowedAnisotropyRatio = 0.25f;
    /// @param nonAnisotropicStretch  The stretch coefficient that should be
    ///   used for points which have no anisotropic neighbourhood (due to
    ///   being isolated or not having enough neighbours to reach the
    ///   specified @sa neighbourThreshold).
    float nonAnisotropicStretch = 1.0;
    /// @param neighbourThreshold  number of points in a given neighbourhood
    ///   that target points must have to be classified as having an elliptical
    ///   distribution. Points with less neighbours than this will end up with
    ///   uniform stretch values of nonAnisotropicStretch and an identity
    ///   transformation.
    /// @note  This number can include the target point itself.
    /// @warning  Changing this value does not change the size or number of
    ///   point neighborhood lookups. As such, increasing this value only
    ///   speeds up the number of calculated covariance matrices. Valid range
    ///   is [0, inf], where 0 effectively disables all covariance calculations
    ///   and 1 enables them for every point.
    size_t neighbourThreshold = 20;
    /// @param The maximum number of source points to gather for contributions
    ///   to each target point, per voxel. When a voxel contains more points
    ///   than this value, source point are trivially stepped over, with the
    ///   step size calculated as max(1, ppv / neighbourThreshold).
    /// @note  There is no prioritisation of nearest points; for example, voxels
    ///   which partially overlap the search radius may end up selecting point
    ///   ids which all lie outside. This is rarely an issue in practice and
    ///   choosing an appropriate value for this setting can significantly
    ///   increase performance without a large impact to the results.
    /// @warning  Valid range is [1, inf). Behaviour is undefined if this value
    ///   is set to zero.
    size_t maxSourcePointsPerVoxel = 8;
    /// @param The maximum number of target points to write anisotropy values
    ///   to. When a voxel contains more points than this value, target point
    ///   are trivially stepped over, with the step size calculated as:
    ///      max(1, ppv / maxTargetPointsPerVoxel).
    ///   default behaviour is to compute for all points. Any points skipped
    ///   will be treated as being isolated and receive an identity transform
    ///   and nonAnisotropicStretch.
    /// @note  When using in conjuction with rasterizeSdf for ellipsoids,
    ///   consider choosing a lower value (e.g. same value as
    ///   maxSourcePointsPerVoxel) to speed up iterations and only
    ///   increase if necessary.
    /// @warning  Valid range is [1, inf). Behaviour is undefined if this value
    ///   is set to zero.
    size_t maxTargetPointsPerVoxel = std::numeric_limits<size_t>::max();
    /// @param averagePositions  the amount (between 0 and 1) to average out
    ///   positions. All points, whether they end up as ellipses or not,
    ///   can have their positions smoothed to account for their neighbourhood
    ///   distribution. This will have no effect for points with no neighbours.
    /// @warning  This options does NOT modify the P attribute - instead, a
    ///   world space position attribute "pws" (default) is created and stores
    ///   the smoothed position.
    /// @warning  Valid range is [0, 1]. Behaviour is undefined when outside
    ///   this range.
    float averagePositions = 1.0f;
    /// @param interrupter  optional interrupter
    util::NullInterrupter* interrupter = nullptr;
};
 
// Suppress spurious warnings on compiler emitted methods (constructors, etc)
// due to deprecated members. Accessing said members still generates the warning.
OPENVDB_NO_DEPRECATION_WARNING_BEGIN
 
/// @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
{
    enum class XformOutput
    {
        // When selected, the "stretch" output attribute won't be created.
        // Instead, only a 3x3 affine unitary "xform" attribute will be
        // created with the combined stretch and rotational transformations
        COMBINED_TRANSFORM,
        // When selected, both the "stretch" vector and "xform" output
        // attributes will be created. The xform will take the form of a
        // quaternion representing the rotation.
        STRETCH_AND_QUATERNION,
        // When selected, both the "stretch" vector and "xform" output
        // attributes will be created. The xform will take the form of an
        // orthogonal 3x3 matrix representing the rotation OR a pure
        // reflection.
        STRETCH_AND_UNITARY_MATRIX
    };
 
    XformOutput xformOutput = XformOutput::STRETCH_AND_UNITARY_MATRIX;
 
    /// @brief  Settings for the "stretch" attribute, a floating point vector
    ///   attribute which represents the scaling components of each points
    ///   ellipse or (1.0,1.0,1.0) for isolated points.
    /// @note   This will only be created if the XformOutput is set to
    ///   STRETCH_AND_QUATERNION or STRETCH_AND_UNITARY_MATRIX.
    using StretchT = math::Vec3<float>;
    std::string stretch = "stretch";
 
    /// @brief  Settings for the "xform" attribute, either a floating point
    ///   matrix or quaternion attribute. See the `xformOutput` setting and
    ///   XformOutput enum above for details.
    using RotationT = math::Mat3<float>;
    using QuatT = math::Quat<float>;
    std::string xform = "xform";
 
    /// @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 in place of the point's actual 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 A point group to create that represents points which have valid
    ///   ellipsoidal neighborhood. Points outside of this group will have
    ///   their stretch and xform attributes set to describe a canonical sphere
    ///   Note however that all points, regardless of this groups membership
    ///   flag, will still contribute to their neighbours and may have their
    ///   world space position deformed in relation to their neighboring points.
    std::string ellipses = "ellipsoids";
 
    /// Old style "xform" attribute which has since been renamed and ONLY
    /// supported rotation. If set, both the above "xform" and this attribute
    /// will be created. This attribute only ever holds rotation, regardless
    /// of the xformOutput specified. This will be removed and should not be
    /// used.
    OPENVDB_DEPRECATED_MESSAGE("Use PcaAttributes::xform") std::string rotation = "";
};
 
OPENVDB_NO_DEPRECATION_WARNING_END
 
}
}
}
 
#include "impl/PrincipalComponentAnalysisImpl.h"
 
#endif // OPENVDB_POINTS_POINT_PRINCIPAL_COMPONENT_ANALYSIS_HAS_BEEN_INCLUDED