/
rect.go
713 lines (625 loc) · 24.3 KB
/
rect.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
// 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(); int8(version) != 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 maximum 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)