#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Polygon_2.h>
#include <CGAL/centroid.h>
#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>

typedef CGAL::Exact_predicates_inexact_constructions_kernel  K;
typedef CGAL::Constrained_Delaunay_triangulation_2<K>        CDT;
typedef CDT::Point                                           Point;
typedef CGAL::Polygon_2<K>                                   Polygon_2;

// Helper function: computes the absolute area of a triangle using the cross product formula.
double triangleArea(const Point& a, const Point& b, const Point& c) {
    double ax = a.x(), ay = a.y();
    double bx = b.x(), by = b.y();
    double cx = c.x(), cy = c.y();
    return std::abs(ax * (by - cy) + bx * (cy - ay) + cx * (ay - by)) / 2.0;
}

int main() {
    // Example: unsorted vertices of a convex polygon.
    std::vector<Point> pts = { Point(-10000, -10000), Point(10000, 10000), Point(-10000, 10000),
    Point(10000,-10000)};

    // Step 1: Sort the vertices in counterclockwise order.
    double cx = 0.0, cy = 0.0;
    for (const auto& p : pts) {
        cx += p.x();
        cy += p.y();
    }
    cx /= pts.size();
    cy /= pts.size();
    Point centroid(cx, cy);

    std::sort(pts.begin(), pts.end(), [&centroid](const Point& a, const Point& b) {
        double angle_a = std::atan2(a.y() - centroid.y(), a.x() - centroid.x());
        double angle_b = std::atan2(b.y() - centroid.y(), b.x() - centroid.x());
        return angle_a < angle_b;
    });

    // Build a CGAL polygon from the sorted vertices.
    Polygon_2 polygon(pts.begin(), pts.end());

    // Step 2: Insert the polygon edges as constraints into a Constrained Delaunay Triangulation (CDT).
    CDT cdt;
    CDT::Vertex_handle prev, first;
    for (auto it = polygon.vertices_begin(); it != polygon.vertices_end(); ++it) {
        CDT::Vertex_handle vh = cdt.insert(*it);
        if (it == polygon.vertices_begin()) {
            first = vh;
        } else {
            cdt.insert_constraint(prev, vh);
        }
        prev = vh;
    }
    // Close the polygon.
    cdt.insert_constraint(prev, first);

    // Step 3: Sum the areas of all triangles whose centroids lie within the polygon.
    double totalArea = 0.0;
    for (auto fit = cdt.finite_faces_begin(); fit != cdt.finite_faces_end(); ++fit) {
        Point a = fit->vertex(0)->point();
        Point b = fit->vertex(1)->point();
        Point c = fit->vertex(2)->point();

        // Compute the triangle's centroid.
        K::FT t_cx = (a.x() + b.x() + c.x()) / 3;
        K::FT t_cy = (a.y() + b.y() + c.y()) / 3;
        Point triCentroid(t_cx, t_cy);

        // Include the triangle if its centroid lies inside (or on) the polygon.
        if (polygon.bounded_side(triCentroid) != CGAL::ON_UNBOUNDED_SIDE) {
            totalArea += triangleArea(a, b, c);
        }
    }

    std::cout << "Area computed by triangulation: " << totalArea << std::endl;
    std::cout << "Polygon area (using Polygon_2::area): " << std::abs(polygon.area()) << std::endl;

    return 0;
}