Prevent polygon corruption from floating point imprecision. (akavel#4)

Floating point imprecision can produce intersections that are infinitesimally close to segment endpoints. Such intersections produce almost-parallel lines that can corrupt the polygon with self-intersections.

Resistance to this corruption is achieved by snapping computed intersection points to segment endpoints within a small floating point tolerance.

bugs_test.go

@@ -237,6 +237,80 @@ func TestResweepingIntersectingEndpoints(t *T) {
 	}.verify(t)
 }
 
+func TestCorruptionResistanceFromFloatingPointImprecision(t *T) {
+	testCases{
+		{
+			op: polyclip.INTERSECTION,
+			subject: polyclip.Polygon{{
+				{2.500054958553072, 20.615428910374025},
+				{42.500054958553065, -19.38457108962598},
+				{82.50005495855308, 20.61542891037402},
+				{42.50005495855307, 60.61542891037402},
+			}},
+			clipping: polyclip.Polygon{{
+				{7.604714313123809, 25.720088264944764},
+				{11.897740920349097, 21.42706165771947},
+				{36.852886624296644, 46.382207361667014},
+				{32.55986001707135, 50.6752339688923},
+			}},
+			result: polyclip.Polygon{{
+				{7.604714313123809, 25.720088264944764},
+				{32.55986001707135, 50.6752339688923},
+				{36.852886624296644, 46.382207361667014},
+				{11.897740920349097, 21.42706165771947},
+			}},
+		},
+		{
+			op: polyclip.DIFFERENCE,
+			subject: polyclip.Polygon{{
+				{2.500054958553072, 20.615428910374025},
+				{42.500054958553065, -19.38457108962598},
+				{82.50005495855308, 20.61542891037402},
+				{42.50005495855307, 60.61542891037402},
+			}},
+			clipping: polyclip.Polygon{{
+				{7.604714313123809, 25.720088264944764},
+				{32.55986001707135, 50.6752339688923},
+				{36.852886624296644, 46.382207361667014},
+				{11.897740920349097, 21.42706165771947},
+			}},
+			result: polyclip.Polygon{{
+				{2.500054958553072, 20.615428910374025},
+				{42.500054958553065, -19.38457108962598},
+				{82.50005495855308, 20.61542891037402},
+				{42.50005495855307, 60.61542891037402},
+				{32.55986001707135, 50.6752339688923},
+				{36.852886624296644, 46.382207361667014},
+				{11.897740920349097, 21.42706165771947},
+				{7.604714313123809, 25.720088264944764},
+			}},
+		},
+		{
+			op: polyclip.UNION,
+			subject: polyclip.Polygon{{
+				{2.500054958553072, 20.615428910374025},
+				{42.500054958553065, -19.38457108962598},
+				{82.50005495855308, 20.61542891037402},
+				{42.50005495855307, 60.61542891037402},
+			}},
+			clipping: polyclip.Polygon{{
+				{7.604714313123809, 25.720088264944764},
+				{32.55986001707135, 50.6752339688923},
+				{36.852886624296644, 46.382207361667014},
+				{11.897740920349097, 21.42706165771947},
+			}},
+			result: polyclip.Polygon{{
+				{2.500054958553072, 20.615428910374025},
+				{42.500054958553065, -19.38457108962598},
+				{82.50005495855308, 20.61542891037402},
+				{42.50005495855307, 60.61542891037402},
+				{32.55986001707135, 50.6752339688923},
+				{7.604714313123809, 25.720088264944764},
+			}},
+		},
+	}.verify(t)
+}
+
 func TestSelfIntersectionAvoidance(t *T) {
 	testCases{
 		{

clipper.go

@@ -389,6 +389,18 @@ func findIntersection2(u0, u1, v0, v1 float64, w *[]float64) int {
 	return 2
 }
 
+// snaps the [pt] to one of [toPts] if they are equal within a tolerance factor.
+// If none of the points are within the tolerance, the original pt is returned.
+func snap(pt Point, toPts ...Point) Point {
+	const tolerance = 3e-14
+	for _, p := range toPts {
+		if pt.equalWithin(p, tolerance) {
+			return p
+		}
+	}
+	return pt
+}
+
 // Returns the endpoints that were divided.
 func (c *clipper) possibleIntersection(e1, e2 *endpoint) []*endpoint {
 	numIntersections, ip1, _ := findIntersection(e1.segment(), e2.segment(), true)
@@ -397,6 +409,10 @@ func (c *clipper) possibleIntersection(e1, e2 *endpoint) []*endpoint {
 		return nil
 	}
 
+	// Adjust for floating point imprecision when intersections are created at endpoints, which
+	// otherwise has the tendency to corrupt the original polygons with new, almost-parallel segments.
+	ip1 = snap(ip1, e1.p, e2.p, e1.other.p, e2.other.p)
+
 	if numIntersections == 1 {
 		switch {
 		case e1.p.Equals(e2.p) || e1.other.p.Equals(e2.other.p):

clipper_test.go

@@ -0,0 +1,11 @@
+package polyclip
+
+import "testing"
+
+func TestSnap(t *testing.T) {
+	p := snap(Point{0, 0}, Point{0, 1e-9}, Point{1e-9, 0}, Point{1e-13, 1e-13})
+	verify(t, p.Equals(Point{0, 0}), "Expected no snapping but snapped to %v", p)
+
+	p = snap(Point{0, 0}, Point{0, 1e-9}, Point{1e-9, 0}, Point{1e-15, 1e-15})
+	verify(t, p.Equals(Point{1e-15, 1e-15}), "Expected snapping to {1e-15, 1e-15}")
+}

geom.go

@@ -27,6 +27,8 @@ package polyclip
 
 import (
 	"math"
+
+	"github.com/gonum/floats"
 )
 
 type Point struct {
@@ -38,6 +40,12 @@ func (p1 Point) Equals(p2 Point) bool {
 	return p1.X == p2.X && p1.Y == p2.Y
 }
 
+// equalWithin returns true if p1 is within tol tolerance of p2
+func (p1 Point) equalWithin(p2 Point, tol float64) bool {
+	return floats.EqualWithinAbsOrRel(p1.X, p2.X, tol, tol) &&
+		floats.EqualWithinAbsOrRel(p1.Y, p2.Y, tol, tol)
+}
+
 // Length returns distance from p to point (0, 0).
 func (p Point) Length() float64 {
 	return math.Sqrt(p.X*p.X + p.Y*p.Y)

geom_test.go

@@ -27,6 +27,7 @@ import (
 )
 
 func verify(t *T, cond bool, format string, args ...interface{}) {
+	t.Helper()
 	if !cond {
 		t.Errorf(format, args...)
 	}
@@ -43,6 +44,12 @@ func TestPoint(t *T) {
 	verify(t, circa(Point{3, 4}.Length(), 5), "Expected length 5")
 }
 
+func TestEqualWithin(t *T) {
+	p := Point{0, 0}
+	verify(t, !p.equalWithin(Point{0, 1e-9}, 1e-10), "Expected not equal")
+	verify(t, p.equalWithin(Point{1e-11, 1e-11}, 1e-10), "Expected equal")
+}
+
 func rect(x, y, w, h float64) Rectangle {
 	return Rectangle{Min: Point{x, y}, Max: Point{x + w, y + h}}
 }