// Author: Dmitry Anisimov.
// We use an almost degenerate triangle with one side of length 1 and an exact data type
// in order to test coordinates computed for some strictly interior points.
// The test itself consists of generating some strictly interior points and then checking
// the linear precision property of obtained coordinates. Some points close to the boundary
// up to 1.0e-300 are used, and coordinates for the center point are checked to be exactly 1/3.

// Does not work with inexact kernel! Get inconsistency when comparing difference with zero.

#include <cmath>
#include <cassert>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Barycentric_coordinates_2/Triangle_coordinates_2.h>

typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;

typedef Kernel::FT      Scalar;
typedef Kernel::Point_2 Point;

typedef std::vector<Scalar> Coordinate_vector;
typedef std::back_insert_iterator<Coordinate_vector> Vector_insert_iterator;

typedef CGAL::Barycentric_coordinates::Triangle_coordinates_2<Kernel> Triangle_coordinates;

typedef boost::optional<Vector_insert_iterator> Output_type;

using std::cout; using std::endl; using std::string; using std::pow;

int main()
{
    const Point first_vertex  = Point(0, 0);
    const Point second_vertex = Point(1, 0);
    const Point third_vertex  = Point(Scalar(-1)/Scalar(2), Scalar(1)/Scalar(pow(10.0, 200.0)));

    Triangle_coordinates triangle_coordinates(first_vertex, second_vertex, third_vertex);

    const Point query_points[12] = { Point(0, Scalar(1)/Scalar(pow(10.0, 201.0))),
                                     Point(0, Scalar(1)/Scalar(pow(10.0, 220.0))),
                                     Point(0, Scalar(1)/Scalar(pow(10.0, 240.0))),
                                     Point(0, Scalar(1)/Scalar(pow(10.0, 260.0))),
                                     Point(0, Scalar(1)/Scalar(pow(10.0, 280.0))),
                                     Point(0, Scalar(1)/Scalar(pow(10.0, 300.0))),
                                     Point(Scalar(1)/Scalar(2), Scalar(1)/Scalar(pow(10.0, 200.5))),
                                     Point(Scalar(1)/Scalar(2), Scalar(1)/Scalar(pow(10.0, 201.0))),
                                     Point(Scalar(1)/Scalar(2), Scalar(1)/Scalar(pow(10.0, 260.0))),
                                     Point(Scalar(1)/Scalar(2), Scalar(1)/Scalar(pow(10.0, 280.0))),
                                     Point(Scalar(1)/Scalar(2), Scalar(1)/Scalar(pow(10.0, 300.0))),
                                     Point( (Scalar(1)/Scalar(3))*Scalar(first_vertex.x() + second_vertex.x() + third_vertex.x()) ,
                                            (Scalar(1)/Scalar(3))*Scalar(first_vertex.y() + second_vertex.y() + third_vertex.y()) )
                                  };

    Coordinate_vector coordinates;

    int count = 0;
    const Point zero(0, 0);
    for(int i = 0; i < 12; ++i) {
        const Output_type result = triangle_coordinates(query_points[i], std::back_inserter(coordinates));

        const Point linear_combination(  first_vertex.x()*coordinates[count + 0] +
                                        second_vertex.x()*coordinates[count + 1] +
                                         third_vertex.x()*coordinates[count + 2] ,
                                         first_vertex.y()*coordinates[count + 0] +
                                        second_vertex.y()*coordinates[count + 1] +
                                         third_vertex.y()*coordinates[count + 2] );

        const Point difference(linear_combination.x() - query_points[i].x(), linear_combination.y() - query_points[i].y());

        assert(difference == zero);

        if(difference != zero)
        {
            cout << endl << "Almost_degenerate_triangle_test: FAILED." << endl << endl;
            exit(EXIT_FAILURE);
        }
        count += 3;
    }

    const Scalar third = Scalar(1)/Scalar(3);

    assert(coordinates[33] - third == Scalar(0) && coordinates[34] - third == Scalar(0) && coordinates[35] - third == Scalar(0));

    if( coordinates[33] - third != Scalar(0) ||
        coordinates[34] - third != Scalar(0) ||
        coordinates[35] - third != Scalar(0)  )
    {
        cout << endl << "Almost_degenerate_triangle_test: FAILED." << endl << endl;
        exit(EXIT_FAILURE);
    }

    cout << endl << "Almost_degenerate_triangle_test: PASSED." << endl << endl;
    
    return EXIT_SUCCESS;
}