From c69e4b183e944d50869d4676d9402ec462539acf Mon Sep 17 00:00:00 2001
From: Peter Stace <peterstace@gmail.com>
Date: Fri, 22 Dec 2023 13:52:53 +1100
Subject: [PATCH] Refactor/simplify MBR algorithm

---
 geom/alg_rotated_mbr.go | 159 +++++++++++++++-------------------------
 1 file changed, 60 insertions(+), 99 deletions(-)

diff --git a/geom/alg_rotated_mbr.go b/geom/alg_rotated_mbr.go
index 97fa82a9..722bbac4 100644
--- a/geom/alg_rotated_mbr.go
+++ b/geom/alg_rotated_mbr.go
@@ -19,132 +19,93 @@ func RotatedMinimumBoundingRectangle(g Geometry) Geometry {
 		panic(fmt.Sprintf("unexpected convex hull geometry type: %s", hull.Type()))
 	}
 
-	// Seek to the first line.
 	lnIdx := 0
-	for {
-		ln, ok := getLine(seq, lnIdx)
-		if ok {
-			break
-		}
-		lnIdx++
+	getLine := func() line {
+		xy0 := seq.GetXY(lnIdx + 0)
+		xy1 := seq.GetXY(lnIdx + 1)
+		return line{a: xy0, b: xy1}
 	}
 
 	// Set up distance calculators.
-	var nIdx, eIdx, wIdx int
-	nDist := func() float64 {
-		pt := seq.GetXY(nIdx)
-		ln, _ := getLine(seq, lnIdx)
-		return -signedPerpendicularDistance(pt, ln)
-	}
-	eDist := func() float64 {
-		pt := seq.GetXY(eIdx)
-		ln, _ := getLine(seq, lnIdx)
+	var rhsIdx, farIdx, lhsIdx int
+	rhsDist := func() float64 {
+		pt := seq.GetXY(rhsIdx)
+		ln := getLine()
 		ln = rotateLineCCW90(ln)
 		return signedPerpendicularDistance(pt, ln)
 	}
-	wDist := func() float64 {
-		pt := seq.GetXY(wIdx)
-		ln, _ := getLine(seq, lnIdx)
+	farDist := func() float64 {
+		pt := seq.GetXY(farIdx)
+		ln := getLine()
+		return -signedPerpendicularDistance(pt, ln)
+	}
+	lhsDist := func() float64 {
+		pt := seq.GetXY(lhsIdx)
+		ln := getLine()
 		ln = rotateLineCCW90(ln)
 		return -signedPerpendicularDistance(pt, ln)
 	}
 
-	// Establish initial antipodal positions.
-	eIdx = 0
-	for {
-		d0 := eDist()
-		eIdx = (eIdx + 1) % seq.Length()
-		d1 := eDist()
-		if d1 < d0 {
-			eIdx = (eIdx + seq.Length() - 1) % seq.Length()
-			break
-		}
-	}
-	nIdx = eIdx
-	for {
-		d0 := nDist()
-		nIdx = (nIdx + 1) % seq.Length()
-		d1 := nDist()
-		if d1 < d0 {
-			nIdx = (nIdx + seq.Length() - 1) % seq.Length()
-			break
-		}
-	}
-	wIdx = nIdx
-	for {
-		d0 := wDist()
-		wIdx = (wIdx + 1) % seq.Length()
-		d1 := wDist()
-		if d1 < d0 {
-			wIdx = (wIdx + seq.Length() - 1) % seq.Length()
-			break
+	advance := func(idx *int, dist func() float64) {
+		for {
+			d0 := dist()
+			*idx = (*idx + 1) % seq.Length()
+			d1 := dist()
+			if d1 < d0 {
+				*idx = (*idx + seq.Length() - 1) % seq.Length()
+				break
+			}
 		}
 	}
 
-	bestArea := nDist() * (eDist() + wDist())
-	bestLine, _ := getLine(seq, lnIdx)
-	bestN := nDist()
-	bestE := eDist()
-	bestW := wDist()
+	// Establish initial antipodal positions.
+	rhsIdx = 0
+	advance(&rhsIdx, rhsDist)
+	farIdx = rhsIdx
+	advance(&farIdx, farDist)
+	lhsIdx = farIdx
+	advance(&lhsIdx, lhsDist)
+
+	bestBaseLine := getLine()
+	bestRHSDist := rhsDist()
+	bestFarDist := farDist()
+	bestLHSDist := lhsDist()
+	bestArea := bestFarDist * (bestRHSDist + bestLHSDist)
 
-	for lnIdx+1 < seq.Length() {
+	for {
 		lnIdx++
-		ln, ok := getLine(seq, lnIdx)
-		if !ok {
-			continue
+		if lnIdx+1 == seq.Length() {
+			break
 		}
+		ln := getLine()
 
 		// Advance calipers:
-		for {
-			d0 := nDist()
-			nIdx = (nIdx + 1) % seq.Length()
-			d1 := nDist()
-			if d1 < d0 {
-				nIdx = (nIdx + seq.Length() - 1) % seq.Length()
-				break
-			}
-		}
-		for {
-			d0 := eDist()
-			eIdx = (eIdx + 1) % seq.Length()
-			d1 := eDist()
-			if d1 < d0 {
-				eIdx = (eIdx + seq.Length() - 1) % seq.Length()
-				break
-			}
-		}
-		for {
-			d0 := wDist()
-			wIdx = (wIdx + 1) % seq.Length()
-			d1 := wDist()
-			if d1 < d0 {
-				wIdx = (wIdx + seq.Length() - 1) % seq.Length()
-				break
-			}
-		}
-		area := nDist() * (eDist() + wDist())
+		advance(&rhsIdx, rhsDist)
+		advance(&farIdx, farDist)
+		advance(&lhsIdx, lhsDist)
+		area := farDist() * (rhsDist() + lhsDist())
 
 		if area < bestArea {
+			bestBaseLine = ln
 			bestArea = area
-			bestLine = ln
-			bestN = nDist()
-			bestE = eDist()
-			bestW = wDist()
+			bestRHSDist = rhsDist()
+			bestFarDist = farDist()
+			bestLHSDist = lhsDist()
 		}
 	}
 
-	eastBasis := bestLine.b.Sub(bestLine.a).Unit()
-	northBasis := eastBasis.rotateCCW90()
+	rhsBasis := bestBaseLine.b.Sub(bestBaseLine.a).Unit()
+	farBasis := rhsBasis.rotateCCW90()
 
-	nVec := northBasis.Scale(bestN)
-	eVec := eastBasis.Scale(bestE)
-	wVec := eastBasis.Scale(-bestW)
-	origin := bestLine.a
+	farVec := farBasis.Scale(bestFarDist)
+	rhsVec := rhsBasis.Scale(bestRHSDist)
+	lhsVec := rhsBasis.Scale(-bestLHSDist)
+	origin := bestBaseLine.a
 
-	pt0 := origin.Add(wVec)
-	pt1 := origin.Add(eVec)
-	pt2 := origin.Add(eVec).Add(nVec)
-	pt3 := origin.Add(wVec).Add(nVec)
+	pt0 := origin.Add(lhsVec)
+	pt1 := origin.Add(rhsVec)
+	pt2 := pt1.Add(farVec)
+	pt3 := pt0.Add(farVec)
 
 	coords := []float64{
 		pt0.X, pt0.Y,