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

using namespace vde::curves;

// ---------------------------------------------------------------------------
//  Helper: create a simple planar NURBS surface (unit square on XY plane)
// ---------------------------------------------------------------------------
static NurbsSurface make_planar_surface() {
    // 2×2 control grid, degree 1×1, spanning [0,1]×[0,1] on XY plane
    std::vector<std::vector<Point3D>> grid = {
        {Point3D(0,0,0), Point3D(0,1,0)},
        {Point3D(1,0,0), Point3D(1,1,0)}
    };
    return NurbsSurface(grid, {0,0,1,1}, {0,0,1,1}, {}, 1, 1);
}

// ---------------------------------------------------------------------------
//  Helper: create a planar NURBS surface from origin + direction vectors
static NurbsSurface make_plane_surface(const Point3D& origin,
                                        const Point3D& u_dir,
                                        const Point3D& v_dir) {
    std::vector<std::vector<Point3D>> grid = {
        {origin, Point3D(origin.x() + v_dir.x(), origin.y() + v_dir.y(), origin.z() + v_dir.z())},
        {Point3D(origin.x() + u_dir.x(), origin.y() + u_dir.y(), origin.z() + u_dir.z()),
         Point3D(origin.x() + u_dir.x() + v_dir.x(),
                 origin.y() + u_dir.y() + v_dir.y(),
                 origin.z() + u_dir.z() + v_dir.z())}
    };
    return NurbsSurface(grid, {0,0,1,1}, {0,0,1,1}, {}, 1, 1);
}

//  Helper: create a simple NURBS line curve
// ---------------------------------------------------------------------------
static NurbsCurve make_line_curve(const Point3D& a, const Point3D& b) {
    return NurbsCurve({a, b}, {0,0,1,1}, {1,1}, 1);
}

// ===========================================================================
//  Test: offset_surface
// ===========================================================================

TEST(NurbsOpsTest, OffsetPlanarSurface) {
    auto plane = make_planar_surface();
    auto offset_plane = offset_surface(plane, 1.0);

    // Center of offset plane should be at (0.5, 0.5, 1.0) — normal of XY plane is +Z
    Point3D center = offset_plane.evaluate(0.5, 0.5);
    EXPECT_NEAR(center.x(), 0.5, 1e-6);
    EXPECT_NEAR(center.y(), 0.5, 1e-6);
    EXPECT_NEAR(center.z(), 1.0, 1e-6);
}

TEST(NurbsOpsTest, OffsetPlanarSurfaceInward) {
    auto plane = make_planar_surface();
    auto offset_plane = offset_surface(plane, -1.0);

    Point3D center = offset_plane.evaluate(0.5, 0.5);
    EXPECT_NEAR(center.x(), 0.5, 1e-6);
    EXPECT_NEAR(center.y(), 0.5, 1e-6);
    EXPECT_NEAR(center.z(), -1.0, 1e-6);
}

TEST(NurbsOpsTest, OffsetZeroDistance) {
    auto plane = make_planar_surface();
    auto same = offset_surface(plane, 0.0);

    // Should return the same surface (no offset)
    Point3D p1 = plane.evaluate(0.3, 0.7);
    Point3D p2 = same.evaluate(0.3, 0.7);
    EXPECT_NEAR(p1.x(), p2.x(), 1e-10);
    EXPECT_NEAR(p1.y(), p2.y(), 1e-10);
    EXPECT_NEAR(p1.z(), p2.z(), 1e-10);
}

// ===========================================================================
//  Test: trim_surface
// ===========================================================================

TEST(NurbsOpsTest, TrimSurfacePreservesGeometry) {
    auto plane = make_planar_surface();

    // Trim to inner half [0.25, 0.75] in both directions
    auto trimmed = trim_surface(plane, 0.25, 0.75, 0.25, 0.75);

    // At parameter 0.5 on trimmed surface →
    // original at 0.25 + 0.5*(0.75-0.25) = 0.5
    Point3D p_orig  = plane.evaluate(0.5, 0.5);
    Point3D p_trim  = trimmed.evaluate(0.5, 0.5);

    EXPECT_NEAR(p_orig.x(), p_trim.x(), 1e-10);
    EXPECT_NEAR(p_orig.y(), p_trim.y(), 1e-10);
    EXPECT_NEAR(p_orig.z(), p_trim.z(), 1e-10);
}

TEST(NurbsOpsTest, TrimSurfaceMapsDomain) {
    auto plane = make_planar_surface();

    // Trim to [0, 0.5] → parameter 1.0 on trimmed should equal 0.5 on original
    auto trimmed = trim_surface(plane, 0.0, 0.5, 0.0, 0.5);

    Point3D p_orig = plane.evaluate(0.5, 0.5);
    Point3D p_trim = trimmed.evaluate(1.0, 1.0);
    EXPECT_NEAR(p_orig.x(), p_trim.x(), 1e-10);
    EXPECT_NEAR(p_orig.y(), p_trim.y(), 1e-10);
    EXPECT_NEAR(p_orig.z(), p_trim.z(), 1e-10);
}

// ===========================================================================
//  Test: ruled_surface
// ===========================================================================

TEST(NurbsOpsTest, RuledSurfaceParallelLines) {
    auto line_a = make_line_curve(Point3D(0,0,0), Point3D(1,0,0));
    auto line_b = make_line_curve(Point3D(0,1,0), Point3D(1,1,0));

    auto ruled = ruled_surface(line_a, line_b);

    // The surface should be planar (z=0 at all points)
    for (double u = 0.0; u <= 1.0; u += 0.25) {
        for (double v = 0.0; v <= 1.0; v += 0.25) {
            Point3D p = ruled.evaluate(u, v);
            EXPECT_NEAR(p.z(), 0.0, 1e-10);
        }
    }

    // Midpoint should be at (0.5, 0.5, 0)
    Point3D mid = ruled.evaluate(0.5, 0.5);
    EXPECT_NEAR(mid.x(), 0.5, 1e-10);
    EXPECT_NEAR(mid.y(), 0.5, 1e-10);
    EXPECT_NEAR(mid.z(), 0.0, 1e-10);
}

TEST(NurbsOpsTest, RuledSurfaceInterpolatesEndCurves) {
    auto line_a = make_line_curve(Point3D(0,0,0), Point3D(1,0,0));
    auto line_b = make_line_curve(Point3D(0,1,1), Point3D(1,1,1)); // Elevated

    auto ruled = ruled_surface(line_a, line_b);

    // v=0 should be on curve_a, v=1 should be on curve_b
    Point3D pa = ruled.evaluate(0.5, 0.0);
    Point3D pb = ruled.evaluate(0.5, 1.0);
    Point3D ca = line_a.evaluate(0.5);
    Point3D cb = line_b.evaluate(0.5);

    EXPECT_NEAR((pa - ca).norm(), 0.0, 1e-10);
    EXPECT_NEAR((pb - cb).norm(), 0.0, 1e-10);
}

// ===========================================================================
//  Test: extrude_curve
// ===========================================================================

TEST(NurbsOpsTest, ExtrudeLineSegment) {
    auto line = make_line_curve(Point3D(0,0,0), Point3D(1,0,0));
    Vector3D dir(0, 0, 1);

    auto extruded = extrude_curve(line, dir, 2.0);

    // The extruded surface should be planar in XZ
    // v=0 along original curve: (0,0,0) → (1,0,0)
    Point3D p00 = extruded.evaluate(0.0, 0.0);
    Point3D p10 = extruded.evaluate(1.0, 0.0);
    EXPECT_NEAR(p00.x(), 0.0, 1e-10);
    EXPECT_NEAR(p00.z(), 0.0, 1e-10);
    EXPECT_NEAR(p10.x(), 1.0, 1e-10);
    EXPECT_NEAR(p10.z(), 0.0, 1e-10);

    // v=1 along extruded curve: (0,0,2) → (1,0,2)
    Point3D p01 = extruded.evaluate(0.0, 1.0);
    Point3D p11 = extruded.evaluate(1.0, 1.0);
    EXPECT_NEAR(p01.x(), 0.0, 1e-10);
    EXPECT_NEAR(p01.z(), 2.0, 1e-10);
    EXPECT_NEAR(p11.x(), 1.0, 1e-10);
    EXPECT_NEAR(p11.z(), 2.0, 1e-10);

    // Midpoint
    Point3D mid = extruded.evaluate(0.5, 0.5);
    EXPECT_NEAR(mid.x(), 0.5, 1e-10);
    EXPECT_NEAR(mid.y(), 0.0, 1e-10);
    EXPECT_NEAR(mid.z(), 1.0, 1e-10);
}

TEST(NurbsOpsTest, ExtrudeCurveYieldsPlanarSurface) {
    auto line = make_line_curve(Point3D(0,0,0), Point3D(2,0,0));
    Vector3D dir(0, 1, 0);

    auto extruded = extrude_curve(line, dir, 3.0);

    // All points should have x+0.5*y on the line, basically flat in XZ at y-dependent values
    for (double u = 0.0; u <= 1.0; u += 0.25) {
        for (double v = 0.0; v <= 1.0; v += 0.25) {
            Point3D p = extruded.evaluate(u, v);
            EXPECT_NEAR(p.x(), u * 2.0, 1e-10);
            EXPECT_NEAR(p.y(), v * 3.0, 1e-10);
            EXPECT_NEAR(p.z(), 0.0, 1e-10);
        }
    }
}

// ===========================================================================
//  Test: coons_patch
// ===========================================================================

TEST(NurbsOpsTest, CoonsPatchPassesThroughCorners) {
    // Four boundary curves forming a square
    auto u0 = make_line_curve(Point3D(0,0,0), Point3D(1,0,0));   // v=0, along u
    auto u1 = make_line_curve(Point3D(0,1,0), Point3D(1,1,0));   // v=1, along u
    auto v0 = make_line_curve(Point3D(0,0,0), Point3D(0,1,0));   // u=0, along v
    auto v1 = make_line_curve(Point3D(1,0,0), Point3D(1,1,0));   // u=1, along v

    auto patch = coons_patch(u0, u1, v0, v1);

    // Corners should match
    Point3D p00 = patch.evaluate(0.0, 0.0);
    Point3D p10 = patch.evaluate(1.0, 0.0);
    Point3D p01 = patch.evaluate(0.0, 1.0);
    Point3D p11 = patch.evaluate(1.0, 1.0);

    EXPECT_NEAR((p00 - Point3D(0,0,0)).norm(), 0.0, 1e-10);
    EXPECT_NEAR((p10 - Point3D(1,0,0)).norm(), 0.0, 1e-10);
    EXPECT_NEAR((p01 - Point3D(0,1,0)).norm(), 0.0, 1e-10);
    EXPECT_NEAR((p11 - Point3D(1,1,0)).norm(), 0.0, 1e-10);
}

TEST(NurbsOpsTest, CoonsPatchMidpoint) {
    auto u0 = make_line_curve(Point3D(0,0,0), Point3D(2,0,0));
    auto u1 = make_line_curve(Point3D(0,2,0), Point3D(2,2,0));
    auto v0 = make_line_curve(Point3D(0,0,0), Point3D(0,2,0));
    auto v1 = make_line_curve(Point3D(2,0,0), Point3D(2,2,0));

    auto patch = coons_patch(u0, u1, v0, v1);

    // Midpoint of planar square → (1,1,0)
    Point3D mid = patch.evaluate(0.5, 0.5);
    EXPECT_NEAR(mid.x(), 1.0, 1e-10);
    EXPECT_NEAR(mid.y(), 1.0, 1e-10);
    EXPECT_NEAR(mid.z(), 0.0, 1e-10);
}

// ===========================================================================
//  Test: extract_boundary_curve
// ===========================================================================

TEST(NurbsOpsTest, ExtractBoundaryUmin) {
    auto plane = make_planar_surface();
    auto curve = extract_boundary_curve(plane, 0); // u=0

    // Should be the line from (0,0,0) to (0,1,0)
    Point3D p0 = curve.evaluate(0.0);
    Point3D p1 = curve.evaluate(1.0);
    EXPECT_NEAR(p0.x(), 0.0, 1e-10);
    EXPECT_NEAR(p0.y(), 0.0, 1e-10);
    EXPECT_NEAR(p1.x(), 0.0, 1e-10);
    EXPECT_NEAR(p1.y(), 1.0, 1e-10);
}

TEST(NurbsOpsTest, ExtractBoundaryUmax) {
    auto plane = make_planar_surface();
    auto curve = extract_boundary_curve(plane, 1); // u=1

    Point3D p0 = curve.evaluate(0.0);
    Point3D p1 = curve.evaluate(1.0);
    EXPECT_NEAR(p0.x(), 1.0, 1e-10);
    EXPECT_NEAR(p0.y(), 0.0, 1e-10);
    EXPECT_NEAR(p1.x(), 1.0, 1e-10);
    EXPECT_NEAR(p1.y(), 1.0, 1e-10);
}

// ===========================================================================
//  Test: blend_surfaces
// ===========================================================================

TEST(NurbsOpsTest, BlendPlanarSurfaces) {
    // Two parallel planes separated in Z
    auto plane1 = make_planar_surface(); // z=0, on XY

    // Plane 2 shifted in Z
    std::vector<std::vector<Point3D>> grid2 = {
        {Point3D(0,0,2), Point3D(0,1,2)},
        {Point3D(1,0,2), Point3D(1,1,2)}
    };
    NurbsSurface plane2(grid2, {0,0,1,1}, {0,0,1,1}, {}, 1, 1);

    // Blend from u=1 edge of plane1 to u=0 edge of plane2
    auto blend = blend_surfaces(plane1, plane2,
                                1,  // edge_a: umax of plane1
                                0,  // edge_b: umin of plane2
                                0.5);

    // The blend should exist and evaluate
    // At v=0.5, u=0 → near plane1's x=1 edge offset inward
    // At v=0.5, u=1 → near plane2's x=0 edge offset inward
    Point3D pa = blend.evaluate(0.0, 0.5);
    Point3D pb = blend.evaluate(1.0, 0.5);

    // pa should be near x=0.5 (1.0 - 0.5 offset inward)
    EXPECT_NEAR(pa.x(), 0.5, 1e-6);
    // pb should be near x=0.5 (0.0 + 0.5 offset inward)
    EXPECT_NEAR(pb.x(), 0.5, 1e-6);
}

// ===========================================================================
//  Test: surface properties preserved
// ===========================================================================

TEST(NurbsOpsTest, ExtrudePreservesDegree) {
    // Degree-3 curve
    NurbsCurve curve(
        {Point3D(0,0,0), Point3D(1,1,0), Point3D(2,-1,0), Point3D(3,0,0)},
        {0,0,0,0,1,1,1,1},
        {1,1,1,1},
        3
    );

    auto extruded = extrude_curve(curve, Vector3D(0,0,1), 1.0);
    EXPECT_EQ(extruded.degree_u(), 3);  // curve degree preserved
    EXPECT_EQ(extruded.degree_v(), 1);  // linear in extrusion direction
}

TEST(NurbsOpsTest, OffsetPreservesDegree) {
    NurbsCurve curve(
        {Point3D(0,0,0), Point3D(1,1,0), Point3D(2,-1,0), Point3D(3,0,0)},
        {0,0,0,0,1,1,1,1},
        {1,1,1,1},
        3
    );
    auto surf = extrude_curve(curve, Vector3D(0,1,0), 1.0);

    auto offset = offset_surface(surf, 0.5);
    EXPECT_EQ(offset.degree_u(), surf.degree_u());
    EXPECT_EQ(offset.degree_v(), surf.degree_v());
}

// ===========================================================================
//  Test: fill_n_sided
// ===========================================================================

TEST(NurbsOpsTest, FillNSided_Triangle) {
    // Create triangular boundary
    auto l1 = make_line_curve(Point3D(0,0,0), Point3D(2,0,0));
    auto l2 = make_line_curve(Point3D(2,0,0), Point3D(1,2,0));
    auto l3 = make_line_curve(Point3D(1,2,0), Point3D(0,0,0));
    std::vector<NurbsCurve> boundaries = {l1, l2, l3};
    auto fill = fill_n_sided(boundaries);
    // Verify fill passes through boundary at midpoint
    Point3D mid = fill.evaluate(0.5, 0.5);
    EXPECT_TRUE(std::isfinite(mid.x()));
}

// ===========================================================================
//  Test: is_g3_continuous
// ===========================================================================

TEST(NurbsOpsTest, G3Continuous) {
    // Two planes should be G3 continuous
    auto p1 = make_plane_surface(Point3D(0,0,0), Point3D(1,0,0), Point3D(0,1,0));
    auto p2 = make_plane_surface(Point3D(1,0,0), Point3D(1,0,0), Point3D(0,1,0));
    EXPECT_TRUE(is_g3_continuous(p1, p2, 1));
}