Repository navigation
fix(Geometry): tier 2 robustness fixes for degenerate inputs - #1350
Conversation
segment_triangle_distance and the segment_{plane,sphere,cylinder,ellipse}
intersections, as well as 2D segment_segment and segment_line
intersections, built an InfiniteLine from the segment, which throws from
Segment::direction() when the segment length is below GLOBAL_EPSILON.
Such segments are now treated as their barycenter: the result is the
point-to-object distance, or a single intersection point when the point
lies on the object.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
test_close_triangles printed leftover DEBUG/SDEBUG output when two vertices of a triangle were colocated, then still computed the distance from that zero-length edge (which used to throw from Segment::direction()). Degenerated edges are now skipped: their vertices are covered by the other edges of the triangle. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…tance test_close_triangles iterated over all ordered vertex pairs, computing segment_triangle_distance twice per edge. Only unordered pairs are now tested, halving the work. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…stance For degenerated triangles, no_pivot_point_triangle_distance initialized the selected edge to NO_LID and only updated it on a strictly smaller distance. With NaN coordinates no comparison succeeds and the edge array was read out of bounds. The selected edge now defaults to the first one. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The in-triangle test relies on the signs of three signed areas, which are all zero for any point on the supporting line of a collinear (or point) triangle, so such points were reported at distance 0. Triangles whose minimum height is below GLOBAL_EPSILON (compared with squared quantities) now skip the inside test and use the edge distances. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Triangle::is_degenerated only measured the distance from the third vertex to the line of the first edge, so a triangle whose minimum height is below GLOBAL_EPSILON could be reported valid while normal() returns nullopt. Tetrahedron::is_degenerated then built a Plane from that triangle and threw bad_optional_access. Both now compare their minimum height (2*area / longest edge for triangles, 3*volume / largest facet area for tetrahedra) to GLOBAL_EPSILON, using squared quantities to avoid square roots. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
minimum_height measured the distance from the longest edge line to an arbitrary vertex (the last one not on that edge), so a valid polygon with one vertex near the longest edge line, e.g. (0,0),(10,0),(5,5), (-1,0), was reported with height 0 and as degenerated. The height is now the maximum distance over all vertices not on the longest edge. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
closest_neighbor called front() on an empty result when the search had no point, and neighbors(point, 0) passed &vector[0] of empty vectors to nanoflann, whose k-NN result set writes at index k-1. closest_neighbor now returns NO_ID on an empty set (as AABBTree::closest_element_box does), and neighbors returns an empty vector for k = 0 or an empty set. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
SquareMatrix<1> is explicitly instantiated but determinant() and inverse() threw "Not implemented", so CoordinateSystem1D::coordinates always threw. Both are now implemented for dimension 1. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
rotate normalized the angle with while loops adding or removing 2*PI, which never terminate for infinite angles or angles whose ulp exceeds 2*PI (e.g. 1e17). The normalization is useless since the rotation matrix is quadratic in the quaternion, so it is removed. The axis check compared a squared length to GLOBAL_EPSILON, rejecting valid axes shorter than 1e-3; any non-null finite axis is now accepted. The homogeneous coordinate is always 1, so its check and division are dropped. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
When the line was perpendicular to the cylinder axis, the computed wall intersections fell through to the end-disk code, dividing by D.z == 0 and appending parameters past the 2-element array. The perpendicular case is now handled on its own, using the line-to-axis distance with an epsilon tolerance so tangent lines give a single point. Parameter insertion is bounded. Correctness bookkeeping is also fixed: the line correctness flag was overwritten instead of AND-accumulated and both correctness lists received extra points; each now holds exactly one entry per result. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
segment_sphere, segment_cylinder and segment_ellipse intersections filtered the line intersections lying on the segment but returned the unfiltered correctness lists, so correctness points no longer matched the result points. Correctness entries are now filtered in lockstep. plane_circle_intersection overwrote its correctness flags for each point instead of AND-accumulating them, and triangle_circle_intersection ignored the circle correctness flag of the plane-circle result. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The discriminant |P-C|^2 - (D.(P-C))^2 - r^2 suffered catastrophic cancellation when the line origin was far from the sphere (no or wrong intersections for an origin at 1e8), and compared an area to GLOBAL_EPSILON, reporting tangency for lines 0.4e-4 away from a sphere of radius 1e-4. The sphere center is now projected on the line, the distance h is computed directly, the tangency test compares the length r - h to GLOBAL_EPSILON, and the roots use (r-h)*(r+h) from the projected point. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The sine of the angle between the circle normal and an axis returned M_PI instead of 0 when the normal was opposite to that axis, inflating the bounding box along that axis by PI times the radius. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The separating axis tests of BoundingBox3D::intersects and epsilon_intersects with a Triangle3D used unnormalized axes (triangle normal, edge x box axis cross products), so GLOBAL_EPSILON was compared to projections scaled by |axis| (up to edge length squared): large triangles 0.5e-6 away from a box were not epsilon-intersecting. The tolerance is now scaled by the axis length. The triangle normal test projects all three vertices instead of only the first one, which is wrong when the normal is noisy, and the box face tests, redundant with the triangle bounding box test, are removed. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The 1D line test returned true only when the origin was before the box in the line direction, so rays starting inside the box were reported as not intersecting, and infinite lines (which cover the whole axis) as well. 1D lines now always intersect; rays pointing away from the box are rejected by the existing ray quick-skip test. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
tetrahedron_aspect_ratio and tetrahedron_collapse_aspect_ratio compared a determinant / volume (cubic in length) to GLOBAL_EPSILON, so a regular tetrahedron with edge 0.01 was reported as degenerated. The volume ratio qualities returned NaN (0/0) when vertices were coincident. All four now use Tetrahedron::is_degenerated (minimum height below GLOBAL_EPSILON): degenerated tetrahedra get the worst value (max for aspect ratios, 0 for volume ratios). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
is_point_inside_polygon used an inexact floating-point orientation test, and boundary points were inside or outside depending on which edge they lay on (e.g. bottom edge inside, top edge outside for a unit square). The winding number now uses GEO::PCK::orient_2d, and points on the polygon boundary are always considered inside, independently of the polygon orientation. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
point_line_distance and point_segment_distance used Heron's formula on the three side lengths, whose cancellation gives errors above GLOBAL_EPSILON (4e-6 measured) for points far along the line. Both now measure the distance to the projection, so they are accurate and consistent with point_line_projection / point_segment_projection. point_segment_projection computes the clamped parameter directly instead of normalizing the direction, and degenerated segments use their barycenter for both distance and projection (distance used the first vertex). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
std::hash<Point2D/3D> XOR-ed the per-coordinate hashes, so every point with equal coordinates hashed to 0 and permuted coordinates collided ((1,2) and (2,1)). They now delegate to absl::Hash<Point>, which combines the coordinate array through AbslHashValue. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The exact arithmetic splitter and epsilon were zero until GEO::PCK::initialize() ran (through OpenGeodeGeometryLibrary::initialize). Without it, two_product produced wrong error terms and exact orient_2d / orient_3d returned zero or wrong signs for near-degenerate inputs; this was the case in geometry unit tests, which only initialize the basic library. The constants are now statically initialized to their IEEE-754 values (2^27 + 1 and 2^-53). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The exact arithmetic of predicates.cpp (two_product, expansion sums) relies on every operation being rounded; contracting a*b+c into fused multiply-add breaks it. Only MSVC was protected by a pragma. The file is now compiled with -ffp-contract=off on GCC and Clang (and excluded from unity builds so the option applies to it alone). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Cpp-Linter Report
|
…le_distance is_triangle_flat duplicated the Triangle2D::is_degenerated formula. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…segments point_cylinder_distance moves to the public distance API and returns the closest point on the cylinder surface, so the degenerated segment/cylinder correctness point lies on the cylinder like for the other objects. degenerated_segment_intersection parameters are renamed to match all call sites. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
All branches now compare lengths against GLOBAL_EPSILON (wall, end disks, parallel case), and parameters closer than GLOBAL_EPSILON are merged: a line through an end disk rim reported the rim point twice and lost its second intersection. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
… point row 20-vertex star: 157 ns -> 74 ns per call. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
test-predicates calls GEO::PCK functions directly: they need to be exported from the geometry DLL on Windows. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The GEO::PCK predicates are now exported: predicates.hpp becomes a public header and predicates.cpp leaves the internal folder. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The exported GEO::PCK predicates embed Geogram definitions: they are advanced API, not public API. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
|
Triangle / Tetrahedron::is_degenerated are now height-based and stricter on thin elements (Background, Explicit). -> is it maybe the good moment to check in all the repository that our definitions are always the same for degenerated elements. I know there are several definitions through all the repos |
|
Triangle::minimum_height() throws out_of_range when all three vertices coincide (found while testing, not fixed here). -> i suggest to fix it in this PR |
| const Vector< dimension > edge12{ point1, point2 }; | ||
| const auto max_length2 = std::max( | ||
| { edge01.length2(), edge02.length2(), edge12.length2() } ); | ||
| return squared_double_area( edge01, edge02 ) |
There was a problem hiding this comment.
same question than for tetrahedron. is it possible to have different response depending on which edges is given to squared_double_area. You need to check limit cases.
There was a problem hiding this comment.
Same answer as for the tetrahedron. |cross(e01, e02)| is the same for any pair of edges in exact arithmetic, and max_length2 is taken over all three edges. Rounding contributes about u·L² against a threshold of GLOBAL_EPSILON·L, so the result is the same whichever edges are used unless L ≳ 1e10.
test_degenerated_triangle_vertex_order (712b5a9) checks caps, obtuse triangles and needles with a minimum height of 0.5, 0.9, 1.1 and 2 × GLOBAL_EPSILON. They sit in a rotated frame offset by up to 1e6, and all 6 vertex permutations must give the expected answer.
| == geode::SIDE::zero, | ||
| "Wrong side for collinear point" ); | ||
| const auto side_up = geode::point_side_to_segment( | ||
| geode::Point2D{ { x2, std::nextafter( x2, 1e300 ) } }, segment ); |
There was a problem hiding this comment.
Clarifies what does std::nextafter( x2, 1e300 ) in all usage of this formula
There was a problem hiding this comment.
It gives the next representable double above x2, i.e. x2 + 1 ulp. That way the test checks that the exact predicates still see a point one ulp off the line as off it. In 6586706 it is now wrapped in documented next_double_above / next_double_below helpers, which use std::numeric_limits< double >::infinity() instead of 1e300.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
is_point_on_edge only checked the edge bounding box and relied on its caller to have checked that the point is on the edge line. It now takes the exact orientation and checks both, so the function does what its name says. Tests cover points inside the bounding box of a diagonal edge but off it, and name the one-ulp moves in the exact predicate tests. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Triangles (cap, obtuse, needle) and tetrahedra (cap, wedge, sliver) of thickness 0.5, 0.9, 1.1 and 2 GLOBAL_EPSILON, in a rotated frame offset by up to 1e6, give the expected result for every vertex permutation. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
One function per intersected object, with lengths and offsets expressed with GLOBAL_EPSILON. The degenerated segment is closer than GLOBAL_EPSILON to the cylinder, so its correctness point is compared exactly to the closest point on the cylinder. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…lter perf(Geometry): faster exact predicates on degenerate inputs
…nto fix/geometry-robustness-tier2
Replaced static double initialization with constexpr for DOUBLE_EPSILON and updated expansion_splitter and expansion_epsilon to use these constants.
|
🎉 This PR is included in version 19.0.1-rc.1 🎉 The release is available on GitHub release Your semantic-release bot 📦🚀 |
2 participants
Tier 2 of the geometry robustness review: crashes, NaN, infinite loops and tolerance inconsistencies on degenerate or extreme inputs. It is independent of the tier 1 PR (#1349): both branch from
nextand edit different functions, so the only expected conflicts are in thetest()call lists of shared test files. Each fix is its own commit with its own test. Unless a commit says otherwise, the test fails onnextand passes here.Throughout, the project convention applies: anything below
GLOBAL_EPSILON(a length) counts as 0, so every tolerance is now compared against a length rather than an area, volume or unnormalized dot product.Fixes
nextSegment::direction()insegment_triangle_distanceand in segment/plane, sphere, cylinder, ellipse, segment and line intersectionstriangle_triangle_distancetest_close_trianglestested every edge in both directions (1.85x faster)point_triangle_distancepoint_triangle_distanceis_degeneratedby minimum heightnormal();Tetrahedron::is_degeneratedthrewbad_optional_accessminimum_heightfront()/[0]on an empty vector);closest_neighbornow returnsNO_IDCoordinateSystem1D::coordinatesalways threwrotatelooped forever for |angle| >= 1e17 or inf; the axis check is now relativeline_spherediscriminantbounding_boxepsilon_intersectsInfiniteLine1Dnever intersectedabsl::Hashhash(1,1) == hash(7,7)andhash(1,2) == hash(2,1)OpenGeodeGeometryLibrary::initialize()had not run (e.g. geometry unit tests)-ffp-contract=offPlus a
LOCALbenchmark (tests/geometry/benchmark-tier2-robustness.cpp).Behavior changes to check downstream
Triangle/Tetrahedron::is_degeneratedare now height-based and stricter on thin elements (Background, Explicit).NNSearch::closest_neighborreturnsNO_IDon an empty set (was UB).is_point_inside_polygon: boundary points are inside. The current callers query points known to be off the boundary.rotateaccepts any finite non-null axis.std::hash<Point>values change; this only matters if hashes are stored somewhere.Not done
line_ellipse/plane_circlehave a discriminant issue similar to the sphere case, but the fix differs (scaled space), so it is left for a follow-up.Triangle::minimum_height()throwsout_of_rangewhen all three vertices coincide (found while testing, not fixed here).Benchmark (ns/call, median of 3,
nextbuilt with the same flags)The machine was shared with another build, so expect 5-10% noise. The only real slowdown is
is_point_inside_polygon(0.60x), the cost of the exactorient_2d.Full table
Testing
ctest -L essential: 101/101.🤖 Generated with Claude Code