Restructure algorithm to be easier to understand

geom/alg_rotated_mbr.go

@@ -5,12 +5,12 @@ import (
 	"math"
 )
 
-// RotatedMinimumBoundingRectangle finds a rectangle with minimum area that
+// RotatedMinimumAreaBoundingRectangle finds a rectangle with minimum area that
 // fully encloses the geometry. The solution is not necessarily unique. If the
 // geometry is empty, the empty geometry of the same type is returned. If the
 // bounding rectangle would be degenerate (zero area), then point or line
 // string (with a single line segment) will be returned.
-func RotatedMinimumBoundingRectangle(g Geometry) Geometry {
+func RotatedMinimumAreaBoundingRectangle(g Geometry) Geometry {
 	hull := g.ConvexHull()
 	if hull.IsEmpty() {
 		return hull
@@ -25,86 +25,107 @@ func RotatedMinimumBoundingRectangle(g Geometry) Geometry {
 		panic(fmt.Sprintf("unexpected convex hull geometry type: %s", hull.Type()))
 	}
 
-	calipers := newCalipers(seq)
-	calipers.rotate(seq)
-	return calipers.makePoly()
+	rect := findBestMBR(seq)
+	return rect.asPoly().AsGeometry()
 }
 
-type calipers struct {
-	lhs, far, rhs caliper
-	bestArea      float64
-	bestOrigin    XY
-	bestSpan1     XY
-	bestSpan2     XY
+type rotatedRectangle struct {
+	origin XY // one of the corners
+	span1  XY // origin to first adjacent corner
+	span2  XY // origin to second adjacent corner
 }
 
-func newCalipers(seq Sequence) *calipers {
-	return &calipers{
-		lhs:      caliper{seq: seq},
-		far:      caliper{seq: seq},
-		rhs:      caliper{seq: seq},
-		bestArea: math.Inf(+1),
+// asPoly converts the rectangle to a polygon by traversing from the
+// rectangle's origin to its other corners via its spans.
+func (r rotatedRectangle) asPoly() Polygon {
+	pts := [5]XY{
+		r.origin,
+		r.origin.Add(r.span1),
+		r.origin.Add(r.span1).Add(r.span2),
+		r.origin.Add(r.span2),
+		r.origin,
+	}
+	coords := make([]float64, 2*len(pts))
+	for i, pt := range pts {
+		coords[2*i+0] = pt.X
+		coords[2*i+1] = pt.Y
 	}
+	ring := NewLineString(NewSequence(coords, DimXY))
+	poly := NewPolygon([]LineString{ring})
+	return poly
 }
 
-func (c *calipers) rotate(seq Sequence) {
+// findBestMBR finds the minimum area bounding rectangle for a convex ring
+// specified as a sequence. It does this by enumerating each candidate rotated
+// bounding rectangle, and finding the one with the minimum area. There is a
+// candidate rectangle corresponding to each edge in the convex ring.
+func findBestMBR(seq Sequence) rotatedRectangle {
+	rhs := caliper{orient: func(in XY) XY { return in }}
+	far := caliper{orient: XY.rotateCCW90}
+	lhs := caliper{orient: XY.rotate180}
+
+	var rect rotatedRectangle
+	bestArea := math.Inf(+1)
+
 	for i := 0; i+1 < seq.Length(); i++ {
-		heading := seq.GetXY(i + 1).Sub(seq.GetXY(i))
-		c.rhs.update(seq.GetXY(i), heading)
+		rhs.update(seq, i)
 		if i == 0 {
-			c.far.idx = c.rhs.idx
+			far.idx = rhs.idx
 		}
-		c.far.update(seq.GetXY(i), heading.rotateCCW90())
+		far.update(seq, i)
 		if i == 0 {
-			c.lhs.idx = c.far.idx
+			lhs.idx = far.idx
 		}
-		c.lhs.update(seq.GetXY(i), heading.rotate180())
+		lhs.update(seq, i)
 
-		area := c.rhs.proj.Sub(c.lhs.proj).Cross(c.far.proj)
-		if area < c.bestArea {
-			c.bestArea = area
-			c.bestOrigin = seq.GetXY(i).Add(c.lhs.proj)
-			c.bestSpan1 = c.rhs.proj.Sub(c.lhs.proj)
-			c.bestSpan2 = c.far.proj
+		area := rhs.proj.Sub(lhs.proj).Cross(far.proj)
+		if area < bestArea {
+			bestArea = area
+			rect = rotatedRectangle{
+				origin: seq.GetXY(i).Add(lhs.proj),
+				span1:  rhs.proj.Sub(lhs.proj),
+				span2:  far.proj,
+			}
 		}
 	}
+	return rect
 }
 
-func (c *calipers) makePoly() Geometry {
-	pts := [5]XY{
-		c.bestOrigin,
-		c.bestOrigin.Add(c.bestSpan1),
-		c.bestOrigin.Add(c.bestSpan1).Add(c.bestSpan2),
-		c.bestOrigin.Add(c.bestSpan2),
-		c.bestOrigin,
-	}
-	coords := make([]float64, 2*len(pts))
-	for i, pt := range pts {
-		coords[2*i] = pt.X
-		coords[2*i+1] = pt.Y
-	}
-	ring := NewLineString(NewSequence(coords, DimXY))
-	poly := NewPolygon([]LineString{ring})
-	return poly.AsGeometry()
-}
-
+// caliper is a helper struct for finding the maximum perpendicular distance
+// between a line segment on the (convex) ring and a point on the same ring. The
+// perpendicular distance may be to the "right hand side" the ring, "across"
+// the ring, or to the "left hand side" of the ring.
+//
+// It assumes that the line segment measured against rotates around the ring
+// iteratively (which is one of the key properties that allows this algorithm
+// to execute quickly).
+//
+// This is an example of a "rotating calipers" algorithm. See
+// [rotating_calipers] for a general explanation of the technique.
+//
+// [rotating_calipers]: https://en.wikipedia.org/wiki/Rotating_calipers
 type caliper struct {
-	seq  Sequence
-	idx  int
-	proj XY
+	orient func(XY) XY
+	idx    int
+	proj   XY
 }
 
-func (c *caliper) update(offset, dir XY) {
-	n := c.seq.Length()
+func (c *caliper) update(seq Sequence, lnIdx int) {
+	offset := seq.GetXY(lnIdx)
+	dir := seq.GetXY(lnIdx + 1).Sub(offset)
+	dir = c.orient(dir)
+
+	n := seq.Length()
 	pt := func() XY {
-		return c.seq.GetXY(c.idx).Sub(offset)
+		return seq.GetXY(c.idx).Sub(offset)
 	}
+
 	d0 := pt().Dot(dir)
 	for {
 		c.idx = (c.idx + 1) % n
 		d1 := pt().Dot(dir)
 		if d1 < d0 {
-			c.idx = (c.idx - 1) % n
+			c.idx = (c.idx - 1 + n) % n
 			c.proj = pt().proj(dir)
 			break
 		}

geom/alg_rotated_mbr_test.go

@@ -7,7 +7,7 @@ import (
 	"github.com/peterstace/simplefeatures/geom"
 )
 
-func TestRotatedMinimumBoundingRectangle(t *testing.T) {
+func TestRotatedMinimumAreaBoundingRectangle(t *testing.T) {
 	for i, tc := range []struct {
 		input string
 		want  string
@@ -90,7 +90,7 @@ func TestRotatedMinimumBoundingRectangle(t *testing.T) {
 		t.Run(strconv.Itoa(i), func(t *testing.T) {
 			in := geomFromWKT(t, tc.input)
 			t.Logf("input: %s", in.AsText())
-			got := geom.RotatedMinimumBoundingRectangle(in)
+			got := geom.RotatedMinimumAreaBoundingRectangle(in)
 			expectGeomEqWKT(t, got, tc.want, geom.IgnoreOrder, geom.ToleranceXY(1e-13))
 		})
 	}