茂之钳
af2ad029b6
- NurbsSurface: fix knot vector validation for half-circle/sphere/empty surfaces - intersect_plane_surface: fix zero-crossing, coplanar detection, analytic plane-plane - intersect_surfaces: analytical plane-plane intersection for degree-1x1 surfaces - G3 continuity: fix parameter reversal, NaN division, domain clamping - Various tolerance relaxations for NURBS approximation
348 lines
12 KiB
C++
348 lines
12 KiB
C++
#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);
|
||
}
|
||
}
|