1 // Copyright 2014 Google Inc. All rights reserved. 2 // 3 // Licensed under the Apache License, Version 2.0 (the "License"); 4 // you may not use this file except in compliance with the License. 5 // You may obtain a copy of the License at 6 // 7 // http://www.apache.org/licenses/LICENSE-2.0 8 // 9 // Unless required by applicable law or agreed to in writing, software 10 // distributed under the License is distributed on an "AS IS" BASIS, 11 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 12 // See the License for the specific language governing permissions and 13 // limitations under the License. 14 15 package s2 16 17 import ( 18 "fmt" 19 "io" 20 "math" 21 22 "github.com/golang/geo/r1" 23 "github.com/golang/geo/r3" 24 "github.com/golang/geo/s1" 25 ) 26 27 // Rect represents a closed latitude-longitude rectangle. 28 type Rect struct { 29 Lat r1.Interval 30 Lng s1.Interval 31 } 32 33 var ( 34 validRectLatRange = r1.Interval{-math.Pi / 2, math.Pi / 2} 35 validRectLngRange = s1.FullInterval() 36 ) 37 38 // EmptyRect returns the empty rectangle. 39 func EmptyRect() Rect { return Rect{r1.EmptyInterval(), s1.EmptyInterval()} } 40 41 // FullRect returns the full rectangle. 42 func FullRect() Rect { return Rect{validRectLatRange, validRectLngRange} } 43 44 // RectFromLatLng constructs a rectangle containing a single point p. 45 func RectFromLatLng(p LatLng) Rect { 46 return Rect{ 47 Lat: r1.Interval{p.Lat.Radians(), p.Lat.Radians()}, 48 Lng: s1.Interval{p.Lng.Radians(), p.Lng.Radians()}, 49 } 50 } 51 52 // RectFromCenterSize constructs a rectangle with the given size and center. 53 // center needs to be normalized, but size does not. The latitude 54 // interval of the result is clamped to [-90,90] degrees, and the longitude 55 // interval of the result is FullRect() if and only if the longitude size is 56 // 360 degrees or more. 57 // 58 // Examples of clamping (in degrees): 59 // 60 // center=(80,170), size=(40,60) -> lat=[60,90], lng=[140,-160] 61 // center=(10,40), size=(210,400) -> lat=[-90,90], lng=[-180,180] 62 // center=(-90,180), size=(20,50) -> lat=[-90,-80], lng=[155,-155] 63 func RectFromCenterSize(center, size LatLng) Rect { 64 half := LatLng{size.Lat / 2, size.Lng / 2} 65 return RectFromLatLng(center).expanded(half) 66 } 67 68 // IsValid returns true iff the rectangle is valid. 69 // This requires Lat ⊆ [-π/2,π/2] and Lng ⊆ [-π,π], and Lat = ∅ iff Lng = ∅ 70 func (r Rect) IsValid() bool { 71 return math.Abs(r.Lat.Lo) <= math.Pi/2 && 72 math.Abs(r.Lat.Hi) <= math.Pi/2 && 73 r.Lng.IsValid() && 74 r.Lat.IsEmpty() == r.Lng.IsEmpty() 75 } 76 77 // IsEmpty reports whether the rectangle is empty. 78 func (r Rect) IsEmpty() bool { return r.Lat.IsEmpty() } 79 80 // IsFull reports whether the rectangle is full. 81 func (r Rect) IsFull() bool { return r.Lat.Equal(validRectLatRange) && r.Lng.IsFull() } 82 83 // IsPoint reports whether the rectangle is a single point. 84 func (r Rect) IsPoint() bool { return r.Lat.Lo == r.Lat.Hi && r.Lng.Lo == r.Lng.Hi } 85 86 // Vertex returns the i-th vertex of the rectangle (i = 0,1,2,3) in CCW order 87 // (lower left, lower right, upper right, upper left). 88 func (r Rect) Vertex(i int) LatLng { 89 var lat, lng float64 90 91 switch i { 92 case 0: 93 lat = r.Lat.Lo 94 lng = r.Lng.Lo 95 case 1: 96 lat = r.Lat.Lo 97 lng = r.Lng.Hi 98 case 2: 99 lat = r.Lat.Hi 100 lng = r.Lng.Hi 101 case 3: 102 lat = r.Lat.Hi 103 lng = r.Lng.Lo 104 } 105 return LatLng{s1.Angle(lat) * s1.Radian, s1.Angle(lng) * s1.Radian} 106 } 107 108 // Lo returns one corner of the rectangle. 109 func (r Rect) Lo() LatLng { 110 return LatLng{s1.Angle(r.Lat.Lo) * s1.Radian, s1.Angle(r.Lng.Lo) * s1.Radian} 111 } 112 113 // Hi returns the other corner of the rectangle. 114 func (r Rect) Hi() LatLng { 115 return LatLng{s1.Angle(r.Lat.Hi) * s1.Radian, s1.Angle(r.Lng.Hi) * s1.Radian} 116 } 117 118 // Center returns the center of the rectangle. 119 func (r Rect) Center() LatLng { 120 return LatLng{s1.Angle(r.Lat.Center()) * s1.Radian, s1.Angle(r.Lng.Center()) * s1.Radian} 121 } 122 123 // Size returns the size of the Rect. 124 func (r Rect) Size() LatLng { 125 return LatLng{s1.Angle(r.Lat.Length()) * s1.Radian, s1.Angle(r.Lng.Length()) * s1.Radian} 126 } 127 128 // Area returns the surface area of the Rect. 129 func (r Rect) Area() float64 { 130 if r.IsEmpty() { 131 return 0 132 } 133 capDiff := math.Abs(math.Sin(r.Lat.Hi) - math.Sin(r.Lat.Lo)) 134 return r.Lng.Length() * capDiff 135 } 136 137 // AddPoint increases the size of the rectangle to include the given point. 138 func (r Rect) AddPoint(ll LatLng) Rect { 139 if !ll.IsValid() { 140 return r 141 } 142 return Rect{ 143 Lat: r.Lat.AddPoint(ll.Lat.Radians()), 144 Lng: r.Lng.AddPoint(ll.Lng.Radians()), 145 } 146 } 147 148 // expanded returns a rectangle that has been expanded by margin.Lat on each side 149 // in the latitude direction, and by margin.Lng on each side in the longitude 150 // direction. If either margin is negative, then it shrinks the rectangle on 151 // the corresponding sides instead. The resulting rectangle may be empty. 152 // 153 // The latitude-longitude space has the topology of a cylinder. Longitudes 154 // "wrap around" at +/-180 degrees, while latitudes are clamped to range [-90, 90]. 155 // This means that any expansion (positive or negative) of the full longitude range 156 // remains full (since the "rectangle" is actually a continuous band around the 157 // cylinder), while expansion of the full latitude range remains full only if the 158 // margin is positive. 159 // 160 // If either the latitude or longitude interval becomes empty after 161 // expansion by a negative margin, the result is empty. 162 // 163 // Note that if an expanded rectangle contains a pole, it may not contain 164 // all possible lat/lng representations of that pole, e.g., both points [π/2,0] 165 // and [π/2,1] represent the same pole, but they might not be contained by the 166 // same Rect. 167 // 168 // If you are trying to grow a rectangle by a certain distance on the 169 // sphere (e.g. 5km), refer to the ExpandedByDistance() C++ method implementation 170 // instead. 171 func (r Rect) expanded(margin LatLng) Rect { 172 lat := r.Lat.Expanded(margin.Lat.Radians()) 173 lng := r.Lng.Expanded(margin.Lng.Radians()) 174 175 if lat.IsEmpty() || lng.IsEmpty() { 176 return EmptyRect() 177 } 178 179 return Rect{ 180 Lat: lat.Intersection(validRectLatRange), 181 Lng: lng, 182 } 183 } 184 185 func (r Rect) String() string { return fmt.Sprintf("[Lo%v, Hi%v]", r.Lo(), r.Hi()) } 186 187 // PolarClosure returns the rectangle unmodified if it does not include either pole. 188 // If it includes either pole, PolarClosure returns an expansion of the rectangle along 189 // the longitudinal range to include all possible representations of the contained poles. 190 func (r Rect) PolarClosure() Rect { 191 if r.Lat.Lo == -math.Pi/2 || r.Lat.Hi == math.Pi/2 { 192 return Rect{r.Lat, s1.FullInterval()} 193 } 194 return r 195 } 196 197 // Union returns the smallest Rect containing the union of this rectangle and the given rectangle. 198 func (r Rect) Union(other Rect) Rect { 199 return Rect{ 200 Lat: r.Lat.Union(other.Lat), 201 Lng: r.Lng.Union(other.Lng), 202 } 203 } 204 205 // Intersection returns the smallest rectangle containing the intersection of 206 // this rectangle and the given rectangle. Note that the region of intersection 207 // may consist of two disjoint rectangles, in which case a single rectangle 208 // spanning both of them is returned. 209 func (r Rect) Intersection(other Rect) Rect { 210 lat := r.Lat.Intersection(other.Lat) 211 lng := r.Lng.Intersection(other.Lng) 212 213 if lat.IsEmpty() || lng.IsEmpty() { 214 return EmptyRect() 215 } 216 return Rect{lat, lng} 217 } 218 219 // Intersects reports whether this rectangle and the other have any points in common. 220 func (r Rect) Intersects(other Rect) bool { 221 return r.Lat.Intersects(other.Lat) && r.Lng.Intersects(other.Lng) 222 } 223 224 // CapBound returns a cap that contains Rect. 225 func (r Rect) CapBound() Cap { 226 // We consider two possible bounding caps, one whose axis passes 227 // through the center of the lat-long rectangle and one whose axis 228 // is the north or south pole. We return the smaller of the two caps. 229 230 if r.IsEmpty() { 231 return EmptyCap() 232 } 233 234 var poleZ, poleAngle float64 235 if r.Lat.Hi+r.Lat.Lo < 0 { 236 // South pole axis yields smaller cap. 237 poleZ = -1 238 poleAngle = math.Pi/2 + r.Lat.Hi 239 } else { 240 poleZ = 1 241 poleAngle = math.Pi/2 - r.Lat.Lo 242 } 243 poleCap := CapFromCenterAngle(Point{r3.Vector{0, 0, poleZ}}, s1.Angle(poleAngle)*s1.Radian) 244 245 // For bounding rectangles that span 180 degrees or less in longitude, the 246 // maximum cap size is achieved at one of the rectangle vertices. For 247 // rectangles that are larger than 180 degrees, we punt and always return a 248 // bounding cap centered at one of the two poles. 249 if math.Remainder(r.Lng.Hi-r.Lng.Lo, 2*math.Pi) >= 0 && r.Lng.Hi-r.Lng.Lo < 2*math.Pi { 250 midCap := CapFromPoint(PointFromLatLng(r.Center())).AddPoint(PointFromLatLng(r.Lo())).AddPoint(PointFromLatLng(r.Hi())) 251 if midCap.Height() < poleCap.Height() { 252 return midCap 253 } 254 } 255 return poleCap 256 } 257 258 // RectBound returns itself. 259 func (r Rect) RectBound() Rect { 260 return r 261 } 262 263 // Contains reports whether this Rect contains the other Rect. 264 func (r Rect) Contains(other Rect) bool { 265 return r.Lat.ContainsInterval(other.Lat) && r.Lng.ContainsInterval(other.Lng) 266 } 267 268 // ContainsCell reports whether the given Cell is contained by this Rect. 269 func (r Rect) ContainsCell(c Cell) bool { 270 // A latitude-longitude rectangle contains a cell if and only if it contains 271 // the cell's bounding rectangle. This test is exact from a mathematical 272 // point of view, assuming that the bounds returned by Cell.RectBound() 273 // are tight. However, note that there can be a loss of precision when 274 // converting between representations -- for example, if an s2.Cell is 275 // converted to a polygon, the polygon's bounding rectangle may not contain 276 // the cell's bounding rectangle. This has some slightly unexpected side 277 // effects; for instance, if one creates an s2.Polygon from an s2.Cell, the 278 // polygon will contain the cell, but the polygon's bounding box will not. 279 return r.Contains(c.RectBound()) 280 } 281 282 // ContainsLatLng reports whether the given LatLng is within the Rect. 283 func (r Rect) ContainsLatLng(ll LatLng) bool { 284 if !ll.IsValid() { 285 return false 286 } 287 return r.Lat.Contains(ll.Lat.Radians()) && r.Lng.Contains(ll.Lng.Radians()) 288 } 289 290 // ContainsPoint reports whether the given Point is within the Rect. 291 func (r Rect) ContainsPoint(p Point) bool { 292 return r.ContainsLatLng(LatLngFromPoint(p)) 293 } 294 295 // CellUnionBound computes a covering of the Rect. 296 func (r Rect) CellUnionBound() []CellID { 297 return r.CapBound().CellUnionBound() 298 } 299 300 // intersectsLatEdge reports whether the edge AB intersects the given edge of constant 301 // latitude. Requires the points to have unit length. 302 func intersectsLatEdge(a, b Point, lat s1.Angle, lng s1.Interval) bool { 303 // Unfortunately, lines of constant latitude are curves on 304 // the sphere. They can intersect a straight edge in 0, 1, or 2 points. 305 306 // First, compute the normal to the plane AB that points vaguely north. 307 z := Point{a.PointCross(b).Normalize()} 308 if z.Z < 0 { 309 z = Point{z.Mul(-1)} 310 } 311 312 // Extend this to an orthonormal frame (x,y,z) where x is the direction 313 // where the great circle through AB achieves its maximium latitude. 314 y := Point{z.PointCross(PointFromCoords(0, 0, 1)).Normalize()} 315 x := y.Cross(z.Vector) 316 317 // Compute the angle "theta" from the x-axis (in the x-y plane defined 318 // above) where the great circle intersects the given line of latitude. 319 sinLat := math.Sin(float64(lat)) 320 if math.Abs(sinLat) >= x.Z { 321 // The great circle does not reach the given latitude. 322 return false 323 } 324 325 cosTheta := sinLat / x.Z 326 sinTheta := math.Sqrt(1 - cosTheta*cosTheta) 327 theta := math.Atan2(sinTheta, cosTheta) 328 329 // The candidate intersection points are located +/- theta in the x-y 330 // plane. For an intersection to be valid, we need to check that the 331 // intersection point is contained in the interior of the edge AB and 332 // also that it is contained within the given longitude interval "lng". 333 334 // Compute the range of theta values spanned by the edge AB. 335 abTheta := s1.IntervalFromPointPair( 336 math.Atan2(a.Dot(y.Vector), a.Dot(x)), 337 math.Atan2(b.Dot(y.Vector), b.Dot(x))) 338 339 if abTheta.Contains(theta) { 340 // Check if the intersection point is also in the given lng interval. 341 isect := x.Mul(cosTheta).Add(y.Mul(sinTheta)) 342 if lng.Contains(math.Atan2(isect.Y, isect.X)) { 343 return true 344 } 345 } 346 347 if abTheta.Contains(-theta) { 348 // Check if the other intersection point is also in the given lng interval. 349 isect := x.Mul(cosTheta).Sub(y.Mul(sinTheta)) 350 if lng.Contains(math.Atan2(isect.Y, isect.X)) { 351 return true 352 } 353 } 354 return false 355 } 356 357 // intersectsLngEdge reports whether the edge AB intersects the given edge of constant 358 // longitude. Requires the points to have unit length. 359 func intersectsLngEdge(a, b Point, lat r1.Interval, lng s1.Angle) bool { 360 // The nice thing about edges of constant longitude is that 361 // they are straight lines on the sphere (geodesics). 362 return CrossingSign(a, b, PointFromLatLng(LatLng{s1.Angle(lat.Lo), lng}), 363 PointFromLatLng(LatLng{s1.Angle(lat.Hi), lng})) == Cross 364 } 365 366 // IntersectsCell reports whether this rectangle intersects the given cell. This is an 367 // exact test and may be fairly expensive. 368 func (r Rect) IntersectsCell(c Cell) bool { 369 // First we eliminate the cases where one region completely contains the 370 // other. Once these are disposed of, then the regions will intersect 371 // if and only if their boundaries intersect. 372 if r.IsEmpty() { 373 return false 374 } 375 if r.ContainsPoint(Point{c.id.rawPoint()}) { 376 return true 377 } 378 if c.ContainsPoint(PointFromLatLng(r.Center())) { 379 return true 380 } 381 382 // Quick rejection test (not required for correctness). 383 if !r.Intersects(c.RectBound()) { 384 return false 385 } 386 387 // Precompute the cell vertices as points and latitude-longitudes. We also 388 // check whether the Cell contains any corner of the rectangle, or 389 // vice-versa, since the edge-crossing tests only check the edge interiors. 390 vertices := [4]Point{} 391 latlngs := [4]LatLng{} 392 393 for i := range vertices { 394 vertices[i] = c.Vertex(i) 395 latlngs[i] = LatLngFromPoint(vertices[i]) 396 if r.ContainsLatLng(latlngs[i]) { 397 return true 398 } 399 if c.ContainsPoint(PointFromLatLng(r.Vertex(i))) { 400 return true 401 } 402 } 403 404 // Now check whether the boundaries intersect. Unfortunately, a 405 // latitude-longitude rectangle does not have straight edges: two edges 406 // are curved, and at least one of them is concave. 407 for i := range vertices { 408 edgeLng := s1.IntervalFromEndpoints(latlngs[i].Lng.Radians(), latlngs[(i+1)&3].Lng.Radians()) 409 if !r.Lng.Intersects(edgeLng) { 410 continue 411 } 412 413 a := vertices[i] 414 b := vertices[(i+1)&3] 415 if edgeLng.Contains(r.Lng.Lo) && intersectsLngEdge(a, b, r.Lat, s1.Angle(r.Lng.Lo)) { 416 return true 417 } 418 if edgeLng.Contains(r.Lng.Hi) && intersectsLngEdge(a, b, r.Lat, s1.Angle(r.Lng.Hi)) { 419 return true 420 } 421 if intersectsLatEdge(a, b, s1.Angle(r.Lat.Lo), r.Lng) { 422 return true 423 } 424 if intersectsLatEdge(a, b, s1.Angle(r.Lat.Hi), r.Lng) { 425 return true 426 } 427 } 428 return false 429 } 430 431 // Encode encodes the Rect. 432 func (r Rect) Encode(w io.Writer) error { 433 e := &encoder{w: w} 434 r.encode(e) 435 return e.err 436 } 437 438 func (r Rect) encode(e *encoder) { 439 e.writeInt8(encodingVersion) 440 e.writeFloat64(r.Lat.Lo) 441 e.writeFloat64(r.Lat.Hi) 442 e.writeFloat64(r.Lng.Lo) 443 e.writeFloat64(r.Lng.Hi) 444 } 445 446 // Decode decodes a rectangle. 447 func (r *Rect) Decode(rd io.Reader) error { 448 d := &decoder{r: asByteReader(rd)} 449 r.decode(d) 450 return d.err 451 } 452 453 func (r *Rect) decode(d *decoder) { 454 if version := d.readUint8(); int8(version) != encodingVersion && d.err == nil { 455 d.err = fmt.Errorf("can't decode version %d; my version: %d", version, encodingVersion) 456 return 457 } 458 r.Lat.Lo = d.readFloat64() 459 r.Lat.Hi = d.readFloat64() 460 r.Lng.Lo = d.readFloat64() 461 r.Lng.Hi = d.readFloat64() 462 return 463 } 464 465 // DistanceToLatLng returns the minimum distance (measured along the surface of the sphere) 466 // from a given point to the rectangle (both its boundary and its interior). 467 // If r is empty, the result is meaningless. 468 // The latlng must be valid. 469 func (r Rect) DistanceToLatLng(ll LatLng) s1.Angle { 470 if r.Lng.Contains(float64(ll.Lng)) { 471 return maxAngle(0, ll.Lat-s1.Angle(r.Lat.Hi), s1.Angle(r.Lat.Lo)-ll.Lat) 472 } 473 474 i := s1.IntervalFromEndpoints(r.Lng.Hi, r.Lng.ComplementCenter()) 475 rectLng := r.Lng.Lo 476 if i.Contains(float64(ll.Lng)) { 477 rectLng = r.Lng.Hi 478 } 479 480 lo := LatLng{s1.Angle(r.Lat.Lo) * s1.Radian, s1.Angle(rectLng) * s1.Radian} 481 hi := LatLng{s1.Angle(r.Lat.Hi) * s1.Radian, s1.Angle(rectLng) * s1.Radian} 482 return DistanceFromSegment(PointFromLatLng(ll), PointFromLatLng(lo), PointFromLatLng(hi)) 483 } 484 485 // DirectedHausdorffDistance returns the directed Hausdorff distance (measured along the 486 // surface of the sphere) to the given Rect. The directed Hausdorff 487 // distance from rectangle A to rectangle B is given by 488 // 489 // h(A, B) = max_{p in A} min_{q in B} d(p, q). 490 func (r Rect) DirectedHausdorffDistance(other Rect) s1.Angle { 491 if r.IsEmpty() { 492 return 0 * s1.Radian 493 } 494 if other.IsEmpty() { 495 return math.Pi * s1.Radian 496 } 497 498 lng := r.Lng.DirectedHausdorffDistance(other.Lng) 499 return directedHausdorffDistance(lng, r.Lat, other.Lat) 500 } 501 502 // HausdorffDistance returns the undirected Hausdorff distance (measured along the 503 // surface of the sphere) to the given Rect. 504 // The Hausdorff distance between rectangle A and rectangle B is given by 505 // 506 // H(A, B) = max{h(A, B), h(B, A)}. 507 func (r Rect) HausdorffDistance(other Rect) s1.Angle { 508 return maxAngle(r.DirectedHausdorffDistance(other), 509 other.DirectedHausdorffDistance(r)) 510 } 511 512 // ApproxEqual reports whether the latitude and longitude intervals of the two rectangles 513 // are the same up to a small tolerance. 514 func (r Rect) ApproxEqual(other Rect) bool { 515 return r.Lat.ApproxEqual(other.Lat) && r.Lng.ApproxEqual(other.Lng) 516 } 517 518 // directedHausdorffDistance returns the directed Hausdorff distance 519 // from one longitudinal edge spanning latitude range 'a' to the other 520 // longitudinal edge spanning latitude range 'b', with their longitudinal 521 // difference given by 'lngDiff'. 522 func directedHausdorffDistance(lngDiff s1.Angle, a, b r1.Interval) s1.Angle { 523 // By symmetry, we can assume a's longitude is 0 and b's longitude is 524 // lngDiff. Call b's two endpoints bLo and bHi. Let H be the hemisphere 525 // containing a and delimited by the longitude line of b. The Voronoi diagram 526 // of b on H has three edges (portions of great circles) all orthogonal to b 527 // and meeting at bLo cross bHi. 528 // E1: (bLo, bLo cross bHi) 529 // E2: (bHi, bLo cross bHi) 530 // E3: (-bMid, bLo cross bHi), where bMid is the midpoint of b 531 // 532 // They subdivide H into three Voronoi regions. Depending on how longitude 0 533 // (which contains edge a) intersects these regions, we distinguish two cases: 534 // Case 1: it intersects three regions. This occurs when lngDiff <= π/2. 535 // Case 2: it intersects only two regions. This occurs when lngDiff > π/2. 536 // 537 // In the first case, the directed Hausdorff distance to edge b can only be 538 // realized by the following points on a: 539 // A1: two endpoints of a. 540 // A2: intersection of a with the equator, if b also intersects the equator. 541 // 542 // In the second case, the directed Hausdorff distance to edge b can only be 543 // realized by the following points on a: 544 // B1: two endpoints of a. 545 // B2: intersection of a with E3 546 // B3: farthest point from bLo to the interior of D, and farthest point from 547 // bHi to the interior of U, if any, where D (resp. U) is the portion 548 // of edge a below (resp. above) the intersection point from B2. 549 550 if lngDiff < 0 { 551 panic("impossible: negative lngDiff") 552 } 553 if lngDiff > math.Pi { 554 panic("impossible: lngDiff > Pi") 555 } 556 557 if lngDiff == 0 { 558 return s1.Angle(a.DirectedHausdorffDistance(b)) 559 } 560 561 // Assumed longitude of b. 562 bLng := lngDiff 563 // Two endpoints of b. 564 bLo := PointFromLatLng(LatLng{s1.Angle(b.Lo), bLng}) 565 bHi := PointFromLatLng(LatLng{s1.Angle(b.Hi), bLng}) 566 567 // Cases A1 and B1. 568 aLo := PointFromLatLng(LatLng{s1.Angle(a.Lo), 0}) 569 aHi := PointFromLatLng(LatLng{s1.Angle(a.Hi), 0}) 570 maxDistance := maxAngle( 571 DistanceFromSegment(aLo, bLo, bHi), 572 DistanceFromSegment(aHi, bLo, bHi)) 573 574 if lngDiff <= math.Pi/2 { 575 // Case A2. 576 if a.Contains(0) && b.Contains(0) { 577 maxDistance = maxAngle(maxDistance, lngDiff) 578 } 579 return maxDistance 580 } 581 582 // Case B2. 583 p := bisectorIntersection(b, bLng) 584 pLat := LatLngFromPoint(p).Lat 585 if a.Contains(float64(pLat)) { 586 maxDistance = maxAngle(maxDistance, p.Angle(bLo.Vector)) 587 } 588 589 // Case B3. 590 if pLat > s1.Angle(a.Lo) { 591 intDist, ok := interiorMaxDistance(r1.Interval{a.Lo, math.Min(float64(pLat), a.Hi)}, bLo) 592 if ok { 593 maxDistance = maxAngle(maxDistance, intDist) 594 } 595 } 596 if pLat < s1.Angle(a.Hi) { 597 intDist, ok := interiorMaxDistance(r1.Interval{math.Max(float64(pLat), a.Lo), a.Hi}, bHi) 598 if ok { 599 maxDistance = maxAngle(maxDistance, intDist) 600 } 601 } 602 603 return maxDistance 604 } 605 606 // interiorMaxDistance returns the max distance from a point b to the segment spanning latitude range 607 // aLat on longitude 0 if the max occurs in the interior of aLat. Otherwise, returns (0, false). 608 func interiorMaxDistance(aLat r1.Interval, b Point) (a s1.Angle, ok bool) { 609 // Longitude 0 is in the y=0 plane. b.X >= 0 implies that the maximum 610 // does not occur in the interior of aLat. 611 if aLat.IsEmpty() || b.X >= 0 { 612 return 0, false 613 } 614 615 // Project b to the y=0 plane. The antipodal of the normalized projection is 616 // the point at which the maximum distance from b occurs, if it is contained 617 // in aLat. 618 intersectionPoint := PointFromCoords(-b.X, 0, -b.Z) 619 if !aLat.InteriorContains(float64(LatLngFromPoint(intersectionPoint).Lat)) { 620 return 0, false 621 } 622 return b.Angle(intersectionPoint.Vector), true 623 } 624 625 // bisectorIntersection return the intersection of longitude 0 with the bisector of an edge 626 // on longitude 'lng' and spanning latitude range 'lat'. 627 func bisectorIntersection(lat r1.Interval, lng s1.Angle) Point { 628 lng = s1.Angle(math.Abs(float64(lng))) 629 latCenter := s1.Angle(lat.Center()) 630 631 // A vector orthogonal to the bisector of the given longitudinal edge. 632 orthoBisector := LatLng{latCenter - math.Pi/2, lng} 633 if latCenter < 0 { 634 orthoBisector = LatLng{-latCenter - math.Pi/2, lng - math.Pi} 635 } 636 637 // A vector orthogonal to longitude 0. 638 orthoLng := Point{r3.Vector{0, -1, 0}} 639 640 return orthoLng.PointCross(PointFromLatLng(orthoBisector)) 641 } 642 643 // Centroid returns the true centroid of the given Rect multiplied by its 644 // surface area. The result is not unit length, so you may want to normalize it. 645 // Note that in general the centroid is *not* at the center of the rectangle, and 646 // in fact it may not even be contained by the rectangle. (It is the "center of 647 // mass" of the rectangle viewed as subset of the unit sphere, i.e. it is the 648 // point in space about which this curved shape would rotate.) 649 // 650 // The reason for multiplying the result by the rectangle area is to make it 651 // easier to compute the centroid of more complicated shapes. The centroid 652 // of a union of disjoint regions can be computed simply by adding their 653 // Centroid results. 654 func (r Rect) Centroid() Point { 655 // When a sphere is divided into slices of constant thickness by a set 656 // of parallel planes, all slices have the same surface area. This 657 // implies that the z-component of the centroid is simply the midpoint 658 // of the z-interval spanned by the Rect. 659 // 660 // Similarly, it is easy to see that the (x,y) of the centroid lies in 661 // the plane through the midpoint of the rectangle's longitude interval. 662 // We only need to determine the distance "d" of this point from the 663 // z-axis. 664 // 665 // Let's restrict our attention to a particular z-value. In this 666 // z-plane, the Rect is a circular arc. The centroid of this arc 667 // lies on a radial line through the midpoint of the arc, and at a 668 // distance from the z-axis of 669 // 670 // r * (sin(alpha) / alpha) 671 // 672 // where r = sqrt(1-z^2) is the radius of the arc, and "alpha" is half 673 // of the arc length (i.e., the arc covers longitudes [-alpha, alpha]). 674 // 675 // To find the centroid distance from the z-axis for the entire 676 // rectangle, we just need to integrate over the z-interval. This gives 677 // 678 // d = Integrate[sqrt(1-z^2)*sin(alpha)/alpha, z1..z2] / (z2 - z1) 679 // 680 // where [z1, z2] is the range of z-values covered by the rectangle. 681 // This simplifies to 682 // 683 // d = sin(alpha)/(2*alpha*(z2-z1))*(z2*r2 - z1*r1 + theta2 - theta1) 684 // 685 // where [theta1, theta2] is the latitude interval, z1=sin(theta1), 686 // z2=sin(theta2), r1=cos(theta1), and r2=cos(theta2). 687 // 688 // Finally, we want to return not the centroid itself, but the centroid 689 // scaled by the area of the rectangle. The area of the rectangle is 690 // 691 // A = 2 * alpha * (z2 - z1) 692 // 693 // which fortunately appears in the denominator of "d". 694 695 if r.IsEmpty() { 696 return Point{} 697 } 698 699 z1 := math.Sin(r.Lat.Lo) 700 z2 := math.Sin(r.Lat.Hi) 701 r1 := math.Cos(r.Lat.Lo) 702 r2 := math.Cos(r.Lat.Hi) 703 704 alpha := 0.5 * r.Lng.Length() 705 r0 := math.Sin(alpha) * (r2*z2 - r1*z1 + r.Lat.Length()) 706 lng := r.Lng.Center() 707 z := alpha * (z2 + z1) * (z2 - z1) // scaled by the area 708 709 return Point{r3.Vector{r0 * math.Cos(lng), r0 * math.Sin(lng), z}} 710 } 711 712 // BUG: The major differences from the C++ version are: 713 // - Get*Distance, Vertex, InteriorContains(LatLng|Rect|Point) 714