Merge pull request #582 from peterstace/rot_mbr

Implement `RotatedMinimumAreaBoundingRectangle`
peterstace · Jan 15, 2024 · 88525b1 · 88525b1
2 parents 03cd665 + 38edcf1
commit 88525b1
Show file tree

Hide file tree

Showing 6 changed files with 322 additions and 0 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -11,6 +11,9 @@ YYYY-MM-DD
   `github.com/peterstace/simplefeatures/geos` package and the package used for
   reference implementation tests.
 
+- Adds a new function `RotatedMinimumAreaBoundingRectangle` that calculates the
+  minimum-area non-axis-aligned bounding rectangle for a geometry.
+
 ## v0.46.0
 
 2023-11-24

diff --git a/geom/alg_rotating_calipers.go b/geom/alg_rotating_calipers.go
@@ -0,0 +1,134 @@
+package geom
+
+import (
+	"fmt"
+	"math"
+)
+
+// RotatedMinimumAreaBoundingRectangle finds a rectangle with minimum area that
+// fully encloses the geometry. 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 RotatedMinimumAreaBoundingRectangle(g Geometry) Geometry {
+	hull := g.ConvexHull()
+	if hull.IsEmpty() {
+		return hull
+	}
+	var seq Sequence
+	switch hull.Type() {
+	case TypePoint, TypeLineString:
+		return hull
+	case TypePolygon:
+		seq = hull.MustAsPolygon().ExteriorRing().Coordinates()
+	default:
+		panic(fmt.Sprintf("unexpected convex hull geometry type: %s", hull.Type()))
+	}
+
+	rect := findBestMBR(seq)
+	return rect.asPoly().AsGeometry()
+}
+
+type rotatedRectangle struct {
+	origin XY // one of the corners
+	span1  XY // origin to first adjacent corner
+	span2  XY // origin to second adjacent corner
+}
+
+// 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
+}
+
+// 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++ {
+		rhs.update(seq, i)
+		if i == 0 {
+			far.idx = rhs.idx
+		}
+		far.update(seq, i)
+		if i == 0 {
+			lhs.idx = far.idx
+		}
+		lhs.update(seq, i)
+
+		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
+}
+
+// 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 {
+	orient func(XY) XY
+	idx    int
+	proj   XY
+}
+
+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 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) % n
+			c.proj = pt().proj(dir)
+			break
+		}
+		d0 = d1
+	}
+}
diff --git a/geom/alg_rotating_calipers_test.go b/geom/alg_rotating_calipers_test.go
@@ -0,0 +1,97 @@
+package geom_test
+
+import (
+	"strconv"
+	"testing"
+
+	"github.com/peterstace/simplefeatures/geom"
+)
+
+func TestRotatedMinimumAreaBoundingRectangle(t *testing.T) {
+	for i, tc := range []struct {
+		input string
+		want  string
+	}{
+		// Convex hull is non-polygon:
+		{
+			"POINT(0 0)",
+			"POINT(0 0)",
+		},
+		{
+			"LINESTRING(0 0,1 1)",
+			"LINESTRING(0 0,1 1)",
+		},
+		{
+			"LINESTRING(0 0,1 1,2 2)",
+			"LINESTRING(0 0,2 2)",
+		},
+
+		// Various polygons:
+		{
+			"POLYGON((0 0,0.5 0.5,0 1,0 0))",
+			"POLYGON((0 0,0.5 0.5,0 1,-0.5 0.5,0 0))",
+		},
+		{
+			"POLYGON((0 0,0.4 0.5,0 1,0 0))",
+			"POLYGON((0 0,0.4 0,0.4 1,0 1,0 0))",
+		},
+		{
+			"POLYGON((0 0,2 0,2 1,3 1,3 2,4 2,4 3,5 3,5 5,3 5,3 4,2 4,2 3,1 3,1 2,0 2,0 0))",
+			"POLYGON((1 -1,6 4,4 6,-1 1,1 -1))",
+		},
+		{
+			"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))",
+		},
+		{
+			"POLYGON((1 0,2 0,3 1,3 2,2 3,1 3,0 2,0 1,1 0))",
+			"POLYGON((1.5 -0.5,3.5 1.5,1.5 3.5,-0.5 1.5,1.5 -0.5))",
+		},
+		{
+			"POLYGON((1 0,2 0,7.5 6.5,6.5 7.5,0 2,0 1,1 0))",
+			"POLYGON((1.5 -0.5,8 6,6 8,-0.5 1.5,1.5 -0.5))",
+		},
+
+		// Randomly generated test cases:
+		{
+			"POLYGON((3 3,0 7,0 3,0 1,4 2,9 2,9 3,5 7,5 7,5 6,3 3))",
+			"POLYGON((9 7,0 7,0 1,9 1,9 7))",
+		},
+		{
+			"POLYGON((9 0,8 6,6 4,4 5,8 7,8 7,9 8,5 8,0 6,2 5,9 0))",
+			"POLYGON((0 6,9 0,12.69230769230769 5.538461538461537,3.692307692307692 11.538461538461538,0 6))",
+		},
+		{
+			"POLYGON((9 2,7 0,3 1,1 3,2 4,2 5,4 3,2 8,6 3,4 6,9 2))",
+			"POLYGON((2 8, -1.0352941176470587 4.45882352941176308, 5.96470588235294308 -1.54117647058823604, 9 2, 2 8))",
+		},
+		{
+			"POLYGON((5 7,4 8,4 5,4 5,7 3,5 7))",
+			"POLYGON((4.4 8.2,2.8 7.4,5.4 2.2,7 3,4.4 8.2))",
+		},
+
+		// Randomly generated "circle" polygons:
+		{
+			// 10 vertices:
+			"POLYGON((0.9556196410375619 0.29460329540458524,0.08050546571016876 0.9967541672803725,-0.015229280363664674 0.999884027785025,-0.7559431898170573 0.6546372230244875,-0.9999999693196223 0.00024771103011236986,-0.999974536263767 -0.00713630324916974,0.14434477669077697 -0.9895274556282356,0.3865654392869984 -0.9222619807564713,0.6109187782623596 -0.791693277959606,0.9966885771071856 -0.08131346914290244,0.9556196410375619 0.29460329540458524))",
+			"POLYGON((0.3151752625729588 -1.13618438607360916,1.29406186013216162 0.0040527140536265,-0.20767112381488262 1.29328132355450642,-1.18655772137408566 0.15304422342727064,0.3151752625729588 -1.13618438607360916))",
+		},
+		{
+			// 20 vertices:
+			"POLYGON((0.368304654881288 0.9297051582048877,0.20377681272133774 0.9790173699159442,0.16140993743170487 0.9868874465197606,-0.4672814384790904 0.8841086229943197,-0.6803702047922049 0.7328685997032572,-0.7454535020291218 0.6665576316512459,-0.8847125071437919 0.46613708252331304,-0.998089433868341 0.06178577506493315,-0.9995816820783242 0.028921632966826333,-0.9935505823458441 -0.1133897716737892,-0.6996146794426619 -0.7145203288279078,-0.48161673621854006 -0.876381948350262,-0.37858354221331025 -0.9255671242893314,-0.19789772311018608 -0.9802226742877376,-0.09045438773456249 -0.9959005993268432,0.38816671081425974 -0.921589173447496,0.5643906565655213 -0.825507835687548,0.668542998447214 -0.7436734896627744,0.7500093415680174 -0.6614272352728674,0.8552707422947223 -0.5181813942767867,0.368304654881288 0.9297051582048877))",
+			"POLYGON((0.29716778709735703 1.14121498796962717,-1.24146472783172235 0.62372841108425137,-0.60893993024358939 -1.25694495038847887,0.92969258468548999 -0.73945837350310351,0.29716778709735703 1.14121498796962717))",
+		},
+		{
+			// 100 vertices:
+			"POLYGON((0.9783225311822975 0.2070869985804575,0.9761682706719957 0.21701499333743116,0.962069932554095 0.27280294147123685,0.9010023908118991 0.43381412119851703,0.7933521058391652 0.6087630377746024,0.7786817019921365 0.6274191637036834,0.7717195720626366 0.6359629722676163,0.7361401192494257 0.676829169607399,0.7292795206122751 0.6842158875790083,0.6417286558283293 0.766931765080027,0.6338590240782793 0.7734486004865038,0.5502189210862213 0.8350204421921147,0.49490779203105156 0.8689454973627226,0.49015309590350964 0.8716363591407857,0.41875035262653315 0.90810139421496,0.3809482902007265 0.9245963444623514,0.36301868578745966 0.9317818595407112,0.3458095838437742 0.9383047115525935,0.3179560664869342 0.9481054476081003,0.29108762276537503 0.95669639691639,0.1200161828343841 0.9927719354705106,0.06375385215969982 0.9979656538853424,-0.10738115946901952 0.9942179271121041,-0.12807960475261765 0.9917638906747983,-0.1686977412889428 0.9856678305007265,-0.19469404249787903 0.9808640220824874,-0.27830909843970725 0.960491564630152,-0.44145277818091067 0.8972844836707897,-0.5654921182653689 0.8247536991003713,-0.592007934933976 0.805932134224222,-0.6073105419456033 0.7944645402041162,-0.610185438117451 0.7922586264026504,-0.6176701365323838 0.7864372844900388,-0.6417064319333313 0.7669503603333093,-0.6433188358486586 0.7655983773770204,-0.6880077135082702 0.725703373392409,-0.7208442488437409 0.693097084764395,-0.7252575559202017 0.6884776522013301,-0.7757075201968318 0.6310925788741949,-0.7998745158990114 0.6001672756951367,-0.807884909676552 0.5893402860121725,-0.8079611864719649 0.5892357093344011,-0.8718520773777028 0.48976928769797745,-0.9197935194772188 0.3924027032624916,-0.9594034748003358 0.2820371829050939,-0.9870802950652542 0.1602263745264516,-0.9876382270676111 0.156750542055669,-0.9901958928833767 0.13968569617857174,-0.9918925987842141 0.12707900092894184,-0.9998852854361413 0.01514648369378918,-0.9973836809900095 -0.07228964583409524,-0.9784320658833705 -0.20656885644065454,-0.9703404749327692 -0.241742347773922,-0.9503244676335602 -0.31126099372872046,-0.9478582515200036 -0.31869222617666965,-0.9337656645004919 -0.35788501477423557,-0.912730324299255 -0.4085625473591243,-0.9100518999868515 -0.4144943176092062,-0.901137798882359 -0.4335327754916082,-0.8581392302174916 -0.5134170444791749,-0.8397755119403577 -0.542933780073878,-0.8298375487293769 -0.558005056176751,-0.7836247646504528 -0.6212344390216324,-0.7792852269269228 -0.6266693985615177,-0.6025945771885426 -0.7980474769980552,-0.5861307068180728 -0.8102165108938144,-0.5408621467397948 -0.8411112519899022,-0.4508468164747215 -0.8926013377060382,-0.33535159222283967 -0.9420930472058513,-0.24360253214673028 -0.9698751498681164,-0.10725751354548622 -0.9942312737929941,-0.07675890649467973 -0.9970496829515273,0.16756073105256494 -0.9858617557290322,0.17005973899098473 -0.9854337548381008,0.17012294131708705 -0.9854228457051434,0.18941051783267693 -0.9818979864193413,0.25642224324558077 -0.9665648623702934,0.41476394903270414 -0.9099290448066797,0.4207576166292053 -0.9071730970705154,0.4648001246587809 -0.8854156335400802,0.5027717836210472 -0.8644191885853245,0.5209436605665887 -0.8535910628137358,0.5243168241487749 -0.8515232632844171,0.5500441677918465 -0.8351355659281762,0.6278936720739633 -0.7782991305208263,0.6867475755237128 -0.726895981218979,0.6880555208695567 -0.7256580463282434,0.6885449696252741 -0.7251936464171004,0.7333043470740315 -0.6799005328445686,0.798157080155061 -0.6024493965457577,0.7998623954804643 -0.6001834288750844,0.8046187456066828 -0.59379177681939,0.936378372445265 -0.35099222728823576,0.9475709274971743 -0.3195455168863192,0.9851065186498259 -0.17194518578204007,0.9933553786070294 -0.11508732246640076,0.9985831015816806 -0.05321455849211797,0.9993065951711081 -0.037233437224182625,0.9995631354820272 -0.029555679409181773,0.9997322331157119 -0.02314005335931624,0.9783225311822975 0.2070869985804575))",
+			"POLYGON((1.04198375442445279 -0.94582002496449491,0.95092601119575626 1.04268011264620686,-1.04425682449689683 0.95131635414775406,-0.9531990812682003 -1.03718378346294759,1.04198375442445279 -0.94582002496449491))",
+		},
+	} {
+		t.Run(strconv.Itoa(i), func(t *testing.T) {
+			in := geomFromWKT(t, tc.input)
+			t.Logf("input: %s", in.AsText())
+			got := geom.RotatedMinimumAreaBoundingRectangle(in)
+			expectGeomEqWKT(t, got, tc.want, geom.IgnoreOrder, geom.ToleranceXY(1e-13))
+		})
+	}
+}
diff --git a/geom/xy.go b/geom/xy.go
@@ -127,3 +127,25 @@ func (w XY) box() rtree.Box {
 		MaxY: w.Y,
 	}
 }
+
+// rotateCCW90 treats the XY as a vector, rotating it 90 degrees in a counter
+// clockwise direction (assuming a right handed/positive orientation).
+func (w XY) rotateCCW90() XY {
+	return XY{
+		X: -w.Y,
+		Y: w.X,
+	}
+}
+
+// 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))
+}
diff --git a/internal/cmprefimpl/cmpgeos/checks.go b/internal/cmprefimpl/cmpgeos/checks.go
@@ -42,6 +42,7 @@ func unaryChecks(g geom.Geometry, lg *log.Logger) error {
 		{"Centroid", checkCentroid},
 		{"PointOnSurface", checkPointOnSurface},
 		{"Simplify", checkSimplify},
+		{"RotatedMinimumAreaBoundingRectangle", checkRotatedMinimumAreaBoundingRectangle},
 	} {
 		lg.Printf("checking %s", check.name)
 		if err := check.fn(g, lg); err != nil {
@@ -577,6 +578,60 @@ func checkSimplify(g geom.Geometry, log *log.Logger) error {
 	return nil
 }
 
+func checkRotatedMinimumAreaBoundingRectangle(g geom.Geometry, log *log.Logger) error {
+	want, err := rawgeos.MinimumRotatedRectangle(g)
+	if err != nil {
+		return err
+	}
+	wantArea := want.Area()
+
+	got := geom.RotatedMinimumAreaBoundingRectangle(g)
+	gotArea := got.Area()
+
+	// The rotated bounding rectangle with minimum area is not always unique
+	// (multiple could exist with the same minimum area). Simplefeatures and
+	// GEOS (both correctly) choose different ones in some cases. To account
+	// for this, the comparison between the GEOS result and the simplefeatures
+	// result is broken into two parts...
+
+	// ...First, the areas are compared.
+	const areaDiffThreshold = 1e-10
+	if math.Abs(wantArea-gotArea) > areaDiffThreshold {
+		log.Printf("areas differ by more than %v", areaDiffThreshold)
+		log.Printf("want: (area %v) %v", wantArea, want.AsText())
+		log.Printf("got:  (area %v) %v", gotArea, got.AsText())
+		return errMismatch
+	}
+
+	// ...Second, we check if any of the input geometry is outside of the
+	// minimum bounding rectangle. Since GEOS cannot compute the difference of
+	// two geometries if one of them is a GeometryCollection, we break into
+	// parts and check the difference of each part individually.
+	var parts []geom.Geometry
+	if gc, ok := g.AsGeometryCollection(); ok {
+		parts = gc.Dump()
+	} else {
+		parts = []geom.Geometry{g}
+	}
+	for i, part := range parts {
+		overhang, err := rawgeos.Difference(part, got)
+		if err != nil {
+			return err
+		}
+		overhangArea := overhang.Area()
+		const overhangAreaThreshold = 1e-14
+		if overhangArea > overhangAreaThreshold {
+			log.Printf("part WKT (%d of %d): %v", i+1, len(parts), part.AsText())
+			log.Printf("part area overhangs MBR by %v (threshold %v)", overhangArea, overhangAreaThreshold)
+			log.Printf("overhang: (area %v) %v", overhangArea, overhang.AsText())
+			log.Printf("want: (area %v) %v", wantArea, want.AsText())
+			log.Printf("got:  (area %v) %v", gotArea, got.AsText())
+			return errMismatch
+		}
+	}
+	return nil
+}
+
 func binaryChecks(g1, g2 geom.Geometry, lg *log.Logger) error {
 	for _, g := range []geom.Geometry{g1, g2} {
 		if valid, err := checkIsValid(g, lg); err != nil {
@@ -620,6 +675,10 @@ func checkIntersects(g1, g2 geom.Geometry, log *log.Logger) error {
 		// https://github.com/peterstace/simplefeatures/issues/274
 		"LINESTRING(0.5 0,0.5000000000000001 0.5)":                              true,
 		"MULTILINESTRING((0 0,2 2.000000000000001),(1 0,-1 2.000000000000001))": true,
+
+		// GEOS gives the wrong result for the intersection of these two inputs:
+		"POLYGON((4.4 8.2,2.8 7.4,5.4 2.2,7 3,4.4 8.2))": true,
+		"POLYGON((1 4,3 4,3 7,1 7,1 4))":                 true,
 	}
 	if skipList[g1.AsText()] || skipList[g2.AsText()] {
 		// Skipping test because GEOS gives the incorrect result for *some*

diff --git a/internal/rawgeos/entrypoints.go b/internal/rawgeos/entrypoints.go
@@ -273,6 +273,13 @@ func Boundary(g geom.Geometry) (geom.Geometry, error) {
 	return result, wrap(err, "executing GEOSBoundary_r")
 }
 
+func MinimumRotatedRectangle(g geom.Geometry) (geom.Geometry, error) {
+	result, err := unaryOpG(g, func(ctx C.GEOSContextHandle_t, g *C.GEOSGeometry) *C.GEOSGeometry {
+		return C.GEOSMinimumRotatedRectangle_r(ctx, g)
+	})
+	return result, wrap(err, "executing GEOSMinimumRotatedRectangle_r")
+}
+
 func AsText(g geom.Geometry) (string, error) {
 	var result string
 	if err := unaryOpE(g, func(h *handle, gh *C.GEOSGeometry) error {