Composite.h

Go to the documentation of this file.

#ifndef PBAT_GEOMETRY_SDF_COMPOSITE_H
#define PBAT_GEOMETRY_SDF_COMPOSITE_H
 
#include "BinaryNode.h"
#include "PhysicsBasedAnimationToolkitExport.h"
#include "Primitive.h"
#include "Transform.h"
#include "TypeDefs.h"
#include "UnaryNode.h"
#include "pbat/HostDevice.h"
 
#include <algorithm>
#include <array>
#include <span>
#include <type_traits>
#include <utility>
#include <variant>
#include <vector>
 
namespace pbat::geometry::sdf {

template <common::CArithmetic TScalar>
using Node = std::variant<
    /* primitives */
    Sphere<TScalar>,
    Box<TScalar>,
    BoxFrame<TScalar>,
    Torus<TScalar>,
    CappedTorus<TScalar>,
    Link<TScalar>,
    InfiniteCylinder<TScalar>,
    Cone<TScalar>,
    InfiniteCone<TScalar>,
    Plane<TScalar>,
    HexagonalPrism<TScalar>,
    Capsule<TScalar>,
    VerticalCapsule<TScalar>,
    CappedCylinder<TScalar>,
    VerticalCappedCylinder<TScalar>,
    RoundedCylinder<TScalar>,
    VerticalCappedCone<TScalar>,
    CutHollowSphere<TScalar>,
    VerticalRoundCone<TScalar>,
    Octahedron<TScalar>,
    Pyramid<TScalar>,
    Triangle<TScalar>,
    Quadrilateral<TScalar>,
    /* unary nodes */
    Scale<TScalar>,
    Elongate<TScalar>,
    Round<TScalar>,
    Onion<TScalar>,
    Symmetrize<TScalar>,
    Repeat<TScalar>,
    RotationalRepeat<TScalar>,
    Bump<TScalar>,
    Twist<TScalar>,
    Bend<TScalar>,
    /* binary nodes */
    Union<TScalar>,
    Difference<TScalar>,
    Intersection<TScalar>,
    ExclusiveOr<TScalar>,
    SmoothUnion<TScalar>,
    SmoothDifference<TScalar>,
    SmoothIntersection<TScalar>>;

enum class ECompositeStatus {
    Valid,             
    InvalidForest,     
    UnexpectedNodeType 
};

template <common::CArithmetic TScalar>
struct Composite
{
    using ScalarType = TScalar; 
    PBAT_HOST_DEVICE
    Composite(
        std::span<Node<ScalarType> const> nodes,
        std::span<Transform<ScalarType> const> transforms,
        std::span<std::pair<int, int> const> children,
        std::span<int const> roots);
    PBAT_HOST_DEVICE
    ScalarType Eval(Vec3<ScalarType> const& p) const;
    PBAT_HOST_DEVICE
    ScalarType Eval(int n, Vec3<ScalarType> const& p) const;
    PBAT_HOST_DEVICE
    Vec3<ScalarType> Grad(Vec3<ScalarType> const& p, ScalarType h) const;
    PBAT_HOST_DEVICE
    ECompositeStatus Status() const { return mStatus; }
 
    std::span<Node<ScalarType> const> mNodes; 
    std::span<Transform<ScalarType> const>
        mTransforms; 
    std::span<std::pair<int, int> const>
        mChildren; 
    std::span<int const> mRoots; 
    ECompositeStatus mStatus;    
};

PBAT_API auto FindRootsAndParents(std::span<std::pair<int, int> const> children)
    -> std::pair<std::vector<int>, std::vector<int>>;
 
template <common::CArithmetic TScalar>
inline Composite<TScalar>::Composite(
    std::span<Node<ScalarType> const> nodes,
    std::span<Transform<ScalarType> const> transforms,
    std::span<std::pair<int, int> const> children,
    std::span<int const> roots)
    : mNodes(nodes),
      mTransforms(transforms),
      mChildren(children),
      mRoots(roots),
      mStatus(ECompositeStatus::Valid)
{
    std::size_t const nNodes = mNodes.size();
    for (std::size_t n = 0ULL; n < nNodes; ++n)
    {
        if (mStatus != ECompositeStatus::Valid)
            break;
        std::visit(
            [&](auto const& node) {
                using NodeType = std::remove_cvref_t<decltype(node)>;
                if constexpr (std::is_base_of_v<Primitive, NodeType>)
                {
                    if (mChildren[n].first >= 0 || mChildren[n].second >= 0)
                        mStatus = ECompositeStatus::InvalidForest;
                }
                else if constexpr (std::is_base_of_v<UnaryNode, NodeType>)
                {
                    if (mChildren[n].first < 0)
                        mStatus = ECompositeStatus::InvalidForest;
                    if (mChildren[n].second >= 0)
                        mStatus = ECompositeStatus::InvalidForest;
                }
                else if constexpr (std::is_base_of_v<BinaryNode, NodeType>)
                {
                    if (mChildren[n].first < 0 || mChildren[n].second < 0)
                        mStatus = ECompositeStatus::InvalidForest;
                }
            },
            mNodes[n]);
    }
}
 
template <common::CArithmetic TScalar>
inline TScalar Composite<TScalar>::Eval(Vec3<ScalarType> const& p) const
{
    using namespace std;
    TScalar sd = numeric_limits<TScalar>::max();
    for (auto root : mRoots)
        sd = min(sd, Eval(root, p));
    return sd;
}
 
template <common::CArithmetic TScalar>
inline TScalar Composite<TScalar>::Eval(int n, Vec3<ScalarType> const& p) const
{
    std::size_t const nStl = static_cast<std::size_t>(n);
    ScalarType sd;
    Vec3<ScalarType> q = mTransforms[nStl] / p;
    std::visit(
        [&](auto const& node) {
            using NodeType = std::remove_cvref_t<decltype(node)>;
            if constexpr (std::is_base_of_v<Primitive, NodeType>)
            {
                sd = node.Eval(q);
            }
            else if constexpr (std::is_base_of_v<UnaryNode, NodeType>)
            {
                int const child = mChildren[nStl].first;
                sd              = node.Eval(q, [this, child](Vec3<ScalarType> const& x) {
                    return Eval(child, x);
                });
            }
            else if constexpr (std::is_base_of_v<BinaryNode, NodeType>)
            {
                auto ci = mChildren[nStl].first;
                auto cj = mChildren[nStl].second;
                sd      = node.Eval(Eval(ci, q), Eval(cj, q));
            }
        },
        mNodes[nStl]);
    return sd;
}
 
template <common::CArithmetic TScalar>
inline Vec3<TScalar> Composite<TScalar>::Grad(Vec3<ScalarType> const& p, ScalarType h) const
{
    ScalarType constexpr one{1};
    ScalarType constexpr none{-1};
    Vec3<ScalarType> g;
    Vec3<ScalarType> k{one, none, none};
    g = k * Eval(p + h * k);
    k = Vec3<ScalarType>{none, none, one};
    g += k * Eval(p + h * k);
    k = Vec3<ScalarType>{none, one, none};
    g += k * Eval(p + h * k);
    k = Vec3<ScalarType>{one, one, one};
    g += k * Eval(p + h * k);
    g *= ScalarType(0.25) / h;
    return g;
}
 
} // namespace pbat::geometry::sdf
 
#endif // PBAT_GEOMETRY_SDF_COMPOSITE_H