geometry/math.hh

#ifndef VTESTBED_GEOMETRY_MATH_HH
#define VTESTBED_GEOMETRY_MATH_HH

#include <ostream>
#include <limits>

#include <vtestbed/geometry/types.hh>

namespace vtb {

    namespace geometry {

        template <class T, int N>
        inline Vertex<T,N>
        unit(const Vertex<T,N>& v) {
            auto l = length(v);
            if (l == T{}) { return v; }
            return v / l;
        }

        template <class T>
        inline Vertex<T,3>
        surface_normal(
            const Vertex<T,3>& v0,
            const Vertex<T,3>& v1,
            const Vertex<T,3>& v2
        ) {
            using blitz::cross;
            Vertex<T,3> a = v1-v0;
            Vertex<T,3> b = v2-v0;
            Vertex<T,3> n = cross(a, b);
            return unit(n);
        }

        template <class T>
        inline Vertex<T,2>
        line_normal(const Vertex<T,2>& v0, const Vertex<T,2>& v1) {
            Vertex<T,2> delta = v1-v0;
            return unit(Vertex<T,2>(delta(1),-delta(0)));
        }

        template <class T>
        inline Vertex<T,1>
        orientation(
            const Vertex<T,2>& origin,
            const Vertex<T,2>& a,
            const Vertex<T,2>& b
        ) {
            using vec = Vertex<T,2>;
            return cross(vec(a-origin), vec(b-origin));
        }

        template <class T>
        inline T
        det(const Vertex<T,2>& a, const Vertex<T,2>& b) {
            return cross(a,b)(0);
        }

        template <class T>
        inline T
        det(const Vertex<T,2>& origin, const Vertex<T,2>& a, const Vertex<T,2>& b) {
            return det(Vertex<T,2>(a-origin), Vertex<T,2>(b-origin));
        }

        template <class T>
        inline T
        det(const Vertex<T,3>& a, const Vertex<T,3>& b, const Vertex<T,3>& c) {
            return dot(a,cross(b,c));
        }

        template <class T>
        inline bool
        clockwise(
            const Vertex<T,2>& origin,
            const Vertex<T,2>& a,
            const Vertex<T,2>& b
        ) { return orientation(origin, a, b)(0) < 0; }

        template <class Container, class T>
        inline int
        winding_number(const Container& vertices, const Vertex<T,2>& origin) {
            int wn = 0;
            const auto n = vertices.size();
            for (size_t i=0,j=n-1; i<n; j=i++) {
                const auto& a = vertices[j];
                const auto& b = vertices[i];
                auto orient = orientation(origin, a, b);
                if (a(1) <= origin(1)) {
                    if (origin(1) < b(1) && all(orient > 0)) { ++wn; }
                } else {
                    if (b(1) <= origin(1) && all(orient < 0)) { --wn; }
                }
            }
            return wn;
        }

        template <class T, int N>
        inline T
        cosine(const Vertex<T,N>& a, const Vertex<T,N>& b) {
            using std::sqrt;
            return dot(a,b) / sqrt(dot(a,a)*dot(b,b));
        }

        template <class Figure>
        typename Figure::value_type
        centroid(const Figure& figure) {
            using value_type = typename Figure::value_type;
            using size_type = typename Figure::size_type;
            const size_type n = figure.size();
            value_type sum; sum = 0;
            if (n == 0) { return sum; }
            for (const auto& v : figure) { sum += v; }
            return sum / n;
        }

        template <class Figure>
        typename Figure::value_type
        perimeter(const Figure& figure) {
            using size_type = typename Figure::size_type;
            using scalar_type = typename Figure::scalar_type;
            const size_type n = figure.size();
            scalar_type sum{};
            if (n <= 1) { return sum; }
            for (size_type i=0, j=n-1; i<n; j=i++) {
                sum += distance(figure[j], figure[i]);
            }
            return sum;
        }

        template <class Figure>
        bool
        border_contains(
            const Figure& figure,
            const typename Figure::value_type& vertex,
            const typename Figure::scalar_type eps
        ) {
            using size_type = typename Figure::size_type;
            using scalar_type = typename Figure::scalar_type;
            const size_type n = figure.size();
            constexpr const int N = Figure::dimensions;
            if (n == 0) { return false; }
            if (n == 1) { return all(abs(figure[0] - vertex) < eps); }
            for (size_type i=0,j=n-1; i<n; j=i++) {
                Line_segment<scalar_type,N> s(figure[j], figure[i]);
                if (s.contains(vertex, eps)) { return true; }
            }
            return false;
        }

        template <class Figure>
        bool
        interior_contains(
            const Figure& figure,
            const typename Figure::value_type& vertex
        ) {
            return winding_number(figure, vertex) != 0;
        }

        template <class Figure>
        bool
        contains(
            const Figure& figure,
            const typename Figure::value_type& vertex,
            const typename Figure::scalar_type eps
        ) {
            return border_contains(figure, vertex, eps) ||
                interior_contains(figure, vertex);
        }

        template <class Figure>
        typename std::enable_if<is_figure<Figure>::value, std::ostream&>::type
        operator<<(std::ostream& out, const Figure& figure) {
            using size_type = typename Figure::size_type;
            const size_type n = figure.size();
            if (n == 0) { return out; }
            out << figure[0];
            for (size_type i=1; i<n; ++i) { out << ' ' << figure[i]; }
            return out;
        }

        template <class Figure>
        auto
        bounding_box(const Figure& figure) ->
        Rectangle<typename Figure::scalar_type,Figure::dimensions> {
            constexpr const int N = Figure::dimensions;
            using T = typename Figure::scalar_type;
            using vertex_type = typename Figure::value_type;
            vertex_type min{std::numeric_limits<T>::max()};
            vertex_type max{std::numeric_limits<T>::lowest()};
            for (const auto& v : figure) {
                for (int i=0; i<N; ++i) {
                    T m = v(i);
                    if (m < min(i)) { min(i) = m; }
                    if (max(i) < m) { max(i) = m; }
                }
            }
            return Rectangle<T,N>{min,max};
        }

        namespace bits {

            template <int dimension, int ndimensions, int degrees, class T>
            struct Rotate;

            #define VTB_MAKE_ROTATE(dim, ndim, deg, ...) \
                template <class T> \
                struct Rotate<dim,ndim,deg,T> { \
                    static inline Vertex<T,ndim> \
                    rotate(const Vertex<T,ndim>& v) { \
                        return {__VA_ARGS__}; \
                    } \
                }

            // 3d x
            VTB_MAKE_ROTATE(0, 3, 90,  v(0), -v(2),  v(1));
            VTB_MAKE_ROTATE(0, 3, 180, v(0), -v(1), -v(2));
            VTB_MAKE_ROTATE(0, 3, 270, v(0),  v(2), -v(1));
            // 3d y
            VTB_MAKE_ROTATE(1, 3, 90,   v(2), v(1), -v(0));
            VTB_MAKE_ROTATE(1, 3, 180, -v(0), v(1), -v(2));
            VTB_MAKE_ROTATE(1, 3, 270, -v(2), v(1),  v(0));
            // 3d z
            VTB_MAKE_ROTATE(2, 3, 90,  -v(1),  v(0), v(2));
            VTB_MAKE_ROTATE(2, 3, 180, -v(0), -v(1), v(2));
            VTB_MAKE_ROTATE(2, 3, 270,  v(1), -v(0), v(2));

            // 2d
            VTB_MAKE_ROTATE(0, 2, 90,  -v(1),  v(0));
            VTB_MAKE_ROTATE(0, 2, 180, -v(0), -v(1));
            VTB_MAKE_ROTATE(0, 2, 270,  v(1), -v(0));

            #undef VTB_MAKE_ROTATE

        }

        template <int dimension, int degrees, class T, int N>
        inline Vertex<T,N>
        rotate(const Vertex<T,N>& v) {
            return bits::Rotate<dimension,N,degrees,T>::rotate(v);
        }

    }

}

#endif // vim:filetype=cpp