Envelope.hpp

#pragma once
 
#include <Eigen/Core>
#include <array>
#include <cstdint>
#include <memory>
#include <vector>
 
// clang-format off
#include <wmtk/utils/DisableWarnings.hpp>
#include <fastenvelope/FastEnvelope.h>
#include <fastenvelope/FastEnvelope2D.h>
#include <SimpleBVH/BVH.hpp>
#include <wmtk/utils/EnableWarnings.hpp>
// clang-format on
 
namespace wmtk {
class Envelope
{
public:
    enum class Kind { Uninitialized, Triangles3d, Edges3d, Edges2d };
 
    Envelope(bool exact = true)
        : use_exact(exact) {};
    virtual ~Envelope() = default;
    double eps2 = 1e-6;
    double eps2_edge = 1e-6;
    double sampling_dist = 1e-3;
 
    void set_max_threads(int n) { m_max_threads = n < 1 ? 1 : n; }
    bool use_exact = true;
 
    bool disabled = false;
    virtual void init(
        const std::vector<Eigen::Vector3d>& m_ver,
        const std::vector<Eigen::Vector3i>& m_faces,
        const double);
    void init(
        const std::vector<Eigen::Vector3d>& m_ver,
        const std::vector<Eigen::Vector2i>& m_edges,
        const double);
    void init(
        const std::vector<Eigen::Vector2d>& m_ver,
        const std::vector<Eigen::Vector2i>& m_edges,
        const double);
 
    void init(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F, const double eps);
    virtual bool is_outside(const std::array<Eigen::Vector3d, 3>& tris) const;
    bool is_outside(const std::array<Eigen::Vector3d, 2>& edge) const;
    virtual bool is_outside(const std::array<Eigen::Vector2d, 2>& edge) const;
    virtual bool is_outside(const Eigen::Vector3d& pts) const;
    virtual bool is_outside(const Eigen::Vector2d& pts) const;
    double nearest_point(const Eigen::Vector3d& pts, Eigen::Vector3d& result) const;
    double nearest_point(const Eigen::Vector2d& pts, Eigen::Vector2d& result) const;
 
    double nearest_point_feature(
        const Eigen::Vector2d& p,
        Eigen::Vector2d& result,
        bool& on_corner,
        Eigen::Vector2d& seg_normal,
        int& feature_id) const;
 
    double nearest_point_feature(
        const Eigen::Vector3d& p,
        Eigen::Vector3d& result,
        int& feature_dim,
        Eigen::Vector3d& dir,
        long long& feature_id) const;
    bool initialized() { return m_bvh != nullptr; };
 
    Kind kind() const { return m_kind; }
 
    double squared_distance(const Eigen::Vector3d& p) const;
    double squared_distance(const Eigen::Vector2d& p) const;
 
    // Per-primitive queries, for a caller that needs the distance d_i to each input primitive i
    // on its own rather than to their union: each d_i is convex (the distance to a convex set),
    // and the min over i is not, so a bound that rests on convexity needs the pieces.
    //
    // A primitive is one entry of the list given to init(), and its id is its index in that
    // list: a triangle for Triangles3d, a segment for Edges3d and Edges2d -- including the
    // degenerate (i, i) entries that stand for isolated points. squared_distance(), and
    // nearest_point(), which returns the same squared distance, is the square of the min over
    // all primitives of primitive_distance(), up to rounding: the BVH computes its own foot point
    // with different arithmetic.
    //
    // Each query comes in the dimension of the init() that built the envelope: the 3D ones on
    // Triangles3d and Edges3d, the 2D ones on Edges2d. Either asked of the other kind, or of an
    // envelope no init() has filled in, throws rather than answering against the wrong copy of
    // the input or dereferencing an empty BVH.
 
    void primitives_in_box(
        const Eigen::Vector3d& lo,
        const Eigen::Vector3d& hi,
        std::vector<int64_t>& ids) const;
    void primitives_in_box(
        const Eigen::Vector2d& lo,
        const Eigen::Vector2d& hi,
        std::vector<int64_t>& ids) const;
 
    double primitive_distance(
        int64_t id,
        const Eigen::Vector3d& p,
        Eigen::Vector3d* grad = nullptr,
        Eigen::Matrix3d* hess = nullptr) const;
    double primitive_distance(
        int64_t id,
        const Eigen::Vector2d& p,
        Eigen::Vector2d* grad = nullptr,
        Eigen::Matrix2d* hess = nullptr) const;
 
    double primitive_triangle_distance(int64_t id, const std::array<Eigen::Vector3d, 3>& tri) const;
 
    double primitive_segment_distance(int64_t id, const std::array<Eigen::Vector2d, 2>& seg) const;
 
private:
    void require_primitive_query(int dim, const char* query) const;
    void require_primitive(int64_t id, int dim, const char* query) const;
    void require_exact_kind(Kind expected, const char* query) const;
    void require_exact_3d(const char* query) const;
    void require_exact_built(const char* query) const;
 
    template <typename VertexList>
    void
    init_exact_edges(const VertexList& V, const std::vector<Eigen::Vector2i>& F, const double _eps);
 
    std::vector<int> geo_vertex_ind;
    std::vector<int> geo_face_ind;
    std::vector<Eigen::Vector2d> m_v2;
    std::vector<Eigen::Vector2i> m_e2;
    std::vector<Eigen::Vector3d> m_v3;
    std::vector<Eigen::Vector3i> m_f3;
    std::vector<Eigen::Vector2i> m_e3;
    std::shared_ptr<SimpleBVH::BVH> m_bvh;
 
private:
    fastEnvelope::FastEnvelope exact_envelope;
    fastEnvelope::FastEnvelope2D exact_envelope_2d;
    Kind m_kind = Kind::Uninitialized;
    bool m_exact_built = false;
 
    int m_max_threads = 1; 
    bool samples_outside_parallel(const std::vector<Eigen::Vector3d>& ps) const;
};
} // namespace wmtk