Repository navigation
fix(Geometry): tier 1 robustness fixes - #1349
Conversation
The global-to-local matrix was cached and never invalidated, so coordinates() kept using the old directions after set_directions() or Frame::set_direction() (used by Grid when rescaling cells). Add CoordinateSystem::set_direction resetting the cache, reset it in set_directions and after deserialization, and validate the frame in set_directions and in the points constructor (the directions constructor no longer checks a moved-from array). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
max/min_elongation_direction used the lexicographic operator< of Vector, returning the direction with the largest/smallest first coordinate instead of the longest/shortest one. Compare squared lengths. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
add_box added the other box's min and max as points; for an empty box (min = DBL_MAX, max = lowest) this blew the box up to the whole space. Take the per-axis min of mins and max of maxs instead, which also fixes AABBTree bounding boxes when some elements have empty boxes. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
When the action returned NaN (or inf) for the hinted element, every comparison against the squared distance failed and the hint was returned with a NaN distance. Restart from an empty best candidate when the hint distance is not finite. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
set_point modified a vertex without invalidating the cached pivot and normal, so normal()/pivot() kept returning values of the old triangle (also on copies, which copy the cache). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The parallel test was fabs(x) < 0., which is never true: parallel lines went through a division by zero and returned an 'incorrect' result with an infinite/NaN point. Test the dot perpendicular against exact zero so parallel lines and segments are reported as parallel. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The sign was computed from the dot product of the normal with the point position vector, i.e. relative to the world origin instead of the circle plane. Use the vector from the circle center to the point. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
segment_triangle_intersection compared the plane parameter against the exact half extent, so a segment ending on the triangle was randomly rejected depending on rounding. Accept parameters up to the half length plus GLOBAL_EPSILON (a length) and clamp the result onto the segment. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
normal() summed unit normals of the fan triangles, so for concave polygons triangles of opposite orientation cancelled out regardless of their area (e.g. a dart gave no normal, and triangulate() then picked an arbitrary orientation). Sum the un-normalized fan cross products (Newell) and return nullopt when the polygon height is below GLOBAL_EPSILON instead of relying on an exception. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
polygon_area_sign returned the sign of the first non-degenerate fan triangle, which is wrong for concave polygons (a dart gave the opposite sign). Use the exact orientation at the lexicographically smallest vertex, which lies on the convex hull, skipping duplicated neighbours, and keep the fan loop as a fallback (now using index_t so polygons with more than 255 vertices are supported). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The exact 3D point/triangle position used determinants with the fixed (1,1,1) vector, which are all zero when the triangle plane contains that direction, giving wrong results (inside points reported as outside). After the exact coplanarity test, drop the coordinate where the triangle normal is the largest (exact) and use the exact 2D predicate; exactly collinear triangles are handled by testing the point against their edges (keeping the 'parallel' convention). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…plane normal segment_segment_intersection_detection in 3D (and the degenerate triangle fallback of segment/triangle detection) dropped the axis with the smallest bounding box extent. For primitives in a plane containing that axis (e.g. a vertical plane) the projection collapsed the plane to a line and crossing segments were reported as parallel or disjoint ones as overlapping. Drop the largest component of the plane normal computed from the longest segment, or, for collinear inputs, the smallest component of its direction. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Heron's formula with an epsilon on p - l_i returned 0 for valid needle triangles (a unit base with a 1e-3 height gave 0, since p - l_i is about h^2/L). Compute the area from the cross product of the two edges adjacent to the vertex opposite to the longest edge, and return 0 only when the triangle height is below GLOBAL_EPSILON (compared on squared quantities). triangle_signed_area in 3D reuses the same cross product for the sign. Re-enable the existing triangle area tests. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
tetrahedron_signed_volume returned 0 as soon as no vertex had two edges making an angle above GLOBAL_ANGULAR_EPSILON, so valid flat tetrahedra (e.g. 2000 long and 0.07 high) got a null volume. Use the plain triple product (same sign convention) and return 0 only when the tetrahedron height relative to its largest facet is below GLOBAL_EPSILON (compared on squared values). Re-enable the existing signed volume test. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…ance Only the facet with the most negative barycentric coordinate was tested, which is not always the closest one for skewed tetrahedra. Return the minimal distance over all facets with a negative coordinate (all facets when the tetrahedron is degenerated). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Several deviations from Eberly's point/hyperellipsoid algorithm made point_ellipse_distance return points off the ellipse for most queries: - the bisection numerator used p * z instead of p^2 * z; - the axes were sorted by index instead of by decreasing extent; - the closest point coordinates were not un-permuted and the reflection signs were applied to the wrong axes; - zeroed coordinates were written at the wrong index. Also compute the axis lengths and unit directions only once. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
… segments Nearly parallel segments bailed out of Eberly's algorithm (log2 > 20 tests) into a dichotomic search returning the same point on both segments with a distance inconsistent with them, and throwing for segments shorter than GLOBAL_EPSILON. The post-processing could also return points that did not match the distance. - compute ac - bb, be - cd and ae - bd with cross products (Binet-Cauchy identity) to avoid cancellations, and use the parallel branch only when sin^2(angle) is below 1e-20; - handle segments shorter than GLOBAL_EPSILON as points; - polish the result once (project P(s) on segment1 then Q(t) on segment0) and compute the distance from the returned points. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Each 'from' direction was greedily matched with its most aligned 'to' direction, so two directions could map to the same one (e.g. on ties or with non-orthogonal frames), and the arrays stayed uninitialized if no dot product was positive. Choose the permutation maximizing the sum of |cos| (first one in lexicographic order on ties). Also remove unused members of NormalFrameTransform::Impl. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Compute the edges, their squared lengths and the cross product of the edges adjacent to the vertex opposite to the longest edge on scalar coordinates with branchless selections, avoiding temporary vectors and store-forwarding stalls of the previous implementation. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Cpp-Linter Report
|
Choosing the edge pair adjacent to the vertex opposite the longest edge only improves precision below GLOBAL_EPSILON, and its data-dependent selects compiled to unpredictable branches: triangle_area was 1.8x slower than the Heron version on random triangles. Use the edges at vertex 0. Benchmark vs next (ns/call): 3D random 11.6 -> 9.3, 2D random 10.5 -> 5.6. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…nto fix/geometry-robustness-tier1
panquez
left a comment
panquez
left a comment
There was a problem hiding this comment.
I think the most important thing is to think about making comparison with epsilon.
When the corner at the lexicographically smallest vertex was flat (spike), polygon_area_sign fell back to the first non-null fan triangle, which is wrong for concave polygons. Such a vertex does not change the polygon signed area: remove it and look at the new smallest vertex. Document the result for self-intersecting polygons. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Replace the dedicated 1e-20 relative threshold on sin^2(angle) by a geometric one: segments are parallel when each one deviates from the direction of the other by less than GLOBAL_EPSILON (det < GLOBAL_EPSILON^2 * min(a, c)). Segment distance tests now check results up to GLOBAL_EPSILON. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Replace the hard-coded 1e-9 / 1e-12 tolerances by GLOBAL_EPSILON, scaled by the longest edge for areas: below GLOBAL_EPSILON, values are considered equal. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
- rename the needle triangle to cap (one angle close to 180 degrees) and test a segment crossing it in its plane; - test the polygon normal on a flat polygon not aligned with the axes; - test segments crossing a triangle close to their first vertex. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* fix(Intersection): handle degenerated segments as points
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>
* fix(Distance): skip degenerated edges in triangle_triangle_distance
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>
* perf(Distance): test each triangle edge once in triangle_triangle_distance
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>
* fix(Distance): avoid out-of-bounds access on NaN in point_triangle_distance
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>
* fix(Distance): handle flat triangles in 2D point_triangle_distance
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>
* fix(Triangle): base is_degenerated on the minimum height
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>
* fix(Polygon): use the farthest vertex in minimum_height
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>
* fix(NNSearch): handle empty point sets and zero neighbor requests
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>
* fix(SquareMatrix): implement 1x1 determinant and inverse
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>
* fix(Rotation): avoid infinite loop on huge angles in rotate
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>
* fix(Intersection): fix line_cylinder_intersection perpendicular case
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>
* fix(Intersection): keep correctness in sync with filtered results
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>
* fix(Intersection): stable discriminant in line_sphere_intersection
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>
* fix(Circle): fix bounding_box for normals opposite to an axis
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>
* fix(BoundingBox): scale-independent triangle epsilon_intersects
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>
* fix(BoundingBox): fix 1D ray and line intersections
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>
* fix(Quality): scale-invariant degeneracy in tetrahedron qualities
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>
* fix(IsPointInside): exact predicates and consistent boundary handling
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>
* fix(Distance): compute point-line/segment distances by projection
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>
* fix(Point): combine coordinates hashes with absl::Hash
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>
* fix(Predicates): statically initialize expansion constants
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>
* fix(Predicates): disable FP contraction for exact predicates
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>
* test(Geometry): add tier2-robustness benchmark
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* refactor(Distance): reuse Triangle::is_degenerated in 2D point_triangle_distance
is_triangle_flat duplicated the Triangle2D::is_degenerated formula.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* fix(Intersection): report the closest cylinder point for degenerated 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>
* fix(Intersection): length tolerances in line_cylinder_intersection
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>
* fix(Rotation): throw on non-finite angle in rotate
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* perf(IsPointInside): run exact orientation only on edges spanning the point row
20-vertex star: 157 ns -> 74 ns per call.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* test(Geometry): remove tier2-robustness benchmark
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* fix(Tetrahedron): include opposite edges in is_degenerated
A tetrahedron is degenerated if its vertices lie in a slab of thickness
GLOBAL_EPSILON. The thinnest slab may be bounded by planes containing two
opposite edges instead of a facet: slivers of thickness below
GLOBAL_EPSILON have heights up to twice as large and were not reported.
Document the strip/slab definitions of Triangle and Tetrahedron.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* perf(Geometry): add zero-filter and missing filters to exact predicates
Exactly degenerate inputs (duplicated points, points sharing coordinates)
always fell back to expansion arithmetic since the floating-point filters
cannot certify a zero sign. A zero-filter (Bartels et al. 2022) now
decides them before the exact computation.
aligned_3d had no filter and dot_2d had none: both always used exact
arithmetic. aligned_3d now concludes as soon as one cross product
component is certified non-zero; dot_2d gets the orient_2d error bound.
Speedups: aligned_3d ~350 -> ~10 ns, dot_2d ~110 -> ~10 ns,
orient_3d on exactly coplanar points ~380 -> ~25 ns. Generic inputs
unchanged, results identical.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* fix(Geometry): exact predicates no longer depend on library initialization
The splitter and epsilon used by the expansion arithmetic were only set
by GEO::PCK::initialize(), called by OpenGeodeGeometryLibrary::initialize().
Code using the geometry library without initializing it (e.g. geometry
tests, which only initialize the basic library) silently got wrong exact
signs: two_product() splits with a zero splitter. They are now constant
initialized with their IEEE 754 double values.
The new consistency test (orient_2d/orient_3d under point permutations on
nearly degenerate inputs) failed before this fix.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* perf(Geometry): add Shewchuk adaptive stages to orient_2d and orient_3d
When the floating-point filter is uncertain, stages B and C of
Shewchuk's adaptive predicates (public domain predicates.c) are now
evaluated before the full exact computation: stage B computes exactly
the determinant of the rounded coordinate differences with short fixed
size expansions, stage C concludes if these differences were exact or
adds the first order corrections of their roundoff errors.
Nearly degenerate and exactly degenerate (without shared coordinates)
configurations are 2.5x to 4x faster, e.g. orient_3d on nearly coplanar
points 660 -> 270 ns, orient_2d on nearly collinear points 128 -> 41 ns.
Configurations still undecided after stage C pay the stages overhead
(~10% slower) as stage D is not ported. Generic inputs unchanged; the
signs of all benchmark calls are identical.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* perf(Geometry): add Shewchuk stage D to orient_2d and orient_3d
When stages B and C cannot decide, stage D now adds exactly all the
remaining roundoff error terms to the stage B expansion, instead of
restarting a full exact evaluation. orient_2d_exact and orient_3d_exact
are no longer used and removed.
This removes the slowdown of configurations undecided after stage C:
orient_3d on nearly collinear points 617 -> 473 ns (574 ns before the
adaptive stages). Other timings unchanged. Signs checked against
Shewchuk's original predicates.c on 20M hard cases (mixed magnitudes,
exact and nearly degenerate configurations): no mismatch.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* test(Geometry): exhaustive non-regression tests for exact predicates
test-predicates now checks every predicate against references that do not
depend on the implementation:
- integer configurations (random, exactly degenerate, one unit away from
degenerate, shared coordinates, duplicated points), at usual scales and
at scales below the floating-point filters underflow guard, compared
with exact 64 bits integer determinants;
- a table of configurations decided by each evaluation stage (filter,
zero-filter, adaptive stages B, C and D, exact evaluation), including
configurations where plain floating-point evaluation gives a wrong sign,
with signs checked with exact rational arithmetic.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* fix(Geometry): adaptive predicates return zero on a zero estimate
sign_of() returned negative for a zero estimate, reachable when all the
products underflow. Also share the double epsilon constant, return a
bool from the aligned_3d filter and document the array expansion
helpers.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* fix(Geometry): export exact predicates for Windows test linking
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>
* refactor(Geometry): move exact predicates out of internal
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>
* refactor(Geometry): move exact predicates to detail
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>
* refactor(Distance): iterate edge vertex pairs with geode::Range
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* refactor(IsPointInside): check collinearity inside is_point_on_edge
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>
* test(Geometry): check is_degenerated against vertex order
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>
* test(Geometry): express tolerances with GLOBAL_EPSILON
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
* test(Intersection): split degenerated segment tests by object
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>
* Replace internal predicates include with detail predicates
* Refactor double initialization in predicates.cpp
Replaced static double initialization with constexpr for DOUBLE_EPSILON and updated expansion_splitter and expansion_epsilon to use these constants.
* Apply prepare changes
* Add test_set_point function to test-triangle.cpp
* Apply prepare changes
* Add test for point skewed tetrahedron distance
* Apply prepare changes
* fix
---------
Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
The merge from next moved predicates.cpp from internal/ to detail/, but SKIP_UNITY and the -ffp-contract=off source property still named the old path, so both silently applied to nothing. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
|
🎉 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 1 of the geometry robustness review: fixes for functions that returned wrong results on ordinary (non-degenerate) inputs. Each fix is its own commit and comes with a test that fails on
nextand passes here (checked for every commit).Fixes (one commit each)
nextCoordinateSystemset_directions/set_directionkept a stale cached inverse matrixFramemax/min_elongation_directioncompared vectors lexicographicallyBoundingBox::add_boxAABBTree::closest_element_boxTriangle::set_pointline_line_intersectionfabs(x) < 0.was never true, so parallel lines gave (-inf, NaN)== 0.point_circle_signed_distancesegment_triangle_intersectionPolygon::normalpolygon_area_signorient_2dat a convex-hull vertexpoint_triangle_position3Dparallelfor planes containing (1,1,1) (e.g. internal faces of a cube cut into tets)segment_segment_intersection_detection3Dparalleltriangle_areap - l_i: needle triangles of height 1000·eps returned 0tetrahedron_signed_volumepoint_tetrahedron_distancepoint_ellipse_distancesegment_segment_distanceFrameTransformPlus two
perf(Mensuration)commits that make thetriangle_areacross product cheaper than Heron, and aLOCALbenchmark (tests/geometry/benchmark-tier1-robustness.cpp).Behavior changes to check downstream
triangle_areaandtetrahedron_signed_volumenow return non-zero values for thin (but not degenerate) elements. Code usingarea == 0as a degeneracy test (Background, Common, Simplex, Conversion, Geosciences) will flag fewer elements.segment_segment_distanceandpoint_tetrahedron_distancenow return correct (sometimes smaller) distances, with consistent closest points.polygon_area_signis correct for concave polygons: Inspector'ssurface_negative_elementsresults may change.line_line_intersection/segment_segment_intersectionreportparallelinstead ofincorrect. All callers found only testoperator bool.CoordinateSystem::set_directionsnow throws on degenerate directions.Deliberately not changed
FrameTransform::apply(Point/Vector)gathers (result[d] = p[dir[d]]) whileapply(Frame)scatters, so the two disagree for 3-cycle permutations. Switching to scatter would silently break Simplex frontal metric coordinates and Hybrid/Numerics PGP, which rely on the gather convention. This needs a design decision, so it is left as is (the test documents the current convention).Benchmark (ns/call, best of 7, pinned to one core,
next670d96d vs this branch, same clang -O3 flags)Significant speedups:
segment_segment_distance(2.3x random, 23x near-parallel),Polygon::normal(14x),triangulate(2.1x),triangle_area2D (1.9x) and 3D (1.25x),triangle_signed_area(1.7x).The slowdowns pay for correctness:
set_directions+coordinatesandset_point+normal:nextreturned a stale cached value (the bug).segment_triangle_intersectionwith an endpoint on the triangle:nextwrongly rejected most of these segments early.triangle_area:nextreturned 0 early.point_tetrahedron_distance: every visible facet is now evaluated.point_ellipse_distance: now converges to the right root.polygon_area_sign: about 3 ns for the minimum-vertex scan.Everything else is within noise.
Full table
Testing
ctest -L essential: 102/102 (newtest-frame).🤖 Generated with Claude Code