Improve triangulation robustness by micycle1 · Pull Request #1212 · locationtech/jts

micycle1 · 2026-07-11T08:37:30Z

isInCircleRobust now computes the in-circle determinant in doubles with a Shewchuk-style error bound: if the magnitude clears the bound, the sign is provably correct and returned at full speed; otherwise it falls back to the existing isInCircleDDNormalized extended-precision evaluation.

Vertex.isCCW now uses Orientation.index.

~12% slower on a 200k random-point triangulation.

Closes #20, #298, #310, #1190 (and therefore #1138), GEOS #1286, and supersedes PR #311.

grootstebozewolf · 2026-08-05T01:01:49Z

Thanks for this — differential check from the NetTopologySuite.Proofs side (mesh / Delaunay lane).

Method

Oracle: extracted b64_inCircle via INCIRCLE_SIGN (same kernel as in-process nts_rocq_in_circle; Shewchuk sign: POS ⇔ CCW triangle and P strictly inside).
Stage A mirror: pure double of this PR’s isInCircleRobust filter (iccerrboundA = (10 + 96·ε)·ε, permanent-scaled; CERTAIN iff |disc| > errbound).
We did not re-run the full Maven suite or claim bit-identity of the DD fallback with the oracle.

Vectors / generator: jts-1212-incircle-lane (table in oracle/jts1212_incircle_vectors.txt).

Results (29 vectors)

Check	Result
Certified pin failures (EXPECTED set)	0
Stage A CERTAIN vs oracle sign conflicts	0
Stage A UNCERTAIN (correct declines)	11

Pins that should hit the DD path (oracle ZERO, Stage A UNCERTAIN): vertex-on-circumcircle (incl. your #1190 constrained triangle sites) and exact cocircular square diagonals. That is expected: |disc| = 0 never clears a positive permanent bound.

Decisive non-zero pins (CERTAIN match): knife-edge ±, flip witness A=(0,0),B=(2,0),C=(1,1),P=(1,-½) → POS, outside-disk witness → NEG, scaled flip (homothety) → POS.

Your regression sites (sample)

Case	Oracle	Stage A
#1190 triangle constraint vertices	ZERO	UNCERTAIN
`testCircle` sample quads	NEG	CERTAIN
GEOS#1040 / JTS#1171 near-cocircular quad	NEG	CERTAIN
JTS#1171 ex.2 (GEOS 955 subset)	POS	ZERO · UNCERTAIN
JTS#20 near-coincident (first 4 WKB sites)	NEG	CERTAIN

The GEOS 955 subset is the useful illustration that Stage A alone is not enough:

A = (18.68285714285716, 100.105)
B = (13.41, 104.82100000000001)
C = (13.41, 107.179)
P = (18.682857142857145, 111.89500000000001)

Stage A double: disc classifies as ZERO, filter UNCERTAIN
Oracle b64_inCircle: POS (strictly inside)

A bare disc > 0 without the uncertain→DD branch would return false here; the fallback path is doing real work. That matches the design of this PR.

Bottom line

Filter shape looks right relative to an independent extracted in-circle kernel: CERTAIN doubles never disagreed with the oracle on this table, and the cases that must use DD (exact zeros + the 955 subset) are correctly marked UNCERTAIN. Happy to re-run if you add more quads to the test suite.

(Affiliation: NetTopologySuite.Proofs — formal mesh/predicate companion, not a JTS committer review.)

grootstebozewolf

LGTM — approve with nits.

The filter is the right Shewchuk Stage A, the fallback is the matching DD form, and Vertex.isCCW now uses the production robust orient. That is the #310 roadmap done as the default (~12% on 200k random points is the cost Martin asked to measure). Vertex.isInCircle and TriDelaunayImprover already call isInCircleRobust.

Do not treat this as a close of #1190. Martin already said a triangulation-as-constraints is outside CDT’s contract; the zero-area face on that one triangle is worth pinning, the extra Steiner splits are not a predicate bug.

Close: #20, #298, #310, GEOS #1286, and the summary #1138 (it lists those three, not #1190). Keep #1190 open, or close it only for the zero-area symptom.

Nits: DOUBLE_EPS javadoc (unit roundoff, not “largest value with 1.0 + ε == 1.0”); testCocircularSitesJTS1171 example 1 duplicates GEOS1040; no direct TrianglePredicate unit test (optional); neither commit has Signed-off-by.

grootstebozewolf · 2026-09-03T04:24:21Z

  /**
-   * Tests if a point is inside the circle defined by 
-   * the triangle with vertices a, b, c (oriented counter-clockwise). 
-   * This method uses more robust computation.
-   * 
+   * The unit roundoff for IEEE 754 double precision (round-to-nearest-even):
+   * the largest value for which {@code 1.0 + DOUBLE_EPS == 1.0}.
+   * This is Shewchuk's <tt>epsilon</tt>, computed by the <tt>exactinit()</tt>
+   * bisection in the original (C) predicates code.
+   */
+  private static final double DOUBLE_EPS = 0x1.0p-53;
+


0x1.0p-53 is Shewchuk’s unit roundoff. It is not “the largest value for which 1.0 + DOUBLE_EPS == 1.0” — larger values still round to 1. Say “Shewchuk epsilon / unit roundoff, 2^-53”.

grootstebozewolf · 2026-09-03T04:24:23Z

  public static boolean isInCircleRobust(
-      Coordinate a, Coordinate b, Coordinate c, 
-      Coordinate p) 
+      Coordinate a, Coordinate b, Coordinate c,
+      Coordinate p)
  {
-    //checkRobustInCircle(a, b, c, p);
-//    return isInCircleNonRobust(a, b, c, p);       
-    return isInCircleNormalized(a, b, c, p);       
+    double adx = a.x - p.x;
+    double ady = a.y - p.y;
+    double bdx = b.x - p.x;
+    double bdy = b.y - p.y;
+    double cdx = c.x - p.x;
+    double cdy = c.y - p.y;
+
+    double bdxcdy = bdx * cdy;
+    double cdxbdy = cdx * bdy;
+    double alift = adx * adx + ady * ady;
+
+    double cdxady = cdx * ady;
+    double adxcdy = adx * cdy;
+    double blift = bdx * bdx + bdy * bdy;
+
+    double adxbdy = adx * bdy;
+    double bdxady = bdx * ady;
+    double clift = cdx * cdx + cdy * cdy;
+
+    double disc = alift * (bdxcdy - cdxbdy)
+                + blift * (cdxady - adxcdy)
+                + clift * (adxbdy - bdxady);
+
+    double permanent = (Math.abs(bdxcdy) + Math.abs(cdxbdy)) * alift
+                     + (Math.abs(cdxady) + Math.abs(adxcdy)) * blift
+                     + (Math.abs(adxbdy) + Math.abs(bdxady)) * clift;
+    double errBound = IN_CIRCLE_ERR_BOUND * permanent;
+
+    if (disc > errBound || -disc > errBound) {
+      return disc > 0;
+    }
+    //-- result is uncertain, so evaluate with extended precision
+    return isInCircleDDNormalized(a, b, c, p);
  }



This is the right Shewchuk Stage A. Algebraically disc is the same as isInCircleNormalized (alift·bcdet + blift·cadet + clift·abdet); the extra products exist so the permanent can be formed. iccerrboundA = (10 + 96ε)ε and disc > err || -disc > err match predicates.c.

Uncertain cases go to isInCircleDDNormalized, which is the same expression in DD — better than #310’s isInCircleDDFast (unnormalized 4-term) as a fallback for this filter. Vertex.isInCircle and TriDelaunayImprover already call this method, so incremental DT and polygon improve pick it up.

An earlier Proofs-side check on this PR (29 quads) had 0 Stage A CERTAIN sign conflicts with an independent b64_inCircle oracle. The GEOS 955 subset is UNCERTAIN with disc classifying as 0 while the oracle is POS — that is why the DD branch is not decorative.

Optional: one direct TrianglePredicate test that a CERTAIN outside-disk is false without needing DT, and that the 955-subset quad is UNCERTAIN and the DD path is true. The suite is integration-only today.

Do not add an IncrementalDelaunayTriangulator “fast unsafe” flag. Making this the default is the #310 outcome.

grootstebozewolf · 2026-09-03T04:24:24Z

+  public final boolean isCCW(Vertex b, Vertex c)
  {
-      /*
-      // test code used to check for robustness of triArea 
-      boolean isCCW = (b.p.x - p.x) * (c.p.y - p.y) 
-      - (b.p.y - p.y) * (c.p.x - p.x) > 0;
-     //boolean isCCW = triArea(this, b, c) > 0;
-     boolean isCCWRobust = CGAlgorithms.orientationIndex(p, b.p, c.p) == CGAlgorithms.COUNTERCLOCKWISE; 
-     if (isCCWRobust != isCCW)
-      System.out.println("CCW failure");
-     //*/
-
-    	// is equal to the signed area of the triangle
-
-      return (b.p.x - p.x) * (c.p.y - p.y) 
-           - (b.p.y - p.y) * (c.p.x - p.x) > 0;
-
-      // original rolled code
-      //boolean isCCW = triArea(this, b, c) > 0;
-      //return isCCW;
-
+      //-- robust orientation test, to avoid triangulation errors
+      //-- caused by predicate failure for nearly-collinear points
+      return Orientation.index(p, b.p, c.p) == Orientation.COUNTERCLOCKWISE;
    }



Right half of the robustness story. Orientation.index → CGAlgorithmsDD.orientationIndex. Collinear stays not-CCW (== COUNTERCLOCKWISE), same as the old > 0. The commented debug that used to live here already knew this was the check that mattered.

grootstebozewolf · 2026-09-03T04:24:25Z

+  /**
+   * Constraints forming a triangle whose vertices are the sites.
+   * The obtuse angle causes the long edge to be split (by design),
+   * but the result must not contain degenerate (zero-area) triangles.
+   *
+   * see https://github.com/locationtech/jts/issues/1190
+   */
+  public void testTriangleConstraints_JTS_1190()
+  throws ParseException
+  {
+    String wkt = "MULTIPOINT ((-221.72957795130824 -26.56505117707799), (-149.72957795130824 -26.56505117707799), (0 -90))";
+    String lineWKT = "MULTILINESTRING ((-221.72957795130824 -26.56505117707799, -149.72957795130824 -26.56505117707799), (-149.72957795130824 -26.56505117707799, 0 -90), (0 -90, -221.72957795130824 -26.56505117707799))";
+    Geometry sites = reader.read(wkt);
+    Geometry constraints = reader.read(lineWKT);
+
+    ConformingDelaunayTriangulationBuilder builder = new ConformingDelaunayTriangulationBuilder();
+    builder.setSites(sites);
+    builder.setConstraints(constraints);
+    Geometry result = builder.getTriangles(new GeometryFactory());
+
+    assertTrue(result.getNumGeometries() > 0);
+    for (int i = 0; i < result.getNumGeometries(); i++) {
+      assertTrue("Zero-area triangle in conforming Delaunay triangulation",
+          result.getGeometryN(i).getArea() > 0);
+    }
+  }
+


Good pin for the zero-area face on the reported triangle. Do not treat this as a close of #1190. Martin already said a triangulation-as-constraints is outside CDT’s contract; the extra Steiner splits are not a predicate bug. Close #20 / #298 / #310 / GEOS #1286 / the summary #1138 (it lists those three, not #1190). Keep #1190 open, or close it only for the zero-area symptom and say so.

grootstebozewolf · 2026-09-03T04:24:26Z

+  public void testCocircularSitesJTS1171() {
+    //-- Example 1 (same sites as GEOS 1040)
+    checkVoronoiValid("MULTIPOINT ((6.6584 53.583000000000006), (6.6576 53.583600000000004), (6.657 53.5848), (6.6572000000000005 53.5842))");
+    //-- Example 2 (subset of the GEOS 955 sites)
+    checkVoronoiValid("MULTIPOINT ((18.68285714285716 100.105), (13.41 104.82100000000001), (13.41 107.179), (18.682857142857145 111.89500000000001))");
+  }


Example 1 is the same WKT as testCocircularSitesGEOS1040. Drop the duplicate; keep example 2 (the 955 subset — that is the quad where Stage A is UNCERTAIN and DD does the work).

grootstebozewolf · 2026-09-03T04:24:27Z

+  /**
+   * The input points are nearly co-circular, so the Delaunay triangulation
+   * is not unique - the expected result is the one chosen by the robust predicates.
+   */
  public void testCircle()
  {
    String wkt = "POLYGON ((42 30, 41.96 29.61, 41.85 29.23, 41.66 28.89, 41.41 28.59, 41.11 28.34, 40.77 28.15, 40.39 28.04, 40 28, 39.61 28.04, 39.23 28.15, 38.89 28.34, 38.59 28.59, 38.34 28.89, 38.15 29.23, 38.04 29.61, 38 30, 38.04 30.39, 38.15 30.77, 38.34 31.11, 38.59 31.41, 38.89 31.66, 39.23 31.85, 39.61 31.96, 40 32, 40.39 31.96, 40.77 31.85, 41.11 31.66, 41.41 31.41, 41.66 31.11, 41.85 30.77, 41.96 30.39, 42 30))";
-    String expected = "MULTILINESTRING ((41.66 31.11, 41.85 30.77), (41.41 31.41, 41.66 31.11), (41.11 31.66, 41.41 31.41), (40.77 31.85, 41.11 31.66), (40.39 31.96, 40.77 31.85), (40 32, 40.39 31.96), (39.61 31.96, 40 32), (39.23 31.85, 39.61 31.96), (38.89 31.66, 39.23 31.85), (38.59 31.41, 38.89 31.66), (38.34 31.11, 38.59 31.41), (38.15 30.77, 38.34 31.11), (38.04 30.39, 38.15 30.77), (38 30, 38.04 30.39), (38 30, 38.04 29.61), (38.04 29.61, 38.15 29.23), (38.15 29.23, 38.34 28.89), (38.34 28.89, 38.59 28.59), (38.59 28.59, 38.89 28.34), (38.89 28.34, 39.23 28.15), (39.23 28.15, 39.61 28.04), (39.61 28.04, 40 28), (40 28, 40.39 28.04), (40.39 28.04, 40.77 28.15), (40.77 28.15, 41.11 28.34), (41.11 28.34, 41.41 28.59), (41.41 28.59, 41.66 28.89), (41.66 28.89, 41.85 29.23), (41.85 29.23, 41.96 29.61), (41.96 29.61, 42 30), (41.96 30.39, 42 30), (41.85 30.77, 41.96 30.39), (41.66 31.11, 41.96 30.39), (41.41 31.41, 41.96 30.39), (41.41 28.59, 41.96 30.39), (41.41 28.59, 41.41 31.41), (38.59 28.59, 41.41 28.59), (38.59 28.59, 41.41 31.41), (38.59 28.59, 38.59 31.41), (38.59 31.41, 41.41 31.41), (38.59 31.41, 39.61 31.96), (39.61 31.96, 41.41 31.41), (39.61 31.96, 40.39 31.96), (40.39 31.96, 41.41 31.41), (40.39 31.96, 41.11 31.66), (38.04 30.39, 38.59 28.59), (38.04 30.39, 38.59 31.41), (38.04 30.39, 38.34 31.11), (38.04 29.61, 38.59 28.59), (38.04 29.61, 38.04 30.39), (39.61 28.04, 41.41 28.59), (38.59 28.59, 39.61 28.04), (38.89 28.34, 39.61 28.04), (40.39 28.04, 41.41 28.59), (39.61 28.04, 40.39 28.04), (41.96 29.61, 41.96 30.39), (41.41 28.59, 41.96 29.61), (41.66 28.89, 41.96 29.61), (40.39 28.04, 41.11 28.34), (38.04 29.61, 38.34 28.89), (38.89 31.66, 39.61 31.96))";
+    String expected = "MULTILINESTRING ((41.96 30.39, 42 30), (41.85 30.77, 41.96 30.39), (41.66 31.11, 41.85 30.77), (41.41 31.41, 41.66 31.11), (41.11 31.66, 41.41 31.41), (40.77 31.85, 41.11 31.66), (40.39 31.96, 40.77 31.85), (40 32, 40.39 31.96), (39.61 31.96, 40 32), (39.23 31.85, 39.61 31.96), (38.89 31.66, 39.23 31.85), (38.59 31.41, 38.89 31.66), (38.34 31.11, 38.59 31.41), (38.15 30.77, 38.34 31.11), (38.04 30.39, 38.15 30.77), (38 30, 38.04 30.39), (38 30, 38.04 29.61), (38.04 29.61, 38.15 29.23), (38.15 29.23, 38.34 28.89), (38.34 28.89, 38.59 28.59), (38.59 28.59, 38.89 28.34), (38.89 28.34, 39.23 28.15), (39.23 28.15, 39.61 28.04), (39.61 28.04, 40 28), (40 28, 40.39 28.04), (40.39 28.04, 40.77 28.15), (40.77 28.15, 41.11 28.34), (41.11 28.34, 41.41 28.59), (41.41 28.59, 41.66 28.89), (41.66 28.89, 41.85 29.23), (41.85 29.23, 41.96 29.61), (41.96 29.61, 42 30), (41.96 29.61, 41.96 30.39), (41.41 28.59, 41.96 29.61), (41.41 28.59, 41.96 30.39), (41.41 28.59, 41.41 31.41), (41.41 31.41, 41.96 30.39), (38.59 28.59, 41.41 28.59), (38.59 28.59, 41.41 31.41), (38.59 28.59, 38.59 31.41), (38.59 31.41, 41.41 31.41), (38.59 31.41, 39.61 31.96), (39.61 31.96, 41.41 31.41), (39.61 31.96, 40.39 31.96), (40.39 31.96, 41.41 31.41), (40.39 31.96, 41.11 31.66), (38.04 29.61, 38.59 28.59), (38.04 29.61, 38.59 31.41), (38.04 29.61, 38.04 30.39), (38.04 30.39, 38.59 31.41), (38.04 30.39, 38.34 31.11), (39.61 28.04, 41.41 28.59), (38.59 28.59, 39.61 28.04), (38.89 28.34, 39.61 28.04), (40.39 28.04, 41.41 28.59), (39.61 28.04, 40.39 28.04), (41.66 28.89, 41.96 29.61), (40.39 28.04, 41.11 28.34), (38.04 29.61, 38.34 28.89), (38.89 31.66, 39.61 31.96), (41.66 31.11, 41.96 30.39))";
    checkDelaunayEdges(wkt, expected);
  }


testCircle’s new edge set is documented as a non-unique Delaunay. That is honest. The #298 / GEOS #1286 pins below are the right integration locks.

grootstebozewolf · 2026-09-03T04:24:29Z

Neither commit has Signed-off-by. Eclipse DCO wants it on every commit. This completes the idea in #311.

micycle1 added 2 commits July 11, 2026 09:03

Improve triangulation robustness

3c258d6

clarify IN_CIRCLE_ERR_BOUND

012465a

This was referenced Jul 11, 2026

VoronoiDiagramBuilder: tolerance-based dedup + safe cell construction (stacked on #1212) #1213

Closed

VoronoiDiagramBuilder: tolerance-based dedup + safe cell construction (stacked on #1212) micycle1/jts#3

Open

grootstebozewolf mentioned this pull request Aug 5, 2026

docs: scout open JTS PRs + #1212 in-circle differential lane grootstebozewolf/NetTopologySuite.Proofs#444

Merged

2 tasks

This was referenced Aug 5, 2026

docs: JTS #1094 in-circle FFI scout + corpus code review grootstebozewolf/NetTopologySuite.Proofs#445

Merged

Make isInCircleRobust robust #1094

Open

isInCircleAdapt sketch Java port #311

Open

grootstebozewolf reviewed Sep 3, 2026

View reviewed changes

Navigation Menu

Uh oh!

Improve triangulation robustness - #1212

Improve triangulation robustness#1212
micycle1 wants to merge 2 commits into
locationtech:masterfrom
micycle1:fixup

micycle1 commented Jul 11, 2026 •

edited

Loading

Uh oh!

grootstebozewolf commented Aug 5, 2026

Uh oh!

grootstebozewolf left a comment

Uh oh!

grootstebozewolf Sep 3, 2026

Uh oh!

grootstebozewolf Sep 3, 2026

Uh oh!

grootstebozewolf Sep 3, 2026

Uh oh!

grootstebozewolf Sep 3, 2026

Uh oh!

grootstebozewolf Sep 3, 2026

Uh oh!

grootstebozewolf Sep 3, 2026

Uh oh!

grootstebozewolf Sep 3, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Uh oh!

Conversation

micycle1 commented Jul 11, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

grootstebozewolf commented Aug 5, 2026

Method

Results (29 vectors)

Your regression sites (sample)

Bottom line

Uh oh!

grootstebozewolf left a comment

Choose a reason for hiding this comment

Uh oh!

grootstebozewolf Sep 3, 2026

Choose a reason for hiding this comment

Uh oh!

grootstebozewolf Sep 3, 2026

Choose a reason for hiding this comment

Uh oh!

grootstebozewolf Sep 3, 2026

Choose a reason for hiding this comment

Uh oh!

grootstebozewolf Sep 3, 2026

Choose a reason for hiding this comment

Uh oh!

grootstebozewolf Sep 3, 2026

Choose a reason for hiding this comment

Uh oh!

grootstebozewolf Sep 3, 2026

Choose a reason for hiding this comment

Uh oh!

grootstebozewolf Sep 3, 2026

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

micycle1 commented Jul 11, 2026 •

edited

Loading