/*---------------------------------------------------------------------------*\

    geophase: A C++ template library for geometrical operations on non-convex
    polyhedra with non-planar faces
     
    Copyright (C) 2019 Tomislav Maric, TU Darmstadt 

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.

Description
    Functions used to either generate common geometrical models, or make one
    model from the other. 

Author
    Tomislav Maric maric@mma.tu-darmstadt.de

\*---------------------------------------------------------------------------*/

#ifndef Make_HPP
#define Make_HPP 

#include "Polygon.hpp"
#include "AxisAlignedBoundingBox.hpp"
#include "Polyhedron.hpp"

#include <cmath>

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace geophase { 

    template<typename Geometry> 
    struct make_it_unit {};

    template<typename FT> 
    struct make_it_unit<ArrayPolygon<Vector<FT, 3>,3>>
    {
        static ArrayPolygon<Vector<FT, 3>, 3> apply()
        {
            return {
                point{FT{0},FT{0},FT{0}}, 
                point{FT{1},FT{0},FT{0}}, 
                point{FT{0},FT{1},FT{0}}
            };
        }
    };

    template<typename FT> 
    struct make_it_unit<ArrayPolygon<Vector<FT,3>,4>>
    {
        static ArrayPolygon<Vector<FT,3>, 4> apply()
        {
            return {
                point{FT{0},FT{0},FT{0}}, 
                point{FT{1},FT{0},FT{0}}, 
                point{FT{1},FT{1},FT{0}}, 
                point{FT{0},FT{1},FT{0}}
            };
        }
    };

    template<>
    struct make_it_unit<ArrayPolygon<Vector<double,3>,6>>
    {
        static ArrayPolygon<Vector<double,3>, 6> apply()
        {
            // Unit-radius hexagon.
            constexpr double deg30rad = (1. / 6.) * M_PI; 
            constexpr double sin30 = std::sin(deg30rad); 
            constexpr double cos30 = std::cos(deg30rad);

            return {
                point{cos30,   sin30,  0.},
                point{0.,      1.,     0.},
                point{-cos30,  sin30,  0.},
                point{-cos30, -sin30,  0.},
                point{0.,      -1.,    0.},
                point{cos30,  -sin30,  0.},
            };
        }
    };

    template<typename Geometry> 
    constexpr Geometry make_unit()
    {
        return make_it_unit<Geometry>::apply();
    }

    template<typename Polygon=vectorPolygon>
    Polygon make_unit_triangle()
    {
        using AT = typename Polygon::value_type::value_type;
        return {
                point{AT{0},AT{0},AT{0}}, 
                point{AT{1},AT{0},AT{0}}, 
                point{AT{0},AT{1},AT{0}}
        };
    }

    template<typename Polygon=vectorPolygon>
    Polygon make_unit_square()
    {
        using AT = typename Polygon::value_type::value_type;
        return {
                point{AT{0},AT{0},AT{0}}, 
                point{AT{1},AT{0},AT{0}}, 
                point{AT{1},AT{1},AT{0}}, 
                point{AT{0},AT{1},AT{0}}
        };
    }

    template<typename Polygon=vectorPolygon>
    Polygon make_unit_hexagon()
    {
        Polygon result(6); 
        auto hexa = make_unit<hexagon>();
        std::copy(hexa.begin(), hexa.end(), result.begin());
        return result;
    }

    //template <typename PointSequence, typename Halfspace>
    //typename std::enable_if
    //<
        //Traits::tag_enabled<Halfspace, halfspace_tag>::value &&
        //Traits::tag_enabled<PointSequence, point_sequence_tag>::value, 
        //PointSequence 
    //>::type 
    //build(Halfspace const& h) 
    //{
        //PointSequence points(4); 

        //const auto& position = h.position(); 
        //const auto& direction = h.direction(); 

        //// Compute a random point that is not lying on the plane.
        //std::random_device device;
        //std::default_random_engine eng(device());
        //// Random coordinate perturbation in [0, ||direction||^2] 
        //std::uniform_real_distribution<double> coord(-1, 1); 

        //// Create a point Q that doesn't lie on the plane.
        //auto Q = position;  

        //while(((Q - position) & direction) == 0)
        //{
            //Q[0] = coord(eng); 
            //Q[1] = coord(eng); 
            //Q[2] = coord(eng); 
        //}

        //// Using the point Q, find its projection R on the plane. 
        //// R = Q + s . direction , s is a real number. 
        //// (R - P) . n  = 0 -> 
        //auto s = ((position - Q) & direction) / (direction & direction); 
        //auto R = Q + (s * direction); 

        //// Compute the point T from R using the orthogonal vector resulting
        //// from the cross product of n and PR 
        //auto displacement = direction ^ (R - position);  
        //auto T = R + displacement;  

        //auto U = position + displacement; 

        //// Keep the point sequence orientation the same as the halfspace
        //// direction. TM.
        //points[0] = position; 
        //points[1] = R; 
        //points[2] = T; 
        //points[3] = U;

        //return points; 
    //}

    template <typename Box, typename Point>
    enable_if_types_tagged_t <Box, Box, axis_box_tag, Point, vector_tag>
    make(Point const& minPoint, Point const& maxPoint)
    {
        return Box(minPoint, maxPoint); 
    } 

    template<typename Box, typename PolygonSequence>
    std::enable_if_t
    <
        tag_base_of<PolygonSequence, sequence_polyhedron_tag>() && 
        type_tagged<Box, axis_box_tag>(), 
        Box
    >
    make(PolygonSequence const& polyhedron)
    {
        Box result; 

        auto extremePointPair = extreme_points(polyhedron); 
        assign_less(result.minPoint(), extremePointPair.first); 
        assign_greater(result.maxPoint(), extremePointPair.second);
        return result; 
    }

    // This constructs a 3D bounding box from a polyhedron. 
    // Dimension checks are not done.
    // This is only used for unit testing. TM.
    template <typename Polyhedron, typename AxisAlignedBoundingBox>
    std::enable_if_t 
    <
        tag_base_of<Polyhedron, sequence_polyhedron_tag>() &&
        type_tagged<AxisAlignedBoundingBox, axis_box_tag>(),
        Polyhedron
    > 
    make(AxisAlignedBoundingBox const & box)
    {
        // Get min point.
        auto min = box.minPoint(); 
        // Get max point.
        auto max = box.maxPoint(); 
        // Get the point difference. 
        auto min_max = max - min;  

        using PointType = typename Polyhedron::value_type::value_type;
        using Float = typename PointType::value_type;

        // Initialize the 8 box points.
        std::vector<PointType> points(8);

        // Build a box result with outward pointing normals. 
        PointType x (Float{1},Float{0},Float{0}); 
        PointType y (Float{0},Float{1},Float{0});
        PointType z (Float{0},Float{0},Float{1}); 

        x *= min_max[0]; 
        y *= min_max[1]; 
        z *= min_max[2]; 

        points[0] = min; 
        points[1] = min + y;  
        points[2] = points[1] + x;  
        points[3] = points[2] - y;

        points[4] = points[0] + z; 
        points[5] = points[1] + z; 
        points[6] = points[2] + z; 
        points[7] = points[3] + z; 

        // Build Polygons.
        typedef typename Polyhedron::value_type PolygonType; 

        PolygonType bottom = {points[0], points[1], points[2], points[3]};
        PolygonType top    = {points[7], points[6], points[5], points[4]}; 
        PolygonType left   = {points[0], points[4], points[5], points[1]}; 
        PolygonType right  = {points[3], points[2], points[6], points[7]}; 
        PolygonType front  = {points[0], points[3], points[7], points[4]}; 
        PolygonType back   = {points[1], points[5], points[6], points[2]}; 

        Polyhedron result = {top, bottom, left, right, front, back};

        return result; 
    }

    template <typename Polyhedron, typename Point>
    std::enable_if_t 
    <
        tag_base_of<Polyhedron, sequence_polyhedron_tag>() &&
        type_tagged<Point, vector_tag>(), 
        Polyhedron 
    >
    make(Point const& p1, Point const& p2)
    {
        return make<Polyhedron>(aabbox(p1, p2));
    }

    template<typename Polyhedron=vectorPolyhedron> 
    Polyhedron make_tetrahedron(
        typename Polyhedron::value_type::value_type const& p0, 
        typename Polyhedron::value_type::value_type const& p1, 
        typename Polyhedron::value_type::value_type const& p2, 
        typename Polyhedron::value_type::value_type const& p3
    )
    {
        using Polygon = typename Polyhedron::value_type;
        using Point = typename Polygon::value_type;
        Polyhedron result(4, Polygon(3, Point{})); 

        auto v1 = p1 - p0; 
        auto v2 = p2 - p0; 
        auto v3 = p3 - p0; 

        Polygon tri0 = {}, tri1 = {}, tri2 = {}, tri3 = {}; 

        if (dot(v3, cross(v1,v2)) < 0)
        {
            tri0 = Polygon{p0, p1, p2}; 
            tri1 = Polygon{p0, p3, p1}; 
            tri2 = Polygon{p0, p2, p3}; 
            tri3 = Polygon{p1, p3, p2}; 
        }
        else
        {
            tri0 = Polygon{p0, p2, p1}; 
            tri1 = Polygon{p0, p1, p3}; 
            tri2 = Polygon{p0, p3, p2}; 
            tri3 = Polygon{p1, p2, p3}; 
        }

        result[0] = tri0; 
        result[1] = tri1; 
        result[2] = tri2; 
        result[3] = tri3; 

        return result; 
    }

    vectorPolyhedron make_unit_cube()
    {
        return make<vectorPolyhedron>(point(0.,0.,0.), point(1.,1.,1.)); 
    }

    vectorPolyhedron make_unit_tetrahedron()
    {
        return make_tetrahedron(
            point(0.,0.,0.), 
            point(1.,0.,0.), 
            point(0.,1.,0.), 
            point(0.,0.,1.)
        );
    }

    std::vector<point> make_dodecahedron_points(
        const double hEdge,
        const double npl 
    )
    {
        // Coordinate values for additional points. 
        const double a = 1. + hEdge;  
        const double b = 1. - (hEdge*hEdge);

        vector nplVec{npl, npl, npl}; 

        // https://en.wikipedia.org/wiki/Dodecahedron

        // Main points (+-1, +-1, +-1)
        return std::vector<point> {

            // Cube points.
            point(-1., -1., -1.) - nplVec, 
            point( 1., -1., -1.) + nplVec, 
            point( 1.,  1., -1.) - nplVec, 
            point(-1.,  1., -1.) + nplVec, 

            point(-1., -1., 1.)  - nplVec, 
            point( 1., -1., 1.)  + nplVec,
            point( 1.,  1., 1.)  - nplVec, 
            point(-1.,  1., 1.)  + nplVec, 

            // Additional edges x = 0: 
            // (0, +-(1 + h), +-(1 − h*h))
            point(0., -a, -b)    - nplVec,
            point(0., +a, -b)    + nplVec,
            point(0., +a, +b)    - nplVec,
            point(0., -a, +b)    + nplVec,

            // Additional edges y = 0: 
            // (±(1 − h*h), 0, ±(1 + h)) 
            point(-b, 0., -a)    - nplVec,
            point(b,  0., -a)    + nplVec,
            point(b,  0.,  a)    - nplVec,
            point(-b, 0.,  a)    + nplVec,
            
            // Additional edges 
            // z = 0: (+-(1 + h), +-(1 − h*2), 0)
            point(-a, -b, 0.)    - nplVec,
            point(+a, -b, 0.)    + nplVec,
            point(+a, +b, 0.)    - nplVec,
            point(-a, +b, 0.)    + nplVec
        };
    }

    vectorPolyhedron make_dodecahedron(
        const double hEdge, 
        const double npl = 0.
    )
    {
        auto points = make_dodecahedron_points(hEdge, npl); 

        // Dodecahedron addressing.
        return vectorPolyhedron{
            {points[15],points[14],points[6],points[10],points[7]},   // Top
            {points[14],points[15],points[4],points[11],points[5]},   // Top
            {points[13],points[12],points[3],points[9],points[2]},    // Bottom
            {points[12],points[13],points[1],points[8],points[0]},    // Bottom 
            {points[8],points[11],points[4],points[16],points[0]},    // Front 
            {points[11],points[8],points[1],points[17],points[5]},    // Front 
            {points[9],points[3],points[19],points[7],points[10]},    // Back 
            {points[9],points[10],points[6],points[18],points[2]},    // Back
            {points[17],points[18],points[6],points[14],points[5]},   // Right 
            {points[18],points[17],points[1],points[13],points[2]},   // Right 
            {points[19],points[16],points[4],points[15],points[7]},   // Left 
            {points[16],points[19],points[3],points[12],points[0]}    // Left 
        };
    }

    vectorPolyhedron make_unit_dodecahedron(const double npl = 0.)
    {
        return make_dodecahedron(0.5 * (-1. + sqrt(5.)), npl);
    }

    vectorPolyhedron make_unit_endo_dodecahedron(const double npl = 0., const double h = -0.25)
    {
        return make_dodecahedron(h, npl);
    }

    //template 
    //<
        //typename MixedPolyhedron, 
        //typename PolygonSequence, 
        //typename Point, 
        //typename IT
    //>
    //typename std::enable_if
    //<
        //Traits::tag_enabled<MixedPolyhedron, mixed_polyhedron_tag>::value &&
        //Traits::tag_enabled<PolygonSequence, polygon_sequence_tag>::value &&
        //Traits::tag_enabled<Point, point_tag>::value, 
        //MixedPolyhedron 
    //>::type 
    //build(PolygonSequence const& polyhedron, Point const& point, IT polygonLabel)
    //{
        //MixedPolyhedron result; 

        //IT count = 0; 

        //for (const auto& polygon : polyhedron)
        //{
            //if (count != polygonLabel)
            //{
                //result.push_back(polygon); 
            //}
            //else
            //{
                //auto curIt = polygon.begin(); 
                //auto nextIt = std::next(curIt); 
                //while(nextIt != polygon.end())
                //{
                    //result.push_back(
                        //build<arrayTriangle>(*curIt, *nextIt, point)
                    //);
                    //++curIt; 
                    //++nextIt; 
                //}
                //nextIt = polygon.begin(); 
                //result.push_back(
                    //build<arrayTriangle>(*curIt, *nextIt, point)
                //);
            //}
            //++count;
        //}

        //return result;
    //}

    //template 
    //<
        //typename MixedPolyhedron, 
        //typename PolygonSequence
    //>
    //typename std::enable_if
    //<
        //Traits::tag_enabled<MixedPolyhedron, mixed_polyhedron_tag>::value &&
        //Traits::tag_enabled<PolygonSequence, polygon_sequence_tag>::value, 
        //MixedPolyhedron 
    //>::type 
    //build(PolygonSequence const& polyhedron)
    //{
        //MixedPolyhedron result; 

        //for (const auto& polygon : polyhedron)
        //{
            //result.push_back(polygon); 
        //}

        //return result; 
    //}

    //template<typename PolygonSequence, typename T1, typename T2> 
    //typename std::enable_if
    //<
        //Traits::tag_enabled<PolygonSequence, polygon_sequence_tag>::value &&  
        //Traits::tag_enabled<T1, point_sequence_tag>::value &&  
        //Traits::tag_enabled<T2, point_collection_tag>::value,
        //PolygonSequence
    //>::type
    //build(T1 const & basePolygon, T2 const & displacements)
    //{
        //PolygonSequence resultPolyhedron (basePolygon.size() + 2); 

        //auto displacedPolygon = displace(basePolygon, displacements);  
        //std::reverse(displacedPolygon.begin(), displacedPolygon.end()); 

        //resultPolyhedron.front() = basePolygon; 
        //resultPolyhedron.back() = displacedPolygon; 

        //auto sideIter = std::next(resultPolyhedron.begin());  
        //const auto sideEndIter = std::prev(resultPolyhedron.end(),2);  
        //auto baseIter = basePolygon.begin(); 
        //auto sweptIter = std::prev(displacedPolygon.end()); 
        //typedef typename PolygonSequence::value_type PolygonType;

        //for(; sideIter != sideEndIter; ++sideIter)
        //{
            //*sideIter = PolygonType{
                //*std::next(baseIter), 
                //*baseIter, 
                //*sweptIter,
                //*std::prev(sweptIter)
            //};
            //--sweptIter; 
            //++baseIter;
        //}

        //*sideIter = PolygonType{
            //basePolygon.front(),  
            //basePolygon.back(),  
            //displacedPolygon.front(), 
            //displacedPolygon.back()
        //};

        //return resultPolyhedron;
    //}

    //template<typename PolygonSequence, typename PointCollection, typename Vector> 
    //typename std::enable_if
    //<
        //Traits::tag_enabled<PolygonSequence, polygon_sequence_tag>::value &&
        //Traits::tag_enabled<PointCollection, point_collection_tag>::value &&
        //Traits::tag_enabled<Vector, vector_tag>::value,
        //PolygonSequence
    //>::type
    //build(PointCollection const & polygon, Vector const & displacement)
    //{
        //return build<PolygonSequence>(
                //polygon, 
                //PointCollection(polygon.size(), displacement)
        //); 
    //} 

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


} // End namespace geophase 

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //
