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

    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/>.

Application
    geophaseTestOrient 

Description
    Orient a set of points using a prescribed orientation vector. 

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

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

#include "Polygon.hpp"
#include "Orient.hpp"

#include <gtest/gtest.h> 
#include "testUtilities.hpp"
#include <random>
#include "Make.hpp"
#include "Area.hpp"
#include "WriteVtkPolyData.hpp"

using namespace geophase;

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:

template<typename PointIterator>
bool pointChainIsNotSelfIntersecting(PointIterator begin, PointIterator end)
{
    // Check if points are chained in a non-self intersecting chain. 
    // Points are not chained if there is an edge (A, B) in the point set S 
    // for which there exist another edge (P, Q) in S \ (A, B) such that 
    // ((B - A) x (P - A)) . ((B - A) x (Q - A)) < 0 , TM

    auto previousPoint = begin; 
    auto nextPoint = std::next(begin); 

    // For all polygon edges
    while(nextPoint != end)
    {
        auto edge = *nextPoint - *previousPoint; 
        auto innerPreviousPoint = std::next(previousPoint); 
        auto innerNextPoint = std::next(nextPoint); 

        while(innerNextPoint != end)
        {
            auto innerPreviousOrientation = cross(
                edge,
                *innerPreviousPoint - *previousPoint
            ); 
            auto innerNextOrientation = cross(
                edge, 
                *innerPreviousPoint - *previousPoint
            ); 

            // If there is such an inner loop edge whose points lie on opposite sides
            // of the outer loop edge, the chain is self-intersecting.
            if (dot(innerPreviousOrientation,innerNextOrientation) < 0) return false; 

            ++innerPreviousPoint; 
            ++innerNextPoint;
        }
        ++previousPoint;
        ++nextPoint;
    }

    return true;
}

template<typename PointIterator, typename Vector>
bool pointChainIsOriented(PointIterator begin, PointIterator end, Vector normal)
{
    auto current = begin; 
    auto next = std::next(current); 
    auto nextNext = std::next(next);

    while(nextNext != end)
    {
        auto firstEdge = *next - *current; 
        auto nextEdge = *nextNext - *next;

        if (dot(cross(firstEdge,nextEdge),normal) < 0) return false;

        ++current; 
        next = std::next(current); 
        nextNext = std::next(next);
    }

    return true; 
}

template<typename PointSequence>
PointSequence generateRandomCirclePoints(unsigned int size)
{
    PointSequence circlePoints(size, point(0.,0.,0.)); 

    std::random_device rd;
    std::default_random_engine e1(rd());
    std::uniform_real_distribution<double> angle_dist(0,3.14159265359); 
    std::bernoulli_distribution sign(0.5);

    for (decltype(size) i = 0; i < size; ++i)
    {
        double angle = angle_dist(e1);
        double x = std::cos(angle); 
        double y = std::sin(angle); 

        if (sign(e1)) {
            x *= -1; 
            y *= -1; 
        }

        circlePoints[i] = point(x, y, 0.); 
    }

    return circlePoints; 
}

template<typename PointSequence>
PointSequence generateRegularOctogon()
{
    PointSequence result = {
        point(1,0,0), 
        point(2,0,0), 
        point(3,1,0), 
        point(3,2,0), 
        point(2,3,0), 
        point(1,3,0), 
        point(0,2,0), 
        point(0,1,0) 
    };

    return result;
}

std::string prependName(std::string name, unsigned int n)
{
    std::stringstream ss; 
    ss << name << "-" << std::setw(4) << std::setfill('0') << n << ".vtk"; 
    return std::string(ss.str());  
}

TEST(orientPoints, trianglePolygon)
{
    auto triangle = make_unit_triangle(); 
    vector reversedNormal = -1.*area_normal_vector(triangle);
    vectorPolygon triangleReversed = orient(triangle, reversedNormal); 

    ASSERT_TRUE(pointChainIsNotSelfIntersecting(triangleReversed.begin(),triangleReversed.end()));
    ASSERT_TRUE(pointChainIsOriented(triangleReversed.begin(),triangleReversed.end(), reversedNormal));
}

TEST(orientPoints, square) 
{
    auto square = make_unit_square<vectorPolygon>();
    vector reversedNormal = -1.*area_normal_vector(square);  
    auto squareReversed = orient(square, reversedNormal); 

    ASSERT_TRUE(pointChainIsNotSelfIntersecting(squareReversed.begin(), squareReversed.end())); 
    ASSERT_TRUE(pointChainIsOriented(squareReversed.begin(), squareReversed.end(), reversedNormal)); 
}

TEST(orientPoints, hexagon)
{
    auto hex = make_unit_hexagon(); 
    vector normal = area_normal_vector(hex); 

    ASSERT_TRUE(pointChainIsNotSelfIntersecting(hex.begin(), hex.end()));
    ASSERT_TRUE(pointChainIsOriented(hex.begin(), hex.end(), normal));

    vector reversedNormal = -1. * normal;
    auto hexReversed = orient(hex, reversedNormal);

    ASSERT_TRUE(pointChainIsNotSelfIntersecting(hexReversed.begin(), hexReversed.end()));
    ASSERT_TRUE(pointChainIsOriented(hexReversed.begin(), hexReversed.end(), reversedNormal));
}

TEST(orientPoints, squareDuplicatePoints)
{
    // Square with repeated points.
    vectorPolygon square = {
        point(0.,0.,0.), 
        point(0.,0.,0.), 
        point(1.,0.,0.), 
        point(1.,0.5,0.), 
        point(1.,0.,0.), 
        point(0.,0.,0.), 
        point(1.,0.25,0.), 
        point(0.,0.,0.), 
        point(1.,0.5,0.), 
        point(1.,0.75,0.), 
        point(1.,1.,0.), 
        point(0.,0.,0.), 
        point(0.,1.,0.),
        point(0.,0.,0.) 
    };

    auto squareNormal = area_normal_vector(square); 
    auto squareOriented = orient(square, squareNormal); 

    ASSERT_TRUE(pointChainIsNotSelfIntersecting(squareOriented.begin(), squareOriented.end()));
    ASSERT_TRUE(pointChainIsOriented(squareOriented.begin(), squareOriented.end(), squareNormal));
}

TEST(orientPoints, heisenBug)
{
    // Encountered orientation heisenbug with random point generation. 
    // Captured input data. 
    vectorPolygon poly = {
       point(0.0493244 ,-0.998783   ,0.),
       point(-0.990411  ,-0.138152  ,0.), 
       point(0.922827  ,0.385214    ,0.), 
       point(-0.976731  ,0.21447    ,0.), 
       point(-0.693507  ,-0.720449  ,0.),
       point(-0.341261  ,-0.939969  ,0.),
       point(0.997168   ,-0.0752104 ,0.),
       point(0.543243   ,0.839575   ,0.),
       point(-0.0643775 ,-0.997926  ,0.),
       point(0.979163   ,-0.203075  ,0.),
       point(-0.24903   ,0.968496   ,0.),
       point(0.996597   ,0.082433   ,0.),
       point(0.904713   ,-0.426022  ,0.),
       point(0.985459   ,-0.169916  ,0.),
       point(-0.320907  ,0.947111   ,0.),
       point(0.440235   ,0.897883   ,0.)
    };
    vector normal(0.,0.,-1.); 
    auto orientedPoly = orient(poly, normal); 

    ASSERT_TRUE(orientedPoly.size() == poly.size()) 
            << "orientedPoly.size() = " << orientedPoly.size() << std::endl
            << "poly.size() = " << poly.size() << std::endl;
    ASSERT_TRUE(pointChainIsNotSelfIntersecting(orientedPoly.begin(), orientedPoly.end()));
    ASSERT_TRUE(pointChainIsOriented(orientedPoly.begin(), orientedPoly.end(), normal));
}

TEST(orientPoints, randomCirclePointsOrientedByProjection)
{
    for (unsigned int nPoints = 3;  nPoints < 100; ++nPoints)
    {
        vectorPolygon circlePolygon = 
            generateRandomCirclePoints<vectorPolygon>(nPoints);  

        write_vtk(circlePolygon, prependName("circlePolygon", nPoints)); 

        vector normal(0.,0.,-1.); 

        auto orientedCirclePolygon = orient(circlePolygon, normal); 

        write_vtk(
            orientedCirclePolygon, 
            prependName("orientedCirclePolygon", nPoints)
        );

        ASSERT_TRUE(orientedCirclePolygon.size() == nPoints) 
            << "orientedCirclePolygon.size() = " << orientedCirclePolygon.size() << std::endl
            << "nPoints = " << nPoints << std::endl;
        ASSERT_TRUE(pointChainIsNotSelfIntersecting(orientedCirclePolygon.begin(), orientedCirclePolygon.end()))
            << "orientedCirclePolygon.size() = " << orientedCirclePolygon.size() << std::endl
            << "nPoints = " << nPoints << std::endl;
        ASSERT_TRUE(pointChainIsOriented(orientedCirclePolygon.begin(), orientedCirclePolygon.end(), normal)) 
            << "orientedCirclePolygon.size() = " << orientedCirclePolygon.size() << std::endl
            << "nPoints = " << nPoints << std::endl;
    }
}

template<typename PointSequence>  
void duplicate_random_polygon_point(PointSequence& polygon)
{
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_int_distribution<> dis(0, polygon.size()-1);

    auto duplicatePointIt = polygon.begin(); 

    int randomPosition = dis(gen); 

    duplicatePointIt = std::next(duplicatePointIt,randomPosition); 

    polygon.push_back(*duplicatePointIt); 
}

TEST(orientPoints, samePoints)
{
    vectorPolygon points = {point(0.,0.,0.),point(0.,0.,0.),point(0.,0.,0.)};

    vector normal (0.,1.,0.); 

    auto orientedPoints = orient(points, normal, EPSILON); 

    ASSERT_TRUE(orientedPoints.size() == 0)
        << "number of points = " << orientedPoints.size() << "\n";
}

int main(int argc, char *argv[])
{
    ::testing::InitGoogleTest(&argc, argv);
    return RUN_ALL_TESTS();

    return 0;
}


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