Bentley-Ottmann: three root-cause fixes for the dense shared-edge failures

path_intersection.go

@@ -294,6 +294,7 @@ type SweepPoint struct {
 	vertical   bool // segment is vertical
 	increasing bool // original direction is left-right (or bottom-top)
 	overlapped bool // segment's overlapping was handled
+	swept      bool // sweep fields have been computed for this left-endpoint
 }
 
 func (s *SweepPoint) InterpolateY(x float64) float64 {
@@ -307,6 +308,14 @@ func (s *SweepPoint) ToleranceEdgeY(xLeft, xRight float64) (float64, float64) {
 	if !s.left {
 		s = s.other
 	}
+	if s.X == s.other.X {
+		// Vertical segments have no x-extent: both endpoints lie in this column and the
+		// y-values at the tolerance edges are simply the endpoints. Interpolating would divide
+		// by zero and yield NaN whenever floating-point rounding puts the segment's x a ulp
+		// outside [xLeft, xRight), which then fails every below/above test and leaves the
+		// segment unbroken across the squares it passes through.
+		return s.Y, s.other.Y
+	}
 
 	y0 := s.Y
 	if s.X < xLeft {
@@ -326,6 +335,7 @@ func (s *SweepPoint) SplitAt(z Point) (*SweepPoint, *SweepPoint) {
 	*r, *l = *s.other, *s
 	r.Point, l.Point = z, z
 	r.end, l.end = false, false
+	r.swept, l.swept = false, false
 
 	// update references
 	r.other, s.other.other = s, l
@@ -401,6 +411,17 @@ func (q *SweepEvents) AddPathEndpoints(p *Path, seg int, clipping bool) int {
 		i += n
 		seg++
 
+		// The sweep decides verticality with an exact comparison, but only snaps coordinates to
+		// the grid after the intersection phase. A segment whose endpoints share a grid column is
+		// therefore driven through the sweep as left-to-right, and turns vertical the moment it
+		// is split; a piece that becomes vertical while pointing downwards has to be reversed to
+		// keep left-endpoints at the bottom, which is impossible once its left-endpoint is in the
+		// sweep status. Make it vertical up front, which moves a point no further than the
+		// snapping that follows the sweep would.
+		if start.X != end.X && snap(start.X, BentleyOttmannEpsilon) == snap(end.X, BentleyOttmannEpsilon) {
+			end.X = start.X
+		}
+
 		if start == end {
 			// skip zero-length lineTo or close command
 			continue
@@ -1659,8 +1680,29 @@ func (a eventSliceV) Swap(i, j int) {
 //	a[i], a[j] = a[j], a[i]
 //}
 
+// computeSweepFieldsInOrder computes the sweep fields of a non-vertical left-endpoint from the
+// segment below it in the sweep status, computing that segment first when it starts in the same
+// point and has not been computed yet. Left-endpoints of a tolerance square are visited in
+// CompareH order, but windings propagate bottom-up along the sweep status, and for overlapping
+// segments starting in the same point the two orders need not agree: a segment would then read
+// the windings of a neighbour that has not been computed yet.
+func (cur *SweepPoint) computeSweepFieldsInOrder(op pathOp, fillRule FillRule) {
+	if cur.swept {
+		return
+	}
+	var prev *SweepPoint
+	if node := cur.node.Prev(); node != nil {
+		prev = node.SweepPoint
+		if !prev.swept && prev.node != nil && prev.Point == cur.Point {
+			prev.computeSweepFieldsInOrder(op, fillRule)
+		}
+	}
+	cur.computeSweepFields(prev, op, fillRule)
+}
+
 func (cur *SweepPoint) computeSweepFields(prev *SweepPoint, op pathOp, fillRule FillRule) {
 	// cur is left-endpoint
+	cur.swept = true
 	if !cur.open {
 		cur.selfWindings = 1
 		if !cur.increasing {
@@ -2281,11 +2323,7 @@ func bentleyOttmann(ps, qs Paths, op pathOp, fillRule FillRule) Paths {
 						event.computeSweepFields(s, op, fillRule)
 					}
 				} else {
-					var s *SweepPoint
-					if event.node.Prev() != nil {
-						s = event.node.Prev().SweepPoint
-					}
-					event.computeSweepFields(s, op, fillRule)
+					event.computeSweepFieldsInOrder(op, fillRule)
 				}
 			}
 		}

path_intersection_coincident_test.go

@@ -0,0 +1,92 @@
+package canvas
+
+import (
+	"math"
+	"math/rand"
+	"testing"
+)
+
+// polygonArea returns the shoelace area of a flattened path (all subpaths).
+func polygonArea(p *Path) float64 {
+	a := 0.0
+	for _, sp := range p.Flatten(1e-3).Split() {
+		var pts []Point
+		d := sp.Data()
+		for i := 0; i < len(d); {
+			cmd := d[i]
+			n := cmdLen(cmd)
+			if cmd != CloseCmd {
+				pts = append(pts, Point{d[i+n-3], d[i+n-2]})
+			}
+			i += n
+		}
+		for i := range pts {
+			j := (i + 1) % len(pts)
+			a += pts[i].X*pts[j].Y - pts[j].X*pts[i].Y
+		}
+	}
+	return a / 2
+}
+
+// TestPathAndNearCoincidentEdge intersects a rounded rectangle with a plain
+// rectangle of the same size. The rounded rectangle is shifted right by less
+// than BentleyOttmannEpsilon, so its right edge nearly coincides with the
+// clip's. The clip's right edge x=50.000000005 lies exactly on a tolerance
+// square boundary (snap(x) - eps/2), which makes the raw-coordinate boundary
+// test in breakupCrossingSegments disagree with snap() about which column the
+// edge belongs to. The polygon walk then finds no next node and closes the
+// contour early: the result has two subpaths and a third of the area missing.
+func TestPathAndNearCoincidentEdge(t *testing.T) {
+	const w, h, r = 50.000000005, 13.0, 2.0
+	for _, tc := range []struct {
+		name    string
+		dx      float64
+		flatten bool
+	}{
+		{"flattened dx=2e-9", 2e-9, true},
+		{"flattened dx=5e-9", 5e-9, true},
+		{"curves dx=2e-9", 2e-9, false},
+		{"control dx=0", 0, true},
+		{"control dx=2e-8", 2e-8, true},
+	} {
+		clip := Rectangle(w, h)
+		p := RoundedRectangle(w, h, r).Translate(tc.dx, 0)
+		if tc.flatten {
+			p = p.Flatten(0.05)
+		}
+		res := p.And(clip)
+		n, a0, a1 := len(res.Split()), polygonArea(p), polygonArea(res)
+		if n != 1 || math.Abs(a1-a0) > 1e-3*a0 {
+			t.Errorf("%s: And() returned %d subpaths with area %.4f, want 1 subpath with area %.4f", tc.name, n, a1, a0)
+		}
+	}
+}
+
+// TestPathAndNearCoincidentEdgeRandom repeats the scenario at random positions
+// and sizes, with the right edge placed on a snap-rounding half-grid line in
+// half of the cases. Without the vertical-segment guard in ToleranceEdgeY
+// about 5% of the cases lose part of the shape.
+func TestPathAndNearCoincidentEdgeRandom(t *testing.T) {
+	rng := rand.New(rand.NewSource(1))
+	fails := 0
+	for i := 0; i < 1000; i++ {
+		x0 := rng.Float64() * 1000
+		y0 := rng.Float64() * 1000
+		w := 5 + rng.Float64()*500
+		h := 5 + rng.Float64()*50
+		if i%2 == 0 {
+			w = math.Round((x0+w)/BentleyOttmannEpsilon)*BentleyOttmannEpsilon + BentleyOttmannEpsilon/2 - x0
+		}
+		dx := (rng.Float64()*2 - 1) * BentleyOttmannEpsilon
+		dy := (rng.Float64()*2 - 1) * BentleyOttmannEpsilon
+		clip := Rectangle(w, h).Translate(x0, y0)
+		p := RoundedRectangle(w, h, 2).Translate(x0+dx, y0+dy).Flatten(0.05)
+		res := p.And(clip)
+		if a0, a1 := polygonArea(p), polygonArea(res); len(res.Split()) != 1 || math.Abs(a1-a0) > 1e-3*a0 {
+			fails++
+		}
+	}
+	if fails != 0 {
+		t.Errorf("%d of 1000 random near-coincident intersections returned a wrong polygon", fails)
+	}
+}

path_intersection_overlap_test.go

@@ -0,0 +1,62 @@
+package canvas
+
+import (
+	"math"
+	"math/rand"
+	"testing"
+)
+
+// noisyGrid builds n×n adjacent cells of size 0.1 whose corners carry random noise of magnitude
+// eps, mimicking floating-point drift in computed geometry such as Voronoi cells.
+func noisyGrid(n int, eps float64, rng *rand.Rand) Paths {
+	var ps Paths
+	nz := func() float64 { return (rng.Float64()*2 - 1) * eps }
+	for i := 0; i < n; i++ {
+		for j := 0; j < n; j++ {
+			x, y := 0.05+0.1*float64(i), 0.05+0.1*float64(j)
+			p := &Path{}
+			p.MoveTo(x+nz(), y+nz())
+			p.LineTo(x+0.1+nz(), y+nz())
+			p.LineTo(x+0.1+nz(), y+0.1+nz())
+			p.LineTo(x+nz(), y+0.1+nz())
+			p.Close()
+			ps = append(ps, p)
+		}
+	}
+	return ps
+}
+
+// TestSettleOverlappingStackOrder strokes a merged grid of nearly-exact cells, which settles the
+// outline with Settle(Positive). The stroke outlines of adjacent cells overlap exactly along the shared edges, so the
+// sweep sees stacks of identical segments starting in the same point. Windings must be computed
+// bottom-up along the sweep status; computing them in the square's CompareH order instead reads
+// a not-yet-computed neighbour, the windings above the stack come out wrong, and the result
+// polygon walk finds no continuation ("next node for result polygon is nil" in debug mode, a
+// silently truncated contour otherwise). The settled outline of the 10×10 grid is the single
+// outer square of side 1.2.
+func TestSettleOverlappingStackOrder(t *testing.T) {
+	DebugPathIntersection = true
+	defer func() { DebugPathIntersection = false }()
+
+	for _, seed := range []int64{6, 9, 12, 13} {
+		rng := rand.New(rand.NewSource(seed))
+		merged := noisyGrid(10, 1e-17, rng).Merge()
+
+		// Stroke settles the outline with Settle(Positive)
+		var res *Path
+		func() {
+			defer func() {
+				if r := recover(); r != nil {
+					t.Errorf("seed %d: panic: %v", seed, r)
+				}
+			}()
+			res = merged.Stroke(0.2, ButtCap, MiterJoin, 0.01)
+		}()
+		if res == nil {
+			continue
+		}
+		if n, a := len(res.Split()), polygonArea(res); n != 1 || math.Abs(a-1.44) > 1e-6 {
+			t.Errorf("seed %d: settled outline has %d subpaths and area %.6f, want 1 subpath of area 1.44", seed, n, a)
+		}
+	}
+}

path_intersection_util.go

@@ -183,6 +183,33 @@ func correctIntersection(z, aMin, aMax, bMin, bMax Point) Point {
 	return z
 }
 
+// correctIntersectionOrder restores the ordering the Bentley-Ottmann sweep relies on between a
+// segment's left-endpoint and an intersection on it: a0 must stay to the left of, or below, z,
+// and likewise b0. correctIntersection constrains x and y independently against each segment's
+// bounding box, which keeps z inside the box but lets it land in the same column as a
+// downwards-sloped segment's left-endpoint and below it. Splitting there makes the piece before
+// the split a vertical segment pointing downwards, which has to be reversed to keep
+// left-endpoints at the bottom, while that left-endpoint may already be in the sweep status where
+// it can no longer be reversed.
+//
+// Collapse z onto the left-endpoint only when the two snap to the same tolerance square, which
+// makes it a split the snapping phase would undo anyway; the segment is then left unsplit by the
+// tangential check in splitAtIntersections and the other segment of the pair is split at that
+// endpoint. A genuinely near-vertical segment keeps its intersections: those sit more than a
+// square away and are split normally, reversing the piece before the split while its
+// left-endpoint is still out of the status. Only that piece is constrained, (z,a1) and (z,b1)
+// are still free to become vertical as documented below.
+func correctIntersectionOrder(z, a0, b0 Point) Point {
+	eps := BentleyOttmannEpsilon
+	if z.X == a0.X && z.Y < a0.Y && snap(z.Y, eps) == snap(a0.Y, eps) {
+		z.Y = a0.Y
+	}
+	if z.X == b0.X && z.Y < b0.Y && snap(z.Y, eps) == snap(b0.Y, eps) {
+		z.Y = b0.Y
+	}
+	return z
+}
+
 // F. Antonio, "Faster Line Segment Intersection", Graphics Gems III, 1992
 func intersectionLineLineBentleyOttmann(zs []Point, a0, a1, b0, b1 Point) []Point {
 	// fast line-line intersection code, with additional constraints for the BentleyOttmann code:
@@ -257,6 +284,7 @@ func intersectionLineLineBentleyOttmann(zs []Point, a0, a1, b0, b1 Point) []Poin
 
 	z := a0.Interpolate(a1, ta)
 	z = correctIntersection(z, aMin, aMax, bMin, bMax)
+	z = correctIntersectionOrder(z, a0, b0)
 	if z != a0 && z != a1 || z != b0 && z != b1 {
 		// not at endpoints for both
 		if a0 != b0 && z != a0 && z != b0 && b0.Sub(z).PerpDot(z.Sub(a0)) == 0.0 {

path_intersection_vertical_test.go

@@ -0,0 +1,45 @@
+package canvas
+
+import (
+	"math"
+	"math/rand"
+	"testing"
+)
+
+// TestSettleVerticalFirstSegment strokes a merged grid of nearly-exact cells whose corners carry
+// noise large enough to leave the shared edges a fraction of the snapping grid apart. The stroke
+// outline then holds segments whose two endpoints fall in the same grid column while differing in
+// x, which the sweep calls left-to-right because it compares verticality exactly and only snaps to
+// the grid after the intersection phase. Splitting such a segment turns the piece before the split
+// vertical, and when it points downwards that piece has to be reversed to keep left-endpoints at
+// the bottom, which is impossible once its left-endpoint sits in the sweep status: "first segment
+// became vertical and needs reversal, but was already in the sweep status". The settled outline of
+// the 10x10 grid is the single outer square of side 1.2.
+func TestSettleVerticalFirstSegment(t *testing.T) {
+	DebugPathIntersection = true
+	defer func() { DebugPathIntersection = false }()
+
+	for _, eps := range []float64{1e-16, 1e-14, 1e-13, 1e-12} {
+		for _, seed := range []int64{1, 2, 3, 4} {
+			rng := rand.New(rand.NewSource(seed))
+			merged := noisyGrid(10, eps, rng).Merge()
+
+			// Stroke settles the outline with Settle(Positive)
+			var res *Path
+			func() {
+				defer func() {
+					if r := recover(); r != nil {
+						t.Errorf("eps %g seed %d: panic: %v", eps, seed, r)
+					}
+				}()
+				res = merged.Stroke(0.2, ButtCap, MiterJoin, 0.01)
+			}()
+			if res == nil {
+				continue
+			}
+			if n, a := len(res.Split()), polygonArea(res); n != 1 || math.Abs(a-1.44) > 1e-6 {
+				t.Errorf("eps %g seed %d: settled outline has %d subpaths and area %.6f, want 1 subpath of area 1.44", eps, seed, n, a)
+			}
+		}
+	}
+}