Minor refactor and add doc comment

geom/alg_rotated_mbr.go

@@ -5,6 +5,11 @@ import (
 	"math"
 )
 
+// RotatedMinimumBoundingRectangle 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 {
 	hull := g.ConvexHull()
 	if hull.IsEmpty() {
@@ -20,64 +25,87 @@ func RotatedMinimumBoundingRectangle(g Geometry) Geometry {
 		panic(fmt.Sprintf("unexpected convex hull geometry type: %s", hull.Type()))
 	}
 
-	heading := seq.GetXY(1).Sub(seq.GetXY(0))
-
-	rhs := caliper{idx: 0}
-	rhs.update(seq, seq.GetXY(0), heading)
-
-	far := caliper{idx: rhs.idx}
-	far.update(seq, seq.GetXY(0), heading.rotateCCW90())
+	calipers := newCalipers(seq)
+	calipers.rotate(seq)
+	return calipers.makePoly()
+}
 
-	lhs := caliper{idx: far.idx}
-	lhs.update(seq, seq.GetXY(0), heading.rotateCCW90().rotateCCW90())
+type calipers struct {
+	lhs, far, rhs caliper
+	bestArea      float64
+	bestOrigin    XY
+	bestSpan1     XY
+	bestSpan2     XY
+}
 
-	var best struct {
-		area float64
-		pts  [5]XY
+func newCalipers(seq Sequence) *calipers {
+	return &calipers{
+		lhs:      caliper{seq: seq},
+		far:      caliper{seq: seq},
+		rhs:      caliper{seq: seq},
+		bestArea: math.Inf(+1),
 	}
-	best.area = math.Inf(+1)
+}
 
+func (c *calipers) rotate(seq Sequence) {
 	for i := 0; i+1 < seq.Length(); i++ {
 		heading := seq.GetXY(i + 1).Sub(seq.GetXY(i))
-		rhs.update(seq, seq.GetXY(i), heading)
-		far.update(seq, seq.GetXY(i), heading.rotateCCW90())
-		lhs.update(seq, seq.GetXY(i), heading.rotateCCW90().rotateCCW90())
-		if area := rhs.proj.Sub(lhs.proj).Cross(far.proj); area < best.area {
-			best.area = area
-			best.pts = [5]XY{
-				seq.GetXY(i).Add(lhs.proj),
-				seq.GetXY(i).Add(rhs.proj),
-				seq.GetXY(i).Add(rhs.proj).Add(far.proj),
-				seq.GetXY(i).Add(lhs.proj).Add(far.proj),
-				seq.GetXY(i).Add(lhs.proj),
-			}
+		c.rhs.update(seq.GetXY(i), heading)
+		if i == 0 {
+			c.far.idx = c.rhs.idx
+		}
+		c.far.update(seq.GetXY(i), heading.rotateCCW90())
+		if i == 0 {
+			c.lhs.idx = c.far.idx
+		}
+		c.lhs.update(seq.GetXY(i), heading.rotate180())
+
+		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
 		}
 	}
+}
 
-	coords := make([]float64, 0, 2*len(best.pts))
-	for _, pt := range best.pts {
-		coords = append(coords, pt.X, pt.Y)
+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))
-	return NewPolygon([]LineString{ring}).AsGeometry()
+	poly := NewPolygon([]LineString{ring})
+	return poly.AsGeometry()
 }
 
 type caliper struct {
+	seq  Sequence
 	idx  int
 	proj XY
 }
 
-func (c *caliper) update(seq Sequence, offset, heading XY) {
-	dot := func() float64 {
-		return seq.GetXY(c.idx).Sub(offset).Dot(heading)
+func (c *caliper) update(offset, dir XY) {
+	n := c.seq.Length()
+	pt := func() XY {
+		return c.seq.GetXY(c.idx).Sub(offset)
 	}
-	d0 := dot()
+	d0 := pt().Dot(dir)
 	for {
-		c.idx = (c.idx + 1) % seq.Length()
-		d1 := dot()
+		c.idx = (c.idx + 1) % n
+		d1 := pt().Dot(dir)
 		if d1 < d0 {
-			c.idx = (c.idx - 1) % seq.Length()
-			c.proj = seq.GetXY(c.idx).Sub(offset).proj(heading)
+			c.idx = (c.idx - 1) % n
+			c.proj = pt().proj(dir)
 			break
 		}
 		d0 = d1

geom/xy.go

@@ -137,6 +137,14 @@ func (w XY) rotateCCW90() XY {
 	}
 }
 
+// rotate180 treats the XY as a vector, rotating it 180 degrees.
+func (w XY) rotate180() XY {
+	return XY{
+		X: -w.X,
+		Y: -w.Y,
+	}
+}
+
 // proj returns the projection of w onto o.
 func (w XY) proj(o XY) XY {
 	return o.Scale(w.Dot(o) / o.Dot(o))