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

    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
    geophaseTestInside

Description
    Unit tests for the is_inside algorithm.

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

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

#include "gtest/gtest.h" 
#include "Perturb.hpp"
#include "Vector.hpp"
#include "Halfspace.hpp"
#include "IsInside.hpp"
#include "Make.hpp"

// Shewchuck's predicates.
#include "predicates.c"

// PNG visualization of the inside test. 
//#include <png++/png.hpp> // TODO: png++ as a CMake dependency

// Simple timing 
#include <chrono>
#include <fstream>

using namespace geophase;

//decltype(auto) color(position pos)
//{
    //if (pos == position::INSIDE)
        //// Red 
        //return png::rgb_pixel(255, 0, 0); 
    //if (pos == position::COINCIDENT)
        //// Yellow
        //return png::rgb_pixel(255, 255, 0); 
    //if (pos == position::OUTSIDE)
        //// Blue 
        //return png::rgb_pixel(0, 0, 255); 
    
    //// Black default.
    //return png::rgb_pixel(0, 0, 0); 
//}

// Function object structs used to simplify testing. 
struct orient2Dtest
{
    position operator()(
        const vector& pPoint, 
        const vector& qPoint, 
        const vector& rPoint
    ) const
    {
        return orient2D(pPoint, qPoint, rPoint); 
    }
};

struct isInsideFloatTest
{
    position operator()(
        const vector& pPoint, 
        const vector& qPoint, 
        const vector& rPoint
    ) const
    {
        halfspace hspace(
            pPoint, 
            vector(
                -qPoint[1] + pPoint[1], 
                qPoint[0] - pPoint[0], 
                rPoint[2]
            )
        );
        //hspace.direction() /= mag(hspace.direction()); 
        return is_inside_float(rPoint, hspace);
    }
};

struct isInsideEpsilonTest
{
    double epsilon_; 

    isInsideEpsilonTest(double epsilon)
        :
            epsilon_(epsilon)
    {}

    position operator()(
        const vector& pPoint, 
        const vector& qPoint, 
        const vector& rPoint
    ) const
    {
        halfspace hspace(
            pPoint,
            vector(
                -qPoint[1] + pPoint[1], 
                qPoint[0] - pPoint[0], 
                rPoint[2]
            )
        );
        return is_inside_epsilon(rPoint, hspace, epsilon_);
    }
};

struct isInsideUlpTest
{
    double NULP_; 

    isInsideUlpTest(std::size_t NULP)
        :
            NULP_(NULP)
    {}

    position operator()(
        const vector& pPoint, 
        const vector& qPoint, 
        const vector& rPoint
    ) const
    {
        halfspace hspace(
            pPoint, 
            vector(
                -qPoint[1] + pPoint[1], 
                qPoint[0] - pPoint[0], 
                rPoint[2]
            )
        );
        return is_inside_ulp(rPoint, hspace, NULP_);
    }
};

struct isInsideErrorTest 
{
    position operator()(
        const vector& pPoint, 
        const vector& qPoint, 
        const vector& rPoint
    ) const
    {
        halfspace hspace(
            pPoint,
            vector(
                -qPoint[1] + pPoint[1], 
                qPoint[0] - pPoint[0], 
                rPoint[2]
            )
        );
        //hspace.orientation() /= mag(hspace.orientation());
        return is_inside_error(rPoint, hspace);
    }
};

//TEST(INSIDE, COLLINEAR_POINTS_BENCHMARK)
//{
    //// Collinearity test for three points from:
    
    //// Kettner, L., Mehlhorn, K., Pion, S., Schirra, S., & Yap, C. (2008).
    //// Classroom examples of robustness problems in geometric computations.
    //// Computational Geometry: Theory and Applications, 40(1), 61–78.
    //// https://doi.org/10.1016/j.comgeo.2007.06.003

    //struct shewchuckOrient2Dtest 
    //{
        //position operator()(
            //const vector& pPoint, 
            //const vector& qPoint, 
            //const vector& rPoint
        //) const
        //{
            //// Cast away constness for Shewchuck's predicates.
            //double* pPtr = const_cast<vector::iterator>(pPoint.begin());
            //double* qPtr = const_cast<vector::iterator>(qPoint.begin());
            //double* rPtr = const_cast<vector::iterator>(rPoint.begin());

            //return static_cast<position>(sign(orient2d(pPtr, qPtr, rPtr)));
        //}
    //};

    //const std::size_t NULP = 16; 
    //const double NEPSILON  = 16*EPSILON;
    //std::ofstream timingFile("inside.Kettner2007.timing.dat"); 
    //timingFile << "orient2D,float,epsilon,error,ulp,adaptive" << std::endl;
    //std::cout  << "orient2D,float,epsilon,error,ulp,adaptive" << std::endl;

    //auto test = [&NULP](
        //const vector& pPoint, 
        //const vector& qPoint, 
        //const vector& rPoint, 
        //const auto& insideTest, 
        //std::string imageName 
    //)
    //{
        //const std::size_t SIZE = 1024; 
        //png::image<png::rgb_pixel> image(SIZE, SIZE);
        //using namespace std::chrono; 
        //high_resolution_clock::time_point start = high_resolution_clock::now();
        //for (size_t nx = 0; nx < SIZE; ++nx)
        //{
            //for (size_t ny = 0; ny < SIZE; ++ny)
            //{
                //vector position (pPoint); 
                //position[0] += nx * ulp(position[0]);
                //position[1] += ny * ulp(position[1]);
                //if (nx != ny) 
                    //image[nx][ny] = color(insideTest(position, rPoint, qPoint));
                //if ((nx == ny) || (nx - ny - 1 == NULP) || (ny - nx - 1 == NULP))
                    //image[nx][ny] = png::rgb_pixel(0,0,0); 
            //}
        //}
        //auto duration = duration_cast<milliseconds>(high_resolution_clock::now() - start);
        //image.write(imageName);
        //return duration.count();
    //};

    //auto testInsideAlgorithms = [&test, &timingFile, &NULP, &NEPSILON](
        //const vector& pPoint, 
        //const vector& qPoint, 
        //const vector& rPoint, 
        //std::string baseName 
    //)
    //{
        //auto durationOrient2D = test(pPoint, qPoint, rPoint, orient2Dtest(), 
                                     //baseName + ".test1.orient2D.png"); 
        //auto durationFloat    = test(pPoint, qPoint, rPoint, isInsideFloatTest(), 
                                     //baseName + ".test2.insideFloat.png"); 
        //auto durationEpsilon  = test(pPoint, qPoint, rPoint, isInsideEpsilonTest(NEPSILON), 
                                     //baseName + ".test3.insideEpsilon.png"); 
        //auto durationError    = test(pPoint, qPoint, rPoint, isInsideErrorTest(), 
                                     //baseName + ".test4.insideError.png"); 
        //auto durationUlp      = test(pPoint, qPoint, rPoint, isInsideUlpTest(NULP), 
                                     //baseName + ".test5.insideUlp.png"); 
        //auto durationAdaptive = test(pPoint, qPoint, rPoint, shewchuckOrient2Dtest(), 
                                     //baseName + ".test6.insideAdaptive.png"); 
        //timingFile << durationOrient2D << "," << durationFloat << "," << durationEpsilon << ","
            //<< durationError << "," << durationUlp << "," << durationAdaptive << std::endl;
        //std::cout << durationOrient2D << "," << durationFloat << "," << durationEpsilon << ","
            //<< durationError << "," << durationUlp << "," << durationAdaptive << std::endl;
    //};


    //// Figure 2a, Kettner et al. 2007, doi:10.1016/j.comgeo.2007.06.003  
    //{
        //const vector pPoint  { 0.5,  0.5, 0.};  
        //const vector qPoint  {12.,  12.,  0.};
        //const vector rPoint  {24.,  24.,  0.}; 
        //testInsideAlgorithms(pPoint, qPoint, rPoint, "Kettner2007.Fig2a"); 
    //}

    //// Figure 2b, Kettner et al. 2007, doi:10.1016/j.comgeo.2007.06.003  
    //{
        //const vector pPoint  { 0.50000000000002531,  0.5000000000000171,   0.};
        //const vector qPoint  {17.300000000000001,   17.300000000000001,     0.};
        //const vector rPoint  {24.00000000000005,    24.0000000000000517765, 0.};
        //testInsideAlgorithms(pPoint, qPoint, rPoint, "Kettner2007.Fig2b"); 
    //}

    //// Figure 2c, Kettner et al. 2007, doi:10.1016/j.comgeo.2007.06.003  
    //{
        //const vector pPoint  { 0.5,                0.5,                    0.};
        //const vector qPoint  { 8.8000000000000007, 8.8000000000000007,     0.};
        //const vector rPoint  {12.1,               12.1,                    0.};
        //testInsideAlgorithms(pPoint, qPoint, rPoint, "Kettner2007.Fig2c"); 
    //}

    //// Overlapping points present in VoF calculations. 
    //// Perturbed (0.,0.,0.) - overlapping with (0.,0.,0.). Normal given by h
    //// (h, h, 0.), h = 2^-8 
    //{
        //const vector pPoint  { 0.,                0.,                    0.};
        //const vector qPoint  { 0.00390625,        0.00390625,            0.};
        //const vector rPoint  { 0.,                0.,                    0.};
        //testInsideAlgorithms(pPoint, qPoint, rPoint, "VoFoverlap.0.2-08.0"); 
    //}
    //// Overlapping points present in VoF calculations. 
    //// Perturbed (0.,0.,0.) - overlapping with (h,h,h). 
    //// h is the discretization length, h = 2^-8
    //{
        //const vector pPoint  { 0.,                0.,                    0.};
        //const vector qPoint  { 0.00390625,        0.00390625,            0.};
        //const vector rPoint  { 0.00390625,        0.00390625,            0.};
        //testInsideAlgorithms(pPoint, qPoint, rPoint, "VoFoverlap.0.2-08.2-08"); 
    //}
//}

TEST(IS_INSIDE_ERROR_UNIT, POINT_HALFSPACE)
{
    point pt {0., 0., 0.}; 

    halfspace hsaboveup (vector{0.,0., 1.}, vector{0.,0.,1.});
    halfspace hsabovedown(hsaboveup);
    hsabovedown.flip(); 

    halfspace hsbelowdown (vector{0.,0., -1.}, vector{0.,0.,-1.});
    halfspace hsbelowup(hsbelowdown); 
    hsbelowup.flip();

    ASSERT_TRUE(is_inside_error(pt, hsaboveup) == position::OUTSIDE);
    ASSERT_TRUE(is_inside_error(pt, hsabovedown) == position::INSIDE);

    ASSERT_TRUE(is_inside_error(pt, hsbelowdown) == position::OUTSIDE);
    ASSERT_TRUE(is_inside_error(pt, hsbelowup) == position::INSIDE);
}

TEST(IS_INSIDE_ERROR_UNIT, TETRAHEDRON_INSIDE_HALFSPACE)
{
    auto tet = make_unit_tetrahedron();

    halfspace hspace (vector{0.,0., 1.}, vector{0.,0.,-1.});

    auto pos = is_inside_error(tet, hspace); 

    ASSERT_TRUE(pos == position::INSIDE);  

    hspace.flip(); 

    pos = is_inside_error(tet, hspace); 

    ASSERT_TRUE(pos != position::INSIDE);  
}

// TODO
//TEST(INSIDE, POINT_IN_POLYGON)
//{
    //pointVector polygon = {point(0,0,0), point(1,0,0), point(1,1,0), point(0,1,0)}; 

    //auto polygon_points_inside = [&]()
    //{
        //// Polygon points are inside the polygon with respect to the normal.
        //vector n = normal_area_vec(polygon);  
        //for (const auto& point : polygon)
            //ASSERT_TRUE(is_inside(polygon, point, n))
                //<< "point = " << point << "\n"
                //<< "normal = " << n << "\n"; 

        //// All points of a polygon are inside the polygon.
        //for (const auto& point : polygon)
            //ASSERT_TRUE(is_inside(polygon, point))
                //<< "point = " << point << "\n"
                //<< "normal = " << n << "\n"; 
    //};

    //polygon_points_inside(); 

    //// Perform a perturbation test. 
    //auto originalPolygon = polygon; 

    //for (unsigned int i = 0; i < 1e5; ++i)
    //{
        //perturb(polygon, EPSILON); 
        //polygon_points_inside(); 
        //polygon = originalPolygon; 
    //}
//}

// TODO
//TEST(INSIDE, OVERLAPPING_POLYGONS)
//{
    //pointVector P = {point(0,0,0), point(1,0,0), point(1,1,0), point(0,1,0)}; 
    //pointVector Q = {point(0.5,0.5,0), point(1.5,0.5,0), point(1.5,1.5,0), point(0.5,1.5,0)}; 

    //pointVector intersectionPoints; 

    //auto iterP = P.begin();
    //auto iterQ = Q.begin(); 

    //for(; iterP != P.end(); ++iterP, ++iterQ)
    //{
        //if (is_inside(Q, *iterP))
            //intersectionPoints.push_back(*iterP); 

        //if (is_inside(P, *iterQ))
            //intersectionPoints.push_back(*iterQ); 
    //}

    //ASSERT_LE(mag(intersectionPoints.front() - Q.front()), EPSILON)
        //<< "Intersection point = " << intersectionPoints.front() << "\n"
        //<< "Q point = " << Q.front() << "\n";
    //ASSERT_LE(mag(intersectionPoints.back() - point(1,1,0)), EPSILON)
        //<< "Intersection point = " << intersectionPoints.back() << "\n"
        //<< "P point = " << point(1,1,0) << "\n";
//}

//TEST(isInside, edgeOverlappingPolygons)
//{
    //const pointVector P = {point(0,0,0), point(1,0,0), point(1,1,0), point(0,1,0)}; 
    //const pointVector Q = {point(1,0,0), point(2,0,0), point(2,1,0), point(1,1,0)}; 

    //auto mean_point_valid = [] (const pointVector& P, const pointVector& Q)
    //{
        //pointVector intersectionPoints; 
        //auto iterP = P.begin();
        //auto iterQ = Q.begin(); 
        //for(; iterP != P.end(); ++iterP, ++iterQ)
        //{
            //if (is_inside(Q, *iterP))
                //intersectionPoints.push_back(*iterP); 

            //if (is_inside(P, *iterQ))
                //intersectionPoints.push_back(*iterQ); 
        //}

        //// Mean of the intersection points equals to point(1, 0.5 , 0)
        //auto meanPoint =  0.5 * (intersectionPoints.front() + intersectionPoints.back());  

        //bool testPassed = (mag(meanPoint - point(1,0.5,0)) <= 2*EPSILON);

        //if (! testPassed)
        //{
            //Info << "First point = " << intersectionPoints.front() << "\n" 
                //<< "Second point = " << intersectionPoints.back() << "\n"
                //<< "Mean point = " << meanPoint << "\n";  
            //for (const auto& point : intersectionPoints)
                //Info << point << "\n";
        //}

        //return testPassed; 
    //};

    //ASSERT_TRUE(mean_point_valid(P,Q));

    //for (unsigned int i = 0; i < MAX_IT; ++i)
    //{
        //auto R = Q;
        //perturb(R, EPSILON); 
        //ASSERT_TRUE(mean_point_valid(P,R));
    //}
//}

int mainArgc; 
char** mainArgv; 

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

int main(int argc, char *argv[])
{
    ::testing::InitGoogleTest(&argc, argv);
   
    mainArgc = argc; 
    mainArgv = argv;

    return RUN_ALL_TESTS();
}


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