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

    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
    Triangulation algorithms.  

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

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

#ifndef Triangulate_HPP
#define Triangulate_HPP

#include "tag.hpp"
#include "Polyhedron.hpp"
#include "Centroid.hpp"
#include "Make.hpp"
#include <numeric>

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

namespace geophase {

    // reserve for std::vector-based centroid triangulations.

    template<typename Polyhedron, typename Polygon>
    void point_reserve(Polyhedron&, Polygon const&) {}

    template
    <
        typename Point,    
        template<typename, typename> class Container,
        template <typename T> class Allocator 
    >
    void point_reserve(
        vectorPolyhedron& triangulation, 
        ContainerPolygon<Point, Container, Allocator>const& polygon 
    ) 
    {
        triangulation.reserve(polygon.size());
    }

    template
    <
        typename Polygon,    
        template<typename, typename> class Container,
        template <typename T> class Allocator 
    >
    void point_reserve(
        vectorPolyhedron& triangulation, 
        ContainerPolyhedron<Polygon, Container, Allocator>const& polyhedron 
    ) 
    {
        using AT = typename vectorPolyhedron::value_type::value_type::value_type; 
        AT sumSizes {0};

        for (const auto& polygon : polyhedron)
            sumSizes += polygon.size(); 

        triangulation.reserve(sumSizes);
    }

    // Reserving functions for std::vector-based oriented triangulations.
    
    template<typename Polyhedron, typename Polygon>
    void oriented_reserve(Polyhedron&, Polygon const&) {}

    template
    <
        typename Point,    
        template<typename, typename> class Container,
        template <typename T> class Allocator 
    >
    void oriented_reserve(
        vectorPolyhedron& triangulation, 
        ContainerPolygon<Point, Container, Allocator>const& polygon 
    ) 
    {
        triangulation.reserve(polygon.size() - 2);
    }

    template
    <
        typename Polygon,    
        template<typename, typename> class Container,
        template <typename T> class Allocator 
    >
    void oriented_reserve(
        vectorPolyhedron& triangulation, 
        ContainerPolyhedron<Polygon, Container, Allocator>const& polyhedron 
    ) 
    {
        using AT = typename vectorPolyhedron::value_type::value_type::value_type; 
        AT sumSizes {0};

        for (const auto& polygon : polyhedron)
            sumSizes += polygon.size() - 2; 

        triangulation.reserve(sumSizes);
    }


    template<typename PolygonSequence=vectorPolyhedron, typename Polygon>
    std::enable_if_t<tag_base_of<Polygon, sequence_polygon_tag>(), PolygonSequence> 
    point_triangulate(
        Polygon const& polygon, 
        typename Polygon::value_type const& triPoint 
    )
    {
        PolygonSequence result;
        point_reserve(result, polygon.size());
        
        auto curr = polygon.begin();  
        auto next = std::next(polygon.begin());
        for (; next != polygon.end(); ++curr, ++next)
            result.push_back({triPoint, *curr, *next});

        result.push_back({triPoint, polygon.back(), polygon.front()});

        return result; 
    }

    template<typename PolygonSequence=vectorPolyhedron, typename Polygon>
    std::enable_if_t<tag_base_of<Polygon, sequence_polygon_tag>(), PolygonSequence> 
    centroid_triangulate(Polygon const& polygon)
    {
        return point_triangulate<PolygonSequence, Polygon>(polygon, centroid(polygon)); 
    }

    template<typename PolygonSequence=vectorPolyhedron, typename Polygon>
    std::enable_if_t<tag_base_of<Polygon, sequence_polygon_tag>(), PolygonSequence> 
    oriented_triangulate(Polygon const& polygon)
    {
        PolygonSequence result; 
        point_reserve(result, polygon);  

        const auto begin = polygon.begin();
        auto first = std::next(begin); 
        auto second = std::next(first); 
        for (; second != polygon.end() ; ++first, ++second)
            result.push_back({*begin, *first, *second});

        return result;
    }

    template<typename PolygonSequence=vectorPolyhedron>
    std::enable_if_t<tag_base_of<PolygonSequence, sequence_polyhedron_tag>(), PolygonSequence>
    centroid_triangulate(PolygonSequence const& polyhedron)
    {
        PolygonSequence result; 
        point_reserve(result, polyhedron); 

        for (const auto& polygon : polyhedron)
        {
            if (polygon.size() > 3)
            {
                auto polygonTriangulation = centroid_triangulate(polygon);  
                for (const auto& triangle : polygonTriangulation)
                    result.push_back(triangle); 
            }
            else if (polygon.size() == 3)
                result.push_back(polygon);
        }
        return result;
    }

    template<typename PolygonSequence=vectorPolyhedron>
    std::enable_if_t<tag_base_of<PolygonSequence, sequence_polyhedron_tag>(), PolygonSequence>
    oriented_triangulate(PolygonSequence const& polyhedron)
    {
        PolygonSequence result; 
        oriented_reserve(result, polyhedron); 

        for (const auto& polygon : polyhedron)
        {
            if (polygon.size() > 3)
            {
                auto polygonTriangulation = oriented_triangulate(polygon);  
                for (const auto& triangle : polygonTriangulation)
                    result.push_back(triangle); 
            }
            else if (polygon.size() == 3)
                result.push_back(polygon);
        }
        return result;
    }

    template<typename PolygonSequence>
    std::enable_if_t<tag_base_of<PolygonSequence, sequence_polyhedron_tag>(), std::vector<PolygonSequence>>
    centroid_vol_triangulate(PolygonSequence const& polyhedron)
    {
        std::vector<PolygonSequence> result;

        if (polyhedron.size() == 4) // It is a tetrahedron.
        {
            result.push_back(polyhedron);
            return result;
        }

        decltype(polyhedron.size()) resultSize = 0; 
        for (const auto& polygon : polyhedron)
            resultSize += polygon.size(); 

        result.reserve(resultSize); 


        auto polyhedronCentroid = centroid(polyhedron); 
        auto surfaceTriangulation = centroid_triangulate(polyhedron);  
        for (const auto& triangle : surfaceTriangulation)
        {
            result.push_back(make_tetrahedron(
                polyhedronCentroid, 
                triangle[0], 
                triangle[1], 
                triangle[2]
                )
            );
        }
        
        return result;
    }

    template<typename PolygonSequence>
    std::enable_if_t<tag_base_of<PolygonSequence, sequence_polyhedron_tag>(), std::vector<aabbox>>
    centroid_vol_boxes(PolygonSequence const& polyhedron)
    {
        std::vector<aabbox> result;

        if (polyhedron.size() == 4) // It is a tetrahedron.
        {
            result.push_back(make<aabbox>(polyhedron));
            return result;
        }

        decltype(polyhedron.size()) resultSize = 0; 
        for (const auto& polygon : polyhedron)
            resultSize += polygon.size(); 

        result.reserve(resultSize); 

        auto polyhedronCentroid = centroid(polyhedron); 
        auto surfaceTriangulation = centroid_triangulate(polyhedron);  
        for (const auto& triangle : surfaceTriangulation)
        {
            aabbox box; 
            box.pointExpand(polyhedronCentroid); 
            box.pointExpand(triangle[0]);
            box.pointExpand(triangle[1]);
            box.pointExpand(triangle[2]);
            result.push_back(box);
        }
        
        return result;
    }
    
    
    /*
    template<typenaem Triangulation=std::vector<vectorPolyhedron>, typename PolygonSequence=vectorPolyhedron>
    std::enable_if_t<tag_base_of<PolygonSequence, sequence_polyhedron_tag>(), Triangulation>
    oriented_vol_triangulate(PolygonSequence const& polyhedron)
    {
        // Volume triangulation is a container of polyhedrons, which one? 
         
        auto surfaceTriangulation = oriented_triangulate(polyhedron);  

        // Compute the centroid of a polyhedron. Centroid.H

        // Create tets from the centroid and the surface triangulation + push them into the result. 
        
        return result;
    }
    */

    //template
    //<
        //typename Triangulation = tetrahedronVector, 
        //typename Box 
    //>
    //typename std::enable_if
    //<
        //Traits::tag_enabled<Triangulation, tetrahedron_sequence_direct_tag>::value &&
        //Traits::tag_enabled<Box, aabbox_tag>::value, 
        //Triangulation
    //>::type 
    //oriented_triangulate(Box const & box)
    //{
        //typedef PointSequence<point, std::vector> polygon;
        //typedef PolygonSequence<polygon, std::vector> polyhedron;

        //auto poly = build<polyhedron, Box>(box);  

        //return oriented_triangulate<Triangulation, polyhedron>(poly);
    //}

    //template
    //<
        //typename Triangulation = tetrahedronVector, 
        //typename Box 
    //>
    //typename std::enable_if
    //<
        //Traits::tag_enabled<Triangulation, tetrahedron_sequence_direct_tag >::value && 
        //Traits::tag_enabled<Box, aabbox_tag>::value,
        //Triangulation
    //>::type 
    //barycentric_triangulate(Box const & box)
    //{
        //typedef PointSequence<point, std::vector> polygon;
        //typedef PolygonSequence<polygon, std::vector> polyhedron;

        //auto poly = build<polyhedron, Box>(box);  

        //return barycentric_triangulate<Triangulation, polyhedron>(poly);
    //}

    //template
    //<
        //typename Triangulation = tetrahedronVector, 
        //typename PolygonSequence
    //>
    //typename std::enable_if
    //<
        //Traits::tag_enabled<Triangulation, tetrahedron_sequence_tag>::value &&
        //Traits::tag_enabled<PolygonSequence, polygon_collection_tag>::value,
        //Triangulation
    //>::type 
    //barycentric_triangulate(PolygonSequence const & polyhedron)
    //{
        //Triangulation result;
        //result.reserve(32); 

        //typedef typename Triangulation::value_type tetrahedron; 
        //typedef typename tetrahedron::value_type triangle; 
        //typedef TriangleSequence<triangle, std::vector> triangulation;

        //auto centerPoint = centroid(polyhedron); 

        //for (const auto & polygon : polyhedron)
        //{
            //auto polygonTriangulation = barycentric_triangulate<triangulation>(polygon); 

            //for (const auto & triangle : polygonTriangulation)
            //{
                //result.push_back(
                    //build<tetrahedron>(
                        //centerPoint, 
                        //triangle[0], 
                        //triangle[1],
                        //triangle[2]
                    //)
                //); 
            //}
        //}

        //return result; 
    //}

    //template
    //<
        //typename Triangulation = tetrahedronVector, 
        //typename MixedPolygonSequence 
    //>
    //typename std::enable_if
    //<
        //Traits::tag_enabled<Triangulation, tetrahedron_sequence_tag>::value &&
        //Traits::tag_enabled<MixedPolygonSequence, mixed_polyhedron_tag>::value,
        //Triangulation
    //>::type 
    //barycentric_triangulate(MixedPolygonSequence const & mixedPolygonSequence)
    //{
        //Triangulation result; 
        
        //typedef typename Triangulation::value_type tetrahedron; 
        //typedef typename tetrahedron::value_type triangle; 
        //typedef TriangleSequence<triangle, std::vector> triangulation;

        //auto centerPoint = centroid(mixedPolygonSequence); 

        //for (const auto & polygon : mixedPolygonSequence.polygons())
        //{
            //auto polygonTriangulation = barycentric_triangulate<triangulation>(polygon); 

            //for (const auto & triangle : polygonTriangulation)
            //{
                //result.push_back(
                    //build<tetrahedron>(
                        //centerPoint, 
                        //triangle[0], 
                        //triangle[1],
                        //triangle[2]
                    //)
                //); 
            //}
        //}

        //for (const auto & triangle : mixedPolygonSequence.triangles())
        //{
            //result.push_back(
                //build<tetrahedron>(
                    //centerPoint, 
                    //triangle[0], 
                    //triangle[1],
                    //triangle[2]
                //)
            //); 
        //}

        //return result; 
    //}

    //template <typename Triangulation> 
    //struct barycentric_triangulation
    //{
        //template<typename Geometry>
        //static Triangulation apply(Geometry const & geo)
        //{
            //return barycentric_triangulate<Triangulation>(geo);
        //}

        //template<typename Geometry>
        //Triangulation operator()(Geometry const & geo) const
        //{
            //return apply(geo);  
        //}
    //};

    //template <typename Triangulation> 
    //struct oriented_triangulation
    //{
        //template <typename Geometry> 
        //static Triangulation apply(Geometry const & geo)
        //{
            //return oriented_triangulate<Triangulation>(geo);
        //}

        //template <typename Geometry> 
        //Triangulation operator()(Geometry const & geo) const
        //{
            //return apply(geo);
        //}
    //};
    //


} // End namespace geophase 

#endif

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