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
|
// Copyright 2005 Google Inc. All Rights Reserved.
#include "s2edgeutil.h"
#include <algorithm>
using std::min;
using std::max;
using std::swap;
using std::reverse;
#include <string>
using std::string;
#include <vector>
using std::vector;
#include "base/commandlineflags.h"
#include "base/logging.h"
#include "base/scoped_ptr.h"
#include "testing/base/public/benchmark.h"
#include "testing/base/public/gunit.h"
#include "s2cap.h"
#include "s2polyline.h"
#include "s2testing.h"
DECLARE_bool(s2debug);
namespace {
const int kDegen = -2;
void CompareResult(int actual, int expected) {
// HACK ALERT: RobustCrossing() is allowed to return 0 or -1 if either edge
// is degenerate. We use the value kDegen to represent this possibility.
if (expected == kDegen) {
EXPECT_LE(actual, 0);
} else {
EXPECT_EQ(expected, actual);
}
}
void TestCrossing(S2Point a, S2Point b, S2Point c, S2Point d,
int robust, bool edge_or_vertex, bool simple) {
a = a.Normalize();
b = b.Normalize();
c = c.Normalize();
d = d.Normalize();
CompareResult(S2EdgeUtil::RobustCrossing(a, b, c, d), robust);
if (simple) {
EXPECT_EQ(robust > 0, S2EdgeUtil::SimpleCrossing(a, b, c, d));
}
S2EdgeUtil::EdgeCrosser crosser(&a, &b, &c);
CompareResult(crosser.RobustCrossing(&d), robust);
CompareResult(crosser.RobustCrossing(&c), robust);
EXPECT_EQ(edge_or_vertex, S2EdgeUtil::EdgeOrVertexCrossing(a, b, c, d));
EXPECT_EQ(edge_or_vertex, crosser.EdgeOrVertexCrossing(&d));
EXPECT_EQ(edge_or_vertex, crosser.EdgeOrVertexCrossing(&c));
}
void TestCrossings(S2Point a, S2Point b, S2Point c, S2Point d,
int robust, bool edge_or_vertex, bool simple) {
TestCrossing(a, b, c, d, robust, edge_or_vertex, simple);
TestCrossing(b, a, c, d, robust, edge_or_vertex, simple);
TestCrossing(a, b, d, c, robust, edge_or_vertex, simple);
TestCrossing(b, a, d, c, robust, edge_or_vertex, simple);
TestCrossing(a, a, c, d, kDegen, 0, false);
TestCrossing(a, b, c, c, kDegen, 0, false);
TestCrossing(a, b, a, b, 0, 1, false);
TestCrossing(c, d, a, b, robust, edge_or_vertex ^ (robust == 0), simple);
}
TEST(S2EdgeUtil, Crossings) {
// The real tests of edge crossings are in s2{loop,polygon}_unittest,
// but we do a few simple tests here.
// Two regular edges that cross.
TestCrossings(S2Point(1, 2, 1), S2Point(1, -3, 0.5),
S2Point(1, -0.5, -3), S2Point(0.1, 0.5, 3), 1, true, true);
// Two regular edges that cross antipodal points.
TestCrossings(S2Point(1, 2, 1), S2Point(1, -3, 0.5),
S2Point(-1, 0.5, 3), S2Point(-0.1, -0.5, -3), -1, false, true);
// Two edges on the same great circle.
TestCrossings(S2Point(0, 0, -1), S2Point(0, 1, 0),
S2Point(0, 1, 1), S2Point(0, 0, 1), -1, false, true);
// Two edges that cross where one vertex is S2::Origin().
TestCrossings(S2Point(1, 0, 0), S2::Origin(),
S2Point(1, -0.1, 1), S2Point(1, 1, -0.1), 1, true, true);
// Two edges that cross antipodal points where one vertex is S2::Origin().
TestCrossings(S2Point(1, 0, 0), S2Point(0, 1, 0),
S2Point(0, 0, -1), S2Point(-1, -1, 1), -1, false, true);
// Two edges that share an endpoint. The Ortho() direction is (-4,0,2),
// and edge CD is further CCW around (2,3,4) than AB.
TestCrossings(S2Point(2, 3, 4), S2Point(-1, 2, 5),
S2Point(7, -2, 3), S2Point(2, 3, 4), 0, false, true);
// Two edges that barely cross each other near the middle of one edge. The
// edge AB is approximately in the x=y plane, while CD is approximately
// perpendicular to it and ends exactly at the x=y plane.
TestCrossings(S2Point(1, 1, 1), S2Point(1, nextafter(1, 0), -1),
S2Point(11, -12, -1), S2Point(10, 10, 1), 1, true, false);
// In this version, the edges are separated by a distance of about 1e-15.
TestCrossings(S2Point(1, 1, 1), S2Point(1, nextafter(1, 2), -1),
S2Point(1, -1, 0), S2Point(1, 1, 0), -1, false, false);
// Two edges that barely cross each other near the end of both edges. This
// example cannot be handled using regular double-precision arithmetic due
// to floating-point underflow.
TestCrossings(S2Point(0, 0, 1), S2Point(2, -1e-323, 1),
S2Point(1, -1, 1), S2Point(1e-323, 0, 1), 1, true, false);
// In this version, the edges are separated by a distance of about 1e-640.
TestCrossings(S2Point(0, 0, 1), S2Point(2, 1e-323, 1),
S2Point(1, -1, 1), S2Point(1e-323, 0, 1), -1, false, false);
// Two edges that barely cross each other near the middle of one edge.
// Computing the exact determinant of some of the triangles in this test
// requires more than 2000 bits of precision.
TestCrossings(S2Point(1, -1e-323, -1e-323), S2Point(1e-323, 1, 1e-323),
S2Point(1, -1, 1e-323), S2Point(1, 1, 0),
1, true, false);
// In this version, the edges are separated by a distance of about 1e-640.
TestCrossings(S2Point(1, 1e-323, -1e-323), S2Point(-1e-323, 1, 1e-323),
S2Point(1, -1, 1e-323), S2Point(1, 1, 0),
-1, false, false);
}
typedef bool CrossingFunction(S2Point const&, S2Point const&,
S2Point const&, S2Point const&);
void BenchmarkCrossing(CrossingFunction crossing, int iters) {
StopBenchmarkTiming();
// We want to avoid cache effects, so kNumPoints should be small enough so
// that the points can be in L1 cache. sizeof(S2Point) == 24, so 400 will
// only take ~9KiB of 64KiB L1 cache.
static const int kNumPoints = 400;
vector<S2Point> p(kNumPoints);
for (int i = 0; i < kNumPoints; ++i)
p[i] = S2Testing::RandomPoint();
StartBenchmarkTiming();
int num_crossings = 0;
for (int i = 0; i < iters; ++i) {
const S2Point& a = p[(i + 0) % kNumPoints];
const S2Point& b = p[(i + 1) % kNumPoints];
const S2Point& c = p[(i + 2) % kNumPoints];
const S2Point& d = p[(i + 3) % kNumPoints];
if (crossing(a, b, c, d))
++num_crossings;
}
VLOG(5) << "Fraction crossing = " << 1.0 * num_crossings / iters;
VLOG(5) << "iters: " << iters;
}
void BM_EdgeOrVertexCrossing(int iters) {
BenchmarkCrossing(&S2EdgeUtil::EdgeOrVertexCrossing, iters);
}
BENCHMARK(BM_EdgeOrVertexCrossing);
bool RobustCrossingBool(S2Point const& a, S2Point const& b,
S2Point const& c, S2Point const& d) {
return S2EdgeUtil::RobustCrossing(a, b, c, d) > 0;
}
void BM_RobustCrossing(int iters) {
BenchmarkCrossing(&RobustCrossingBool, iters);
}
BENCHMARK(BM_RobustCrossing);
S2LatLngRect GetEdgeBound(double x1, double y1, double z1,
double x2, double y2, double z2) {
S2EdgeUtil::RectBounder bounder;
S2Point p1 = S2Point(x1, y1, z1).Normalize();
S2Point p2 = S2Point(x2, y2, z2).Normalize();
bounder.AddPoint(&p1);
bounder.AddPoint(&p2);
return bounder.GetBound();
}
TEST(S2EdgeUtil, RectBounder) {
// Check cases where min/max latitude is not at a vertex.
EXPECT_DOUBLE_EQ(M_PI_4,
GetEdgeBound(1,1,1, 1,-1,1).lat().hi()); // Max, CW
EXPECT_DOUBLE_EQ(M_PI_4,
GetEdgeBound(1,-1,1, 1,1,1).lat().hi()); // Max, CCW
EXPECT_DOUBLE_EQ(-M_PI_4,
GetEdgeBound(1,-1,-1, -1,-1,-1).lat().lo()); // Min, CW
EXPECT_DOUBLE_EQ(-M_PI_4,
GetEdgeBound(-1,1,-1, -1,-1,-1).lat().lo()); // Min, CCW
// Check cases where the edge passes through one of the poles.
EXPECT_DOUBLE_EQ(M_PI_2, GetEdgeBound(.3,.4,1, -.3,-.4,1).lat().hi());
EXPECT_DOUBLE_EQ(-M_PI_2, GetEdgeBound(.3,.4,-1, -.3,-.4,-1).lat().lo());
// Check cases where the min/max latitude is attained at a vertex.
static double const kCubeLat = asin(sqrt(1./3)); // 35.26 degrees
EXPECT_TRUE(GetEdgeBound(1,1,1, 1,-1,-1).lat().
ApproxEquals(R1Interval(-kCubeLat, kCubeLat)));
EXPECT_TRUE(GetEdgeBound(1,-1,1, 1,1,-1).lat().
ApproxEquals(R1Interval(-kCubeLat, kCubeLat)));
}
TEST(S2EdgeUtil, LongitudePruner) {
S2EdgeUtil::LongitudePruner pruner1(S1Interval(0.75 * M_PI, -0.75 * M_PI),
S2Point(0, 1, 2));
EXPECT_FALSE(pruner1.Intersects(S2Point(1, 1, 3)));
EXPECT_TRUE(pruner1.Intersects(S2Point(-1 - 1e-15, -1, 0)));
EXPECT_TRUE(pruner1.Intersects(S2Point(-1, 0, 0)));
EXPECT_TRUE(pruner1.Intersects(S2Point(-1, 0, 0)));
EXPECT_TRUE(pruner1.Intersects(S2Point(1, -1, 8)));
EXPECT_FALSE(pruner1.Intersects(S2Point(1, 0, -2)));
EXPECT_TRUE(pruner1.Intersects(S2Point(-1, -1e-15, 0)));
S2EdgeUtil::LongitudePruner pruner2(S1Interval(0.25 * M_PI, 0.25 * M_PI),
S2Point(1, 0, 0));
EXPECT_FALSE(pruner2.Intersects(S2Point(2, 1, 2)));
EXPECT_TRUE(pruner2.Intersects(S2Point(1, 2, 3)));
EXPECT_FALSE(pruner2.Intersects(S2Point(0, 1, 4)));
EXPECT_FALSE(pruner2.Intersects(S2Point(-1e-15, -1, -1)));
}
void TestWedge(S2Point a0, S2Point ab1, S2Point a2, S2Point b0, S2Point b2,
bool contains, bool intersects,
S2EdgeUtil::WedgeRelation wedge_relation) {
a0 = a0.Normalize();
ab1 = ab1.Normalize();
a2 = a2.Normalize();
b0 = b0.Normalize();
b2 = b2.Normalize();
EXPECT_EQ(contains, S2EdgeUtil::WedgeContains(a0, ab1, a2, b0, b2));
EXPECT_EQ(intersects, S2EdgeUtil::WedgeIntersects(a0, ab1, a2, b0, b2));
EXPECT_EQ(wedge_relation, S2EdgeUtil::GetWedgeRelation(a0, ab1, a2, b0, b2));
}
TEST(S2EdgeUtil, Wedges) {
// For simplicity, all of these tests use an origin of (0, 0, 1).
// This shouldn't matter as long as the lower-level primitives are
// implemented correctly.
// Intersection in one wedge.
TestWedge(S2Point(-1, 0, 10), S2Point(0, 0, 1), S2Point(1, 2, 10),
S2Point(0, 1, 10), S2Point(1, -2, 10),
false, true, S2EdgeUtil::WEDGE_PROPERLY_OVERLAPS);
// Intersection in two wedges.
TestWedge(S2Point(-1, -1, 10), S2Point(0, 0, 1), S2Point(1, -1, 10),
S2Point(1, 0, 10), S2Point(-1, 1, 10),
false, true, S2EdgeUtil::WEDGE_PROPERLY_OVERLAPS);
// Normal containment.
TestWedge(S2Point(-1, -1, 10), S2Point(0, 0, 1), S2Point(1, -1, 10),
S2Point(-1, 0, 10), S2Point(1, 0, 10),
true, true, S2EdgeUtil::WEDGE_PROPERLY_CONTAINS);
// Containment with equality on one side.
TestWedge(S2Point(2, 1, 10), S2Point(0, 0, 1), S2Point(-1, -1, 10),
S2Point(2, 1, 10), S2Point(1, -5, 10),
true, true, S2EdgeUtil::WEDGE_PROPERLY_CONTAINS);
// Containment with equality on the other side.
TestWedge(S2Point(2, 1, 10), S2Point(0, 0, 1), S2Point(-1, -1, 10),
S2Point(1, -2, 10), S2Point(-1, -1, 10),
true, true, S2EdgeUtil::WEDGE_PROPERLY_CONTAINS);
// Containment with equality on both sides.
TestWedge(S2Point(-2, 3, 10), S2Point(0, 0, 1), S2Point(4, -5, 10),
S2Point(-2, 3, 10), S2Point(4, -5, 10),
true, true, S2EdgeUtil::WEDGE_EQUALS);
// Disjoint with equality on one side.
TestWedge(S2Point(-2, 3, 10), S2Point(0, 0, 1), S2Point(4, -5, 10),
S2Point(4, -5, 10), S2Point(-2, -3, 10),
false, false, S2EdgeUtil::WEDGE_IS_DISJOINT);
// Disjoint with equality on the other side.
TestWedge(S2Point(-2, 3, 10), S2Point(0, 0, 1), S2Point(0, 5, 10),
S2Point(4, -5, 10), S2Point(-2, 3, 10),
false, false, S2EdgeUtil::WEDGE_IS_DISJOINT);
// Disjoint with equality on both sides.
TestWedge(S2Point(-2, 3, 10), S2Point(0, 0, 1), S2Point(4, -5, 10),
S2Point(4, -5, 10), S2Point(-2, 3, 10),
false, false, S2EdgeUtil::WEDGE_IS_DISJOINT);
// B contains A with equality on one side.
TestWedge(S2Point(2, 1, 10), S2Point(0, 0, 1), S2Point(1, -5, 10),
S2Point(2, 1, 10), S2Point(-1, -1, 10),
false, true, S2EdgeUtil::WEDGE_IS_PROPERLY_CONTAINED);
// B contains A with equality on the other side.
TestWedge(S2Point(2, 1, 10), S2Point(0, 0, 1), S2Point(1, -5, 10),
S2Point(-2, 1, 10), S2Point(1, -5, 10),
false, true, S2EdgeUtil::WEDGE_IS_PROPERLY_CONTAINED);
}
// Given a point X and an edge AB, check that the distance from X to AB is
// "distance_radians" and the closest point on AB is "expected_closest".
void CheckDistance(S2Point x, S2Point a, S2Point b,
double distance_radians, S2Point expected_closest) {
x = x.Normalize();
a = a.Normalize();
b = b.Normalize();
expected_closest = expected_closest.Normalize();
EXPECT_DOUBLE_EQ(distance_radians,
S2EdgeUtil::GetDistance(x, a, b).radians());
S2Point closest = S2EdgeUtil::GetClosestPoint(x, a, b);
if (expected_closest == S2Point(0, 0, 0)) {
// This special value says that the result should be A or B.
EXPECT_TRUE(closest == a || closest == b);
} else {
EXPECT_TRUE(S2::ApproxEquals(closest, expected_closest));
}
}
TEST(S2EdgeUtil, Distance) {
CheckDistance(S2Point(1, 0, 0), S2Point(1, 0, 0), S2Point(0, 1, 0),
0, S2Point(1, 0, 0));
CheckDistance(S2Point(0, 1, 0), S2Point(1, 0, 0), S2Point(0, 1, 0),
0, S2Point(0, 1, 0));
CheckDistance(S2Point(1, 3, 0), S2Point(1, 0, 0), S2Point(0, 1, 0),
0, S2Point(1, 3, 0));
CheckDistance(S2Point(0, 0, 1), S2Point(1, 0, 0), S2Point(0, 1, 0),
M_PI_2, S2Point(1, 0, 0));
CheckDistance(S2Point(0, 0, -1), S2Point(1, 0, 0), S2Point(0, 1, 0),
M_PI_2, S2Point(1, 0, 0));
CheckDistance(S2Point(-1, -1, 0), S2Point(1, 0, 0), S2Point(0, 1, 0),
0.75 * M_PI, S2Point(0, 0, 0));
CheckDistance(S2Point(0, 1, 0), S2Point(1, 0, 0), S2Point(1, 1, 0),
M_PI_4, S2Point(1, 1, 0));
CheckDistance(S2Point(0, -1, 0), S2Point(1, 0, 0), S2Point(1, 1, 0),
M_PI_2, S2Point(1, 0, 0));
CheckDistance(S2Point(0, -1, 0), S2Point(1, 0, 0), S2Point(-1, 1, 0),
M_PI_2, S2Point(1, 0, 0));
CheckDistance(S2Point(-1, -1, 0), S2Point(1, 0, 0), S2Point(-1, 1, 0),
M_PI_2, S2Point(-1, 1, 0));
CheckDistance(S2Point(1, 1, 1), S2Point(1, 0, 0), S2Point(0, 1, 0),
asin(sqrt(1./3)), S2Point(1, 1, 0));
CheckDistance(S2Point(1, 1, -1), S2Point(1, 0, 0), S2Point(0, 1, 0),
asin(sqrt(1./3)), S2Point(1, 1, 0));
CheckDistance(S2Point(-1, 0, 0), S2Point(1, 1, 0), S2Point(1, 1, 0),
0.75 * M_PI, S2Point(1, 1, 0));
CheckDistance(S2Point(0, 0, -1), S2Point(1, 1, 0), S2Point(1, 1, 0),
M_PI_2, S2Point(1, 1, 0));
CheckDistance(S2Point(-1, 0, 0), S2Point(1, 0, 0), S2Point(1, 0, 0),
M_PI, S2Point(1, 0, 0));
}
void CheckInterpolate(double t, S2Point a, S2Point b, S2Point expected) {
a = a.Normalize();
b = b.Normalize();
expected = expected.Normalize();
S2Point actual = S2EdgeUtil::Interpolate(t, a, b);
// We allow a bit more than the usual 1e-15 error tolerance because
// Interpolate() uses trig functions.
EXPECT_TRUE(S2::ApproxEquals(expected, actual, 3e-15))
<< "Expected: " << expected << ", actual: " << actual;
}
TEST(S2EdgeUtil, Interpolate) {
// A zero-length edge.
CheckInterpolate(0, S2Point(1, 0, 0), S2Point(1, 0, 0), S2Point(1, 0, 0));
CheckInterpolate(1, S2Point(1, 0, 0), S2Point(1, 0, 0), S2Point(1, 0, 0));
// Start, end, and middle of a medium-length edge.
CheckInterpolate(0, S2Point(1, 0, 0), S2Point(0, 1, 0), S2Point(1, 0, 0));
CheckInterpolate(1, S2Point(1, 0, 0), S2Point(0, 1, 0), S2Point(0, 1, 0));
CheckInterpolate(0.5, S2Point(1, 0, 0), S2Point(0, 1, 0), S2Point(1, 1, 0));
// Test that interpolation is done using distances on the sphere rather than
// linear distances.
CheckInterpolate(1./3, S2Point(1, 0, 0), S2Point(0, 1, 0),
S2Point(sqrt(3), 1, 0));
CheckInterpolate(2./3, S2Point(1, 0, 0), S2Point(0, 1, 0),
S2Point(1, sqrt(3), 0));
// Test that interpolation is accurate on a long edge (but not so long that
// the definition of the edge itself becomes too unstable).
{
double const kLng = M_PI - 1e-2;
S2Point a = S2LatLng::FromRadians(0, 0).ToPoint();
S2Point b = S2LatLng::FromRadians(0, kLng).ToPoint();
for (double f = 0.4; f > 1e-15; f *= 0.1) {
CheckInterpolate(f, a, b,
S2LatLng::FromRadians(0, f * kLng).ToPoint());
CheckInterpolate(1 - f, a, b,
S2LatLng::FromRadians(0, (1 - f) * kLng).ToPoint());
}
}
}
TEST(S2EdgeUtil, InterpolateCanExtrapolate) {
const S2Point i(1, 0, 0);
const S2Point j(0, 1, 0);
// Initial vectors at 90 degrees.
CheckInterpolate(0, i, j, S2Point(1, 0, 0));
CheckInterpolate(1, i, j ,S2Point(0, 1, 0));
CheckInterpolate(1.5, i, j ,S2Point(-1, 1, 0));
CheckInterpolate(2, i, j, S2Point(-1, 0, 0));
CheckInterpolate(3, i, j, S2Point(0, -1, 0));
CheckInterpolate(4, i, j, S2Point(1, 0, 0));
// Negative values of t.
CheckInterpolate(-1, i, j, S2Point(0, -1, 0));
CheckInterpolate(-2, i, j, S2Point(-1, 0, 0));
CheckInterpolate(-3, i, j ,S2Point(0, 1, 0));
CheckInterpolate(-4, i, j, S2Point(1, 0, 0));
// Initial vectors at 45 degrees.
CheckInterpolate(2, i, S2Point(1, 1, 0), S2Point(0, 1, 0));
CheckInterpolate(3, i, S2Point(1, 1, 0), S2Point(-1, 1, 0));
CheckInterpolate(4, i, S2Point(1, 1, 0), S2Point(-1, 0, 0));
// Initial vectors at 135 degrees.
CheckInterpolate(2, i, S2Point(-1, 1, 0), S2Point(0, -1, 0));
// Take a small fraction along the curve.
S2Point p(S2EdgeUtil::Interpolate(0.001, i, j));
// We should get back where we started.
CheckInterpolate(1000, i, p, j);
}
TEST(S2EdgeUtil, RepeatedInterpolation) {
// Check that points do not drift away from unit length when repeated
// interpolations are done.
for (int i = 0; i < 100; ++i) {
S2Point a = S2Testing::RandomPoint();
S2Point b = S2Testing::RandomPoint();
for (int j = 0; j < 1000; ++j) {
a = S2EdgeUtil::Interpolate(0.01, a, b);
}
EXPECT_TRUE(S2::IsUnitLength(a));
}
}
TEST(S2EdgeUtil, IntersectionTolerance) {
// We repeatedly construct two edges that cross near a random point "p",
// and measure the distance from the actual intersection point "x" to the
// the expected intersection point "p" and also to the edges that cross
// near "p".
//
// Note that GetIntersection() does not guarantee that "x" and "p" will be
// close together (since the intersection point is numerically unstable
// when the edges cross at a very small angle), but it does guarantee that
// "x" will be close to both of the edges that cross.
S1Angle max_point_dist, max_edge_dist;
S2Testing::Random* rnd = &S2Testing::rnd;
for (int i = 0; i < 1000; ++i) {
// We construct two edges AB and CD that intersect near "p". The angle
// between AB and CD (expressed as a slope) is chosen randomly between
// 1e-15 and 1.0 such that its logarithm is uniformly distributed. This
// implies that small values are much more likely to be chosen.
//
// Once the slope is chosen, the four points ABCD must be offset from P
// by at least (1e-15 / slope) so that the points are guaranteed to have
// the correct circular ordering around P. This is the distance from P
// at which the two edges are separated by about 1e-15, which is
// approximately the minimum distance at which we can expect computed
// points on the two lines to be distinct and have the correct ordering.
//
// The actual offset distance from P is chosen randomly in the range
// [1e-15 / slope, 1.0], again uniformly distributing the logarithm.
// This ensures that we test both long and very short segments that
// intersect at both large and very small angles.
S2Point p, d1, d2;
S2Testing::GetRandomFrame(&p, &d1, &d2);
double slope = pow(1e-15, rnd->RandDouble());
d2 = d1 + slope * d2;
S2Point a = (p + pow(1e-15 / slope, rnd->RandDouble()) * d1).Normalize();
S2Point b = (p - pow(1e-15 / slope, rnd->RandDouble()) * d1).Normalize();
S2Point c = (p + pow(1e-15 / slope, rnd->RandDouble()) * d2).Normalize();
S2Point d = (p - pow(1e-15 / slope, rnd->RandDouble()) * d2).Normalize();
S2Point x = S2EdgeUtil::GetIntersection(a, b, c, d);
S1Angle dist_ab = S2EdgeUtil::GetDistance(x, a, b);
S1Angle dist_cd = S2EdgeUtil::GetDistance(x, c, d);
EXPECT_LE(dist_ab, S2EdgeUtil::kIntersectionTolerance);
EXPECT_LE(dist_cd, S2EdgeUtil::kIntersectionTolerance);
max_edge_dist = max(max_edge_dist, max(dist_ab, dist_cd));
max_point_dist = max(max_point_dist, S1Angle(p, x));
}
LOG(INFO) << "Max distance to either edge being intersected: "
<< max_edge_dist.radians();
LOG(INFO) << "Maximum distance to expected intersection point: "
<< max_point_dist.radians();
}
bool IsEdgeBNearEdgeA(string const& a_str, const string& b_str,
double max_error_degrees) {
scoped_ptr<S2Polyline> a(S2Testing::MakePolyline(a_str));
EXPECT_EQ(2, a->num_vertices());
scoped_ptr<S2Polyline> b(S2Testing::MakePolyline(b_str));
EXPECT_EQ(2, b->num_vertices());
return S2EdgeUtil::IsEdgeBNearEdgeA(a->vertex(0), a->vertex(1),
b->vertex(0), b->vertex(1),
S1Angle::Degrees(max_error_degrees));
}
TEST(S2EdgeUtil, EdgeBNearEdgeA) {
// Edge is near itself.
EXPECT_TRUE(IsEdgeBNearEdgeA("5:5, 10:-5", "5:5, 10:-5", 1e-6));
// Edge is near its reverse
EXPECT_TRUE(IsEdgeBNearEdgeA("5:5, 10:-5", "10:-5, 5:5", 1e-6));
// Short edge is near long edge.
EXPECT_TRUE(IsEdgeBNearEdgeA("10:0, -10:0", "2:1, -2:1", 1.0));
// Long edges cannot be near shorter edges.
EXPECT_FALSE(IsEdgeBNearEdgeA("2:1, -2:1", "10:0, -10:0", 1.0));
// Orthogonal crossing edges are not near each other...
EXPECT_FALSE(IsEdgeBNearEdgeA("10:0, -10:0", "0:1.5, 0:-1.5", 1.0));
// ... unless all points on B are within tolerance of A.
EXPECT_TRUE(IsEdgeBNearEdgeA("10:0, -10:0", "0:1.5, 0:-1.5", 2.0));
// Very long edges whose endpoints are close may have interior points that are
// far apart. An implementation that only considers the vertices of polylines
// will incorrectly consider such edges as "close" when they are not.
// Consider, for example, two consecutive lines of longitude. As they
// approach the poles, they become arbitrarily close together, but along the
// equator they bow apart.
EXPECT_FALSE(IsEdgeBNearEdgeA("89:1, -89:1", "89:2, -89:2", 0.5));
EXPECT_TRUE(IsEdgeBNearEdgeA("89:1, -89:1", "89:2, -89:2", 1.5));
// The two arcs here are nearly as long as S2 edges can be (just shy of 180
// degrees), and their endpoints are less than 1 degree apart. Their
// midpoints, however, are at opposite ends of the sphere along its equator.
EXPECT_FALSE(IsEdgeBNearEdgeA(
"0:-179.75, 0:-0.25", "0:179.75, 0:0.25", 1.0));
// At the equator, the second arc here is 9.75 degrees from the first, and
// closer at all other points. However, the southern point of the second arc
// (-1, 9.75) is too far from the first arc for the short-circuiting logic in
// IsEdgeBNearEdgeA to apply.
EXPECT_TRUE(IsEdgeBNearEdgeA("40:0, -5:0", "39:0.975, -1:0.975", 1.0));
// Same as above, but B's orientation is reversed, causing the angle between
// the normal vectors of circ(B) and circ(A) to be (180-9.75) = 170.5 degrees,
// not 9.75 degrees. The greatest separation between the planes is still 9.75
// degrees.
EXPECT_TRUE(IsEdgeBNearEdgeA("10:0, -10:0", "-.4:0.975, 0.4:0.975", 1.0));
// A and B are on the same great circle, A and B partially overlap, but the
// only part of B that does not overlap A is shorter than tolerance.
EXPECT_TRUE(IsEdgeBNearEdgeA("0:0, 1:0", "0.9:0, 1.1:0", 0.25));
// A and B are on the same great circle, all points on B are close to A at its
// second endpoint, (1,0).
EXPECT_TRUE(IsEdgeBNearEdgeA("0:0, 1:0", "1.1:0, 1.2:0", 0.25));
// Same as above, but B's orientation is reversed. This case is special
// because the projection of the normal defining A onto the plane containing B
// is the null vector, and must be handled by a special case.
EXPECT_TRUE(IsEdgeBNearEdgeA("0:0, 1:0", "1.2:0, 1.1:0", 0.25));
}
TEST(S2EdgeUtil, CollinearEdgesThatDontTouch) {
const int kIters = 1000;
for (int iter = 0; iter < kIters; ++iter) {
S2Point a = S2Testing::RandomPoint();
S2Point d = S2Testing::RandomPoint();
S2Point b = S2EdgeUtil::Interpolate(0.05, a, d);
S2Point c = S2EdgeUtil::Interpolate(0.95, a, d);
EXPECT_GT(0, S2EdgeUtil::RobustCrossing(a, b, c, d));
EXPECT_GT(0, S2EdgeUtil::RobustCrossing(a, b, c, d));
S2EdgeUtil::EdgeCrosser crosser(&a, &b, &c);
EXPECT_GT(0, crosser.RobustCrossing(&d));
EXPECT_GT(0, crosser.RobustCrossing(&c));
}
}
TEST(S2EdgeUtil, CoincidentZeroLengthEdgesThatDontTouch) {
// Since normalization is not perfect, and vertices are not always perfectly
// normalized anyway for various reasons, it is important that the edge
// primitives can handle vertices that exactly coincident when projected
// onto the unit sphere but that are not identical.
//
// This test checks pairs of edges AB and CD where A,B,C,D are exactly
// coincident on the sphere and the norms of A,B,C,D are monotonically
// increasing. Such edge pairs should never intersect. (This is not
// obvious, since it depends on the particular symbolic perturbations used
// by S2::RobustCCW(). It would be better to replace this with a test that
// says that the CCW results must be consistent with each other.)
const int kIters = 1000;
for (int iter = 0; iter < kIters; ++iter) {
// Construct a point P where every component is zero or a power of 2.
S2Point p;
for (int i = 0; i < 2; ++i) {
int binary_exp = S2Testing::rnd.Skewed(11);
p[i] = (binary_exp > 1022) ? 0 : pow(2, -binary_exp);
}
// If all components were zero, try again. Note that normalization may
// convert a non-zero point into a zero one due to underflow (!)
p = p.Normalize();
if (p == S2Point(0, 0, 0)) { --iter; continue; }
// Now every non-zero component should have exactly the same mantissa.
// This implies that if we scale the point by an arbitrary factor, every
// non-zero component will still have the same mantissa. Scale the points
// so that they are all distinct and yet still satisfy S2::IsNormalized().
S2Point a = (1-3e-16) * p;
S2Point b = (1-1e-16) * p;
S2Point c = p;
S2Point d = (1+2e-16) * p;
// Verify that the expected edges do not cross.
EXPECT_GT(0, S2EdgeUtil::RobustCrossing(a, b, c, d));
S2EdgeUtil::EdgeCrosser crosser(&a, &b, &c);
EXPECT_GT(0, crosser.RobustCrossing(&d));
EXPECT_GT(0, crosser.RobustCrossing(&c));
}
}
static void BM_RobustCrosserEdgesCross(int iters) {
// Copied from BenchmarkCrossing() above.
StopBenchmarkTiming();
static const int kNumPoints = 400;
vector<S2Point> p(kNumPoints);
for (int i = 0; i < kNumPoints; ++i) p[i] = S2Testing::RandomPoint();
// 1/4th of the points will cross ab.
const S2Point& a = S2Testing::RandomPoint();
const S2Point& b = (-a + S2Point(0.1, 0.1, 0.1)).Normalize();
StartBenchmarkTiming();
S2EdgeUtil::EdgeCrosser crosser(&a, &b, &p[0]);
int num_crossings = 0;
for (int i = 0; i < iters; ++i) {
const S2Point& d = p[i % kNumPoints];
if (crosser.RobustCrossing(&d) > 0)
++num_crossings;
}
VLOG(5) << "Fraction crossing = " << 1.0 * num_crossings / iters;
VLOG(5) << "num_crossings = " << num_crossings;
}
BENCHMARK(BM_RobustCrosserEdgesCross);
} // namespace
|