← Back to branch summary

~evarlast/ubuntu/utopic/mongodb/upstart-workaround-debian-bug-718702

« back to all changes in this revision

Viewing changes to src/third_party/s2/s2latlngrect.cc

  • Committer: Package Import Robot
  • Author(s): James Page, James Page, Robie Basak
  • Date: 2013-05-29 17:44:42 UTC
  • mfrom: (44.1.7 sid)
  • Revision ID: package-import@ubuntu.com-20130529174442-z0a4qmoww4y0t458
Tags: 1:2.4.3-1ubuntu1
[ James Page ]
* Merge from Debian unstable, remaining changes:
  - Enable SSL support:
    + d/control: Add libssl-dev to BD's.
    + d/rules: Enabled --ssl option.
    + d/mongodb.conf: Add example SSL configuration options.
  - d/mongodb-server.mongodb.upstart: Add upstart configuration.
  - d/rules: Don't strip binaries during scons build for Ubuntu.
  - d/control: Add armhf to target archs.
  - d/p/SConscript.client.patch: fixup install of client libraries.
  - d/p/0010-install-libs-to-usr-lib-not-usr-lib64-Closes-588557.patch:
    Install libraries to lib not lib64.
* Dropped changes:
  - d/p/arm-support.patch: Included in Debian.
  - d/p/double-alignment.patch: Included in Debian.
  - d/rules,control: Debian also builds with avaliable system libraries
    now.
* Fix FTBFS due to gcc and boost upgrades in saucy:
  - d/p/0008-ignore-unused-local-typedefs.patch: Add -Wno-unused-typedefs
    to unbreak building with g++-4.8.
  - d/p/0009-boost-1.53.patch: Fixup signed/unsigned casting issue.

[ Robie Basak ]
* d/p/0011-Use-a-signed-char-to-store-BSONType-enumerations.patch: Fixup
  build failure on ARM due to missing signed'ness of char cast.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/third_party/s2/s2latlngrect.cc
 
1
// Copyright 2005 Google Inc. All Rights Reserved.
 
2
 
 
3
#include <algorithm>
 
4
using std::min;
 
5
using std::max;
 
6
using std::swap;
 
7
using 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
 
 
19
static const unsigned char kCurrentEncodingVersionNumber = 1;
 
20
 
 
21
S2LatLngRect S2LatLngRect::FromCenterSize(S2LatLng const& center,
 
22
                                          S2LatLng const& size) {
 
23
  return FromPoint(center).Expanded(0.5 * size);
 
24
}
 
25
 
 
26
S2LatLngRect S2LatLngRect::FromPoint(S2LatLng const& p) {
 
27
  DCHECK(p.is_valid());
 
28
  return S2LatLngRect(p, p);
 
29
}
 
30
 
 
31
S2LatLngRect 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
 
 
41
S2LatLngRect* S2LatLngRect::Clone() const {
 
42
  return new S2LatLngRect(*this);
 
43
}
 
44
 
 
45
S2LatLng 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
 
 
50
S2LatLng S2LatLngRect::GetCenter() const {
 
51
  return S2LatLng::FromRadians(lat_.GetCenter(), lng_.GetCenter());
 
52
}
 
53
 
 
54
S2LatLng S2LatLngRect::GetSize() const {
 
55
  return S2LatLng::FromRadians(lat_.GetLength(), lng_.GetLength());
 
56
}
 
57
 
 
58
double 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
 
 
66
bool 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
 
 
72
bool S2LatLngRect::InteriorContains(S2Point const& p) const {
 
73
  return InteriorContains(S2LatLng(p));
 
74
}
 
75
 
 
76
bool 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
 
 
82
bool S2LatLngRect::Contains(S2LatLngRect const& other) const {
 
83
  return lat_.Contains(other.lat_) && lng_.Contains(other.lng_);
 
84
}
 
85
 
 
86
bool S2LatLngRect::InteriorContains(S2LatLngRect const& other) const {
 
87
  return (lat_.InteriorContains(other.lat_) &&
 
88
          lng_.InteriorContains(other.lng_));
 
89
}
 
90
 
 
91
bool S2LatLngRect::Intersects(S2LatLngRect const& other) const {
 
92
  return lat_.Intersects(other.lat_) && lng_.Intersects(other.lng_);
 
93
}
 
94
 
 
95
bool S2LatLngRect::InteriorIntersects(S2LatLngRect const& other) const {
 
96
  return (lat_.InteriorIntersects(other.lat_) &&
 
97
          lng_.InteriorIntersects(other.lng_));
 
98
}
 
99
 
 
100
void S2LatLngRect::AddPoint(S2Point const& p) {
 
101
  AddPoint(S2LatLng(p));
 
102
}
 
103
 
 
104
void S2LatLngRect::AddPoint(S2LatLng const& ll) {
 
105
  DCHECK(ll.is_valid());
 
106
  lat_.AddPoint(ll.lat().radians());
 
107
  lng_.AddPoint(ll.lng().radians());
 
108
}
 
109
 
 
110
S2LatLngRect 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
 
 
118
S2LatLngRect S2LatLngRect::Union(S2LatLngRect const& other) const {
 
119
  return S2LatLngRect(lat_.Union(other.lat_),
 
120
                      lng_.Union(other.lng_));
 
121
}
 
122
 
 
123
S2LatLngRect 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
 
 
133
S2LatLngRect 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
 
 
150
S2Cap 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
 
 
188
S2LatLngRect S2LatLngRect::GetRectBound() const {
 
189
  return *this;
 
190
}
 
191
 
 
192
bool 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
 
 
205
bool 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
 
 
210
void 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
 
 
222
bool 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
 
 
238
bool 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
 
 
249
bool 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
 
 
302
bool 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
 
 
350
S1Angle 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
 
 
409
S1Angle 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
 
 
435
S1Angle S2LatLngRect::GetHausdorffDistance(S2LatLngRect const& other) const {
 
436
  return max(GetDirectedHausdorffDistance(other),
 
437
             other.GetDirectedHausdorffDistance(*this));
 
438
}
 
439
 
 
440
S1Angle 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'.
 
457
S1Angle 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'.
 
542
S2Point 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.
 
561
S1Angle 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
 
 
579
bool S2LatLngRect::Contains(S2Point const& p) const {
 
580
  return Contains(S2LatLng(p));
 
581
}
 
582
 
 
583
bool 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
 
 
589
ostream& operator<<(ostream& os, S2LatLngRect const& r) {
 
590
  return os << "[Lo" << r.lo() << ", Hi" << r.hi() << "]";
 
591
}