Skip to content

Commit

Permalink
Discard non-reductive segment divisions that otherwise recurse infini…
Browse files Browse the repository at this point in the history 
…tely. (#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.
  • Loading branch information
@oceanful
oceanful authored May 9, 2019
1 parent 2cfdb71 commit fcff711
Show file tree
Hide file tree
Showing 5 changed files with 403 additions and 27 deletions.
194 changes: 170 additions & 24 deletions bugs_test.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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) {
Expand Down
17 changes: 14 additions & 3 deletions clipper.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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
Expand Down
49 changes: 49 additions & 0 deletions endpoint.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
}
Loading

0 comments on commit fcff711

Please sign in to comment.