1// Copyright 2005 Google Inc. All Rights Reserved.
2
3#include <algorithm>
4using std::min;
5using std::max;
6using std::swap;
7using std::reverse;
8
9
10#include "s2latlngrect.h"
11
12#include "base/logging.h"
13#include "util/coding/coder.h"
14#include "s2cap.h"
15#include "s2cell.h"
16#include "s2edgeutil.h"
17#include "util/math/mathutil.h"
18
19static const unsigned char kCurrentEncodingVersionNumber = 1;
20
21S2LatLngRect S2LatLngRect::FromCenterSize(S2LatLng const& center,
22 S2LatLng const& size) {
23 return FromPoint(center).Expanded(0.5 * size);
24}
25
26S2LatLngRect S2LatLngRect::FromPoint(S2LatLng const& p) {
27 DCHECK(p.is_valid());
28 return S2LatLngRect(p, p);
29}
30
31S2LatLngRect S2LatLngRect::FromPointPair(S2LatLng const& p1,
32 S2LatLng const& p2) {
33 DCHECK(p1.is_valid()) << p1;
34 DCHECK(p2.is_valid()) << p2;
35 return S2LatLngRect(R1Interval::FromPointPair(p1.lat().radians(),
36 p2.lat().radians()),
37 S1Interval::FromPointPair(p1.lng().radians(),
38 p2.lng().radians()));
39}
40
41S2LatLngRect* S2LatLngRect::Clone() const {
42 return new S2LatLngRect(*this);
43}
44
45S2LatLng S2LatLngRect::GetVertex(int k) const {
46 // Twiddle bits to return the points in CCW order (SW, SE, NE, NW).
47 return S2LatLng::FromRadians(lat_.bound(k>>1), lng_.bound((k>>1) ^ (k&1)));
48}
49
50S2LatLng S2LatLngRect::GetCenter() const {
51 return S2LatLng::FromRadians(lat_.GetCenter(), lng_.GetCenter());
52}
53
54S2LatLng S2LatLngRect::GetSize() const {
55 return S2LatLng::FromRadians(lat_.GetLength(), lng_.GetLength());
56}
57
58double S2LatLngRect::Area() const {
59 if (is_empty()) return 0.0;
60 // This is the size difference of the two spherical caps, multiplied by
61 // the longitude ratio.
62 return lng().GetLength()* fabs(sin(lat_hi().radians()) -
63 sin(lat_lo().radians()));
64}
65
66bool S2LatLngRect::Contains(S2LatLng const& ll) const {
67 DCHECK(ll.is_valid());
68 return (lat_.Contains(ll.lat().radians()) &&
69 lng_.Contains(ll.lng().radians()));
70}
71
72bool S2LatLngRect::InteriorContains(S2Point const& p) const {
73 return InteriorContains(S2LatLng(p));
74}
75
76bool S2LatLngRect::InteriorContains(S2LatLng const& ll) const {
77 DCHECK(ll.is_valid());
78 return (lat_.InteriorContains(ll.lat().radians()) &&
79 lng_.InteriorContains(ll.lng().radians()));
80}
81
82bool S2LatLngRect::Contains(S2LatLngRect const& other) const {
83 return lat_.Contains(other.lat_) && lng_.Contains(other.lng_);
84}
85
86bool S2LatLngRect::InteriorContains(S2LatLngRect const& other) const {
87 return (lat_.InteriorContains(other.lat_) &&
88 lng_.InteriorContains(other.lng_));
89}
90
91bool S2LatLngRect::Intersects(S2LatLngRect const& other) const {
92 return lat_.Intersects(other.lat_) && lng_.Intersects(other.lng_);
93}
94
95bool S2LatLngRect::InteriorIntersects(S2LatLngRect const& other) const {
96 return (lat_.InteriorIntersects(other.lat_) &&
97 lng_.InteriorIntersects(other.lng_));
98}
99
100void S2LatLngRect::AddPoint(S2Point const& p) {
101 AddPoint(S2LatLng(p));
102}
103
104void S2LatLngRect::AddPoint(S2LatLng const& ll) {
105 DCHECK(ll.is_valid());
106 lat_.AddPoint(ll.lat().radians());
107 lng_.AddPoint(ll.lng().radians());
108}
109
110S2LatLngRect S2LatLngRect::Expanded(S2LatLng const& margin) const {
111 DCHECK_GE(margin.lat().radians(), 0);
112 DCHECK_GE(margin.lng().radians(), 0);
113 return S2LatLngRect(
114 lat_.Expanded(margin.lat().radians()).Intersection(FullLat()),
115 lng_.Expanded(margin.lng().radians()));
116}
117
118S2LatLngRect S2LatLngRect::Union(S2LatLngRect const& other) const {
119 return S2LatLngRect(lat_.Union(other.lat_),
120 lng_.Union(other.lng_));
121}
122
123S2LatLngRect S2LatLngRect::Intersection(S2LatLngRect const& other) const {
124 R1Interval lat = lat_.Intersection(other.lat_);
125 S1Interval lng = lng_.Intersection(other.lng_);
126 if (lat.is_empty() || lng.is_empty()) {
127 // The lat/lng ranges must either be both empty or both non-empty.
128 return Empty();
129 }
130 return S2LatLngRect(lat, lng);
131}
132
133S2LatLngRect S2LatLngRect::ConvolveWithCap(S1Angle const& angle) const {
134 // The most straightforward approach is to build a cap centered on each
135 // vertex and take the union of all the bounding rectangles (including the
136 // original rectangle; this is necessary for very large rectangles).
137
138 // Optimization: convert the angle to a height exactly once.
139 S2Cap cap = S2Cap::FromAxisAngle(S2Point(1, 0, 0), angle);
140
141 S2LatLngRect r = *this;
142 for (int k = 0; k < 4; ++k) {
143 S2Cap vertex_cap = S2Cap::FromAxisHeight(GetVertex(k).ToPoint(),
144 cap.height());
145 r = r.Union(vertex_cap.GetRectBound());
146 }
147 return r;
148}
149
150S2Cap S2LatLngRect::GetCapBound() const {
151 // We consider two possible bounding caps, one whose axis passes
152 // through the center of the lat-long rectangle and one whose axis
153 // is the north or south pole. We return the smaller of the two caps.
154
155 if (is_empty()) return S2Cap::Empty();
156
157 double pole_z, pole_angle;
158 if (lat_.lo() + lat_.hi() < 0) {
159 // South pole axis yields smaller cap.
160 pole_z = -1;
161 pole_angle = M_PI_2 + lat_.hi();
162 } else {
163 pole_z = 1;
164 pole_angle = M_PI_2 - lat_.lo();
165 }
166 S2Cap pole_cap = S2Cap::FromAxisAngle(S2Point(0, 0, pole_z),
167 S1Angle::Radians(pole_angle));
168
169 // For bounding rectangles that span 180 degrees or less in longitude, the
170 // maximum cap size is achieved at one of the rectangle vertices. For
171 // rectangles that are larger than 180 degrees, we punt and always return a
172 // bounding cap centered at one of the two poles.
173 double lng_span = lng_.hi() - lng_.lo();
174 if (drem(lng_span, 2 * M_PI) >= 0) {
175 if (lng_span < 2 * M_PI) {
176 S2Cap mid_cap = S2Cap::FromAxisAngle(GetCenter().ToPoint(),
177 S1Angle::Radians(0));
178 for (int k = 0; k < 4; ++k) {
179 mid_cap.AddPoint(GetVertex(k).ToPoint());
180 }
181 if (mid_cap.height() < pole_cap.height())
182 return mid_cap;
183 }
184 }
185 return pole_cap;
186}
187
188S2LatLngRect S2LatLngRect::GetRectBound() const {
189 return *this;
190}
191
192bool S2LatLngRect::Contains(S2Cell const& cell) const {
193 // A latitude-longitude rectangle contains a cell if and only if it contains
194 // the cell's bounding rectangle. This test is exact from a mathematical
195 // point of view, assuming that the bounds returned by S2Cell::GetRectBound()
196 // are tight. However, note that there can be a loss of precision when
197 // converting between representations -- for example, if an S2Cell is
198 // converted to a polygon, the polygon's bounding rectangle may not contain
199 // the cell's bounding rectangle. This has some slightly unexpected side
200 // effects; for instance, if one creates an S2Polygon from an S2Cell, the
201 // polygon will contain the cell, but the polygon's bounding box will not.
202 return Contains(cell.GetRectBound());
203}
204
205bool S2LatLngRect::MayIntersect(S2Cell const& cell) const {
206 // This test is cheap but is NOT exact (see s2latlngrect.h).
207 return Intersects(cell.GetRectBound());
208}
209
210void S2LatLngRect::Encode(Encoder* encoder) const {
211 encoder->Ensure(40); // sufficient
212
213 encoder->put8(kCurrentEncodingVersionNumber);
214 encoder->putdouble(lat_.lo());
215 encoder->putdouble(lat_.hi());
216 encoder->putdouble(lng_.lo());
217 encoder->putdouble(lng_.hi());
218
219 DCHECK_GE(encoder->avail(), 0);
220}
221
222bool S2LatLngRect::Decode(Decoder* decoder) {
223 unsigned char version = decoder->get8();
224 if (version > kCurrentEncodingVersionNumber) return false;
225
226 double lat_lo = decoder->getdouble();
227 double lat_hi = decoder->getdouble();
228 lat_ = R1Interval(lat_lo, lat_hi);
229 double lng_lo = decoder->getdouble();
230 double lng_hi = decoder->getdouble();
231 lng_ = S1Interval(lng_lo, lng_hi);
232
233 DCHECK(is_valid());
234
235 return decoder->avail() >= 0;
236}
237
238bool S2LatLngRect::IntersectsLngEdge(S2Point const& a, S2Point const& b,
239 R1Interval const& lat, double lng) {
240 // Return true if the segment AB intersects the given edge of constant
241 // longitude. The nice thing about edges of constant longitude is that
242 // they are straight lines on the sphere (geodesics).
243
244 return S2EdgeUtil::SimpleCrossing(
245 a, b, S2LatLng::FromRadians(lat.lo(), lng).ToPoint(),
246 S2LatLng::FromRadians(lat.hi(), lng).ToPoint());
247}
248
249bool S2LatLngRect::IntersectsLatEdge(S2Point const& a, S2Point const& b,
250 double lat, S1Interval const& lng) {
251 // Return true if the segment AB intersects the given edge of constant
252 // latitude. Unfortunately, lines of constant latitude are curves on
253 // the sphere. They can intersect a straight edge in 0, 1, or 2 points.
254 DCHECK(S2::IsUnitLength(a));
255 DCHECK(S2::IsUnitLength(b));
256
257 // First, compute the normal to the plane AB that points vaguely north.
258 S2Point z = S2::RobustCrossProd(a, b).Normalize();
259 if (z[2] < 0) z = -z;
260
261 // Extend this to an orthonormal frame (x,y,z) where x is the direction
262 // where the great circle through AB achieves its maximium latitude.
263 S2Point y = S2::RobustCrossProd(z, S2Point(0, 0, 1)).Normalize();
264 S2Point x = y.CrossProd(z);
265 DCHECK(S2::IsUnitLength(x));
266 DCHECK_GE(x[2], 0);
267
268 // Compute the angle "theta" from the x-axis (in the x-y plane defined
269 // above) where the great circle intersects the given line of latitude.
270 double sin_lat = sin(lat);
271 if (fabs(sin_lat) >= x[2]) {
272 return false; // The great circle does not reach the given latitude.
273 }
274 DCHECK_GT(x[2], 0);
275 double cos_theta = sin_lat / x[2];
276 double sin_theta = sqrt(1 - cos_theta * cos_theta);
277 double theta = atan2(sin_theta, cos_theta);
278
279 // The candidate intersection points are located +/- theta in the x-y
280 // plane. For an intersection to be valid, we need to check that the
281 // intersection point is contained in the interior of the edge AB and
282 // also that it is contained within the given longitude interval "lng".
283
284 // Compute the range of theta values spanned by the edge AB.
285 S1Interval ab_theta = S1Interval::FromPointPair(
286 atan2(a.DotProd(y), a.DotProd(x)),
287 atan2(b.DotProd(y), b.DotProd(x)));
288
289 if (ab_theta.Contains(theta)) {
290 // Check if the intersection point is also in the given "lng" interval.
291 S2Point isect = x * cos_theta + y * sin_theta;
292 if (lng.Contains(atan2(isect[1], isect[0]))) return true;
293 }
294 if (ab_theta.Contains(-theta)) {
295 // Check if the intersection point is also in the given "lng" interval.
296 S2Point isect = x * cos_theta - y * sin_theta;
297 if (lng.Contains(atan2(isect[1], isect[0]))) return true;
298 }
299 return false;
300}
301
302bool S2LatLngRect::Intersects(S2Cell const& cell) const {
303 // First we eliminate the cases where one region completely contains the
304 // other. Once these are disposed of, then the regions will intersect
305 // if and only if their boundaries intersect.
306
307 if (is_empty()) return false;
308 if (Contains(cell.GetCenterRaw())) return true;
309 if (cell.Contains(GetCenter().ToPoint())) return true;
310
311 // Quick rejection test (not required for correctness).
312 if (!Intersects(cell.GetRectBound())) return false;
313
314 // Precompute the cell vertices as points and latitude-longitudes. We also
315 // check whether the S2Cell contains any corner of the rectangle, or
316 // vice-versa, since the edge-crossing tests only check the edge interiors.
317
318 S2Point cell_v[4];
319 S2LatLng cell_ll[4];
320 for (int i = 0; i < 4; ++i) {
321 cell_v[i] = cell.GetVertex(i); // Must be normalized.
322 cell_ll[i] = S2LatLng(cell_v[i]);
323 if (Contains(cell_ll[i])) return true;
324 if (cell.Contains(GetVertex(i).ToPoint())) return true;
325 }
326
327 // Now check whether the boundaries intersect. Unfortunately, a
328 // latitude-longitude rectangle does not have straight edges -- two edges
329 // are curved, and at least one of them is concave.
330
331 for (int i = 0; i < 4; ++i) {
332 S1Interval edge_lng = S1Interval::FromPointPair(
333 cell_ll[i].lng().radians(), cell_ll[(i+1)&3].lng().radians());
334 if (!lng_.Intersects(edge_lng)) continue;
335
336 S2Point const& a = cell_v[i];
337 S2Point const& b = cell_v[(i+1)&3];
338 if (edge_lng.Contains(lng_.lo())) {
339 if (IntersectsLngEdge(a, b, lat_, lng_.lo())) return true;
340 }
341 if (edge_lng.Contains(lng_.hi())) {
342 if (IntersectsLngEdge(a, b, lat_, lng_.hi())) return true;
343 }
344 if (IntersectsLatEdge(a, b, lat_.lo(), lng_)) return true;
345 if (IntersectsLatEdge(a, b, lat_.hi(), lng_)) return true;
346 }
347 return false;
348}
349
350S1Angle S2LatLngRect::GetDistance(S2LatLngRect const& other) const {
351 S2LatLngRect const& a = *this;
352 S2LatLngRect const& b = other;
353 DCHECK(!a.is_empty());
354 DCHECK(!b.is_empty());
355
356 // First, handle the trivial cases where the longitude intervals overlap.
357 if (a.lng().Intersects(b.lng())) {
358 if (a.lat().Intersects(b.lat()))
359 return S1Angle::Radians(0); // Intersection between a and b.
360
361 // We found an overlap in the longitude interval, but not in the latitude
362 // interval. This means the shortest path travels along some line of
363 // longitude connecting the high-latitude of the lower rect with the
364 // low-latitude of the higher rect.
365 S1Angle lo, hi;
366 if (a.lat().lo() > b.lat().hi()) {
367 lo = b.lat_hi();
368 hi = a.lat_lo();
369 } else {
370 lo = a.lat_hi();
371 hi = b.lat_lo();
372 }
373 return hi - lo;
374 }
375
376 // The longitude intervals don't overlap. In this case, the closest points
377 // occur somewhere on the pair of longitudinal edges which are nearest in
378 // longitude-space.
379 S1Angle a_lng, b_lng;
380 S1Interval lo_hi = S1Interval::FromPointPair(a.lng().lo(), b.lng().hi());
381 S1Interval hi_lo = S1Interval::FromPointPair(a.lng().hi(), b.lng().lo());
382 if (lo_hi.GetLength() < hi_lo.GetLength()) {
383 a_lng = a.lng_lo();
384 b_lng = b.lng_hi();
385 } else {
386 a_lng = a.lng_hi();
387 b_lng = b.lng_lo();
388 }
389
390 // The shortest distance between the two longitudinal segments will include at
391 // least one segment endpoint. We could probably narrow this down further to a
392 // single point-edge distance by comparing the relative latitudes of the
393 // endpoints, but for the sake of clarity, we'll do all four point-edge
394 // distance tests.
395 S2Point a_lo = S2LatLng(a.lat_lo(), a_lng).ToPoint();
396 S2Point a_hi = S2LatLng(a.lat_hi(), a_lng).ToPoint();
397 S2Point a_lo_cross_hi =
398 S2LatLng::FromRadians(0, a_lng.radians() - M_PI_2).Normalized().ToPoint();
399 S2Point b_lo = S2LatLng(b.lat_lo(), b_lng).ToPoint();
400 S2Point b_hi = S2LatLng(b.lat_hi(), b_lng).ToPoint();
401 S2Point b_lo_cross_hi =
402 S2LatLng::FromRadians(0, b_lng.radians() - M_PI_2).Normalized().ToPoint();
403 return min(S2EdgeUtil::GetDistance(a_lo, b_lo, b_hi, b_lo_cross_hi),
404 min(S2EdgeUtil::GetDistance(a_hi, b_lo, b_hi, b_lo_cross_hi),
405 min(S2EdgeUtil::GetDistance(b_lo, a_lo, a_hi, a_lo_cross_hi),
406 S2EdgeUtil::GetDistance(b_hi, a_lo, a_hi, a_lo_cross_hi))));
407}
408
409S1Angle S2LatLngRect::GetDistance(S2LatLng const& p) const {
410 // The algorithm here is the same as in GetDistance(S2LagLngRect), only
411 // with simplified calculations.
412 S2LatLngRect const& a = *this;
413 DCHECK(!a.is_empty());
414 DCHECK(p.is_valid());
415
416 if (a.lng().Contains(p.lng().radians())) {
417 return S1Angle::Radians(max(0.0, max(p.lat().radians() - a.lat().hi(),
418 a.lat().lo() - p.lat().radians())));
419 }
420
421 S1Interval interval(a.lng().hi(), a.lng().GetComplementCenter());
422 double a_lng;
423 if (interval.Contains(p.lng().radians())) {
424 a_lng = a.lng().hi();
425 } else {
426 a_lng = a.lng().lo();
427 }
428 S2Point lo = S2LatLng::FromRadians(a.lat().lo(), a_lng).ToPoint();
429 S2Point hi = S2LatLng::FromRadians(a.lat().hi(), a_lng).ToPoint();
430 S2Point lo_cross_hi =
431 S2LatLng::FromRadians(0, a_lng - M_PI_2).Normalized().ToPoint();
432 return S2EdgeUtil::GetDistance(p.ToPoint(), lo, hi, lo_cross_hi);
433}
434
435S1Angle S2LatLngRect::GetHausdorffDistance(S2LatLngRect const& other) const {
436 return max(GetDirectedHausdorffDistance(other),
437 other.GetDirectedHausdorffDistance(*this));
438}
439
440S1Angle S2LatLngRect::GetDirectedHausdorffDistance(
441 S2LatLngRect const& other) const {
442 if (is_empty()) {
443 return S1Angle::Radians(0);
444 }
445 if (other.is_empty()) {
446 return S1Angle::Radians(M_PI); // maximum possible distance on S2
447 }
448
449 double lng_distance = lng().GetDirectedHausdorffDistance(other.lng());
450 DCHECK_GE(lng_distance, 0);
451 return GetDirectedHausdorffDistance(lng_distance, lat(), other.lat());
452}
453
454// Return the directed Hausdorff distance from one longitudinal edge spanning
455// latitude range 'a_lat' to the other longitudinal edge spanning latitude
456// range 'b_lat', with their longitudinal difference given by 'lng_diff'.
457S1Angle S2LatLngRect::GetDirectedHausdorffDistance(
458 double lng_diff, R1Interval const& a, R1Interval const& b) {
459 // By symmetry, we can assume a's longtitude is 0 and b's longtitude is
460 // lng_diff. Call b's two endpoints b_lo and b_hi. Let H be the hemisphere
461 // containing a and delimited by the longitude line of b. The Voronoi diagram
462 // of b on H has three edges (portions of great circles) all orthogonal to b
463 // and meeting at b_lo_cross_b_hi.
464 // E1: (b_lo, b_lo_cross_b_hi)
465 // E2: (b_hi, b_lo_cross_b_hi)
466 // E3: (-b_mid, b_lo_cross_b_hi), where b_mid is the midpoint of b
467 //
468 // They subdivide H into three Voronoi regions. Depending on how longitude 0
469 // (which contains edge a) intersects these regions, we distinguish two cases:
470 // Case 1: it intersects three regions. This occurs when lng_diff <= M_PI_2.
471 // Case 2: it intersects only two regions. This occurs when lng_diff > M_PI_2.
472 //
473 // In the first case, the directed Hausdorff distance to edge b can only be
474 // realized by the following points on a:
475 // A1: two endpoints of a.
476 // A2: intersection of a with the equator, if b also intersects the equator.
477 //
478 // In the second case, the directed Hausdorff distance to edge b can only be
479 // realized by the following points on a:
480 // B1: two endpoints of a.
481 // B2: intersection of a with E3
482 // B3: farthest point from b_lo to the interior of D, and farthest point from
483 // b_hi to the interior of U, if any, where D (resp. U) is the portion
484 // of edge a below (resp. above) the intersection point from B2.
485
486 DCHECK_GE(lng_diff, 0);
487 DCHECK_LE(lng_diff, M_PI);
488
489 if (lng_diff == 0) {
490 return S1Angle::Radians(a.GetDirectedHausdorffDistance(b));
491 }
492
493 // Assumed longtitude of b.
494 double b_lng = lng_diff;
495 // Two endpoints of b.
496 S2Point b_lo = S2LatLng::FromRadians(b.lo(), b_lng).ToPoint();
497 S2Point b_hi = S2LatLng::FromRadians(b.hi(), b_lng).ToPoint();
498 // Cross product of b_lo and b_hi.
499 const S2Point& b_lo_cross_b_hi =
500 S2LatLng::FromRadians(0, b_lng - M_PI_2).ToPoint();
501
502 // Handling of each case outlined at the top of the function starts here.
503 // This is initialized a few lines below.
504 S1Angle max_distance;
505
506 // Cases A1 and B1.
507 S2Point a_lo = S2LatLng::FromRadians(a.lo(), 0).ToPoint();
508 S2Point a_hi = S2LatLng::FromRadians(a.hi(), 0).ToPoint();
509 max_distance = S2EdgeUtil::GetDistance(a_lo, b_lo, b_hi, b_lo_cross_b_hi);
510 max_distance = max(
511 max_distance, S2EdgeUtil::GetDistance(a_hi, b_lo, b_hi, b_lo_cross_b_hi));
512
513 if (lng_diff <= M_PI_2) {
514 // Case A2.
515 if (a.Contains(0) && b.Contains(0)) {
516 max_distance = max(max_distance, S1Angle::Radians(lng_diff));
517 }
518 } else {
519 // Case B2.
520 const S2Point& p = GetBisectorIntersection(b, b_lng);
521 double p_lat = S2LatLng::Latitude(p).radians();
522 if (a.Contains(p_lat)) {
523 max_distance = max(max_distance, S1Angle::Radians(p.Angle(b_lo)));
524 }
525
526 // Case B3.
527 if (p_lat > a.lo()) {
528 max_distance = max(max_distance, GetInteriorMaxDistance(
529 R1Interval(a.lo(), min(p_lat, a.hi())), b_lo));
530 }
531 if (p_lat < a.hi()) {
532 max_distance = max(max_distance, GetInteriorMaxDistance(
533 R1Interval(max(p_lat, a.lo()), a.hi()), b_hi));
534 }
535 }
536
537 return max_distance;
538}
539
540// Return the intersection of longitude 0 with the bisector of an edge
541// on longitude 'lng' and spanning latitude range 'lat'.
542S2Point S2LatLngRect::GetBisectorIntersection(R1Interval const& lat,
543 double lng) {
544 lng = fabs(lng);
545 double lat_center = lat.GetCenter();
546 // A vector orthogonal to the bisector of the given longitudinal edge.
547 S2LatLng ortho_bisector;
548 if (lat_center >= 0) {
549 ortho_bisector = S2LatLng::FromRadians(lat_center - M_PI_2, lng);
550 } else {
551 ortho_bisector = S2LatLng::FromRadians(-lat_center - M_PI_2, lng - M_PI);
552 }
553 // A vector orthogonal to longitude 0.
554 static const S2Point ortho_lng = S2Point(0, -1, 0);
555 return S2::RobustCrossProd(ortho_lng, ortho_bisector.ToPoint());
556}
557
558// Return max distance from a point b to the segment spanning latitude range
559// a_lat on longitude 0, if the max occurs in the interior of a_lat. Otherwise
560// return -1.
561S1Angle S2LatLngRect::GetInteriorMaxDistance(R1Interval const& a_lat,
562 S2Point const& b) {
563 // Longitude 0 is in the y=0 plane. b.x() >= 0 implies that the maximum
564 // does not occur in the interior of a_lat.
565 if (a_lat.is_empty() || b.x() >= 0) return S1Angle::Radians(-1);
566
567 // Project b to the y=0 plane. The antipodal of the normalized projection is
568 // the point at which the maxium distance from b occurs, if it is contained
569 // in a_lat.
570 S2Point intersection_point = S2Point(-b.x(), 0, -b.z()).Normalize();
571 if (a_lat.InteriorContains(
572 S2LatLng::Latitude(intersection_point).radians())) {
573 return S1Angle::Radians(b.Angle(intersection_point));
574 } else {
575 return S1Angle::Radians(-1);
576 }
577}
578
579bool S2LatLngRect::Contains(S2Point const& p) const {
580 return Contains(S2LatLng(p));
581}
582
583bool S2LatLngRect::ApproxEquals(S2LatLngRect const& other,
584 double max_error) const {
585 return (lat_.ApproxEquals(other.lat_, max_error) &&
586 lng_.ApproxEquals(other.lng_, max_error));
587}
588
589ostream& operator<<(ostream& os, S2LatLngRect const& r) {
590 return os << "[Lo" << r.lo() << ", Hi" << r.hi() << "]";
591}
592