triangle.hh

#ifndef VTESTBED_GEOMETRY_TRIANGLE_HH
#define VTESTBED_GEOMETRY_TRIANGLE_HH

#include <iosfwd>

#include <vtestbed/geometry/math.hh>
#include <vtestbed/geometry/plane.hh>
#include <vtestbed/geometry/ray.hh>
#include <vtestbed/geometry/types.hh>

namespace vtb {

    namespace geometry {

        template <class T, int N>
        class Triangle: public std::array<Vertex<T,N>,3> {

        private:
            using base_type = std::array<Vertex<T,N>,3>;

        public:
            static constexpr const int dimensions = N;
            using scalar_type = T;
            using value_type = typename base_type::value_type;
            using pointer = typename base_type::pointer;
            using const_pointer = typename base_type::const_pointer;
            using reference = typename base_type::reference;
            using const_reference = typename base_type::const_reference;
            using iterator = typename base_type::iterator;
            using const_iterator = typename base_type::const_iterator;
            using size_type = typename base_type::size_type;
            using difference_type = typename base_type::difference_type;
            using reverse_iterator = typename base_type::reverse_iterator;
            using const_reverse_iterator = typename base_type::const_reverse_iterator;

        public:
            using base_type::base_type;
            using base_type::operator=;
            using base_type::operator[];

        public:

            inline
            Triangle(const_reference p0, const_reference p1, const_reference p2):
            base_type{p0,p1,p2} {}

            inline void flip() { std::swap(this->front(), this->back()); }
            inline void reverse() { this->flip(); }
            inline void invert() { this->flip(); }

            Plane<T,N> plane() const;
            Polyline<T,N> polyline() const;
            Polygon<T,N> polygon() const;
            bool degenerate(scalar_type eps) const;

            void
            gnuplot(std::ostream& out) const;

        };

        template <class T, int N>
        struct is_figure<Triangle<T,N>>: public std::true_type {};

        template <class T>
        T tetrahedron_signed_volume(const Triangle<T,3>& tr) {
            return dot(tr[2], cross(tr[0], tr[1])) / T{6};
        }

        template <class T>
        T tetrahedron_signed_volume(const Triangle<T,2>& tr) {
            return T{};
        }

        template <class T, int N>
        inline bool
        intersect(const Ray<T,N>& ray, const Triangle<T,N>& triangle, T eps) {
            using vec = Vertex<T,N>;
            vec e1{triangle[1] - triangle[0]}, e2{triangle[2] - triangle[0]};
            vec n{cross(ray.direction(), e2)};
            T det = dot(e1, n);
            if (std::abs(det) < eps) {
                return false;
            }
            T rdet{T{1}/det};
            vec tvec{ray.origin() - triangle[0]};
            T u = dot(tvec, n) * rdet;
            if (u < T{0} || u > T{1}) {
                return false;
            }
            vec qvec{cross(tvec, e1)};
            T v = dot(ray.direction(), qvec) * rdet;
            if (v < T{0} || u+v > T{1}) {
                return false;
            }
            return true;
        }

        template <class T, int N>
        inline Vertex<Plane_position,3>
        compare(const Plane<T,N> plane, const Triangle<T,N>& triangle, T eps) {
            Vertex<Plane_position,3> order;
            for (int i=0; i<3; ++i) {
                order(i) = compare(plane, triangle[i], eps);
            }
            return order;
        }

        template <class T, int N>
        void
        intersect(
            const Plane<T,N>& plane,
            const Triangle<T,N>& triangle,
            const Vertex<Plane_position,3>& positions,
            Polygon<T,N>& front,
            Polygon<T,N>& back
        );

        template <class T, int N>
        inline Vertex<T,N>
        centroid(const Triangle<T,N>& t) {
            return (t[0]+t[1]+t[2])/T{3};
        }

        template <class T, int N>
        inline T
        area(const Triangle<T,N>& t) {
            Vertex<T,N> e1 = t[1] - t[0];
            Vertex<T,N> e2 = t[2] - t[0];
            return T{0.5} * length(cross(e1,e2));
        }

        template <class T>
        inline Vertex<T,3>
        surface_normal(const Triangle<T,3>& t) {
            return surface_normal(t[0],t[1],t[2]);
        }

    }

}

#endif // vim:filetype=cpp