tetrahedron.h

#ifndef TETRAHEDRON_H
#define TETRAHEDRON_H

// GLM
#include <glm/glm.hpp>

// CGAL
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/Bbox_3.h>

// Local
#include "libvcad/utils/math_utils.h"
#include "libvcad/geometry/primitives.h"

//! \class Tetrahedron
//! \brief A class that represents a tetrahedron in 3D space
//! \details This class is used to represent a tetrahedron in 3D space. The tetrahedron is defined by four points in 3D space.
//!          The class also stores a value (template) associated with each point in the tetrahedron. Methods are provided to interpolate given
//!          a point in space within the tetrahedron.
template <typename T>
class Tetrahedron
{
    public:
        //! \brief Constructor
        Tetrahedron() = default;

        //! \brief Constructor
        //! \param a The first point of the tetrahedron
        //! \param b The second point of the tetrahedron
        //! \param c The third point of the tetrahedron
        //! \param d The fourth point of the tetrahedron
        //! \param a_val The value associated with the first point
        //! \param b_val The value associated with the second point
        //! \param c_val The value associated with the third point
        //! \param d_val The value associated with the fourth point
        Tetrahedron(const glm::vec3& a, const glm::vec3& b, const glm::vec3& c, const glm::vec3& d,
                    const T& a_val, const T& b_val, const T& c_val, const T& d_val) :
            m_a(a), m_b(b), m_c(c), m_d(d), m_a_val(a_val), m_b_val(b_val), m_c_val(c_val), m_d_val(d_val) {}

        //! \brief Compute the volume of the tetrahedron
        //! \return The volume of the tetrahedron in mm^3
        double volume() const
        {
            return MathUtils::tet_volume(m_a, m_b, m_c, m_d);
        }

        //! \brief Interpolate a value at a point in space within the tetrahedron
        //! \param point The point in space to interpolate the value at
        //! \return The interpolated value at the point
        T sample(const glm::vec3& point) const
        {
            glm::vec4 bary_coords = barycentric_coordinates(point);
            return MathUtils::tet_interpolate(bary_coords, m_a_val, m_b_val, m_c_val, m_d_val);
        }

        //! \brief Compute the barycentric coordinates of a point in space within the tetrahedron
        //! \param point The point in space to compute the barycentric coordinates for
        //! \return The barycentric coordinates of the point
        glm::vec4 barycentric_coordinates(const glm::vec3& point) const
        {
            float V1234 =   MathUtils::tet_volume(m_a, m_b, m_c, m_d);
            float lambda1 = MathUtils::tet_volume(point, m_b, m_c, m_d) / V1234;
            float lambda2 = MathUtils::tet_volume(m_a, point, m_c, m_d) / V1234;
            float lambda3 = MathUtils::tet_volume(m_a, m_b, point, m_d) / V1234;
            float lambda4 = MathUtils::tet_volume(m_a, m_b, m_c, point) / V1234;

            return {lambda1, lambda2, lambda3, lambda4};
        }

        //! \brief Check if a point is inside the tetrahedron
        //! \note depends on barycentric_coordinates()
        //! \param point The point to check
        //! \return True if the point is inside the tetrahedron, false otherwise
        bool is_point_inside(const glm::vec3& point) const
        {
            glm::vec4 bary_coords = barycentric_coordinates(point);
            double sum = (bary_coords.x + bary_coords.y + bary_coords.z + bary_coords.w);
            double epsilon = std::numeric_limits<float>::epsilon() * 10;
            auto inside = (bary_coords.x >= -epsilon
                        && bary_coords.y >= -epsilon
                        && bary_coords.z >= -epsilon
                        && bary_coords.w >= -epsilon
                        && glm::abs(sum - 1.0) <= epsilon);
            return inside;
        }

        //! \brief Compute the minimum point of the tetrahedron
        //! \return The minimum point of the tetrahedron
        glm::vec3 min_point() const
        {
            return glm::vec3(
                    std::min({m_a.x, m_b.x, m_c.x, m_d.x}),
                    std::min({m_a.y, m_b.y, m_c.y, m_d.y}),
                    std::min({m_a.z, m_b.z, m_c.z, m_d.z})
            );
        }

        //! \brief Compute the maximum point of the tetrahedron
        //! \return The maximum point of the tetrahedron
        glm::vec3 max_point() const
        {
            return glm::vec3(
                    std::max({m_a.x, m_b.x, m_c.x, m_d.x}),
                    std::max({m_a.y, m_b.y, m_c.y, m_d.y}),
                    std::max({m_a.z, m_b.z, m_c.z, m_d.z})
            );
        }

        //! \brief Compute the bounding box of the tetrahedron
        //! \return A pair of the minimum and maximum points of the bounding box
        std::pair<glm::vec3, glm::vec3> bounding_box() const
        {
            return std::make_pair(min_point(), max_point());
        }

        //! \brief Get the first point of the tetrahedron
        glm::vec3 a() const { return m_a; }

        //! \brief Get the second point of the tetrahedron
        glm::vec3 b() const { return m_b; }

        //! \brief Get the third point of the tetrahedron
        glm::vec3 c() const { return m_c; }

        //! \brief Get the fourth point of the tetrahedron
        glm::vec3 d() const { return m_d; }

        //! \brief Set the first point of the tetrahedron
        void set_a(const glm::vec3& a) { m_a = a; }

        //! \brief Set the second point of the tetrahedron
        void set_b(const glm::vec3& b) { m_b = b; }

        //! \brief Set the third point of the tetrahedron
        void set_c(const glm::vec3& c) { m_c = c; }

        //! \brief Set the fourth point of the tetrahedron
        void set_d(const glm::vec3& d) { m_d = d; }

        //! \brief Get the value associated with the first point
        T a_val() const { return m_a_val; }

        //! \brief Get the value associated with the second point
        T b_val() const { return m_b_val; }

        //! \brief Get the value associated with the third point
        T c_val() const { return m_c_val; }

        //! \brief Get the value associated with the fourth point
        T d_val() const { return m_d_val; }

        //! \brief Set the value associated with the first point
        void set_a_val(const T& val) { m_a_val = val; }

        //! \brief Set the value associated with the second point
        void set_b_val(const T& val) { m_b_val = val; }

        //! \brief Set the value associated with the third point
        void set_c_val(const T& val) { m_c_val = val; }

        //! \brief Set the value associated with the fourth point
        void set_d_val(const T& val) { m_d_val = val; }

        //! \brief Compute the barycenter of the tetrahedron
        //! \return The barycenter or center of mass of the tetrahedron
        glm::vec3 barycenter() const
        {
            return (m_a + m_b + m_c + m_d) / 4.0f;
        }

    private:
        glm::vec3 m_a;
        glm::vec3 m_b;
        glm::vec3 m_c;
        glm::vec3 m_d;

        T m_a_val;
        T m_b_val;
        T m_c_val;
        T m_d_val;
};

//! \class TetrahedronPrimitive
//! \brief A class that wraps a Tetrahedron object to be used in the AABB trees
template <typename T>
struct TetrahedronPrimitive
{
    typedef const Tetrahedron<T>* Id;
    typedef SC_Kernel::Iso_cuboid_3 Datum;
    typedef SC_Point Point;

    Tetrahedron<T>* m_ptr;

    //! \brief Constructor that wraps a Tetrahedron object
    //! \param tet The pointer to the Tetrahedron object
    TetrahedronPrimitive(Tetrahedron<T>* tet) : m_ptr(tet) {}

    //! \brief Get the ID of the Tetrahedron object. The ID is simply the pointer memory address
    //! \note Required by CGAL AABBTraits
    //! \return The ID of the Tetrahedron object
    Id id() const { return m_ptr; }

    //! \brief Get the bounding box of the Tetrahedron object
    //! \note Required by CGAL AABBTraits
    //! \return The bounding box of the Tetrahedron object
    Datum datum() const
    {
        auto bbox_coords = m_ptr->bounding_box();
        return Datum(Point(bbox_coords.first.x, bbox_coords.first.y, bbox_coords.first.z),
                     Point(bbox_coords.second.x, bbox_coords.second.y, bbox_coords.second.z));
    }

    //! \brief Get the reference point of the Tetrahedron object. This can be any point that represents the object,
    //!        in this case, the barycenter of the Tetrahedron object is used.
    //! \note Required by CGAL AABBTraits
    //! \return The reference point of the Tetrahedron object
    Point reference_point() const
    {
        auto center_of_mass = m_ptr->barycenter();
        return Point(center_of_mass.x, center_of_mass.y, center_of_mass.z);
    }
};

//! \brief Define the AABB traits for the TetrahedronPrimitive
typedef CGAL::AABB_traits<SC_Kernel, TetrahedronPrimitive<glm::vec3>> AABB_tetrahedron_traits;

#endif // TETRAHEDRON_H