#include <gtest/gtest.h>
#include "vde/curves/surface_intersection.h"
#include "vde/curves/nurbs_surface.h"
#include "vde/curves/nurbs_curve.h"
#include "vde/curves/nurbs_operations.h"
#include "vde/core/plane.h"
#include <cmath>

using namespace vde::curves;
using vde::core::Point3D;
using vde::core::Vector3D;

// ===========================================================================
//  Helper: unit-square planar surface on z=0
// ===========================================================================
static NurbsSurface make_plane_xy() {
    return NurbsSurface(
        { {Point3D(0,0,0), Point3D(0,1,0)},
          {Point3D(1,0,0), Point3D(1,1,0)} },
        {0,0,1,1}, {0,0,1,1}, {}, 1, 1);
}

// ===========================================================================
//  Helper: NURBS cylinder (radius 1, height 2, along Z, centered at origin)
//  Degree 2 in u (circular), degree 1 in v (axial)
// ===========================================================================
static NurbsSurface make_cylinder(double radius = 1.0, double height = 2.0) {
    double hh = height * 0.5;
    double r  = radius;
    double w  = std::sqrt(2.0) / 2.0;

    // 9 control points around the full circle (degree 2)
    std::vector<Point3D> circle_bottom = {
        { r, 0, -hh}, { r, r, -hh}, { 0, r, -hh}, {-r, r, -hh},
        {-r, 0, -hh}, {-r,-r, -hh}, { 0,-r, -hh}, { r,-r, -hh},
        { r, 0, -hh}   // wrap closing point
    };
    std::vector<Point3D> circle_top = {
        { r, 0,  hh}, { r, r,  hh}, { 0, r,  hh}, {-r, r,  hh},
        {-r, 0,  hh}, {-r,-r,  hh}, { 0,-r,  hh}, { r,-r,  hh},
        { r, 0,  hh}
    };

    std::vector<std::vector<Point3D>> grid;
    std::vector<std::vector<double>> weights;
    std::vector<double> u_weights = {1, w, 1, w, 1, w, 1, w, 1};
    for (int i = 0; i < 9; ++i) {
        grid.push_back({circle_bottom[i], circle_top[i]});
        weights.push_back({u_weights[i], u_weights[i]});
    }

    return NurbsSurface(
        grid,
        {0,0,0, 0.25,0.25, 0.5,0.5, 0.75,0.75, 1,1,1},
        {0,0,1,1},
        weights,
        2, 1);
}

// ===========================================================================
//  Helper: tilted planar surface (degree 1×1, spans [0,1]×[0,1])
//  Rotates the standard XY plane by angle around X axis
// ===========================================================================
static NurbsSurface make_tilted_plane(double angle_rad) {
    double c = std::cos(angle_rad), s = std::sin(angle_rad);
    return NurbsSurface(
        { {Point3D(0, 0, 0),         Point3D(0, c, s)},
          {Point3D(1, 0, 0),         Point3D(1, c, s)} },
        {0,0,1,1}, {0,0,1,1}, {}, 1, 1);
}

// ===========================================================================
//  Test: Plane–surface intersection
// ===========================================================================

TEST(SurfaceIntersection, PlaneCylinderHorizontal) {
    auto cyl = make_cylinder(1.0, 4.0);

    // Cut with z=0 plane (horizontal mid-plane)
    auto curves = intersect_plane_surface(cyl, Point3D(0,0,0), Vector3D(0,0,1), 32);

    ASSERT_GE(curves.size(), 1u);
    EXPECT_GE(curves[0].points.size(), 8u);

    // All points must lie in the z=0 plane
    for (const auto& p : curves[0].points) {
        EXPECT_NEAR(p.z(), 0.0, 1e-4);
    }

    // All points must be on the cylinder surface: x² + y² ≈ radius²
    for (const auto& p : curves[0].points) {
        double r_sq = p.x()*p.x() + p.y()*p.y();
        EXPECT_NEAR(r_sq, 1.0, 0.15);
    }
}

TEST(SurfaceIntersection, PlaneCylinderVertical) {
    auto cyl = make_cylinder(1.0, 4.0);

    // Cut with x=0 plane (vertical slice)
    auto curves = intersect_plane_surface(cyl, Point3D(0,0,0), Vector3D(1,0,0), 32);

    ASSERT_GE(curves.size(), 1u);

    // Should get two curves (two lines along the cylinder walls at x=0)
    // The curves should have points with x ≈ 0
    for (const auto& curve : curves) {
        for (const auto& p : curve.points) {
            EXPECT_NEAR(p.x(), 0.0, 1e-3);
        }
    }
}

TEST(SurfaceIntersection, PlaneCylinderOblique) {
    auto cyl = make_cylinder(1.0, 4.0);

    // Cut with 45° plane: z = y
    Vector3D normal(0, -1, 1);
    normal.normalize();
    auto curves = intersect_plane_surface(cyl, Point3D(0,0,0), normal, 32);

    ASSERT_GE(curves.size(), 1u);

    // All points should satisfy: z ≈ y (i.e., y - z ≈ 0)
    for (const auto& p : curves[0].points) {
        EXPECT_NEAR(p.y(), p.z(), 1e-3);
    }
}

TEST(SurfaceIntersection, PlanePlanarSurface) {
    auto plane = make_plane_xy();

    // Cut xy-plane with z=0 plane → intersection is the entire surface
    auto curves = intersect_plane_surface(plane, Point3D(0,0,0), Vector3D(0,0,1), 16);

    ASSERT_GE(curves.size(), 1u);

    // All points on the intersection curve should have z ≈ 0
    for (const auto& p : curves[0].points) {
        EXPECT_NEAR(p.z(), 0.0, 1e-4);
    }
}

TEST(SurfaceIntersection, PlaneNoIntersection) {
    auto plane = make_plane_xy();

    // Cut with plane z=5 (no intersection)
    auto curves = intersect_plane_surface(plane, Point3D(0,0,5), Vector3D(0,0,1), 16);

    EXPECT_EQ(curves.size(), 0u);
}

// ===========================================================================
//  Test: Surface–surface intersection
// ===========================================================================

TEST(SurfaceIntersection, ParallelPlanesNoIntersection) {
    auto p1 = make_plane_xy();
    // Plane shifted to z=2
    NurbsSurface p2(
        { {Point3D(0,0,2), Point3D(0,1,2)},
          {Point3D(1,0,2), Point3D(1,1,2)} },
        {0,0,1,1}, {0,0,1,1}, {}, 1, 1);

    auto curves = intersect_surfaces(p1, p2, 4, 1e-4);
    EXPECT_EQ(curves.size(), 0u);
}

TEST(SurfaceIntersection, NonIntersectingDistant) {
    auto p1 = make_plane_xy();
    // Plane shifted to z=10
    NurbsSurface p2(
        { {Point3D(0,0,10), Point3D(0,1,10)},
          {Point3D(1,0,10), Point3D(1,1,10)} },
        {0,0,1,1}, {0,0,1,1}, {}, 1, 1);

    auto curves = intersect_surfaces(p1, p2, 3, 1e-4);
    EXPECT_EQ(curves.size(), 0u);
}

TEST(SurfaceIntersection, TwoIntersectingPlanes) {
    auto p1 = make_plane_xy();  // z = 0

    // Tilted plane: rotate around X axis by 45° → line of intersection at y=0, z=0
    auto p2 = make_tilted_plane(M_PI / 4.0);  // z = y

    auto curves = intersect_surfaces(p1, p2, 5, 1e-4);

    // Two planes intersect in a line
    ASSERT_GE(curves.size(), 1u);
    ASSERT_GE(curves[0].points.size(), 3u);

    // Intersection line: y = 0, z = 0 (along X axis)
    for (const auto& p : curves[0].points) {
        EXPECT_NEAR(p.z(), 0.0, 1e-3);
        EXPECT_NEAR(p.y(), 0.0, 1e-3);
    }
}

TEST(SurfaceIntersection, SelfIntersectionSameSurface) {
    auto p1 = make_plane_xy();
    auto p2 = make_plane_xy();  // identical copy

    auto curves = intersect_surfaces(p1, p2, 4, 1e-3);
    // Two identical coplanar surfaces → should find intersection
    ASSERT_GE(curves.size(), 1u);
    EXPECT_GE(curves[0].points.size(), 3u);

    // All intersection points should be coplanar (z ≈ 0)
    for (const auto& p : curves[0].points) {
        EXPECT_NEAR(p.z(), 0.0, 1e-3);
    }
}

TEST(SurfaceIntersection, CylinderAndPlane) {
    auto cyl = make_cylinder(2.0, 6.0);
    auto plane = make_plane_xy();  // z=0 plane

    auto curves = intersect_surfaces(cyl, plane, 5, 1e-3);

    // Should find intersection curves or at least fail gracefully
    // A cylinder intersected with its mid-plane should yield a circle
    if (curves.size() >= 1 && curves[0].points.size() >= 4) {
        for (const auto& p : curves[0].points) {
            EXPECT_NEAR(p.z(), 0.0, 1e-2);
        }
    }
    // Note: surface-surface intersection with flat surface may produce partial results
    SUCCEED();  // at minimum the algorithm should not crash
}

TEST(SurfaceIntersection, ReturnsEmptyForNonIntersecting) {
    auto p1 = make_plane_xy();

    // Plane far away at z=100
    NurbsSurface p2(
        { {Point3D(100,100,100), Point3D(100,101,100)},
          {Point3D(101,100,100), Point3D(101,101,100)} },
        {0,0,1,1}, {0,0,1,1}, {}, 1, 1);

    auto curves = intersect_surfaces(p1, p2, 2, 1e-4);
    EXPECT_EQ(curves.size(), 0u);
}

// ===========================================================================
//  Test: Sphere–cylinder intersection (using revolve for sphere)
// ===========================================================================

TEST(SurfaceIntersection, SphereCylinder) {
    // Create a NURBS sphere via revolution of a half-circle
    // Half-circle in XZ plane, centered at origin, radius 2
    double r = 2.0;
    double w = std::sqrt(2.0) / 2.0;

    NurbsCurve half_circle(
        {Point3D(r, 0, 0), Point3D(r, 0, r), Point3D(0, 0, r),
         Point3D(-r, 0, r), Point3D(-r, 0, 0)},
        {0,0,0, 0.5, 0.5, 1, 1,1},
        {1, w, 1, w, 1},
        2);

    auto sphere = NurbsSurface::revolve(half_circle, Point3D(0,0,0),
                                         Vector3D(0,0,1), 2.0 * M_PI);

    // Cylinder along Z, radius 1.5, centered at origin
    auto cyl = make_cylinder(1.5, 6.0);

    auto curves = intersect_surfaces(sphere, cyl, 4, 1e-3);

    // If intersection is found, verify points lie on both surfaces
    for (const auto& curve : curves) {
        for (const auto& p : curve.points) {
            // On sphere: x² + y² + z² ≈ r²
            double d_sphere = std::abs(p.x()*p.x() + p.y()*p.y() + p.z()*p.z() - r*r);
            EXPECT_LT(d_sphere, 0.5);  // generous tolerance for NURBS approximation

            // On cylinder: x² + y² ≈ R²
            double d_cyl = std::abs(p.x()*p.x() + p.y()*p.y() - 1.5*1.5);
            EXPECT_LT(d_cyl, 0.5);
        }
    }
    SUCCEED();
}

// ===========================================================================
//  Test: Curve fitting
// ===========================================================================

TEST(SurfaceIntersection, FittedCurveMatchesPoints) {
    auto p1 = make_plane_xy();
    auto p2 = make_tilted_plane(M_PI / 6.0);

    auto curves = intersect_surfaces(p1, p2, 5, 1e-3);

    if (curves.size() >= 1) {
        // Verify intersection points lie near the Y=0, Z=0 line
        for (const auto& p : curves[0].points) {
            EXPECT_NEAR(p.z(), 0.0, 1e-2);
            EXPECT_NEAR(p.y(), 0.0, 5e-2);  // looser due to angle
        }
    }
    SUCCEED();
}

// ===========================================================================
//  Test: Algorithm robustness
// ===========================================================================

TEST(SurfaceIntersection, HandlesEmptySurface) {
    // This just validates the algorithm doesn't crash on edge cases
    auto p1 = make_plane_xy();
    // Degenerate: single-point surface with valid NURBS (1 CP, degree 1 needs 3 knots)
    NurbsSurface tiny(
        { {Point3D(0,0,0)} },
        {0,0,1}, {0,0,1}, {}, 1, 1);

    auto curves = intersect_surfaces(p1, tiny, 2, 1e-3);
    // May or may not find intersection; just shouldn't crash
    SUCCEED();
}

TEST(SurfaceIntersection, MaxDepthLimitsSubdivision) {
    auto p1 = make_plane_xy();
    auto p2 = make_tilted_plane(M_PI / 4.0);

    // Max depth = 1 should still work (coarser but functional)
    auto curves = intersect_surfaces(p1, p2, 1, 1e-2);
    // Either finds intersection or returns empty; shouldn't crash
    SUCCEED();
}

// ===========================================================================
//  Test: intersect_plane_surface boundary cases
// ===========================================================================

TEST(SurfaceIntersection, PlaneSurfaceBoundaryTouching) {
    auto plane = make_plane_xy();

    // Plane z=0 just touches the surface (which IS at z=0)
    auto curves = intersect_plane_surface(plane, Point3D(0.5,0.5,0), Vector3D(0,0,1), 16);

    // The entire surface is on the plane → should find boundary contours
    ASSERT_GE(curves.size(), 1u);
    for (const auto& p : curves[0].points) {
        EXPECT_NEAR(p.z(), 0.0, 1e-4);
    }
}