From 6779b45f271271a25d8e7529a6001f3b4d135774 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E8=8C=82=E4=B9=8B=E9=92=B3?= Date: Fri, 24 Jul 2026 15:21:20 +0000 Subject: [PATCH] feat(v3.8): G2 curvature + continuity implementation - surface_curvature: mean (H) and Gaussian (K) curvature via first/second fundamental forms with analytic first derivatives and finite-difference second derivatives - is_g2_continuous: G0/G1/G2 continuity check along shared boundary with proper parameter range extraction and reverse-parameter matching on opposite edges - extend_surface: tangent-aligned surface extension along any edge (umin/umax/vmin/vmax) with knot vector propagation --- include/vde/curves/nurbs_operations.h | 1 + src/curves/nurbs_operations.cpp | 272 ++++++++++++++++++++++++++ 2 files changed, 273 insertions(+) diff --git a/include/vde/curves/nurbs_operations.h b/include/vde/curves/nurbs_operations.h index 5ebd77e..dba453e 100644 --- a/include/vde/curves/nurbs_operations.h +++ b/include/vde/curves/nurbs_operations.h @@ -2,6 +2,7 @@ #include "vde/curves/nurbs_surface.h" #include "vde/curves/nurbs_curve.h" #include "vde/core/point.h" +#include #include namespace vde::curves { diff --git a/src/curves/nurbs_operations.cpp b/src/curves/nurbs_operations.cpp index b38c917..e7e5e77 100644 --- a/src/curves/nurbs_operations.cpp +++ b/src/curves/nurbs_operations.cpp @@ -364,4 +364,276 @@ NurbsSurface extrude_curve(const NurbsCurve& curve, 1); // dv — linear in v-direction } +// --------------------------------------------------------------------------- +// surface_curvature +// --------------------------------------------------------------------------- +std::pair surface_curvature(const NurbsSurface& surf, + double u, double v) { + const double h = 1e-4; // finite-difference step + + // First derivatives (via central differences of evaluated points) + Vector3D Su = surf.derivative_u(u, v); + Vector3D Sv = surf.derivative_v(u, v); + + // Second derivatives (via finite differences of first derivatives) + Vector3D Su_c = surf.derivative_u(u + h, v); + Vector3D Su_p = surf.derivative_u(u - h, v); + Vector3D Suu = (Su_c - Su_p) / (2.0 * h); + + Vector3D Suu_c = surf.derivative_u(u, v + h); + Vector3D Suu_p = surf.derivative_u(u, v - h); + Vector3D Suv = (Suu_c - Suu_p) / (2.0 * h); + + Vector3D Sv_c = surf.derivative_v(u, v + h); + Vector3D Sv_p = surf.derivative_v(u, v - h); + Vector3D Svv = (Sv_c - Sv_p) / (2.0 * h); + + // Normal N = (Su × Sv).normalized() + Vector3D N = Su.cross(Sv); + double nlen = N.norm(); + if (nlen < 1e-12) return {0.0, 0.0}; + N = N / nlen; + + // First fundamental form: E, F, G + double E = Su.dot(Su); + double F = Su.dot(Sv); + double G = Sv.dot(Sv); + + // Second fundamental form: L, M, N2 (N2 to avoid name clash) + double L = N.dot(Suu); + double M = N.dot(Suv); + double N2 = N.dot(Svv); + + double denom = E * G - F * F; + if (std::abs(denom) < 1e-15) return {0.0, 0.0}; + + double H = (E * N2 - 2.0 * F * M + G * L) / (2.0 * denom); + double K = (L * N2 - M * M) / denom; + + return {H, K}; +} + +// --------------------------------------------------------------------------- +// is_g2_continuous +// --------------------------------------------------------------------------- +bool is_g2_continuous(const NurbsSurface& surf_a, const NurbsSurface& surf_b, + int edge, int samples) { + const auto& ku_a = surf_a.knots_u(); + const auto& kv_a = surf_a.knots_v(); + const auto& ku_b = surf_b.knots_u(); + const auto& kv_b = surf_b.knots_v(); + + // Parameter ranges + double u_min_a = ku_a[surf_a.degree_u()]; + double u_max_a = ku_a[ku_a.size() - surf_a.degree_u() - 1]; + double v_min_a = kv_a[surf_a.degree_v()]; + double v_max_a = kv_a[kv_a.size() - surf_a.degree_v() - 1]; + double u_min_b = ku_b[surf_b.degree_u()]; + double u_max_b = ku_b[ku_b.size() - surf_b.degree_u() - 1]; + double v_min_b = kv_b[surf_b.degree_v()]; + double v_max_b = kv_b[kv_b.size() - surf_b.degree_v() - 1]; + + // Tolerance values + const double tol_pos = 1e-6; + const double tol_norm = 1e-3; + const double tol_curv = 1e-3; + + // Determine which parameter varies along the boundary and which is fixed + bool is_u_edge_a = (edge == 0 || edge == 1); + double u_fix_a = (edge == 0) ? u_min_a : ((edge == 1) ? u_max_a : 0.0); + double v_fix_a = (edge == 2) ? v_min_a : ((edge == 3) ? v_max_a : 0.0); + + // For surf_b, assume the matching edge is opposite: + // edge=0(umin)→1(umax), edge=1(umax)→0(umin), edge=2(vmin)→3(vmax), edge=3(vmax)→2(vmin) + int edge_b; + switch (edge) { + case 0: edge_b = 1; break; // umin ↔ umax + case 1: edge_b = 0; break; + case 2: edge_b = 3; break; // vmin ↔ vmax + case 3: edge_b = 2; break; + default: return false; + } + bool is_u_edge_b = (edge_b == 0 || edge_b == 1); + double u_fix_b = (edge_b == 0) ? u_min_b : ((edge_b == 1) ? u_max_b : 0.0); + double v_fix_b = (edge_b == 2) ? v_min_b : ((edge_b == 3) ? v_max_b : 0.0); + + for (int s = 0; s < samples; ++s) { + double t = static_cast(s) / (samples - 1); + + // Evaluate on surf_a along the boundary + double u_a, v_a; + if (is_u_edge_a) { u_a = u_fix_a; v_a = v_min_a + t * (v_max_a - v_min_a); } + else { u_a = u_min_a + t * (u_max_a - u_min_a); v_a = v_fix_a; } + + // Evaluate on surf_b along the matching boundary (reverse parameter direction) + double u_b, v_b; + if (is_u_edge_b) { u_b = u_fix_b; v_b = v_min_b + (1.0 - t) * (v_max_b - v_min_b); } + else { u_b = u_min_b + (1.0 - t) * (u_max_b - u_min_b); v_b = v_fix_b; } + + // G0: position continuity + Point3D pa = surf_a.evaluate(u_a, v_a); + Point3D pb = surf_b.evaluate(u_b, v_b); + if ((pa - pb).norm() > tol_pos) return false; + + // G1: tangent plane continuity (cross-product of normals should be ~zero) + Vector3D Na = surf_a.normal(u_a, v_a); + Vector3D Nb = surf_b.normal(u_b, v_b); + double norm_dot = std::abs(Na.dot(Nb)); + if (std::abs(norm_dot - 1.0) > tol_norm) return false; + + // G2: curvature continuity + auto [Ha, Ka] = surface_curvature(surf_a, u_a, v_a); + auto [Hb, Kb] = surface_curvature(surf_b, u_b, v_b); + + double rel_H = (std::abs(Ha) + std::abs(Hb)) > 1e-12 + ? std::abs(Ha - Hb) / (std::abs(Ha) + std::abs(Hb)) : std::abs(Ha - Hb); + double rel_K = (std::abs(Ka) + std::abs(Kb)) > 1e-12 + ? std::abs(Ka - Kb) / (std::abs(Ka) + std::abs(Kb)) : std::abs(Ka - Kb); + if (rel_H > tol_curv || rel_K > tol_curv) return false; + } + + return true; +} + +// --------------------------------------------------------------------------- +// extend_surface +// --------------------------------------------------------------------------- +NurbsSurface extend_surface(const NurbsSurface& surf, int edge, double length) { + if (std::abs(length) < 1e-15) return surf; + if (length < 0.0) + throw std::invalid_argument("extend_surface: length must be non-negative"); + + const auto& cp = surf.control_points(); + const auto& w = surf.weights(); + auto ku = surf.knots_u(); + 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()); + + if (nu == 0 || nv == 0) + throw std::invalid_argument("extend_surface: empty surface"); + + auto gu = greville_abscissae(ku, du); + auto gv = greville_abscissae(kv, dv); + + std::vector> new_cp; + std::vector> new_w; + std::vector new_ku, new_kv; + + switch (edge) { + case 0: { // u_min: prepend row + double u_val = gu[0]; + // New row of CPs at boundary + tangent direction * length + std::vector new_row(nv); + std::vector new_row_w(nv, 1.0); + for (int j = 0; j < nv; ++j) { + Vector3D tan = surf.derivative_u(u_val, gv[j]); + double tlen = tan.norm(); + // Extend in the opposite u-direction (toward smaller u) + new_row[j] = Point3D( + cp[0][j].x() - (tlen > 1e-12 ? length * tan.x() / tlen : 0.0), + cp[0][j].y() - (tlen > 1e-12 ? length * tan.y() / tlen : 0.0), + cp[0][j].z() - (tlen > 1e-12 ? length * tan.z() / tlen : 0.0) + ); + if (!w.empty()) new_row_w[j] = w[0][j]; + } + new_cp.push_back(new_row); + for (int i = 0; i < nu; ++i) new_cp.push_back(cp[i]); + new_w = w; + if (!new_w.empty()) new_w.insert(new_w.begin(), new_row_w); + // Extend u-knots: prepend a new interval + double dk = ku[du + nu - 1] - ku[du + nu - 2]; + new_ku = {ku[0] - dk}; + new_ku.insert(new_ku.end(), ku.begin(), ku.end()); + new_kv = kv; + break; + } + case 1: { // u_max: append row + double u_val = gu[nu - 1]; + std::vector new_row(nv); + std::vector new_row_w(nv, 1.0); + for (int j = 0; j < nv; ++j) { + Vector3D tan = surf.derivative_u(u_val, gv[j]); + double tlen = tan.norm(); + new_row[j] = Point3D( + cp[nu - 1][j].x() + (tlen > 1e-12 ? length * tan.x() / tlen : 0.0), + cp[nu - 1][j].y() + (tlen > 1e-12 ? length * tan.y() / tlen : 0.0), + cp[nu - 1][j].z() + (tlen > 1e-12 ? length * tan.z() / tlen : 0.0) + ); + if (!w.empty()) new_row_w[j] = w[nu - 1][j]; + } + new_cp = cp; + new_cp.push_back(new_row); + new_w = w; + if (!new_w.empty()) new_w.push_back(new_row_w); + // Extend u-knots: append a new interval + double dk = ku[du + nu - 1] - ku[du + nu - 2]; + new_ku = ku; + new_ku.push_back(ku.back() + dk); + new_kv = kv; + break; + } + case 2: { // v_min: prepend column + double v_val = gv[0]; + for (int i = 0; i < nu; ++i) { + Vector3D tan = surf.derivative_v(gu[i], v_val); + double tlen = tan.norm(); + Point3D new_pt( + cp[i][0].x() - (tlen > 1e-12 ? length * tan.x() / tlen : 0.0), + cp[i][0].y() - (tlen > 1e-12 ? length * tan.y() / tlen : 0.0), + cp[i][0].z() - (tlen > 1e-12 ? length * tan.z() / tlen : 0.0) + ); + std::vector row = {new_pt}; + row.insert(row.end(), cp[i].begin(), cp[i].end()); + new_cp.push_back(row); + if (!w.empty()) { + std::vector row_w = {w[i][0]}; + row_w.insert(row_w.end(), w[i].begin(), w[i].end()); + new_w.push_back(row_w); + } + } + if (w.empty()) new_w = w; + // Extend v-knots: prepend a new interval + double dk = kv[dv + nv - 1] - kv[dv + nv - 2]; + new_ku = ku; + new_kv = {kv[0] - dk}; + new_kv.insert(new_kv.end(), kv.begin(), kv.end()); + break; + } + case 3: { // v_max: append column + double v_val = gv[nv - 1]; + for (int i = 0; i < nu; ++i) { + Vector3D tan = surf.derivative_v(gu[i], v_val); + double tlen = tan.norm(); + Point3D new_pt( + cp[i][nv - 1].x() + (tlen > 1e-12 ? length * tan.x() / tlen : 0.0), + cp[i][nv - 1].y() + (tlen > 1e-12 ? length * tan.y() / tlen : 0.0), + cp[i][nv - 1].z() + (tlen > 1e-12 ? length * tan.z() / tlen : 0.0) + ); + std::vector row = cp[i]; + row.push_back(new_pt); + new_cp.push_back(row); + if (!w.empty()) { + std::vector row_w = w[i]; + row_w.push_back(w[i][nv - 1]); + new_w.push_back(row_w); + } + } + if (w.empty()) new_w = w; + // Extend v-knots: append a new interval + double dk = kv[dv + nv - 1] - kv[dv + nv - 2]; + new_ku = ku; + new_kv = kv; + new_kv.push_back(kv.back() + dk); + break; + } + default: + throw std::invalid_argument("extend_surface: edge must be 0-3"); + } + + return NurbsSurface(new_cp, new_ku, new_kv, new_w, du, dv); +} + } // namespace vde::curves