Robust geometry: epsilons, snapping and exact predicates
Why geometry code that is right on paper fails on real models: rounding, tolerances that are not transitive, exact orientation tests, and welding and snapping vertices.
Doubles have gaps
double a = 0.1 + 0.2;Console.WriteLine($"0.1 + 0.2 == 0.3: {a == 0.3}, off by {a - 0.3:0.0E0}"); double x = 0;foreach (double step in new[] { 0.1, 0.2, 0.3 }) x += step;x -= 0.6;Console.WriteLine($"0.1 + 0.2 + 0.3 - 0.6 = {x:0.0E0}, back at zero: {x == 0}"); Console.WriteLine("gap to the next number up:");foreach (double v in new[] { 1.0, 1_000, 512_340, 6_000_000 }){ double gd = Math.BitIncrement(v) - v; float f = (float)v; double gf = MathF.BitIncrement(f) - f; Console.WriteLine($" {v,9:0}: double {gd:0.0E0} m, float {gf:0.0E0} m");}0.1 + 0.2 == 0.3: False, off by 5.6E-17
0.1 + 0.2 + 0.3 - 0.6 = 1.1E-16, back at zero: False
gap to the next number up:
1: double 2.2E-16 m, float 1.2E-7 m
1000: double 1.1E-13 m, float 6.1E-5 m
512340: double 5.8E-11 m, float 3.1E-2 m
6000000: double 9.3E-10 m, float 5.0E-1 mTolerances help, but “near” is not equality
Console.WriteLine($"Near(0, 0.0008) = {Near(0, 0.0008)}");Console.WriteLine($"Near(0.0008, 0.0016) = {Near(0.0008, 0.0016)}");Console.WriteLine($"Near(0, 0.0016) = {Near(0, 0.0016)}"); // Merging points one by one: the answer depends on the order they arrive inConsole.WriteLine($"merge 0, 0.0008, 0.0016: {Merge(0, 0.0008, 0.0016)} kept");Console.WriteLine($"merge 0.0008, 0, 0.0016: {Merge(0.0008, 0, 0.0016)} kept"); static bool Near(double a, double b) => Math.Abs(a - b) <= 0.001; // 1 mm static int Merge(params double[] xs){ var kept = new List<double>(); foreach (var x in xs) if (!kept.Any(k => Near(k, x))) kept.Add(x); return kept.Count;}Near(0, 0.0008) = True Near(0.0008, 0.0016) = True Near(0, 0.0016) = False merge 0, 0.0008, 0.0016: 2 kept merge 0.0008, 0, 0.0016: 1 kept
The orientation test, done properly
using System.Numerics; const double H = 1.0 / (1L << 53); // 2^-53, the gap between doubles near 0.5double qx = 17.3, qy = 0.5 + 0.7 * 16.8 + 20 * H, rx = 30.1, ry = 0.5 + 0.7 * 29.6 + 20 * H; int wrong = 0, zero = 0, fallbacks = 0, adaptiveWrong = 0;for (int j = 0; j < 40; j++) for (int i = 0; i < 40; i++) { double px = 0.5 + i * H, py = 0.5 + j * H; int truth = OrientExact(px, py, qx, qy, rx, ry); int plain = Math.Sign(OrientPlain(px, py, qx, qy, rx, ry)); if (plain == 0 && truth != 0) zero++; else if (plain != truth) wrong++; if (OrientAdaptive(px, py, qx, qy, rx, ry, ref fallbacks) != truth) adaptiveWrong++; }Console.WriteLine($"plain doubles: {wrong} on the wrong side, {zero} wrongly on the line, of 1600");Console.WriteLine($"adaptive: {adaptiveWrong} wrong, exact arithmetic needed {fallbacks} times"); int f2 = 0;Console.WriteLine($"(3, 9) against the same line: {OrientAdaptive(3, 9, qx, qy, rx, ry, ref f2)}, exact needed {f2} times"); static double OrientPlain(double px, double py, double qx, double qy, double rx, double ry) => (qx - px) * (ry - py) - (qy - py) * (rx - px); // Fast double test first; only when the result is smaller than its possible error, do it exactly.static int OrientAdaptive(double px, double py, double qx, double qy, double rx, double ry, ref int fallbacks){ double l = (qx - px) * (ry - py), r = (qy - py) * (rx - px), det = l - r; double bound = 3.3306690738754716e-16 * (Math.Abs(l) + Math.Abs(r)); // Shewchuk's error bound if (Math.Abs(det) > bound) return Math.Sign(det); fallbacks++; return OrientExact(px, py, qx, qy, rx, ry);} // A double is a whole number times a power of two, so scaling turns these into exact integers.static int OrientExact(double px, double py, double qx, double qy, double rx, double ry){ static BigInteger X(double v) => new BigInteger(v * Math.Pow(2, 60)); var d = (X(qx) - X(px)) * (X(ry) - X(py)) - (X(qy) - X(py)) * (X(rx) - X(px)); return d.Sign;}plain doubles: 277 on the wrong side, 686 wrongly on the line, of 1600 adaptive: 0 wrong, exact arithmetic needed 1600 times (3, 9) against the same line: 1, exact needed 0 times
Welding and snapping
Snapping rounds every coordinate to a grid (1 mm, or 2⁻²⁰ m). Afterwards equality is exact equality, and integer coordinates make exact arithmetic easy. The risk: snapping can create new problems, such as duplicate vertices, zero-length edges and three points that become collinear. Welding merges vertices closer than a tolerance into one shared vertex, so triangles that should share an edge actually reference the same indices. A plain “round to the grid and compare” misses pairs that straddle a grid line; a hash grid that also checks the neighbouring cells does not.
// Two triangles from different sources that should share the edge (0,0)-(4,3)var soup = new V3[]{ new(0, 0, 0), new(4, 0, 0), new(4, 3, 0), new(0.0004, 0.0002, 0), new(4.0006, 3.0001, 0), new(0, 3, 0),};Console.WriteLine($"rounding to the nearest mm: {WeldByRounding(soup, 0.001)} vertices");var (verts, map) = Weld(soup, 0.001);Console.WriteLine($"welding within 1 mm: {verts.Count} vertices");Console.WriteLine($"triangles: ({map[0]}, {map[1]}, {map[2]}) and ({map[3]}, {map[4]}, {map[5]})"); static int WeldByRounding(V3[] pts, double tol) => pts.Select(p => (Math.Round(p.X / tol), Math.Round(p.Y / tol), Math.Round(p.Z / tol))).Distinct().Count(); // Hash grid with cells one tolerance wide: a match can only be in the same cell or a neighbour.static (List<V3> Verts, int[] Map) Weld(V3[] pts, double tol){ var verts = new List<V3>(); var grid = new Dictionary<(long, long, long), List<int>>(); var map = new int[pts.Length]; for (int k = 0; k < pts.Length; k++) { var p = pts[k]; long cx = (long)Math.Floor(p.X / tol), cy = (long)Math.Floor(p.Y / tol), cz = (long)Math.Floor(p.Z / tol); int found = -1; for (long dx = -1; dx <= 1; dx++) for (long dy = -1; dy <= 1; dy++) for (long dz = -1; dz <= 1; dz++) if (found < 0 && grid.TryGetValue((cx + dx, cy + dy, cz + dz), out var list)) foreach (int i in list) if ((verts[i] - p).Length <= tol) { found = i; break; } if (found < 0) { found = verts.Count; verts.Add(p); if (!grid.TryGetValue((cx, cy, cz), out var cell)) grid[(cx, cy, cz)] = cell = new(); cell.Add(found); } map[k] = found; } return (verts, map);} readonly record struct V3(double X, double Y, double Z){ public static V3 operator +(V3 a, V3 b) => new(a.X + b.X, a.Y + b.Y, a.Z + b.Z); public static V3 operator -(V3 a, V3 b) => new(a.X - b.X, a.Y - b.Y, a.Z - b.Z); public static V3 operator *(V3 v, double k) => new(v.X * k, v.Y * k, v.Z * k); public static double Dot(V3 a, V3 b) => a.X * b.X + a.Y * b.Y + a.Z * b.Z; public static V3 Cross(V3 a, V3 b) => new(a.Y * b.Z - a.Z * b.Y, a.Z * b.X - a.X * b.Z, a.X * b.Y - a.Y * b.X); public double Length => Math.Sqrt(Dot(this, this)); public V3 Normalized() => this * (1 / Length); static string F(double v) => Math.Abs(v) < 5e-4 ? "0" : v.ToString("0.###"); public override string ToString() => $"({F(X)}, {F(Y)}, {F(Z)})";}rounding to the nearest mm: 5 vertices welding within 1 mm: 4 vertices triangles: (0, 1, 2) and (0, 2, 3)
A practical policy
Work near the origin. Subtract a project origin before doing geometry; keep large survey coordinates for input and output only. Clean the input once. Weld or snap vertices, remove zero-length edges and repeated points, and merge nearly collinear vertices when data comes in, not halfway through an algorithm. Use exact or adaptive predicates for decisions. Orientation, in-circle and segment-crossing tests should give consistent signs; quantities like lengths and areas can stay approximate. Choose one tolerance and use it everywhere. Two parts of a program with different ε values will disagree about whether two points are the same. Keep topology separate from coordinates. Decide which edges meet by index, not by comparing positions again later.
Show a hint
Show the solution
var poly = new V2[] { new(0, 0), new(5.0012, 0.0011), new(5.0018, 0.0009), new(5, 2.5004), new(4.9989, 4.0013), new(0.0021, 3.9992) };const double Grid = 0.005;var snapped = poly.Select(p => new V2(Math.Round(p.X / Grid) * Grid, Math.Round(p.Y / Grid) * Grid)); var clean = new List<V2>();foreach (var p in snapped) if (clean.Count == 0 || (p - clean[^1]).Length > 1e-9) clean.Add(p);if ((clean[0] - clean[^1]).Length <= 1e-9) clean.RemoveAt(clean.Count - 1); for (bool changed = true; changed;){ changed = false; for (int i = 0; i < clean.Count && !changed; i++) { V2 a = clean[(i + clean.Count - 1) % clean.Count], b = clean[i], c = clean[(i + 1) % clean.Count]; double cross = (b.X - a.X) * (c.Y - a.Y) - (b.Y - a.Y) * (c.X - a.X); if (Math.Abs(cross) < 1e-9) { clean.RemoveAt(i); changed = true; } }} Console.WriteLine($"before: {poly.Length} vertices, area {Area(poly):0.000000} m2");Console.WriteLine($"after: {clean.Count} vertices, area {Area(clean):0.000000} m2");Console.WriteLine(string.Join(" ", clean)); static double Area(IList<V2> p){ double a = 0; for (int i = 0; i < p.Count; i++) { var q = p[(i + 1) % p.Count]; a += p[i].X * q.Y - q.X * p[i].Y; } return a / 2;} readonly record struct V2(double X, double Y){ public static V2 operator +(V2 a, V2 b) => new(a.X + b.X, a.Y + b.Y); public static V2 operator -(V2 a, V2 b) => new(a.X - b.X, a.Y - b.Y); public static V2 operator *(V2 v, double k) => new(v.X * k, v.Y * k); public double Length => Math.Sqrt(X * X + Y * Y); static string F(double v) => Math.Abs(v) < 5e-4 ? "0" : v.ToString("0.###"); public override string ToString() => $"({F(X)}, {F(Y)})";}before: 6 vertices, area 19.995724 m2 after: 4 vertices, area 20.000000 m2 (0, 0) (5, 0) (5, 4) (0, 4)
Remember
Doubles have gaps that grow with the size of the number; geometry far from the origin loses precision. Geometric decisions (left or right, inside or outside, same point or not) can flip because of rounding, and one flipped decision breaks topology. Tolerances are not transitive, so results can depend on processing order; use one tolerance consistently. Exact and adaptive predicates give correct signs: a fast double test, with exact arithmetic only when the result is within its error bound. Clean input once by welding or snapping vertices, then decide topology by index, not by comparing coordinates again.