diff options
Diffstat (limited to 'vendor/github.com/golang/geo/s2/rect.go')
| -rw-r--r-- | vendor/github.com/golang/geo/s2/rect.go | 710 |
1 files changed, 710 insertions, 0 deletions
diff --git a/vendor/github.com/golang/geo/s2/rect.go b/vendor/github.com/golang/geo/s2/rect.go new file mode 100644 index 000000000..cb4c93180 --- /dev/null +++ b/vendor/github.com/golang/geo/s2/rect.go @@ -0,0 +1,710 @@ +// Copyright 2014 Google Inc. All rights reserved. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +package s2 + +import ( + "fmt" + "io" + "math" + + "github.com/golang/geo/r1" + "github.com/golang/geo/r3" + "github.com/golang/geo/s1" +) + +// Rect represents a closed latitude-longitude rectangle. +type Rect struct { + Lat r1.Interval + Lng s1.Interval +} + +var ( + validRectLatRange = r1.Interval{-math.Pi / 2, math.Pi / 2} + validRectLngRange = s1.FullInterval() +) + +// EmptyRect returns the empty rectangle. +func EmptyRect() Rect { return Rect{r1.EmptyInterval(), s1.EmptyInterval()} } + +// FullRect returns the full rectangle. +func FullRect() Rect { return Rect{validRectLatRange, validRectLngRange} } + +// RectFromLatLng constructs a rectangle containing a single point p. +func RectFromLatLng(p LatLng) Rect { + return Rect{ + Lat: r1.Interval{p.Lat.Radians(), p.Lat.Radians()}, + Lng: s1.Interval{p.Lng.Radians(), p.Lng.Radians()}, + } +} + +// RectFromCenterSize constructs a rectangle with the given size and center. +// center needs to be normalized, but size does not. The latitude +// interval of the result is clamped to [-90,90] degrees, and the longitude +// interval of the result is FullRect() if and only if the longitude size is +// 360 degrees or more. +// +// Examples of clamping (in degrees): +// center=(80,170), size=(40,60) -> lat=[60,90], lng=[140,-160] +// center=(10,40), size=(210,400) -> lat=[-90,90], lng=[-180,180] +// center=(-90,180), size=(20,50) -> lat=[-90,-80], lng=[155,-155] +func RectFromCenterSize(center, size LatLng) Rect { + half := LatLng{size.Lat / 2, size.Lng / 2} + return RectFromLatLng(center).expanded(half) +} + +// IsValid returns true iff the rectangle is valid. +// This requires Lat ⊆ [-π/2,π/2] and Lng ⊆ [-π,π], and Lat = ∅ iff Lng = ∅ +func (r Rect) IsValid() bool { + return math.Abs(r.Lat.Lo) <= math.Pi/2 && + math.Abs(r.Lat.Hi) <= math.Pi/2 && + r.Lng.IsValid() && + r.Lat.IsEmpty() == r.Lng.IsEmpty() +} + +// IsEmpty reports whether the rectangle is empty. +func (r Rect) IsEmpty() bool { return r.Lat.IsEmpty() } + +// IsFull reports whether the rectangle is full. +func (r Rect) IsFull() bool { return r.Lat.Equal(validRectLatRange) && r.Lng.IsFull() } + +// IsPoint reports whether the rectangle is a single point. +func (r Rect) IsPoint() bool { return r.Lat.Lo == r.Lat.Hi && r.Lng.Lo == r.Lng.Hi } + +// Vertex returns the i-th vertex of the rectangle (i = 0,1,2,3) in CCW order +// (lower left, lower right, upper right, upper left). +func (r Rect) Vertex(i int) LatLng { + var lat, lng float64 + + switch i { + case 0: + lat = r.Lat.Lo + lng = r.Lng.Lo + case 1: + lat = r.Lat.Lo + lng = r.Lng.Hi + case 2: + lat = r.Lat.Hi + lng = r.Lng.Hi + case 3: + lat = r.Lat.Hi + lng = r.Lng.Lo + } + return LatLng{s1.Angle(lat) * s1.Radian, s1.Angle(lng) * s1.Radian} +} + +// Lo returns one corner of the rectangle. +func (r Rect) Lo() LatLng { + return LatLng{s1.Angle(r.Lat.Lo) * s1.Radian, s1.Angle(r.Lng.Lo) * s1.Radian} +} + +// Hi returns the other corner of the rectangle. +func (r Rect) Hi() LatLng { + return LatLng{s1.Angle(r.Lat.Hi) * s1.Radian, s1.Angle(r.Lng.Hi) * s1.Radian} +} + +// Center returns the center of the rectangle. +func (r Rect) Center() LatLng { + return LatLng{s1.Angle(r.Lat.Center()) * s1.Radian, s1.Angle(r.Lng.Center()) * s1.Radian} +} + +// Size returns the size of the Rect. +func (r Rect) Size() LatLng { + return LatLng{s1.Angle(r.Lat.Length()) * s1.Radian, s1.Angle(r.Lng.Length()) * s1.Radian} +} + +// Area returns the surface area of the Rect. +func (r Rect) Area() float64 { + if r.IsEmpty() { + return 0 + } + capDiff := math.Abs(math.Sin(r.Lat.Hi) - math.Sin(r.Lat.Lo)) + return r.Lng.Length() * capDiff +} + +// AddPoint increases the size of the rectangle to include the given point. +func (r Rect) AddPoint(ll LatLng) Rect { + if !ll.IsValid() { + return r + } + return Rect{ + Lat: r.Lat.AddPoint(ll.Lat.Radians()), + Lng: r.Lng.AddPoint(ll.Lng.Radians()), + } +} + +// expanded returns a rectangle that has been expanded by margin.Lat on each side +// in the latitude direction, and by margin.Lng on each side in the longitude +// direction. If either margin is negative, then it shrinks the rectangle on +// the corresponding sides instead. The resulting rectangle may be empty. +// +// The latitude-longitude space has the topology of a cylinder. Longitudes +// "wrap around" at +/-180 degrees, while latitudes are clamped to range [-90, 90]. +// This means that any expansion (positive or negative) of the full longitude range +// remains full (since the "rectangle" is actually a continuous band around the +// cylinder), while expansion of the full latitude range remains full only if the +// margin is positive. +// +// If either the latitude or longitude interval becomes empty after +// expansion by a negative margin, the result is empty. +// +// Note that if an expanded rectangle contains a pole, it may not contain +// all possible lat/lng representations of that pole, e.g., both points [π/2,0] +// and [π/2,1] represent the same pole, but they might not be contained by the +// same Rect. +// +// If you are trying to grow a rectangle by a certain distance on the +// sphere (e.g. 5km), refer to the ExpandedByDistance() C++ method implementation +// instead. +func (r Rect) expanded(margin LatLng) Rect { + lat := r.Lat.Expanded(margin.Lat.Radians()) + lng := r.Lng.Expanded(margin.Lng.Radians()) + + if lat.IsEmpty() || lng.IsEmpty() { + return EmptyRect() + } + + return Rect{ + Lat: lat.Intersection(validRectLatRange), + Lng: lng, + } +} + +func (r Rect) String() string { return fmt.Sprintf("[Lo%v, Hi%v]", r.Lo(), r.Hi()) } + +// PolarClosure returns the rectangle unmodified if it does not include either pole. +// If it includes either pole, PolarClosure returns an expansion of the rectangle along +// the longitudinal range to include all possible representations of the contained poles. +func (r Rect) PolarClosure() Rect { + if r.Lat.Lo == -math.Pi/2 || r.Lat.Hi == math.Pi/2 { + return Rect{r.Lat, s1.FullInterval()} + } + return r +} + +// Union returns the smallest Rect containing the union of this rectangle and the given rectangle. +func (r Rect) Union(other Rect) Rect { + return Rect{ + Lat: r.Lat.Union(other.Lat), + Lng: r.Lng.Union(other.Lng), + } +} + +// Intersection returns the smallest rectangle containing the intersection of +// this rectangle and the given rectangle. Note that the region of intersection +// may consist of two disjoint rectangles, in which case a single rectangle +// spanning both of them is returned. +func (r Rect) Intersection(other Rect) Rect { + lat := r.Lat.Intersection(other.Lat) + lng := r.Lng.Intersection(other.Lng) + + if lat.IsEmpty() || lng.IsEmpty() { + return EmptyRect() + } + return Rect{lat, lng} +} + +// Intersects reports whether this rectangle and the other have any points in common. +func (r Rect) Intersects(other Rect) bool { + return r.Lat.Intersects(other.Lat) && r.Lng.Intersects(other.Lng) +} + +// CapBound returns a cap that contains Rect. +func (r Rect) CapBound() Cap { + // We consider two possible bounding caps, one whose axis passes + // through the center of the lat-long rectangle and one whose axis + // is the north or south pole. We return the smaller of the two caps. + + if r.IsEmpty() { + return EmptyCap() + } + + var poleZ, poleAngle float64 + if r.Lat.Hi+r.Lat.Lo < 0 { + // South pole axis yields smaller cap. + poleZ = -1 + poleAngle = math.Pi/2 + r.Lat.Hi + } else { + poleZ = 1 + poleAngle = math.Pi/2 - r.Lat.Lo + } + poleCap := CapFromCenterAngle(Point{r3.Vector{0, 0, poleZ}}, s1.Angle(poleAngle)*s1.Radian) + + // For bounding rectangles that span 180 degrees or less in longitude, the + // maximum cap size is achieved at one of the rectangle vertices. For + // rectangles that are larger than 180 degrees, we punt and always return a + // bounding cap centered at one of the two poles. + if math.Remainder(r.Lng.Hi-r.Lng.Lo, 2*math.Pi) >= 0 && r.Lng.Hi-r.Lng.Lo < 2*math.Pi { + midCap := CapFromPoint(PointFromLatLng(r.Center())).AddPoint(PointFromLatLng(r.Lo())).AddPoint(PointFromLatLng(r.Hi())) + if midCap.Height() < poleCap.Height() { + return midCap + } + } + return poleCap +} + +// RectBound returns itself. +func (r Rect) RectBound() Rect { + return r +} + +// Contains reports whether this Rect contains the other Rect. +func (r Rect) Contains(other Rect) bool { + return r.Lat.ContainsInterval(other.Lat) && r.Lng.ContainsInterval(other.Lng) +} + +// ContainsCell reports whether the given Cell is contained by this Rect. +func (r Rect) ContainsCell(c Cell) bool { + // A latitude-longitude rectangle contains a cell if and only if it contains + // the cell's bounding rectangle. This test is exact from a mathematical + // point of view, assuming that the bounds returned by Cell.RectBound() + // are tight. However, note that there can be a loss of precision when + // converting between representations -- for example, if an s2.Cell is + // converted to a polygon, the polygon's bounding rectangle may not contain + // the cell's bounding rectangle. This has some slightly unexpected side + // effects; for instance, if one creates an s2.Polygon from an s2.Cell, the + // polygon will contain the cell, but the polygon's bounding box will not. + return r.Contains(c.RectBound()) +} + +// ContainsLatLng reports whether the given LatLng is within the Rect. +func (r Rect) ContainsLatLng(ll LatLng) bool { + if !ll.IsValid() { + return false + } + return r.Lat.Contains(ll.Lat.Radians()) && r.Lng.Contains(ll.Lng.Radians()) +} + +// ContainsPoint reports whether the given Point is within the Rect. +func (r Rect) ContainsPoint(p Point) bool { + return r.ContainsLatLng(LatLngFromPoint(p)) +} + +// CellUnionBound computes a covering of the Rect. +func (r Rect) CellUnionBound() []CellID { + return r.CapBound().CellUnionBound() +} + +// intersectsLatEdge reports whether the edge AB intersects the given edge of constant +// latitude. Requires the points to have unit length. +func intersectsLatEdge(a, b Point, lat s1.Angle, lng s1.Interval) bool { + // Unfortunately, lines of constant latitude are curves on + // the sphere. They can intersect a straight edge in 0, 1, or 2 points. + + // First, compute the normal to the plane AB that points vaguely north. + z := Point{a.PointCross(b).Normalize()} + if z.Z < 0 { + z = Point{z.Mul(-1)} + } + + // Extend this to an orthonormal frame (x,y,z) where x is the direction + // where the great circle through AB achieves its maximium latitude. + y := Point{z.PointCross(PointFromCoords(0, 0, 1)).Normalize()} + x := y.Cross(z.Vector) + + // Compute the angle "theta" from the x-axis (in the x-y plane defined + // above) where the great circle intersects the given line of latitude. + sinLat := math.Sin(float64(lat)) + if math.Abs(sinLat) >= x.Z { + // The great circle does not reach the given latitude. + return false + } + + cosTheta := sinLat / x.Z + sinTheta := math.Sqrt(1 - cosTheta*cosTheta) + theta := math.Atan2(sinTheta, cosTheta) + + // The candidate intersection points are located +/- theta in the x-y + // plane. For an intersection to be valid, we need to check that the + // intersection point is contained in the interior of the edge AB and + // also that it is contained within the given longitude interval "lng". + + // Compute the range of theta values spanned by the edge AB. + abTheta := s1.IntervalFromPointPair( + math.Atan2(a.Dot(y.Vector), a.Dot(x)), + math.Atan2(b.Dot(y.Vector), b.Dot(x))) + + if abTheta.Contains(theta) { + // Check if the intersection point is also in the given lng interval. + isect := x.Mul(cosTheta).Add(y.Mul(sinTheta)) + if lng.Contains(math.Atan2(isect.Y, isect.X)) { + return true + } + } + + if abTheta.Contains(-theta) { + // Check if the other intersection point is also in the given lng interval. + isect := x.Mul(cosTheta).Sub(y.Mul(sinTheta)) + if lng.Contains(math.Atan2(isect.Y, isect.X)) { + return true + } + } + return false +} + +// intersectsLngEdge reports whether the edge AB intersects the given edge of constant +// longitude. Requires the points to have unit length. +func intersectsLngEdge(a, b Point, lat r1.Interval, lng s1.Angle) bool { + // The nice thing about edges of constant longitude is that + // they are straight lines on the sphere (geodesics). + return CrossingSign(a, b, PointFromLatLng(LatLng{s1.Angle(lat.Lo), lng}), + PointFromLatLng(LatLng{s1.Angle(lat.Hi), lng})) == Cross +} + +// IntersectsCell reports whether this rectangle intersects the given cell. This is an +// exact test and may be fairly expensive. +func (r Rect) IntersectsCell(c Cell) bool { + // First we eliminate the cases where one region completely contains the + // other. Once these are disposed of, then the regions will intersect + // if and only if their boundaries intersect. + if r.IsEmpty() { + return false + } + if r.ContainsPoint(Point{c.id.rawPoint()}) { + return true + } + if c.ContainsPoint(PointFromLatLng(r.Center())) { + return true + } + + // Quick rejection test (not required for correctness). + if !r.Intersects(c.RectBound()) { + return false + } + + // Precompute the cell vertices as points and latitude-longitudes. We also + // check whether the Cell contains any corner of the rectangle, or + // vice-versa, since the edge-crossing tests only check the edge interiors. + vertices := [4]Point{} + latlngs := [4]LatLng{} + + for i := range vertices { + vertices[i] = c.Vertex(i) + latlngs[i] = LatLngFromPoint(vertices[i]) + if r.ContainsLatLng(latlngs[i]) { + return true + } + if c.ContainsPoint(PointFromLatLng(r.Vertex(i))) { + return true + } + } + + // Now check whether the boundaries intersect. Unfortunately, a + // latitude-longitude rectangle does not have straight edges: two edges + // are curved, and at least one of them is concave. + for i := range vertices { + edgeLng := s1.IntervalFromEndpoints(latlngs[i].Lng.Radians(), latlngs[(i+1)&3].Lng.Radians()) + if !r.Lng.Intersects(edgeLng) { + continue + } + + a := vertices[i] + b := vertices[(i+1)&3] + if edgeLng.Contains(r.Lng.Lo) && intersectsLngEdge(a, b, r.Lat, s1.Angle(r.Lng.Lo)) { + return true + } + if edgeLng.Contains(r.Lng.Hi) && intersectsLngEdge(a, b, r.Lat, s1.Angle(r.Lng.Hi)) { + return true + } + if intersectsLatEdge(a, b, s1.Angle(r.Lat.Lo), r.Lng) { + return true + } + if intersectsLatEdge(a, b, s1.Angle(r.Lat.Hi), r.Lng) { + return true + } + } + return false +} + +// Encode encodes the Rect. +func (r Rect) Encode(w io.Writer) error { + e := &encoder{w: w} + r.encode(e) + return e.err +} + +func (r Rect) encode(e *encoder) { + e.writeInt8(encodingVersion) + e.writeFloat64(r.Lat.Lo) + e.writeFloat64(r.Lat.Hi) + e.writeFloat64(r.Lng.Lo) + e.writeFloat64(r.Lng.Hi) +} + +// Decode decodes a rectangle. +func (r *Rect) Decode(rd io.Reader) error { + d := &decoder{r: asByteReader(rd)} + r.decode(d) + return d.err +} + +func (r *Rect) decode(d *decoder) { + if version := d.readUint8(); int(version) != int(encodingVersion) && d.err == nil { + d.err = fmt.Errorf("can't decode version %d; my version: %d", version, encodingVersion) + return + } + r.Lat.Lo = d.readFloat64() + r.Lat.Hi = d.readFloat64() + r.Lng.Lo = d.readFloat64() + r.Lng.Hi = d.readFloat64() + return +} + +// DistanceToLatLng returns the minimum distance (measured along the surface of the sphere) +// from a given point to the rectangle (both its boundary and its interior). +// If r is empty, the result is meaningless. +// The latlng must be valid. +func (r Rect) DistanceToLatLng(ll LatLng) s1.Angle { + if r.Lng.Contains(float64(ll.Lng)) { + return maxAngle(0, ll.Lat-s1.Angle(r.Lat.Hi), s1.Angle(r.Lat.Lo)-ll.Lat) + } + + i := s1.IntervalFromEndpoints(r.Lng.Hi, r.Lng.ComplementCenter()) + rectLng := r.Lng.Lo + if i.Contains(float64(ll.Lng)) { + rectLng = r.Lng.Hi + } + + lo := LatLng{s1.Angle(r.Lat.Lo) * s1.Radian, s1.Angle(rectLng) * s1.Radian} + hi := LatLng{s1.Angle(r.Lat.Hi) * s1.Radian, s1.Angle(rectLng) * s1.Radian} + return DistanceFromSegment(PointFromLatLng(ll), PointFromLatLng(lo), PointFromLatLng(hi)) +} + +// DirectedHausdorffDistance returns the directed Hausdorff distance (measured along the +// surface of the sphere) to the given Rect. The directed Hausdorff +// distance from rectangle A to rectangle B is given by +// h(A, B) = max_{p in A} min_{q in B} d(p, q). +func (r Rect) DirectedHausdorffDistance(other Rect) s1.Angle { + if r.IsEmpty() { + return 0 * s1.Radian + } + if other.IsEmpty() { + return math.Pi * s1.Radian + } + + lng := r.Lng.DirectedHausdorffDistance(other.Lng) + return directedHausdorffDistance(lng, r.Lat, other.Lat) +} + +// HausdorffDistance returns the undirected Hausdorff distance (measured along the +// surface of the sphere) to the given Rect. +// The Hausdorff distance between rectangle A and rectangle B is given by +// H(A, B) = max{h(A, B), h(B, A)}. +func (r Rect) HausdorffDistance(other Rect) s1.Angle { + return maxAngle(r.DirectedHausdorffDistance(other), + other.DirectedHausdorffDistance(r)) +} + +// ApproxEqual reports whether the latitude and longitude intervals of the two rectangles +// are the same up to a small tolerance. +func (r Rect) ApproxEqual(other Rect) bool { + return r.Lat.ApproxEqual(other.Lat) && r.Lng.ApproxEqual(other.Lng) +} + +// directedHausdorffDistance returns the directed Hausdorff distance +// from one longitudinal edge spanning latitude range 'a' to the other +// longitudinal edge spanning latitude range 'b', with their longitudinal +// difference given by 'lngDiff'. +func directedHausdorffDistance(lngDiff s1.Angle, a, b r1.Interval) s1.Angle { + // By symmetry, we can assume a's longitude is 0 and b's longitude is + // lngDiff. Call b's two endpoints bLo and bHi. Let H be the hemisphere + // containing a and delimited by the longitude line of b. The Voronoi diagram + // of b on H has three edges (portions of great circles) all orthogonal to b + // and meeting at bLo cross bHi. + // E1: (bLo, bLo cross bHi) + // E2: (bHi, bLo cross bHi) + // E3: (-bMid, bLo cross bHi), where bMid is the midpoint of b + // + // They subdivide H into three Voronoi regions. Depending on how longitude 0 + // (which contains edge a) intersects these regions, we distinguish two cases: + // Case 1: it intersects three regions. This occurs when lngDiff <= π/2. + // Case 2: it intersects only two regions. This occurs when lngDiff > π/2. + // + // In the first case, the directed Hausdorff distance to edge b can only be + // realized by the following points on a: + // A1: two endpoints of a. + // A2: intersection of a with the equator, if b also intersects the equator. + // + // In the second case, the directed Hausdorff distance to edge b can only be + // realized by the following points on a: + // B1: two endpoints of a. + // B2: intersection of a with E3 + // B3: farthest point from bLo to the interior of D, and farthest point from + // bHi to the interior of U, if any, where D (resp. U) is the portion + // of edge a below (resp. above) the intersection point from B2. + + if lngDiff < 0 { + panic("impossible: negative lngDiff") + } + if lngDiff > math.Pi { + panic("impossible: lngDiff > Pi") + } + + if lngDiff == 0 { + return s1.Angle(a.DirectedHausdorffDistance(b)) + } + + // Assumed longitude of b. + bLng := lngDiff + // Two endpoints of b. + bLo := PointFromLatLng(LatLng{s1.Angle(b.Lo), bLng}) + bHi := PointFromLatLng(LatLng{s1.Angle(b.Hi), bLng}) + + // Cases A1 and B1. + aLo := PointFromLatLng(LatLng{s1.Angle(a.Lo), 0}) + aHi := PointFromLatLng(LatLng{s1.Angle(a.Hi), 0}) + maxDistance := maxAngle( + DistanceFromSegment(aLo, bLo, bHi), + DistanceFromSegment(aHi, bLo, bHi)) + + if lngDiff <= math.Pi/2 { + // Case A2. + if a.Contains(0) && b.Contains(0) { + maxDistance = maxAngle(maxDistance, lngDiff) + } + return maxDistance + } + + // Case B2. + p := bisectorIntersection(b, bLng) + pLat := LatLngFromPoint(p).Lat + if a.Contains(float64(pLat)) { + maxDistance = maxAngle(maxDistance, p.Angle(bLo.Vector)) + } + + // Case B3. + if pLat > s1.Angle(a.Lo) { + intDist, ok := interiorMaxDistance(r1.Interval{a.Lo, math.Min(float64(pLat), a.Hi)}, bLo) + if ok { + maxDistance = maxAngle(maxDistance, intDist) + } + } + if pLat < s1.Angle(a.Hi) { + intDist, ok := interiorMaxDistance(r1.Interval{math.Max(float64(pLat), a.Lo), a.Hi}, bHi) + if ok { + maxDistance = maxAngle(maxDistance, intDist) + } + } + + return maxDistance +} + +// interiorMaxDistance returns the max distance from a point b to the segment spanning latitude range +// aLat on longitude 0 if the max occurs in the interior of aLat. Otherwise, returns (0, false). +func interiorMaxDistance(aLat r1.Interval, b Point) (a s1.Angle, ok bool) { + // Longitude 0 is in the y=0 plane. b.X >= 0 implies that the maximum + // does not occur in the interior of aLat. + if aLat.IsEmpty() || b.X >= 0 { + return 0, false + } + + // Project b to the y=0 plane. The antipodal of the normalized projection is + // the point at which the maxium distance from b occurs, if it is contained + // in aLat. + intersectionPoint := PointFromCoords(-b.X, 0, -b.Z) + if !aLat.InteriorContains(float64(LatLngFromPoint(intersectionPoint).Lat)) { + return 0, false + } + return b.Angle(intersectionPoint.Vector), true +} + +// bisectorIntersection return the intersection of longitude 0 with the bisector of an edge +// on longitude 'lng' and spanning latitude range 'lat'. +func bisectorIntersection(lat r1.Interval, lng s1.Angle) Point { + lng = s1.Angle(math.Abs(float64(lng))) + latCenter := s1.Angle(lat.Center()) + + // A vector orthogonal to the bisector of the given longitudinal edge. + orthoBisector := LatLng{latCenter - math.Pi/2, lng} + if latCenter < 0 { + orthoBisector = LatLng{-latCenter - math.Pi/2, lng - math.Pi} + } + + // A vector orthogonal to longitude 0. + orthoLng := Point{r3.Vector{0, -1, 0}} + + 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: +// - Get*Distance, Vertex, InteriorContains(LatLng|Rect|Point) |