diff --git a/include/vde/curves/nurbs_operations.h b/include/vde/curves/nurbs_operations.h new file mode 100644 index 0000000..6a9d0fe --- /dev/null +++ b/include/vde/curves/nurbs_operations.h @@ -0,0 +1,87 @@ +#pragma once +#include "vde/curves/nurbs_surface.h" +#include "vde/curves/nurbs_curve.h" +#include "vde/core/point.h" +#include + +namespace vde::curves { + +using core::Point3D; +using core::Vector3D; + +/// Offset a NURBS surface by a constant distance along its normal. +/// Computes the normal at each control point's Greville abscissa and +/// offsets the control point along that normal (Piegl & Tiller approximation). +/// @param surf Input NURBS surface +/// @param distance Offset distance (positive = outward along normal) +/// @return Offset NURBS surface (approximation) +[[nodiscard]] NurbsSurface offset_surface(const NurbsSurface& surf, double distance); + +/// Trim a NURBS surface to a parameter sub-domain via reparameterization. +/// Maps the old domain [u_min,u_max]×[v_min,v_max] to [0,1]×[0,1] by +/// shifting and scaling knot vectors. Control points and weights are unchanged; +/// the surface geometry is preserved, only the parameter range changes. +/// @param surf Input surface +/// @param u_min, u_max Trimmed u range (must be within original domain) +/// @param v_min, v_max Trimmed v range (must be within original domain) +/// @return Trimmed NURBS surface +[[nodiscard]] NurbsSurface trim_surface(const NurbsSurface& surf, + double u_min, double u_max, + double v_min, double v_max); + +/// Blend between two surfaces along their boundary edges. +/// Extracts boundary curves, offsets them inward by blend_radius along +/// the surface normal, and creates a ruled (lofted) surface between them. +/// @param surf_a First surface +/// @param surf_b Second surface +/// @param edge_a Edge of surf_a to blend from (0=umin, 1=umax, 2=vmin, 3=vmax) +/// @param edge_b Edge of surf_b to blend to (0=umin, 1=umax, 2=vmin, 3=vmax) +/// @param blend_radius Radius of the blend transition +/// @return Blend surface (ruled surface between offset boundary curves) +[[nodiscard]] NurbsSurface blend_surfaces(const NurbsSurface& surf_a, + const NurbsSurface& surf_b, + int edge_a, int edge_b, + double blend_radius); + +/// Create a ruled surface between two NURBS curves. +/// The surface is linear (degree=1) in the ruling direction and inherits +/// the curve degree in the longitudinal direction. +/// @param curve_a First curve (v=0 edge of ruled surface) +/// @param curve_b Second curve (v=1 edge of ruled surface) +/// @return Ruled NURBS surface +[[nodiscard]] NurbsSurface ruled_surface(const NurbsCurve& curve_a, + const NurbsCurve& curve_b); + +/// Create a Coons patch from four boundary curves. +/// Bilinear interpolation: S(u,v) = (1-u)*C_v0(v) + u*C_v1(v) +/// + (1-v)*C_u0(u) + v*C_u1(u) +/// - corner correction terms. +/// All four curves should have compatible degree and knot structure. +/// @param curve_u0 Boundary at v=0 (runs along u) +/// @param curve_u1 Boundary at v=1 (runs along u) +/// @param curve_v0 Boundary at u=0 (runs along v) +/// @param curve_v1 Boundary at u=1 (runs along v) +/// @return Coons patch NURBS surface +[[nodiscard]] NurbsSurface coons_patch(const NurbsCurve& curve_u0, + const NurbsCurve& curve_u1, + const NurbsCurve& curve_v0, + const NurbsCurve& curve_v1); + +/// Extrude a NURBS curve along a direction vector. +/// Creates a degree(d)×1 NURBS surface where d is the curve degree. +/// First row of control points = curve CPs, last row = offset CPs. +/// @param curve Profile curve +/// @param direction Extrusion direction (will be normalized internally) +/// @param length Extrusion length +/// @return Extruded NURBS surface +[[nodiscard]] NurbsSurface extrude_curve(const NurbsCurve& curve, + const Vector3D& direction, + double length); + +/// Extract a boundary isoparametric curve from a NURBS surface. +/// @param surf Input surface +/// @param edge Which boundary: 0=umin, 1=umax, 2=vmin, 3=vmax +/// @return NURBS curve along the specified boundary +[[nodiscard]] NurbsCurve extract_boundary_curve(const NurbsSurface& surf, int edge); + +} // namespace vde::curves diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 2c4d7b6..34e2bbe 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -42,6 +42,7 @@ add_library(vde_curves STATIC curves/bezier_surface.cpp curves/bspline_surface.cpp curves/nurbs_surface.cpp + curves/nurbs_operations.cpp curves/tessellation.cpp ) target_include_directories(vde_curves diff --git a/src/brep/modeling.cpp b/src/brep/modeling.cpp index 5c96fa3..aeb9292 100644 --- a/src/brep/modeling.cpp +++ b/src/brep/modeling.cpp @@ -1,5 +1,6 @@ #include "vde/brep/modeling.h" #include +#include #include #include @@ -746,14 +747,16 @@ BrepModel shell(const BrepModel& body, int open_face, double thickness) { BrepModel result; // ── 1. Offset all vertices inward ── - std::vector outer_vid(body.num_vertices()); // map old id → new outer id - std::vector inner_vid(body.num_vertices()); // map old id → new inner id + std::map outer_vid; // vertex ID → new outer vertex ID + std::map inner_vid; // vertex ID → new inner vertex ID for (size_t i = 0; i < body.num_vertices(); ++i) { - Point3D p_outer = body.vertex(static_cast(i)).point; + const auto& v = body.vertex(static_cast(i)); + int vid = v.id; + Point3D p_outer = v.point; Point3D p_inner = offset_vertex(body, static_cast(i), -thickness); - outer_vid[i] = result.add_vertex(p_outer); - inner_vid[i] = result.add_vertex(p_inner); + outer_vid[vid] = result.add_vertex(p_outer); + inner_vid[vid] = result.add_vertex(p_inner); } std::vector all_faces; diff --git a/src/curves/nurbs_operations.cpp b/src/curves/nurbs_operations.cpp new file mode 100644 index 0000000..4c1a564 --- /dev/null +++ b/src/curves/nurbs_operations.cpp @@ -0,0 +1,389 @@ +#include "vde/curves/nurbs_operations.h" +#include "vde/curves/bspline_curve.h" +#include +#include +#include + +namespace vde::curves { + +// --------------------------------------------------------------------------- +// Helper: Greville abscissae for knot vector +// --------------------------------------------------------------------------- +static std::vector greville_abscissae(const std::vector& knots, int degree) { + int n = static_cast(knots.size()) - degree - 1; // number of control points + std::vector greville(n); + for (int i = 0; i < n; ++i) { + double sum = 0.0; + for (int k = 1; k <= degree; ++k) + sum += knots[i + k]; + greville[i] = sum / degree; + } + return greville; +} + +// --------------------------------------------------------------------------- +// offset_surface +// --------------------------------------------------------------------------- +NurbsSurface offset_surface(const NurbsSurface& surf, double distance) { + if (std::abs(distance) < 1e-15) + return surf; + + const auto& cp = surf.control_points(); + const auto& w = surf.weights(); + const auto& ku = surf.knots_u(); + const auto& kv = surf.knots_v(); + int du = surf.degree_u(); + int dv = surf.degree_v(); + + auto gu = greville_abscissae(ku, du); + auto gv = greville_abscissae(kv, dv); + + int nu = static_cast(cp.size()); + int nv = static_cast(cp.empty() ? 0 : cp[0].size()); + + std::vector> new_cp(nu, std::vector(nv)); + for (int i = 0; i < nu; ++i) { + for (int j = 0; j < nv; ++j) { + Vector3D n = surf.normal(gu[i], gv[j]); + new_cp[i][j] = Point3D( + cp[i][j].x() + distance * n.x(), + cp[i][j].y() + distance * n.y(), + cp[i][j].z() + distance * n.z() + ); + } + } + + return NurbsSurface(new_cp, ku, kv, w, du, dv); +} + +// --------------------------------------------------------------------------- +// trim_surface +// --------------------------------------------------------------------------- +NurbsSurface trim_surface(const NurbsSurface& surf, + double u_min, double u_max, + double v_min, double v_max) { + if (u_min >= u_max || v_min >= v_max) + throw std::invalid_argument("trim_surface: invalid trim domain"); + + const auto& old_ku = surf.knots_u(); + const auto& old_kv = surf.knots_v(); + + double u_scale = u_max - u_min; + double v_scale = v_max - v_min; + + std::vector new_ku = old_ku; + std::vector new_kv = old_kv; + + for (auto& k : new_ku) k = (k - u_min) / u_scale; + for (auto& k : new_kv) k = (k - v_min) / v_scale; + + return NurbsSurface(surf.control_points(), new_ku, new_kv, + surf.weights(), surf.degree_u(), surf.degree_v()); +} + +// --------------------------------------------------------------------------- +// extract_boundary_curve +// --------------------------------------------------------------------------- +NurbsCurve extract_boundary_curve(const NurbsSurface& surf, int edge) { + const auto& cp = surf.control_points(); + const auto& w = surf.weights(); + int nu = static_cast(cp.size()); + int nv = static_cast(cp.empty() ? 0 : cp[0].size()); + + if (nu == 0 || nv == 0) + throw std::invalid_argument("extract_boundary_curve: empty surface"); + + std::vector curve_cp; + std::vector curve_w; + std::vector curve_knots; + int curve_deg; + + switch (edge) { + case 0: { // u = u_min → row 0, knots_v, degree_v + curve_cp = cp[0]; + curve_w = w.empty() ? std::vector(nv, 1.0) : w[0]; + curve_knots = surf.knots_v(); + curve_deg = surf.degree_v(); + break; + } + case 1: { // u = u_max → last row, knots_v, degree_v + curve_cp = cp[nu - 1]; + curve_w = w.empty() ? std::vector(nv, 1.0) : w[nu - 1]; + curve_knots = surf.knots_v(); + curve_deg = surf.degree_v(); + break; + } + case 2: { // v = v_min → column 0, knots_u, degree_u + curve_cp.reserve(nu); + curve_w.reserve(nu); + for (int i = 0; i < nu; ++i) { + curve_cp.push_back(cp[i][0]); + curve_w.push_back(w.empty() ? 1.0 : w[i][0]); + } + curve_knots = surf.knots_u(); + curve_deg = surf.degree_u(); + break; + } + case 3: { // v = v_max → last column, knots_u, degree_u + curve_cp.reserve(nu); + curve_w.reserve(nu); + for (int i = 0; i < nu; ++i) { + curve_cp.push_back(cp[i][nv - 1]); + curve_w.push_back(w.empty() ? 1.0 : w[i][nv - 1]); + } + curve_knots = surf.knots_u(); + curve_deg = surf.degree_u(); + break; + } + default: + throw std::invalid_argument("extract_boundary_curve: edge must be 0-3"); + } + + return NurbsCurve(curve_cp, curve_knots, curve_w, curve_deg); +} + +// --------------------------------------------------------------------------- +// offset a curve by moving its control points along surface normals +// --------------------------------------------------------------------------- +static NurbsCurve offset_boundary_curve(const NurbsSurface& surf, + int edge, double dist) { + const auto& cp = surf.control_points(); + const auto& w = surf.weights(); + const auto& ku = surf.knots_u(); + const auto& kv = surf.knots_v(); + int du = surf.degree_u(); + int dv = surf.degree_v(); + int nu = static_cast(cp.size()); + int nv = static_cast(cp.empty() ? 0 : cp[0].size()); + + auto gu = greville_abscissae(ku, du); + auto gv = greville_abscissae(kv, dv); + + // Negative normal = inward offset + double sign = -1.0; + + std::vector curve_cp; + std::vector curve_w; + std::vector curve_knots; + int curve_deg; + + if (edge == 0) { // u=0 + curve_cp = cp[0]; + curve_w = w.empty() ? std::vector(nv, 1.0) : w[0]; + curve_knots = kv; + curve_deg = dv; + for (int j = 0; j < nv; ++j) { + Vector3D n = surf.normal(gu[0], gv[j]); + curve_cp[j] = Point3D(curve_cp[j].x() + sign * dist * n.x(), + curve_cp[j].y() + sign * dist * n.y(), + curve_cp[j].z() + sign * dist * n.z()); + } + } else if (edge == 1) { // u=1 + curve_cp = cp[nu - 1]; + curve_w = w.empty() ? std::vector(nv, 1.0) : w[nu - 1]; + curve_knots = kv; + curve_deg = dv; + for (int j = 0; j < nv; ++j) { + Vector3D n = surf.normal(gu[nu - 1], gv[j]); + curve_cp[j] = Point3D(curve_cp[j].x() + sign * dist * n.x(), + curve_cp[j].y() + sign * dist * n.y(), + curve_cp[j].z() + sign * dist * n.z()); + } + } else if (edge == 2) { // v=0 + curve_cp.reserve(nu); + curve_w.reserve(nu); + for (int i = 0; i < nu; ++i) { + Vector3D n = surf.normal(gu[i], gv[0]); + Point3D pt = cp[i][0]; + curve_cp.push_back(Point3D(pt.x() + sign * dist * n.x(), + pt.y() + sign * dist * n.y(), + pt.z() + sign * dist * n.z())); + curve_w.push_back(w.empty() ? 1.0 : w[i][0]); + } + curve_knots = ku; + curve_deg = du; + } else { // v=1 (edge == 3) + curve_cp.reserve(nu); + curve_w.reserve(nu); + for (int i = 0; i < nu; ++i) { + Vector3D n = surf.normal(gu[i], gv[nv - 1]); + Point3D pt = cp[i][nv - 1]; + curve_cp.push_back(Point3D(pt.x() + sign * dist * n.x(), + pt.y() + sign * dist * n.y(), + pt.z() + sign * dist * n.z())); + curve_w.push_back(w.empty() ? 1.0 : w[i][nv - 1]); + } + curve_knots = ku; + curve_deg = du; + } + + return NurbsCurve(curve_cp, curve_knots, curve_w, curve_deg); +} + +// --------------------------------------------------------------------------- +// blend_surfaces +// --------------------------------------------------------------------------- +NurbsSurface blend_surfaces(const NurbsSurface& surf_a, + const NurbsSurface& surf_b, + int edge_a, int edge_b, + double blend_radius) { + auto curve_a = offset_boundary_curve(surf_a, edge_a, blend_radius); + auto curve_b = offset_boundary_curve(surf_b, edge_b, blend_radius); + + // Build a ruled surface between the two offset boundary curves. + // Use the curve with more control points as the longitudinal direction. + return ruled_surface(curve_a, curve_b); +} + +// --------------------------------------------------------------------------- +// ruled_surface +// --------------------------------------------------------------------------- +NurbsSurface ruled_surface(const NurbsCurve& curve_a, + const NurbsCurve& curve_b) { + // Delegate to the existing static factory + return NurbsSurface::ruled(curve_a, curve_b); +} + +// --------------------------------------------------------------------------- +// coons_patch +// --------------------------------------------------------------------------- +NurbsSurface coons_patch(const NurbsCurve& curve_u0, + const NurbsCurve& curve_u1, + const NurbsCurve& curve_v0, + const NurbsCurve& curve_v1) { + // u-direction curves: curve_u0 (v=0), curve_u1 (v=1) — both parameterized by u + // v-direction curves: curve_v0 (u=0), curve_v1 (u=1) — both parameterized by v + // + // For simplicity we require compatible structures. + // Degree: degree_u = degree of v-boundaries, degree_v = degree of u-boundaries. + int deg_u = curve_v0.degree(); + int deg_v = curve_u0.degree(); + + const auto& cp_v0 = curve_v0.control_points(); + const auto& cp_v1 = curve_v1.control_points(); + const auto& cp_u0 = curve_u0.control_points(); + const auto& cp_u1 = curve_u1.control_points(); + + const auto& knots_u = curve_v0.knots(); // u-knots from v-boundary + const auto& knots_v = curve_u0.knots(); // v-knots from u-boundary + + int nu = static_cast(cp_v0.size()); + int nv = static_cast(cp_u0.size()); + + auto gu = greville_abscissae(knots_u, deg_u); + auto gv = greville_abscissae(knots_v, deg_v); + + // Pre-compute weights + const auto& w_u0 = curve_u0.weights(); + const auto& w_u1 = curve_u1.weights(); + const auto& w_v0 = curve_v0.weights(); + const auto& w_v1 = curve_v1.weights(); + + // Build control grid and weights + std::vector> grid(nu, std::vector(nv)); + std::vector> w(nu, std::vector(nv, 1.0)); + + // Corners + Point3D P00 = curve_u0.evaluate(gu[0]); // u=0, v=0 + Point3D P10 = curve_u0.evaluate(gu[nu - 1]); // u=1, v=0 + Point3D P01 = curve_u1.evaluate(gu[0]); // u=0, v=1 + Point3D P11 = curve_u1.evaluate(gu[nu - 1]); // u=1, v=1 + + // Edge control points (from boundary curves) + // u=0 edge: v0 curve CPs + // u=1 edge: v1 curve CPs + // v=0 edge: u0 curve CPs + // v=1 edge: u1 curve CPs + + for (int i = 0; i < nu; ++i) { + double u_i = gu[i]; + // v=0 edge + grid[i][0] = cp_u0[i]; + w[i][0] = i < static_cast(w_u0.size()) ? w_u0[i] : 1.0; + // v=1 edge + grid[i][nv - 1] = cp_u1[i]; + w[i][nv - 1] = i < static_cast(w_u1.size()) ? w_u1[i] : 1.0; + } + for (int j = 0; j < nv; ++j) { + double v_j = gv[j]; + // u=0 edge + grid[0][j] = cp_v0[j]; + w[0][j] = j < static_cast(w_v0.size()) ? w_v0[j] : 1.0; + // u=1 edge + grid[nu - 1][j] = cp_v1[j]; + w[nu - 1][j] = j < static_cast(w_v1.size()) ? w_v1[j] : 1.0; + } + + // Corners from boundary CPs + grid[0][0] = cp_u0[0]; // P00 + grid[nu-1][0] = cp_u0[nu-1]; // P10 + grid[0][nv-1] = cp_u1[0]; // P01 + grid[nu-1][nv-1] = cp_u1[nu-1]; // P11 + + // Interior points via bilinear blending (Coons formula) + for (int i = 1; i < nu - 1; ++i) { + double u_i = gu[i]; + for (int j = 1; j < nv - 1; ++j) { + double v_j = gv[j]; + + // S(u,v) = (1-u)*C_v0(v) + u*C_v1(v) + // + (1-v)*C_u0(u) + v*C_u1(u) + // - corner terms + Point3D term1 = (1.0 - u_i) * grid[0][j] + u_i * grid[nu - 1][j]; + Point3D term2 = (1.0 - v_j) * grid[i][0] + v_j * grid[i][nv - 1]; + Point3D term3 = (1.0 - u_i) * (1.0 - v_j) * grid[0][0] + + u_i * (1.0 - v_j) * grid[nu - 1][0] + + (1.0 - u_i) * v_j * grid[0][nv - 1] + + u_i * v_j * grid[nu - 1][nv - 1]; + + grid[i][j] = Point3D( + term1.x() + term2.x() - term3.x(), + term1.y() + term2.y() - term3.y(), + term1.z() + term2.z() - term3.z() + ); + } + } + + return NurbsSurface(grid, knots_u, knots_v, w, deg_u, deg_v); +} + +// --------------------------------------------------------------------------- +// extrude_curve +// --------------------------------------------------------------------------- +NurbsSurface extrude_curve(const NurbsCurve& curve, + const Vector3D& direction, + double length) { + const auto& cp = curve.control_points(); + const auto& weights = curve.weights(); + int n = static_cast(cp.size()); + + Vector3D dir_n = direction.normalized(); + Vector3D offset = dir_n * length; + + std::vector> grid(2); + std::vector> w(2); + grid[0] = cp; + grid[1].reserve(n); + w[0] = weights; + w[1].reserve(n); + + for (int i = 0; i < n; ++i) { + grid[1].push_back(Point3D( + cp[i].x() + offset.x(), + cp[i].y() + offset.y(), + cp[i].z() + offset.z() + )); + w[1].push_back(i < static_cast(weights.size()) ? weights[i] : 1.0); + } + + // v-knots: linear degree-1 in extrusion direction + // u-knots: inherit from curve + return NurbsSurface(grid, + {0.0, 0.0, 1.0, 1.0}, // v-direction knots (degree 1) + curve.knots(), // u-direction knots + w, + 1, // degree_v = 1 (linear extrusion) + curve.degree()); // degree_u = curve's degree +} + +} // namespace vde::curves diff --git a/tests/curves/CMakeLists.txt b/tests/curves/CMakeLists.txt index 6fa47dd..0050782 100644 --- a/tests/curves/CMakeLists.txt +++ b/tests/curves/CMakeLists.txt @@ -1,2 +1,3 @@ add_vde_test(test_bezier) add_vde_test(test_nurbs) +add_vde_test(test_nurbs_operations) diff --git a/tests/curves/test_nurbs_operations.cpp b/tests/curves/test_nurbs_operations.cpp new file mode 100644 index 0000000..cf47922 --- /dev/null +++ b/tests/curves/test_nurbs_operations.cpp @@ -0,0 +1,322 @@ +#include +#include "vde/curves/nurbs_operations.h" +#include "vde/curves/nurbs_curve.h" +#include "vde/curves/nurbs_surface.h" +#include + +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> 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 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> 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()); +}