s2: Add Rect.Centroid.

Signed-off-by: David Symonds <dsymonds@golang.org>

Author: rsned

s2/rect.go

@@ -637,5 +637,74 @@ func bisectorIntersection(lat r1.Interval, lng s1.Angle) Point {
 	return orthoLng.PointCross(PointFromLatLng(orthoBisector))
 }
 
+// Centroid returns the true centroid of the given Rect multiplied by its
+// surface area. The result is not unit length, so you may want to normalize it.
+// Note that in general the centroid is *not* at the center of the rectangle, and
+// in fact it may not even be contained by the rectangle. (It is the "center of
+// mass" of the rectangle viewed as subset of the unit sphere, i.e. it is the
+// point in space about which this curved shape would rotate.)
+//
+// The reason for multiplying the result by the rectangle area is to make it
+// easier to compute the centroid of more complicated shapes. The centroid
+// of a union of disjoint regions can be computed simply by adding their
+// Centroid results.
+func (r Rect) Centroid() Point {
+	// When a sphere is divided into slices of constant thickness by a set
+	// of parallel planes, all slices have the same surface area. This
+	// implies that the z-component of the centroid is simply the midpoint
+	// of the z-interval spanned by the Rect.
+	//
+	// Similarly, it is easy to see that the (x,y) of the centroid lies in
+	// the plane through the midpoint of the rectangle's longitude interval.
+	// We only need to determine the distance "d" of this point from the
+	// z-axis.
+	//
+	// Let's restrict our attention to a particular z-value. In this
+	// z-plane, the Rect is a circular arc. The centroid of this arc
+	// lies on a radial line through the midpoint of the arc, and at a
+	// distance from the z-axis of
+	//
+	//     r * (sin(alpha) / alpha)
+	//
+	// where r = sqrt(1-z^2) is the radius of the arc, and "alpha" is half
+	// of the arc length (i.e., the arc covers longitudes [-alpha, alpha]).
+	//
+	// To find the centroid distance from the z-axis for the entire
+	// rectangle, we just need to integrate over the z-interval. This gives
+	//
+	//    d = Integrate[sqrt(1-z^2)*sin(alpha)/alpha, z1..z2] / (z2 - z1)
+	//
+	// where [z1, z2] is the range of z-values covered by the rectangle.
+	// This simplifies to
+	//
+	//    d = sin(alpha)/(2*alpha*(z2-z1))*(z2*r2 - z1*r1 + theta2 - theta1)
+	//
+	// where [theta1, theta2] is the latitude interval, z1=sin(theta1),
+	// z2=sin(theta2), r1=cos(theta1), and r2=cos(theta2).
+	//
+	// Finally, we want to return not the centroid itself, but the centroid
+	// scaled by the area of the rectangle. The area of the rectangle is
+	//
+	//    A = 2 * alpha * (z2 - z1)
+	//
+	// which fortunately appears in the denominator of "d".
+
+	if r.IsEmpty() {
+		return Point{}
+	}
+
+	z1 := math.Sin(r.Lat.Lo)
+	z2 := math.Sin(r.Lat.Hi)
+	r1 := math.Cos(r.Lat.Lo)
+	r2 := math.Cos(r.Lat.Hi)
+
+	alpha := 0.5 * r.Lng.Length()
+	r0 := math.Sin(alpha) * (r2*z2 - r1*z1 + r.Lat.Length())
+	lng := r.Lng.Center()
+	z := alpha * (z2 + z1) * (z2 - z1) // scaled by the area
+
+	return Point{r3.Vector{r0 * math.Cos(lng), r0 * math.Sin(lng), z}}
+}
+
 // BUG: The major differences from the C++ version are:
-//   - GetCentroid, Get*Distance, Vertex, InteriorContains(LatLng|Rect|Point)
+//  - Get*Distance, Vertex, InteriorContains(LatLng|Rect|Point)

s2/rect_test.go

@@ -19,6 +19,7 @@ import (
 	"testing"
 
 	"github.com/golang/geo/r1"
+	"github.com/golang/geo/r2"
 	"github.com/golang/geo/r3"
 	"github.com/golang/geo/s1"
 )
@@ -1129,3 +1130,79 @@ func TestRectApproxEqual(t *testing.T) {
 		}
 	}
 }
+
+func TestRectCentroidEmptyFull(t *testing.T) {
+	// Empty and full rectangles.
+	if got, want := EmptyRect().Centroid(), (Point{}); !got.ApproxEqual(want) {
+		t.Errorf("%v.Centroid() = %v, want %v", EmptyRect(), got, want)
+	}
+	if got, want := FullRect().Centroid().Norm(), epsilon; !float64Eq(got, want) {
+		t.Errorf("%v.Centroid().Norm() = %v, want %v", FullRect(), got, want)
+	}
+}
+
+func testRectCentroidSplitting(t *testing.T, r Rect, leftSplits int) {
+	// Recursively verify that when a rectangle is split into two pieces, the
+	// centroids of the children sum to give the centroid of the parent.
+	var child0, child1 Rect
+	if oneIn(2) {
+		lat := randomUniformFloat64(r.Lat.Lo, r.Lat.Hi)
+		child0 = Rect{r1.Interval{r.Lat.Lo, lat}, r.Lng}
+		child1 = Rect{r1.Interval{lat, r.Lat.Hi}, r.Lng}
+	} else {
+		lng := randomUniformFloat64(r.Lng.Lo, r.Lng.Hi)
+		child0 = Rect{r.Lat, s1.Interval{r.Lng.Lo, lng}}
+		child1 = Rect{r.Lat, s1.Interval{lng, r.Lng.Hi}}
+	}
+
+	if got, want := r.Centroid().Sub(child0.Centroid().Vector).Sub(child1.Centroid().Vector).Norm(), 1e-15; got > want {
+		t.Errorf("%v.Centroid() - %v.Centroid() - %v.Centroid = %v, want ~0", r, child0, child1, got)
+	}
+	if leftSplits > 0 {
+		testRectCentroidSplitting(t, child0, leftSplits-1)
+		testRectCentroidSplitting(t, child1, leftSplits-1)
+	}
+}
+
+func TestRectCentroidFullRange(t *testing.T) {
+	// Rectangles that cover the full longitude range.
+	for i := 0; i < 100; i++ {
+		lat1 := randomUniformFloat64(-math.Pi/2, math.Pi/2)
+		lat2 := randomUniformFloat64(-math.Pi/2, math.Pi/2)
+		r := Rect{r1.Interval{lat1, lat2}, s1.FullInterval()}
+		centroid := r.Centroid()
+		if want := 0.5 * (math.Sin(lat1) + math.Sin(lat2)) * r.Area(); !float64Near(want, centroid.Z, epsilon) {
+			t.Errorf("%v.Centroid().Z was %v, want %v", r, centroid.Z, want)
+		}
+		if got := (r2.Point{centroid.X, centroid.Y}.Norm()); got > epsilon {
+			t.Errorf("%v.Centroid().Norm() was %v, want > %v ", r, got, epsilon)
+		}
+	}
+
+	// Rectangles that cover the full latitude range.
+	for i := 0; i < 100; i++ {
+		lat1 := randomUniformFloat64(-math.Pi, math.Pi)
+		lat2 := randomUniformFloat64(-math.Pi, math.Pi)
+		r := Rect{r1.Interval{-math.Pi / 2, math.Pi / 2}, s1.Interval{lat1, lat2}}
+		centroid := r.Centroid()
+
+		if got, want := math.Abs(centroid.Z), epsilon; got > want {
+			t.Errorf("math.Abs(%v.Centroid().Z) = %v, want <= %v", r, got, want)
+		}
+
+		if got, want := LatLngFromPoint(centroid).Lng.Radians(), r.Lng.Center(); !float64Near(got, want, epsilon) {
+			t.Errorf("%v.Lng.Radians() = %v, want %v", centroid, got, want)
+		}
+
+		alpha := 0.5 * r.Lng.Length()
+		if got, want := (r2.Point{centroid.X, centroid.Y}.Norm()), (0.25 * math.Pi * math.Sin(alpha) / alpha * r.Area()); !float64Near(got, want, epsilon) {
+			t.Errorf("%v.Centroid().Norm() = %v, want ~%v", got, want, epsilon)
+		}
+	}
+
+	// Finally, verify that when a rectangle is recursively split into pieces,
+	// the centroids of the pieces add to give the centroid of their parent.
+	// To make the code simpler we avoid rectangles that cross the 180 degree
+	// line of longitude.
+	testRectCentroidSplitting(t, Rect{r1.Interval{-math.Pi / 2, math.Pi / 2}, s1.Interval{-math.Pi, math.Pi}}, 10)
+}