Refactor to use dot product and vector projection

geom/alg_rotated_mbr.go

@@ -25,32 +25,50 @@ func RotatedMinimumBoundingRectangle(g Geometry) Geometry {
 		xy1 := seq.GetXY(lnIdx + 1)
 		return line{a: xy0, b: xy1}
 	}
+	getVect := func() XY {
+		ln := getLine()
+		return ln.b.Sub(ln.a)
+	}
 
 	// Set up distance calculators.
 	var rhsIdx, farIdx, lhsIdx int
-	rhsDist := func() float64 {
-		pt := seq.GetXY(rhsIdx)
-		ln := getLine()
-		ln = rotateLineCCW90(ln)
-		return signedPerpendicularDistance(pt, ln)
+	rhsProj := func() XY {
+		pt := seq.GetXY(rhsIdx).Sub(seq.GetXY(lnIdx))
+		onto := getVect()
+		return pt.proj(onto)
 	}
-	farDist := func() float64 {
-		pt := seq.GetXY(farIdx)
-		ln := getLine()
-		return -signedPerpendicularDistance(pt, ln)
+	farProj := func() XY {
+		pt := seq.GetXY(farIdx).Sub(seq.GetXY(lnIdx))
+		onto := getVect().rotateCCW90()
+		return pt.proj(onto)
 	}
-	lhsDist := func() float64 {
-		pt := seq.GetXY(lhsIdx)
-		ln := getLine()
-		ln = rotateLineCCW90(ln)
-		return -signedPerpendicularDistance(pt, ln)
+	lhsProj := func() XY {
+		pt := seq.GetXY(lhsIdx).Sub(seq.GetXY(lnIdx))
+		onto := getVect().rotateCCW90().rotateCCW90()
+		return pt.proj(onto)
+	}
+
+	rhsDot := func() float64 {
+		pt := seq.GetXY(rhsIdx).Sub(seq.GetXY(lnIdx))
+		onto := getVect()
+		return pt.Dot(onto)
+	}
+	farDot := func() float64 {
+		pt := seq.GetXY(farIdx).Sub(seq.GetXY(lnIdx))
+		onto := getVect().rotateCCW90()
+		return pt.Dot(onto)
+	}
+	lhsDot := func() float64 {
+		pt := seq.GetXY(lhsIdx).Sub(seq.GetXY(lnIdx))
+		onto := getVect().rotateCCW90().rotateCCW90()
+		return pt.Dot(onto)
 	}
 
-	advance := func(idx *int, dist func() float64) {
-		d0 := dist()
+	advance := func(idx *int, dot func() float64) {
+		d0 := dot()
 		for {
 			*idx = (*idx + 1) % seq.Length()
-			d1 := dist()
+			d1 := dot()
 			if d1 < d0 {
 				*idx = (*idx - 1) % seq.Length()
 				break
@@ -61,11 +79,11 @@ func RotatedMinimumBoundingRectangle(g Geometry) Geometry {
 
 	// Establish initial antipodal positions.
 	rhsIdx = 0
-	advance(&rhsIdx, rhsDist)
+	advance(&rhsIdx, rhsDot)
 	farIdx = rhsIdx
-	advance(&farIdx, farDist)
+	advance(&farIdx, farDot)
 	lhsIdx = farIdx
-	advance(&lhsIdx, lhsDist)
+	advance(&lhsIdx, lhsDot)
 
 	best := struct {
 		base int
@@ -75,7 +93,7 @@ func RotatedMinimumBoundingRectangle(g Geometry) Geometry {
 		area float64
 	}{
 		lnIdx, rhsIdx, farIdx, lhsIdx,
-		farDist() * (rhsDist() + lhsDist()),
+		rhsProj().Sub(lhsProj()).Cross(farProj()),
 	}
 
 	for {
@@ -85,10 +103,10 @@ func RotatedMinimumBoundingRectangle(g Geometry) Geometry {
 		}
 
 		// Advance calipers:
-		advance(&rhsIdx, rhsDist)
-		advance(&farIdx, farDist)
-		advance(&lhsIdx, lhsDist)
-		area := farDist() * (rhsDist() + lhsDist())
+		advance(&rhsIdx, rhsDot)
+		advance(&farIdx, farDot)
+		advance(&lhsIdx, lhsDot)
+		area := rhsProj().Sub(lhsProj()).Cross(farProj())
 
 		if area < best.area {
 			best.base = lnIdx
@@ -104,20 +122,10 @@ func RotatedMinimumBoundingRectangle(g Geometry) Geometry {
 	farIdx = best.far
 	lhsIdx = best.lhs
 
-	baseLine := getLine()
-
-	rhsBasis := baseLine.b.Sub(baseLine.a).Unit()
-	farBasis := rhsBasis.rotateCCW90()
-
-	rhsVec := rhsBasis.Scale(rhsDist())
-	farVec := farBasis.Scale(farDist())
-	lhsVec := rhsBasis.Scale(-lhsDist())
-	origin := baseLine.a
-
-	pt0 := origin.Add(lhsVec)
-	pt1 := origin.Add(rhsVec)
-	pt2 := pt1.Add(farVec)
-	pt3 := pt0.Add(farVec)
+	pt0 := seq.GetXY(lnIdx).Add(lhsProj())
+	pt1 := seq.GetXY(lnIdx).Add(rhsProj())
+	pt2 := pt1.Add(farProj())
+	pt3 := pt0.Add(farProj())
 
 	coords := []float64{
 		pt0.X, pt0.Y,
@@ -129,23 +137,3 @@ func RotatedMinimumBoundingRectangle(g Geometry) Geometry {
 	ring := NewLineString(NewSequence(coords, DimXY))
 	return NewPolygon([]LineString{ring}).AsGeometry()
 }
-
-func rotateLineCCW90(ln line) line {
-	return line{
-		a: ln.a,
-		b: ln.b.Sub(ln.a).rotateCCW90().Add(ln.a),
-	}
-}
-
-func signedPerpendicularDistance(pt XY, ln line) float64 {
-	x0 := pt.X
-	x1 := ln.a.X
-	x2 := ln.b.X
-	y0 := pt.Y
-	y1 := ln.a.Y
-	y2 := ln.b.Y
-
-	numer := (x2-x1)*(y1-y0) - (x1-x0)*(y2-y1) // TODO: determinant?
-	denom := ln.a.distanceTo(ln.b)
-	return numer / denom
-}

geom/alg_rotated_mbr_test.go

@@ -40,7 +40,7 @@ func TestRotatedMinimumBoundingRectangle(t *testing.T) {
 			"POLYGON((1 -1,6 4,4 6,-1 1,1 -1))",
 		},
 		{
-			"POLYGON((0 0,2 3,3 2,0 0))",
+			"POLYGON((0 0,-0.25 0.25,2 3,3 2,0 0))",
 			"POLYGON((0.5 -0.5,3 2,2 3,-0.5 0.5,0.5 -0.5))",
 		},
 		{
@@ -89,6 +89,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)
 			expectGeomEqWKT(t, got, tc.want, geom.IgnoreOrder, geom.ToleranceXY(1e-13))
 		})

geom/xy.go

@@ -136,3 +136,8 @@ func (w XY) rotateCCW90() XY {
 		Y: w.X,
 	}
 }
+
+// proj returns the projection of w onto o.
+func (w XY) proj(o XY) XY {
+	return o.Scale(w.Dot(o) / o.Dot(o))
+}