From fcff711bc788861605e677d13c3f58ccd9fd2fbe Mon Sep 17 00:00:00 2001 From: oceanful <35951082+oceanful@users.noreply.github.com> Date: Thu, 9 May 2019 09:57:22 -0700 Subject: [PATCH] Discard non-reductive segment divisions that otherwise recurse infinitely. (#1) Floating point imprecision in findIntersection() can create "non-reductive" divisions that result in infinite recursion. One class of non-reductive divisions can be detected at segment division time; divisions that result in segments going in the wrong direction (including zero-length segments, which have no direction) are non-reductive. This is detected by checking if endpoint.isValidDirection() in divideSegment(). The other class of non-reductive divisions does create "valid" segments, but one being infinitesimally small, and the other which recurses into a similar non-reductive division. This happens when the left or right endpoints of the two segments are very close but not equal, the problematic division for which manifests as an intersection point falling outside of the endpoints on a horizontal or vertical line. Note that theoretically, both classes of non-reductive division could be detected by comparing the length of the original segment against the length of the resulting segments, the latter of which should always be less than the former. Unfortunately, computation of segment length is subject to floating point imprecision and can introduce false positives. The rationale behind the current approach (isValidDirection and isValidSingleIntersection) is to rely only on boolean comparisons. --- bugs_test.go | 194 +++++++++++++++++++++++++++++++++++++------ clipper.go | 17 +++- endpoint.go | 49 +++++++++++ endpoint_test.go | 105 +++++++++++++++++++++++ intersection_test.go | 65 +++++++++++++++ 5 files changed, 403 insertions(+), 27 deletions(-) create mode 100644 intersection_test.go diff --git a/bugs_test.go b/bugs_test.go index 7dade6d..bcef72b 100644 --- a/bugs_test.go +++ b/bugs_test.go @@ -2,11 +2,12 @@ package polyclip_test import ( "fmt" + "math" "sort" . "testing" "time" - "github.com/akavel/polyclip-go" + polyclip "github.com/akavel/polyclip-go" ) type sorter polyclip.Polygon @@ -117,41 +118,186 @@ func TestBug3(t *T) { } } -func TestBug4(t *T) { +func TestNonReductiveSegmentDivisions(t *T) { if Short() { return } - cases := []struct{ subject, clipping, result polyclip.Polygon }{ - // original reported github issue #4, resulting in infinite loop + cases := []struct{ subject, clipping polyclip.Polygon }{ + { + // original reported github issue #4, resulting in infinite loop + subject: polyclip.Polygon{{ + {X: 1.427255375e+06, Y: -2.3283064365386963e-10}, + {X: 1.4271285e+06, Y: 134.7111358642578}, + {X: 1.427109e+06, Y: 178.30108642578125}}}, + clipping: polyclip.Polygon{{ + {X: 1.416e+06, Y: -12000}, + {X: 1.428e+06, Y: -12000}, + {X: 1.428e+06, Y: 0}, + {X: 1.416e+06, Y: 0}, + {X: 1.416e+06, Y: -12000}}}, + }, + // Test cases from https://github.com/ctessum/polyclip-go/blob/master/bugs_test.go + { + subject: polyclip.Polygon{{ + {X: 1.7714672107465276e+06, Y: -102506.68254093888}, + {X: 1.7713768917571804e+06, Y: -102000.75485953009}, + {X: 1.7717109214841307e+06, Y: -101912.19625031832}}}, + clipping: polyclip.Polygon{{ + {X: 1.7714593229229522e+06, Y: -102470.35230830211}, + {X: 1.7714672107465276e+06, Y: -102506.68254093867}, + {X: 1.771439738086082e+06, Y: -102512.92027456204}}}, + }, + { + subject: polyclip.Polygon{{ + {X: -1.8280000000000012e+06, Y: -492999.99999999953}, + {X: -1.8289999999999995e+06, Y: -494000.0000000006}, + {X: -1.828e+06, Y: -493999.9999999991}, + {X: -1.8280000000000012e+06, Y: -492999.99999999953}}}, + clipping: polyclip.Polygon{{ + {X: -1.8280000000000005e+06, Y: -495999.99999999977}, + {X: -1.8280000000000007e+06, Y: -492000.0000000014}, + {X: -1.8240000000000007e+06, Y: -492000.0000000014}, + {X: -1.8280000000000005e+06, Y: -495999.99999999977}}}, + }, + { + subject: polyclip.Polygon{{ + {X: -2.0199999999999988e+06, Y: -394999.99999999825}, + {X: -2.0199999999999988e+06, Y: -392000.0000000009}, + {X: -2.0240000000000012e+06, Y: -395999.9999999993}, + {X: -2.0199999999999988e+06, Y: -394999.99999999825}}}, + clipping: polyclip.Polygon{{ + {X: -2.0199999999999988e+06, Y: -394999.99999999825}, + {X: -2.020000000000001e+06, Y: -394000.0000000001}, + {X: -2.0190000000000005e+06, Y: -394999.9999999997}, + {X: -2.0199999999999988e+06, Y: -394999.99999999825}}}, + }, + { + subject: polyclip.Polygon{{ + {X: -47999.99999999992, Y: -23999.999999998756}, + {X: 0, Y: -24000.00000000017}, + {X: 0, Y: 24000.00000000017}, + {X: -48000.00000000014, Y: 24000.00000000017}, + {X: -47999.99999999992, Y: -23999.999999998756}}}, + clipping: polyclip.Polygon{{ + {X: -48000, Y: -24000}, + {X: 0, Y: -24000}, + {X: 0, Y: 24000}, + {X: -48000, Y: 24000}, + {X: -48000, Y: -24000}}}, + }, + { + subject: polyclip.Polygon{{ + {X: -2.137000000000001e+06, Y: -122000.00000000093}, + {X: -2.1360000000000005e+06, Y: -121999.99999999907}, + {X: -2.1360000000000014e+06, Y: -121000.00000000186}}}, + clipping: polyclip.Polygon{{ + {X: -2.1120000000000005e+06, Y: -120000}, + {X: -2.136000000000001e+06, Y: -120000.00000000093}, + {X: -2.1360000000000005e+06, Y: -144000}}}, + }, + { + subject: polyclip.Polygon{{ + {X: 1.556e+06, Y: -1.139999999999999e+06}, + {X: 1.5600000000000002e+06, Y: -1.140000000000001e+06}, + {X: 1.56e+06, Y: -1.136000000000001e+06}}}, + clipping: polyclip.Polygon{{ + {X: 1.56e+06, Y: -1.127999999999999e+06}, + {X: 1.5600000000000002e+06, Y: -1.151999999999999e+06}}}, + }, + { + subject: polyclip.Polygon{{ + {X: 1.0958876176594219e+06, Y: -567467.5197556159}, + {X: 1.0956330600760083e+06, Y: -567223.72588934}, + {X: 1.0958876176594219e+06, Y: -567467.5197556159}}}, + clipping: polyclip.Polygon{{ + {X: 1.0953516248896217e+06, Y: -564135.1861293605}, + {X: 1.0959085007300845e+06, Y: -568241.1879245406}, + {X: 1.0955136237022132e+06, Y: -581389.3748769956}}}, + }, + { + subject: polyclip.Polygon{{ + {X: 608000, Y: -113151.36476426799}, + {X: 608000, Y: -114660.04962779157}, + {X: 612000, Y: -115414.39205955336}, + {X: 1.616e+06, Y: -300000}, + {X: 1.608e+06, Y: -303245.6575682382}, + {X: 0, Y: 0}}}, + clipping: polyclip.Polygon{{ + {X: 1.612e+06, Y: -296000}}}, + }, + { + subject: polyclip.Polygon{{ + {X: 1.1458356382266793e+06, Y: -251939.4635597784}, + {X: 1.1460824662209095e+06, Y: -251687.86194535438}, + {X: 1.1458356382266793e+06, Y: -251939.4635597784}}}, + clipping: polyclip.Polygon{{ + {X: 1.1486683769211173e+06, Y: -251759.06331944838}, + {X: 1.1468807511323579e+06, Y: -251379.90576799586}, + {X: 1.1457914974731328e+06, Y: -251816.31287551578}}}, + }, { + // From https://github.com/ctessum/polyclip-go/commit/6614925d6d7087b7afcd4c55571554f67efd2ec3 subject: polyclip.Polygon{{ - {1.427255375e+06, -2.3283064365386963e-10}, - {1.4271285e+06, 134.7111358642578}, - {1.427109e+06, 178.30108642578125}}}, + {X: 426694.6365274183, Y: -668547.1611580737}, + {X: 426714.57523030025, Y: -668548.9238652373}, + {X: 426745.39648089616, Y: -668550.4651249861}}}, clipping: polyclip.Polygon{{ - {1.416e+06, -12000}, - {1.428e+06, -12000}, - {1.428e+06, 0}, - {1.416e+06, 0}, - {1.416e+06, -12000}}}, + {X: 426714.5752302991, Y: -668548.9238652373}, + {X: 426744.63718662335, Y: -668550.0591896093}, + {X: 426745.3964821229, Y: -668550.4652243527}}}, + }, + { + // Produces invalid divisions that would otherwise continually generate new segments. + subject: polyclip.Polygon{{ + {X: 99.67054939325573, Y: 23.50752393246498}, + {X: 99.88993946188153, Y: 20.999883973365655}, + {X: 100.01468418889, Y: 20.53433031419374}}}, + clipping: polyclip.Polygon{{ + {X: 100.15374164547939, Y: 20.015360821030836}, + {X: 95.64222842284941, Y: 36.85255738690467}, + {X: 100.15374164547939, Y: -14.714274712355238}}}, }, } + for _, c := range cases { - // check that we get a result in finite time - - ch := make(chan polyclip.Polygon) - go func() { - ch <- c.subject.Construct(polyclip.UNION, c.clipping) - }() - - select { - case <-ch: - case <-time.After(1 * time.Second): - // panicking in attempt to get full stacktrace - panic(fmt.Sprintf("case UNION:\nsubject: %v\nclipping: %v\ntimed out.", c.subject, c.clipping)) + const rotations = 360 + // Test multiple rotations of each case to catch any orientation assumptions. + for i := 0; i < rotations; i++ { + angle := 2 * math.Pi * float64(i) / float64(rotations) + subject := rotate(c.subject, angle) + clipping := rotate(c.clipping, angle) + + for _, op := range []polyclip.Op{polyclip.UNION, polyclip.INTERSECTION, polyclip.DIFFERENCE} { + ch := make(chan polyclip.Polygon) + go func() { + ch <- subject.Construct(op, clipping) + }() + + select { + case <-ch: + // check that we get a result in finite time + case <-time.After(1 * time.Second): + // panicking in attempt to get full stacktrace + panic(fmt.Sprintf("case %v:\nsubject: %v\nclipping: %v\ntimed out.", op, subject, clipping)) + } + } + } + } +} + +func rotate(p polyclip.Polygon, radians float64) polyclip.Polygon { + result := p.Clone() + for i, contour := range p { + result[i] = make(polyclip.Contour, len(contour)) + for j, point := range contour { + result[i][j] = polyclip.Point{ + X: point.X*math.Cos(radians) - point.Y*math.Sin(radians), + Y: point.Y*math.Cos(radians) + point.X*math.Sin(radians), + } } } + return result } func TestBug5(t *T) { diff --git a/clipper.go b/clipper.go index 057fe9f..a3c7d8b 100644 --- a/clipper.go +++ b/clipper.go @@ -284,7 +284,7 @@ func findIntersection(seg0, seg1 segment) (int, Point, Point) { d0 := Point{seg0.end.X - p0.X, seg0.end.Y - p0.Y} p1 := seg1.start d1 := Point{seg1.end.X - p1.X, seg1.end.Y - p1.Y} - sqrEpsilon := 1e-7 // was 1e-3 earlier + sqrEpsilon := 1e-15 // was originally 1e-3, which is very prone to false positives E := Point{p1.X - p0.X, p1.Y - p0.Y} kross := d0.X*d1.Y - d0.Y*d1.X sqrKross := kross * kross @@ -384,8 +384,13 @@ func (c *clipper) possibleIntersection(e1, e2 *endpoint) { return } - if numIntersections == 1 && (e1.p.Equals(e2.p) || e1.other.p.Equals(e2.other.p)) { - return // the line segments intersect at an endpoint of both line segments + if numIntersections == 1 { + if e1.p.Equals(e2.p) || e1.other.p.Equals(e2.other.p) { + return // the line segments intersect at an endpoint of both line segments + } else if !isValidSingleIntersection(e1, e2, ip1) { + _DBG(func() { fmt.Printf("Dropping invalid intersection %v between %v and %v\n", ip1, e1, e2) }) + return + } } //if numIntersections == 2 && e1.p.Equals(e2.p) { @@ -487,6 +492,12 @@ func (c *clipper) divideSegment(e *endpoint, p Point) { // "Left event" of the "right line segment" resulting from dividing e (the line segment associated to e) l := &endpoint{p: p, left: true, polygonType: e.polygonType, other: e.other, edgeType: e.other.edgeType} + // Discard segments of the wrong-direction (including zero-length). See isValidSingleIntersection() for reasoning. + if !l.isValidDirection() || !r.isValidDirection() { + _DBG(func() { fmt.Printf("Dropping invalid division of %v at %v:\n - %v\n - %v\n", *e, p, l, r) }) + return + } + if endpointLess(l, e.other) { // avoid a rounding error. The left event would be processed after the right event // println("Oops") e.other.left = true diff --git a/endpoint.go b/endpoint.go index ea907ab..ffdd882 100644 --- a/endpoint.go +++ b/endpoint.go @@ -78,3 +78,52 @@ func (se *endpoint) below(x Point) bool { func (se *endpoint) above(x Point) bool { return !se.below(x) } + +// leftRight() returns the left and right endpoints, in that order. +func (se *endpoint) leftRight() (Point, Point) { + if se.left { + return se.p, se.other.p + } + return se.other.p, se.p +} + +// isValid() is true if the segment has the correct direction. +// Note that segments of zero length have no direction and are thus not considered valid. +func (se *endpoint) isValidDirection() bool { + lp, rp := se.leftRight() + return lp.X < rp.X || (lp.X == rp.X && lp.Y < rp.Y) +} + +// Floating point imprecision in findIntersection() can create "non-reductive" +// divisions that result in infinite recursion. +// +// One class of non-reductive divisions can be detected at segment division time; +// divisions that result in segments going in the wrong direction +// (including zero-length segments, which have no direction) are non-reductive. +// This is detected by checking if endpoint.isValidDirection() in divideSegment(). +// +// The other class of non-reductive divisions does create "valid" segments; one +// being infinitesimally small, and the other which recurses into a similar +// non-reductive division. This happens when the left or right endpoints of the +// two segments are very close but not equal, the problematic division for which manifests +// as an intersection point falling outside of the endpoints on a horizontal or vertical line. +// +// Note that theoretically, both classes of non-reductive division could be detected by +// comparing the length of the original segment against the length of the resulting segments, +// the latter of which should always be less than the former. Unfortunately, computation of +// segment length is subject to floating point imprecision and can introduce false positives. +// The rationale behind the current approach (isValidDirection and isValidSingleIntersection) +// is to rely only on boolean comparisons. +func isValidSingleIntersection(e1, e2 *endpoint, ip Point) bool { + switch { + case e1.p.X == ip.X && e2.p.X == ip.X: // e1.p, ip, e2.p on a vertical line + return (ip.Y-e1.p.Y > 0) != (ip.Y-e2.p.Y > 0) // ip is above (or below) both e1.p and e2.p + case e1.p.Y == ip.Y && e2.p.Y == ip.Y: // e1.p, ip, e2.p on a horizontal line + return (ip.X-e1.p.X > 0) != (ip.X-e2.p.X > 0) // ip is to the left (or right) of both e1.p and e2.p + case e1.other.p.X == ip.X && e2.other.p.X == ip.X: // e1.other.p, ip, e2.other.p on a vertical line + return (ip.Y-e1.other.p.Y > 0) != (ip.Y-e2.other.p.Y > 0) // ip is above (or below) both e1.other.p and e2.other.p + case e1.other.p.Y == ip.Y && e2.other.p.Y == ip.Y: // e1.other.p, ip, e2.other.p on a horizontal line + return (ip.X-e1.other.p.X > 0) != (ip.X-e2.other.p.X > 0) // ip is to the left (or right) of both e1.other.p and e2.other.p + } + return true +} diff --git a/endpoint_test.go b/endpoint_test.go index 4d81ab4..3076bc5 100644 --- a/endpoint_test.go +++ b/endpoint_test.go @@ -39,3 +39,108 @@ func TestAbove(t *T) { verify(t, e.above(v.x) == v.result, "Expected %v above %v (case %d/b)", e, v.x, i) } } + +const ( + left = true + right = false + valid = true + invalid = false +) + +func TestEndpointIsValidDirection(t *T) { + cases := []struct { + left, right Point + dir bool + isValid bool + }{ + {Point{0, 1}, Point{0, 1}, left, invalid}, // Zero-length + {Point{0, 1}, Point{0, 1}, right, invalid}, // Zero-length + {Point{0, 1}, Point{1, 1}, left, valid}, // Horizontally valid + {Point{0, 1}, Point{-1, 1}, right, valid}, // Horizontally valid + {Point{0, 1}, Point{-1, 1}, left, invalid}, // Horizontally invalid + {Point{0, 1}, Point{1, 1}, right, invalid}, // Horizontally invalid + {Point{0, 1}, Point{0, 2}, left, valid}, // Vertically valid + {Point{0, 1}, Point{0, -1}, right, valid}, // Vertically valid + {Point{0, 1}, Point{0, -1}, left, invalid}, // Vertically invalid + {Point{0, 1}, Point{0, 2}, right, invalid}, // Vertically invalid + } + for i, v := range cases { + e := &endpoint{p: v.left, left: v.dir == left, other: &endpoint{p: v.right, left: v.dir != left}} + verify(t, e.isValidDirection() == v.isValid, "Expected %v isValidDirection()=%v (case %d)", e, v.isValid, i) + } +} + +func TestInvalidSingleIntersection(t *T) { + cases := []struct { + l1, r1, l2, r2, intersection Point + isValid bool + }{ + { + Point{0, 1.00000000000000}, Point{1, 1}, + Point{0, 1.00000000000001}, Point{3, 2}, + Point{0, 1.00000000000002}, + invalid, + }, + { + Point{0, 1.00000000000001}, Point{1, 1}, + Point{0, 1.00000000000002}, Point{3, 2}, + Point{0, 1.00000000000000}, + invalid, + }, + { + Point{0, 1.00000000000000}, Point{1, 1}, + Point{0, 1.00000000000001}, Point{3, 2}, + Point{0, 1.00000000000002}, + invalid, + }, + { + Point{0, 1.00000000000000}, Point{1, 1}, + Point{0, 1.00000000000002}, Point{3, 2}, + Point{0, 1.00000000000001}, + valid, + }, + { + Point{1.00000000000000, 0}, Point{1, 1}, + Point{1.00000000000001, 0}, Point{3, 2}, + Point{1.00000000000002, 0}, + invalid, + }, + { + Point{1.00000000000001, 0}, Point{1, 1}, + Point{1.00000000000002, 0}, Point{3, 2}, + Point{1.00000000000000, 0}, + invalid, + }, + { + Point{1.00000000000000, 0}, Point{1, 1}, + Point{1.00000000000001, 0}, Point{3, 2}, + Point{1.00000000000002, 0}, + invalid, + }, + { + Point{1.00000000000000, 0}, Point{1, 1}, + Point{1.00000000000002, 0}, Point{3, 2}, + Point{1.00000000000001, 0}, + valid, + }, + } + for i, v := range cases { + e1 := &endpoint{p: v.l1, left: true, other: &endpoint{p: v.r1, left: false}} + e2 := &endpoint{p: v.l2, left: true, other: &endpoint{p: v.r2, left: false}} + verify(t, isValidSingleIntersection(e1, e2, v.intersection) == v.isValid, + "Case %d: Expected intersection at %v of (%v, %v) to be isValidSingleIntersection()=%v", i, v.intersection, e1, e2, v.isValid) + e3 := &endpoint{p: v.r1, left: false, other: &endpoint{p: v.l1, left: true}} + e4 := &endpoint{p: v.r2, left: false, other: &endpoint{p: v.l2, left: true}} + verify(t, isValidSingleIntersection(e3, e4, v.intersection) == v.isValid, + "Case %d: Expected intersection at %v of (%v, %v) to be isValidSingleIntersection()=%v", i, v.intersection, e3, e4, v.isValid) + e5 := &endpoint{p: v.l1, left: false, other: &endpoint{p: v.r1, left: true}} + e6 := &endpoint{p: v.l2, left: false, other: &endpoint{p: v.r2, left: true}} + verify(t, isValidSingleIntersection(e5, e6, v.intersection) == v.isValid, + "Case %d: Expected intersection at %v of (%v, %v) to be isValidSingleIntersection()=%v", i, v.intersection, e5, e6, v.isValid) + e7 := &endpoint{p: v.r1, left: true, other: &endpoint{p: v.l1, left: false}} + e8 := &endpoint{p: v.r2, left: true, other: &endpoint{p: v.l2, left: false}} + verify(t, isValidSingleIntersection(e7, e8, v.intersection) == v.isValid, + "Case %d: Expected intersection at %v of (%v, %v) to be isValidSingleIntersection()=%v", i, v.intersection, e7, e8, v.isValid) + } + +} diff --git a/intersection_test.go b/intersection_test.go new file mode 100644 index 0000000..9fe00e5 --- /dev/null +++ b/intersection_test.go @@ -0,0 +1,65 @@ +package polyclip + +import "testing" + +func TestFindIntersection(t *testing.T) { + cases := []struct { + s1, s2 segment + numIntersections int + ip1, ip2 Point + }{ + { + // Almost (but not) parallel lines + segment{Point{0, 0}, Point{100, 0.0001}}, + segment{Point{1, 0}, Point{100, 0}}, + 0, Point{}, Point{}, + }, + { + // Almost (but not) parallel lines + segment{Point{0, 0}, Point{100, 0.0000001}}, + segment{Point{1, 0}, Point{100, 0}}, + 0, Point{}, Point{}, + }, + { + // Cross + segment{Point{1, 0}, Point{1, 3}}, + segment{Point{0, 1}, Point{3, 1}}, + 1, Point{1, 1}, Point{}, + }, + { + // Rays + segment{Point{0, 1}, Point{1, 3}}, + segment{Point{0, 1}, Point{3, 1}}, + 1, Point{0, 1}, Point{}, + }, + { + // Colinear rays + segment{Point{2, 1}, Point{0, 1}}, + segment{Point{2, 1}, Point{1, 1}}, + 2, Point{2, 1}, Point{}, // Why isn't this 2 intersections at {1,1} {2,1}? + }, + { + // Colinear rays + segment{Point{0, 3}, Point{0, 1}}, + segment{Point{0, 3}, Point{0, 2}}, + 2, Point{0, 3}, Point{}, // Why isn't this 2 intersections at {0,2}, {0,3}? + }, + { + // Overlapping segments + segment{Point{0, 1}, Point{3, 1}}, + segment{Point{1, 1}, Point{2, 1}}, + 1, Point{3, 1}, Point{}, // Why isn't this 2 intersections at {1,1} {2,1}? + }, + { + // Overlapping segments + segment{Point{0, 1}, Point{0, 4}}, + segment{Point{0, 2}, Point{0, 3}}, + 1, Point{0, 4}, Point{}, // Why isn't this 2 intersections at {0,2}, {0,3}? + }, + } + for i, v := range cases { + num, ip1, _ := findIntersection(v.s1, v.s2) + verify(t, num == v.numIntersections, "Case %d: Expected numIntersections to be %d, but got %d", i, v.numIntersections, num) + verify(t, ip1.Equals(v.ip1), "Case %d: Expected ip1 to be %v, but got %v", i, v.ip1, ip1) + } +}