1 // Copyright 2005 Google Inc. All Rights Reserved.
2 //
3 // To run the benchmarks, use:
4
5 #include "s2polygon.h"
6
7 #include <algorithm>
8 using std::min;
9 using std::max;
10 using std::swap;
11 using std::reverse;
12
13 #include <cstdio>
14 #include <string>
15 using std::string;
16
17 #include <vector>
18 using std::vector;
19
20
21 #include "base/commandlineflags.h"
22 #include "base/logging.h"
23 #include "base/macros.h"
24 #include "base/scoped_ptr.h"
25 #include "strings/stringprintf.h"
26 #include "testing/base/public/benchmark.h"
27 #include "testing/base/public/gunit.h"
28 #include "util/coding/coder.h"
29 #include "s2.h"
30 #include "s2cap.h"
31 #include "s2cellunion.h"
32 #include "s2latlng.h"
33 #include "s2loop.h"
34 #include "s2polygonbuilder.h"
35 #include "s2polyline.h"
36 #include "s2regioncoverer.h"
37 #include "s2testing.h"
38 #include "util/math/matrix3x3.h"
39 #include "util/math/matrix3x3-inl.h"
40
41 DEFINE_int32(num_loops_per_polygon_for_bm,
42 10,
43 "Number of loops per polygon to use for an s2polygon "
44 "encode/decode benchmark. Can be a maximum of 90.");
45
46 // A set of nested loops around the point 0:0 (lat:lng).
47 // Every vertex of kNear0 is a vertex of kNear1.
48 char const kNearPoint[] = "0:0";
49 string const kNear0 = "-1:0, 0:1, 1:0, 0:-1;";
50 string const kNear1 = "-1:-1, -1:0, -1:1, 0:1, 1:1, 1:0, 1:-1, 0:-1;";
51 string const kNear2 = "5:-2, -2:5, -1:-2;";
52 string const kNear3 = "6:-3, -3:6, -2:-2;";
53 string const kNearHemi = "0:-90, -90:0, 0:90, 90:0;";
54
55 // A set of nested loops around the point 0:180 (lat:lng).
56 // Every vertex of kFar0 and kFar2 belongs to kFar1, and all
57 // the loops except kFar2 are non-convex.
58 string const kFar0 = "0:179, 1:180, 0:-179, 2:-180;";
59 string const kFar1 =
60 "0:179, -1:179, 1:180, -1:-179, 0:-179, 3:-178, 2:-180, 3:178;";
61 string const kFar2 = "-1:-179, -1:179, 3:178, 3:-178;"; // opposite direction
62 string const kFar3 = "-3:-178, -2:179, -3:178, 4:177, 4:-177;";
63 string const kFarHemi = "0:-90, 60:90, -60:90;";
64
65 // A set of nested loops around the point -90:0 (lat:lng).
66 string const kSouthPoint = "-89.9999:0.001";
67 string const kSouth0a = "-90:0, -89.99:0, -89.99:0.01;";
68 string const kSouth0b = "-90:0, -89.99:0.02, -89.99:0.03;";
69 string const kSouth0c = "-90:0, -89.99:0.04, -89.99:0.05;";
70 string const kSouth1 = "-90:0, -89.9:-0.1, -89.9:0.1;";
71 string const kSouth2 = "-90:0, -89.8:-0.2, -89.8:0.2;";
72 string const kSouthHemi = "0:-180, 0:60, 0:-60;";
73
74 // Two different loops that surround all the Near and Far loops except
75 // for the hemispheres.
76 string const kNearFar1 = "-1:-9, -9:-9, -9:9, 9:9, 9:-9, 1:-9, "
77 "1:-175, 9:-175, 9:175, -9:175, -9:-175, -1:-175;";
78 string const kNearFar2 = "-8:-4, 8:-4, 2:15, 2:170, "
79 "8:-175, -8:-175, -2:170, -2:15;";
80
81 // Loops that result from intersection of other loops.
82 string const kFarHSouthH = "0:-180, 0:90, -60:90, 0:-90;";
83
84 // Rectangles that form a cross, with only shared vertices, no crossing edges.
85 // Optional holes outside the intersecting region.
86 string const kCross1 = "-2:1, -1:1, 1:1, 2:1, 2:-1, 1:-1, -1:-1, -2:-1;";
87 string const kCross1SideHole = "-1.5:0.5, -1.2:0.5, -1.2:-0.5, -1.5:-0.5;";
88 string const kCross2 = "1:-2, 1:-1, 1:1, 1:2, -1:2, -1:1, -1:-1, -1:-2;";
89 string const kCross2SideHole = "0.5:-1.5, 0.5:-1.2, -0.5:-1.2, -0.5:-1.5;";
90 string const kCrossCenterHole = "-0.5:0.5, 0.5:0.5, 0.5:-0.5, -0.5:-0.5;";
91
92 // Two rectangles that intersect, but no edges cross and there's always
93 // local containment (rather than crossing) at each shared vertex.
94 // In this ugly ASCII art, 1 is A+B, 2 is B+C:
95 // +---+---+---+
96 // | A | B | C |
97 // +---+---+---+
98 string const kOverlap1 = "0:1, 1:1, 2:1, 2:0, 1:0, 0:0;";
99 string const kOverlap1SideHole = "0.2:0.8, 0.8:0.8, 0.8:0.2, 0.2:0.2;";
100 string const kOverlap2 = "1:1, 2:1, 3:1, 3:0, 2:0, 1:0;";
101 string const kOverlap2SideHole = "2.2:0.8, 2.8:0.8, 2.8:0.2, 2.2:0.2;";
102 string const kOverlapCenterHole = "1.2:0.8, 1.8:0.8, 1.8:0.2, 1.2:0.2;";
103
104 class S2PolygonTestBase : public testing::Test {
105 public:
106 S2PolygonTestBase();
107 ~S2PolygonTestBase();
108 protected:
109 // Some standard polygons to use in the tests.
110 S2Polygon const* const near_0;
111 S2Polygon const* const near_10;
112 S2Polygon const* const near_30;
113 S2Polygon const* const near_32;
114 S2Polygon const* const near_3210;
115 S2Polygon const* const near_H3210;
116
117 S2Polygon const* const far_10;
118 S2Polygon const* const far_21;
119 S2Polygon const* const far_321;
120 S2Polygon const* const far_H20;
121 S2Polygon const* const far_H3210;
122
123 S2Polygon const* const south_0ab;
124 S2Polygon const* const south_2;
125 S2Polygon const* const south_210b;
126 S2Polygon const* const south_H21;
127 S2Polygon const* const south_H20abc;
128
129 S2Polygon const* const nf1_n10_f2_s10abc;
130
131 S2Polygon const* const nf2_n2_f210_s210ab;
132
133 S2Polygon const* const f32_n0;
134 S2Polygon const* const n32_s0b;
135
136 S2Polygon const* const cross1;
137 S2Polygon const* const cross1_side_hole;
138 S2Polygon const* const cross1_center_hole;
139 S2Polygon const* const cross2;
140 S2Polygon const* const cross2_side_hole;
141 S2Polygon const* const cross2_center_hole;
142
143 S2Polygon const* const overlap1;
144 S2Polygon const* const overlap1_side_hole;
145 S2Polygon const* const overlap1_center_hole;
146 S2Polygon const* const overlap2;
147 S2Polygon const* const overlap2_side_hole;
148 S2Polygon const* const overlap2_center_hole;
149
150 S2Polygon const* const far_H;
151 S2Polygon const* const south_H;
152 S2Polygon const* const far_H_south_H;
153 };
154
MakePolygon(string const & str)155 static S2Polygon* MakePolygon(string const& str) {
156 scoped_ptr<S2Polygon> polygon(S2Testing::MakePolygon(str));
157 Encoder encoder;
158 polygon->Encode(&encoder);
159 Decoder decoder(encoder.base(), encoder.length());
160 scoped_ptr<S2Polygon> decoded_polygon(new S2Polygon);
161 decoded_polygon->Decode(&decoder);
162 return decoded_polygon.release();
163 }
164
CheckContains(string const & a_str,string const & b_str)165 static void CheckContains(string const& a_str, string const& b_str) {
166 S2Polygon* a = MakePolygon(a_str);
167 S2Polygon* b = MakePolygon(b_str);
168 scoped_ptr<S2Polygon> delete_a(a);
169 scoped_ptr<S2Polygon> delete_b(b);
170 EXPECT_TRUE(a->Contains(b));
171 EXPECT_TRUE(a->ApproxContains(b, S1Angle::Radians(1e-15)));
172 }
173
CheckContainsPoint(string const & a_str,string const & b_str)174 static void CheckContainsPoint(string const& a_str, string const& b_str) {
175 scoped_ptr<S2Polygon> a(S2Testing::MakePolygon(a_str));
176 EXPECT_TRUE(a->VirtualContainsPoint(S2Testing::MakePoint(b_str)))
177 << " " << a_str << " did not contain " << b_str;
178 }
179
TEST(S2Polygon,Init)180 TEST(S2Polygon, Init) {
181 CheckContains(kNear1, kNear0);
182 CheckContains(kNear2, kNear1);
183 CheckContains(kNear3, kNear2);
184 CheckContains(kNearHemi, kNear3);
185 CheckContains(kFar1, kFar0);
186 CheckContains(kFar2, kFar1);
187 CheckContains(kFar3, kFar2);
188 CheckContains(kFarHemi, kFar3);
189 CheckContains(kSouth1, kSouth0a);
190 CheckContains(kSouth1, kSouth0b);
191 CheckContains(kSouth1, kSouth0c);
192 CheckContains(kSouthHemi, kSouth2);
193 CheckContains(kNearFar1, kNear3);
194 CheckContains(kNearFar1, kFar3);
195 CheckContains(kNearFar2, kNear3);
196 CheckContains(kNearFar2, kFar3);
197
198 CheckContainsPoint(kNear0, kNearPoint);
199 CheckContainsPoint(kNear1, kNearPoint);
200 CheckContainsPoint(kNear2, kNearPoint);
201 CheckContainsPoint(kNear3, kNearPoint);
202 CheckContainsPoint(kNearHemi, kNearPoint);
203 CheckContainsPoint(kSouth0a, kSouthPoint);
204 CheckContainsPoint(kSouth1, kSouthPoint);
205 CheckContainsPoint(kSouth2, kSouthPoint);
206 CheckContainsPoint(kSouthHemi, kSouthPoint);
207 }
208
S2PolygonTestBase()209 S2PolygonTestBase::S2PolygonTestBase():
210 near_0(MakePolygon(kNear0)),
211 near_10(MakePolygon(kNear0 + kNear1)),
212 near_30(MakePolygon(kNear3 + kNear0)),
213 near_32(MakePolygon(kNear2 + kNear3)),
214 near_3210(MakePolygon(kNear0 + kNear2 + kNear3 + kNear1)),
215 near_H3210(MakePolygon(kNear0 + kNear2 + kNear3 + kNearHemi + kNear1)),
216
217 far_10(MakePolygon(kFar0 + kFar1)),
218 far_21(MakePolygon(kFar2 + kFar1)),
219 far_321(MakePolygon(kFar2 + kFar3 + kFar1)),
220 far_H20(MakePolygon(kFar2 + kFarHemi + kFar0)),
221 far_H3210(MakePolygon(kFar2 + kFarHemi + kFar0 + kFar1 + kFar3)),
222
223 south_0ab(MakePolygon(kSouth0a + kSouth0b)),
224 south_2(MakePolygon(kSouth2)),
225 south_210b(MakePolygon(kSouth2 + kSouth0b + kSouth1)),
226 south_H21(MakePolygon(kSouth2 + kSouthHemi + kSouth1)),
227 south_H20abc(MakePolygon(
228 kSouth2 + kSouth0b + kSouthHemi + kSouth0a + kSouth0c)),
229
230 nf1_n10_f2_s10abc(MakePolygon(kSouth0c + kFar2 + kNear1 + kNearFar1 +
231 kNear0 + kSouth1 + kSouth0b + kSouth0a)),
232
233 nf2_n2_f210_s210ab(MakePolygon(kFar2 + kSouth0a + kFar1 + kSouth1 + kFar0 +
234 kSouth0b + kNearFar2 + kSouth2 + kNear2)),
235
236 f32_n0(MakePolygon(kFar2 + kNear0 + kFar3)),
237 n32_s0b(MakePolygon(kNear3 + kSouth0b + kNear2)),
238
239 cross1(MakePolygon(kCross1)),
240 cross1_side_hole(MakePolygon(kCross1 + kCross1SideHole)),
241 cross1_center_hole(MakePolygon(kCross1 + kCrossCenterHole)),
242 cross2(MakePolygon(kCross2)),
243 cross2_side_hole(MakePolygon(kCross2 + kCross2SideHole)),
244 cross2_center_hole(MakePolygon(kCross2 + kCrossCenterHole)),
245
246 overlap1(MakePolygon(kOverlap1)),
247 overlap1_side_hole(MakePolygon(kOverlap1 + kOverlap1SideHole)),
248 overlap1_center_hole(MakePolygon(kOverlap1 + kOverlapCenterHole)),
249 overlap2(MakePolygon(kOverlap2)),
250 overlap2_side_hole(MakePolygon(kOverlap2 + kOverlap2SideHole)),
251 overlap2_center_hole(MakePolygon(kOverlap2 + kOverlapCenterHole)),
252
253 far_H(MakePolygon(kFarHemi)),
254 south_H(MakePolygon(kSouthHemi)),
255 far_H_south_H(MakePolygon(kFarHSouthH))
256 {}
257
~S2PolygonTestBase()258 S2PolygonTestBase::~S2PolygonTestBase() {
259 delete near_0;
260 delete near_10;
261 delete near_30;
262 delete near_32;
263 delete near_3210;
264 delete near_H3210;
265
266 delete far_10;
267 delete far_21;
268 delete far_321;
269 delete far_H20;
270 delete far_H3210;
271
272 delete south_0ab;
273 delete south_2;
274 delete south_210b;
275 delete south_H21;
276 delete south_H20abc;
277
278 delete nf1_n10_f2_s10abc;
279
280 delete nf2_n2_f210_s210ab;
281
282 delete f32_n0;
283 delete n32_s0b;
284
285 delete cross1;
286 delete cross1_side_hole;
287 delete cross1_center_hole;
288 delete cross2;
289 delete cross2_side_hole;
290 delete cross2_center_hole;
291
292 delete overlap1;
293 delete overlap1_side_hole;
294 delete overlap1_center_hole;
295 delete overlap2;
296 delete overlap2_side_hole;
297 delete overlap2_center_hole;
298
299 delete far_H;
300 delete south_H;
301 delete far_H_south_H;
302 }
303
CheckEqual(S2Polygon const * a,S2Polygon const * b,double max_error=1e-31)304 static void CheckEqual(S2Polygon const* a, S2Polygon const* b,
305 double max_error = 1e-31) {
306 if (a->IsNormalized() && b->IsNormalized()) {
307 ASSERT_TRUE(a->BoundaryApproxEquals(b, max_error));
308 } else {
309 S2PolygonBuilder builder(S2PolygonBuilderOptions::DIRECTED_XOR());
310 S2Polygon a2, b2;
311 builder.AddPolygon(a);
312 ASSERT_TRUE(builder.AssemblePolygon(&a2, NULL));
313 builder.AddPolygon(b);
314 ASSERT_TRUE(builder.AssemblePolygon(&b2, NULL));
315 ASSERT_TRUE(a2.BoundaryApproxEquals(&b2, max_error));
316 }
317 }
318
TestUnion(S2Polygon const * a,S2Polygon const * b)319 static void TestUnion(S2Polygon const* a, S2Polygon const* b) {
320 S2Polygon c_union;
321 c_union.InitToUnion(a, b);
322
323 vector<S2Polygon*> polygons;
324 polygons.push_back(new S2Polygon);
325 polygons.back()->Copy(a);
326 polygons.push_back(new S2Polygon);
327 polygons.back()->Copy(b);
328 scoped_ptr<S2Polygon> c_destructive_union(
329 S2Polygon::DestructiveUnion(&polygons));
330
331 CheckEqual(&c_union, c_destructive_union.get());
332 }
333
TestContains(S2Polygon const * a,S2Polygon const * b)334 static void TestContains(S2Polygon const* a, S2Polygon const* b) {
335 S2Polygon c, d, e;
336 c.InitToUnion(a, b);
337 CheckEqual(&c, a);
338 TestUnion(a, b);
339
340 d.InitToIntersection(a, b);
341 CheckEqual(&d, b);
342
343 e.InitToDifference(b, a);
344 EXPECT_EQ(0, e.num_loops());
345 }
346
TEST(S2Polygon,TestApproxContains)347 TEST(S2Polygon, TestApproxContains) {
348 // Get a random S2Cell as a polygon.
349 S2CellId id = S2CellId::FromLatLng(S2LatLng::FromE6(40565459, -74645276));
350 S2Cell cell(id.parent(10));
351 S2Polygon cell_as_polygon(cell);
352
353 // We want to roughly bisect the polygon, so we make a rectangle that is the
354 // top half of the current polygon's bounding rectangle.
355 S2LatLngRect const& bounds = cell_as_polygon.GetRectBound();
356 S2LatLngRect upper_half = bounds;
357 upper_half.mutable_lat()->set_lo(bounds.lat().GetCenter());
358
359 // Turn the S2LatLngRect into an S2Polygon
360 vector<S2Point> points;
361 for (int i = 0; i < 4; i++)
362 points.push_back(upper_half.GetVertex(i).ToPoint());
363 vector<S2Loop*> loops;
364 loops.push_back(new S2Loop(points));
365 S2Polygon upper_half_polygon(&loops);
366
367 // Get the intersection. There is no guarantee that the intersection will be
368 // contained by A or B.
369 S2Polygon intersection;
370 intersection.InitToIntersection(&cell_as_polygon, &upper_half_polygon);
371 EXPECT_FALSE(cell_as_polygon.Contains(&intersection));
372
373 EXPECT_TRUE(
374 cell_as_polygon.ApproxContains(&intersection,
375 S2EdgeUtil::kIntersectionTolerance));
376 }
377
TestDisjoint(S2Polygon const * a,S2Polygon const * b)378 static void TestDisjoint(S2Polygon const* a, S2Polygon const* b) {
379 S2PolygonBuilder builder(S2PolygonBuilderOptions::DIRECTED_XOR());
380 builder.AddPolygon(a);
381 builder.AddPolygon(b);
382 S2Polygon ab;
383 ASSERT_TRUE(builder.AssemblePolygon(&ab, NULL));
384
385 S2Polygon c, d, e, f;
386 c.InitToUnion(a, b);
387 CheckEqual(&c, &ab);
388 TestUnion(a, b);
389
390 d.InitToIntersection(a, b);
391 EXPECT_EQ(0, d.num_loops());
392
393 e.InitToDifference(a, b);
394 CheckEqual(&e, a);
395
396 f.InitToDifference(b, a);
397 CheckEqual(&f, b);
398 }
399
TestRelationWithDesc(S2Polygon const * a,S2Polygon const * b,int contains,bool intersects,const char * test_description)400 static void TestRelationWithDesc(S2Polygon const* a, S2Polygon const* b,
401 int contains, bool intersects,
402 const char *test_description) {
403 SCOPED_TRACE(test_description);
404 EXPECT_EQ(contains > 0, a->Contains(b));
405 EXPECT_EQ(contains < 0, b->Contains(a));
406 EXPECT_EQ(intersects, a->Intersects(b));
407 if (contains > 0) {
408 TestContains(a, b);
409 } else if (contains < 0) {
410 TestContains(b, a);
411 }
412 if (!intersects) {
413 TestDisjoint(a, b);
414 }
415 }
416
TEST_F(S2PolygonTestBase,Relations)417 TEST_F(S2PolygonTestBase, Relations) {
418 #define TestRelation(a, b, contains, intersects) \
419 TestRelationWithDesc(a, b, contains, intersects, "args " #a ", " #b)
420 TestRelation(near_10, near_30, -1, true);
421 TestRelation(near_10, near_32, 0, false);
422 TestRelation(near_10, near_3210, -1, true);
423 TestRelation(near_10, near_H3210, 0, false);
424 TestRelation(near_30, near_32, 1, true);
425 TestRelation(near_30, near_3210, 1, true);
426 TestRelation(near_30, near_H3210, 0, true);
427 TestRelation(near_32, near_3210, -1, true);
428 TestRelation(near_32, near_H3210, 0, false);
429 TestRelation(near_3210, near_H3210, 0, false);
430
431 TestRelation(far_10, far_21, 0, false);
432 TestRelation(far_10, far_321, -1, true);
433 TestRelation(far_10, far_H20, 0, false);
434 TestRelation(far_10, far_H3210, 0, false);
435 TestRelation(far_21, far_321, 0, false);
436 TestRelation(far_21, far_H20, 0, false);
437 TestRelation(far_21, far_H3210, -1, true);
438 TestRelation(far_321, far_H20, 0, true);
439 TestRelation(far_321, far_H3210, 0, true);
440 TestRelation(far_H20, far_H3210, 0, true);
441
442 TestRelation(south_0ab, south_2, -1, true);
443 TestRelation(south_0ab, south_210b, 0, true);
444 TestRelation(south_0ab, south_H21, -1, true);
445 TestRelation(south_0ab, south_H20abc, -1, true);
446 TestRelation(south_2, south_210b, 1, true);
447 TestRelation(south_2, south_H21, 0, true);
448 TestRelation(south_2, south_H20abc, 0, true);
449 TestRelation(south_210b, south_H21, 0, true);
450 TestRelation(south_210b, south_H20abc, 0, true);
451 TestRelation(south_H21, south_H20abc, 1, true);
452
453 TestRelation(nf1_n10_f2_s10abc, nf2_n2_f210_s210ab, 0, true);
454 TestRelation(nf1_n10_f2_s10abc, near_32, 1, true);
455 TestRelation(nf1_n10_f2_s10abc, far_21, 0, false);
456 TestRelation(nf1_n10_f2_s10abc, south_0ab, 0, false);
457 TestRelation(nf1_n10_f2_s10abc, f32_n0, 1, true);
458
459 TestRelation(nf2_n2_f210_s210ab, near_10, 0, false);
460 TestRelation(nf2_n2_f210_s210ab, far_10, 1, true);
461 TestRelation(nf2_n2_f210_s210ab, south_210b, 1, true);
462 TestRelation(nf2_n2_f210_s210ab, south_0ab, 1, true);
463 TestRelation(nf2_n2_f210_s210ab, n32_s0b, 1, true);
464
465 TestRelation(cross1, cross2, 0, true);
466 TestRelation(cross1_side_hole, cross2, 0, true);
467 TestRelation(cross1_center_hole, cross2, 0, true);
468 TestRelation(cross1, cross2_side_hole, 0, true);
469 TestRelation(cross1, cross2_center_hole, 0, true);
470 TestRelation(cross1_side_hole, cross2_side_hole, 0, true);
471 TestRelation(cross1_center_hole, cross2_side_hole, 0, true);
472 TestRelation(cross1_side_hole, cross2_center_hole, 0, true);
473 TestRelation(cross1_center_hole, cross2_center_hole, 0, true);
474
475 // These cases, when either polygon has a hole, test a different code path
476 // from the other cases.
477 TestRelation(overlap1, overlap2, 0, true);
478 TestRelation(overlap1_side_hole, overlap2, 0, true);
479 TestRelation(overlap1_center_hole, overlap2, 0, true);
480 TestRelation(overlap1, overlap2_side_hole, 0, true);
481 TestRelation(overlap1, overlap2_center_hole, 0, true);
482 TestRelation(overlap1_side_hole, overlap2_side_hole, 0, true);
483 TestRelation(overlap1_center_hole, overlap2_side_hole, 0, true);
484 TestRelation(overlap1_side_hole, overlap2_center_hole, 0, true);
485 TestRelation(overlap1_center_hole, overlap2_center_hole, 0, true);
486 #undef TestRelation
487 }
488
489 struct TestCase {
490 char const* a;
491 char const* b;
492 char const* a_and_b;
493 char const* a_or_b;
494 char const* a_minus_b;
495 };
496
497 TestCase test_cases[] = {
498 // Two triangles that share an edge.
499 { "4:2, 3:1, 3:3;",
500
501 "3:1, 2:2, 3:3;",
502
503 "",
504
505 "4:2, 3:1, 2:2, 3:3;",
506
507 "4:2, 3:1, 3:3;"
508 },
509
510 // Two vertical bars and a horizontal bar connecting them.
511 { "0:0, 0:2, 3:2, 3:0; 0:3, 0:5, 3:5, 3:3;",
512
513 "1:1, 1:4, 2:4, 2:1;",
514
515 "1:1, 1:2, 2:2, 2:1; 1:3, 1:4, 2:4, 2:3;",
516
517 "0:0, 0:2, 1:2, 1:3, 0:3, 0:5, 3:5, 3:3, 2:3, 2:2, 3:2, 3:0;",
518
519 "0:0, 0:2, 1:2, 1:1, 2:1, 2:2, 3:2, 3:0; "
520 "0:3, 0:5, 3:5, 3:3, 2:3, 2:4, 1:4, 1:3;"
521 },
522
523 // Two vertical bars and two horizontal bars centered around S2::Origin().
524 { "1:88, 1:93, 2:93, 2:88; -1:88, -1:93, 0:93, 0:88;",
525
526 "-2:89, -2:90, 3:90, 3:89; -2:91, -2:92, 3:92, 3:91;",
527
528 "1:89, 1:90, 2:90, 2:89; 1:91, 1:92, 2:92, 2:91; "
529 "-1:89, -1:90, 0:90, 0:89; -1:91, -1:92, 0:92, 0:91;",
530
531 "-1:88, -1:89, -2:89, -2:90, -1:90, -1:91, -2:91, -2:92, -1:92, -1:93,"
532 "0:93, 0:92, 1:92, 1:93, 2:93, 2:92, 3:92, 3:91, 2:91, 2:90, 3:90,"
533 "3:89, 2:89, 2:88, 1:88, 1:89, 0:89, 0:88; "
534 "0:90, 0:91, 1:91, 1:90;",
535
536 "1:88, 1:89, 2:89, 2:88; 1:90, 1:91, 2:91, 2:90; "
537 "1:92, 1:93, 2:93, 2:92; -1:88, -1:89, 0:89, 0:88; "
538 "-1:90, -1:91, 0:91, 0:90; -1:92, -1:93, 0:93, 0:92;"
539 },
540
541 // Two interlocking square doughnuts centered around -S2::Origin().
542 { "-1:-93, -1:-89, 3:-89, 3:-93; 0:-92, 0:-90, 2:-90, 2:-92;",
543
544 "-3:-91, -3:-87, 1:-87, 1:-91; -2:-90, -2:-88, 0:-88, 0:-90;",
545
546 "-1:-91, -1:-90, 0:-90, 0:-91; 0:-90, 0:-89, 1:-89, 1:-90;",
547
548 "-1:-93, -1:-91, -3:-91, -3:-87, 1:-87, 1:-89, 3:-89, 3:-93; "
549 "0:-92, 0:-91, 1:-91, 1:-90, 2:-90, 2:-92; "
550 "-2:-90, -2:-88, 0:-88, 0:-89, -1:-89, -1:-90;",
551
552 "-1:-93, -1:-91, 0:-91, 0:-92, 2:-92, 2:-90, 1:-90, 1:-89, 3:-89, 3:-93; "
553 "-1:-90, -1:-89, 0:-89, 0:-90;"
554 },
555
556 // An incredibly thin triangle intersecting a square, such that the two
557 // intersection points of the triangle with the square are identical.
558 // This results in a degenerate loop that needs to be handled correctly.
559 { "10:44, 10:46, 12:46, 12:44;",
560
561 "11:45, 89:45.00000000000001, 90:45;",
562
563 "", // Empty intersection!
564
565 // Original square with extra vertex, and triangle disappears (due to
566 // default vertex_merge_radius of S2EdgeUtil::kIntersectionTolerance).
567 "10:44, 10:46, 12:46, 12:45, 12:44;",
568
569 "10:44, 10:46, 12:46, 12:45, 12:44;"
570 },
571 };
572
TEST_F(S2PolygonTestBase,Operations)573 TEST_F(S2PolygonTestBase, Operations) {
574 S2Polygon far_south;
575 far_south.InitToIntersection(far_H, south_H);
576 CheckEqual(&far_south, far_H_south_H);
577
578 for (int i = 0; i < arraysize(test_cases); ++i) {
579 SCOPED_TRACE(StringPrintf("Polygon operation test case %d", i));
580 TestCase* test = test_cases + i;
581 scoped_ptr<S2Polygon> a(MakePolygon(test->a));
582 scoped_ptr<S2Polygon> b(MakePolygon(test->b));
583 scoped_ptr<S2Polygon> expected_a_and_b(MakePolygon(test->a_and_b));
584 scoped_ptr<S2Polygon> expected_a_or_b(MakePolygon(test->a_or_b));
585 scoped_ptr<S2Polygon> expected_a_minus_b(MakePolygon(test->a_minus_b));
586
587 // The intersections in the "expected" data were computed in lat-lng
588 // space, while the actual intersections are computed using geodesics.
589 // The error due to this depends on the length and direction of the line
590 // segment being intersected, and how close the intersection is to the
591 // endpoints of the segment. The worst case is for a line segment between
592 // two points at the same latitude, where the intersection point is in the
593 // middle of the segment. In this case the error is approximately
594 // (p * t^2) / 8, where "p" is the absolute latitude in radians, "t" is
595 // the longitude difference in radians, and both "p" and "t" are small.
596 // The test cases all have small latitude and longitude differences.
597 // If "p" and "t" are converted to degrees, the following error bound is
598 // valid as long as (p * t^2 < 150).
599
600 static double const kMaxError = 1e-4;
601
602 S2Polygon a_and_b, a_or_b, a_minus_b;
603 a_and_b.InitToIntersection(a.get(), b.get());
604 CheckEqual(&a_and_b, expected_a_and_b.get(), kMaxError);
605 a_or_b.InitToUnion(a.get(), b.get());
606 TestUnion(a.get(), b.get());
607 CheckEqual(&a_or_b, expected_a_or_b.get(), kMaxError);
608 a_minus_b.InitToDifference(a.get(), b.get());
609 CheckEqual(&a_minus_b, expected_a_minus_b.get(), kMaxError);
610 }
611 }
612
ClearPolylineVector(vector<S2Polyline * > * polylines)613 void ClearPolylineVector(vector<S2Polyline*>* polylines) {
614 for (vector<S2Polyline*>::const_iterator it = polylines->begin();
615 it != polylines->end(); ++it) {
616 delete *it;
617 }
618 polylines->clear();
619 }
620
PolylineIntersectionSharedEdgeTest(const S2Polygon * p,int start_vertex,int direction)621 static void PolylineIntersectionSharedEdgeTest(const S2Polygon *p,
622 int start_vertex,
623 int direction) {
624 SCOPED_TRACE(StringPrintf("Polyline intersection shared edge test "
625 " start=%d direction=%d",
626 start_vertex, direction));
627 vector<S2Point> points;
628 points.push_back(p->loop(0)->vertex(start_vertex));
629 points.push_back(p->loop(0)->vertex(start_vertex + direction));
630 S2Polyline polyline(points);
631 vector<S2Polyline*> polylines;
632 if (direction < 0) {
633 p->IntersectWithPolyline(&polyline, &polylines);
634 EXPECT_EQ(0, polylines.size());
635 ClearPolylineVector(&polylines);
636 p->SubtractFromPolyline(&polyline, &polylines);
637 ASSERT_EQ(1, polylines.size());
638 ASSERT_EQ(2, polylines[0]->num_vertices());
639 EXPECT_EQ(points[0], polylines[0]->vertex(0));
640 EXPECT_EQ(points[1], polylines[0]->vertex(1));
641 } else {
642 p->IntersectWithPolyline(&polyline, &polylines);
643 ASSERT_EQ(1, polylines.size());
644 ASSERT_EQ(2, polylines[0]->num_vertices());
645 EXPECT_EQ(points[0], polylines[0]->vertex(0));
646 EXPECT_EQ(points[1], polylines[0]->vertex(1));
647 ClearPolylineVector(&polylines);
648 p->SubtractFromPolyline(&polyline, &polylines);
649 EXPECT_EQ(0, polylines.size());
650 }
651 ClearPolylineVector(&polylines);
652 }
653
654 // This tests polyline-polyline intersections.
655 // It covers the same edge cases as TestOperations and also adds some
656 // extra tests for shared edges.
TEST_F(S2PolygonTestBase,PolylineIntersection)657 TEST_F(S2PolygonTestBase, PolylineIntersection) {
658 for (int v = 0; v < 3; ++v) {
659 PolylineIntersectionSharedEdgeTest(cross1, v, 1);
660 PolylineIntersectionSharedEdgeTest(cross1, v + 1, -1);
661 PolylineIntersectionSharedEdgeTest(cross1_side_hole, v, 1);
662 PolylineIntersectionSharedEdgeTest(cross1_side_hole, v + 1, -1);
663 }
664
665 // See comments in TestOperations about the vlue of this constant.
666 static double const kMaxError = 1e-4;
667
668 // This duplicates some of the tests in TestOperations by
669 // converting the outline of polygon A to a polyline then intersecting
670 // it with the polygon B. It then converts B to a polyline and intersects
671 // it with A. It then feeds all of the results into a polygon builder and
672 // tests that the output is equal to doing an intersection between A and B.
673 for (int i = 0; i < arraysize(test_cases); ++i) {
674 SCOPED_TRACE(StringPrintf("Polyline intersection test case %d", i));
675 TestCase* test = test_cases + i;
676 scoped_ptr<S2Polygon> a(MakePolygon(test->a));
677 scoped_ptr<S2Polygon> b(MakePolygon(test->b));
678 scoped_ptr<S2Polygon> expected_a_and_b(MakePolygon(test->a_and_b));
679
680 vector<S2Point> points;
681 vector<S2Polyline *> polylines;
682 for (int ab = 0; ab < 2; ab++) {
683 S2Polygon *tmp = ab ? a.get() : b.get();
684 S2Polygon *tmp2 = ab ? b.get() : a.get();
685 for (int l = 0; l < tmp->num_loops(); l++) {
686 points.clear();
687 if (tmp->loop(l)->is_hole()) {
688 for (int v = tmp->loop(l)->num_vertices(); v >=0 ; v--) {
689 points.push_back(tmp->loop(l)->vertex(v));
690 }
691 } else {
692 for (int v = 0; v <= tmp->loop(l)->num_vertices(); v++) {
693 points.push_back(tmp->loop(l)->vertex(v));
694 }
695 }
696 S2Polyline polyline(points);
697 vector<S2Polyline *> tmp;
698 tmp2->IntersectWithPolyline(&polyline, &tmp);
699 polylines.insert(polylines.end(), tmp.begin(), tmp.end());
700 }
701 }
702
703 S2PolygonBuilder builder(S2PolygonBuilderOptions::DIRECTED_XOR());
704 for (int i = 0; i < polylines.size(); i++) {
705 for (int j = 0; j < polylines[i]->num_vertices() - 1; j++) {
706 builder.AddEdge(polylines[i]->vertex(j), polylines[i]->vertex(j + 1));
707 VLOG(3) << " ... Adding edge: " << polylines[i]->vertex(j) << " - " <<
708 polylines[i]->vertex(j + 1);
709 }
710 }
711 ClearPolylineVector(&polylines);
712
713 S2Polygon a_and_b;
714 ASSERT_TRUE(builder.AssemblePolygon(&a_and_b, NULL));
715 CheckEqual(&a_and_b, expected_a_and_b.get(), kMaxError);
716 }
717 }
718
719 // Remove a random polygon from "pieces" and return it.
ChoosePiece(vector<S2Polygon * > * pieces)720 static S2Polygon* ChoosePiece(vector<S2Polygon*> *pieces) {
721 int i = S2Testing::rnd.Uniform(pieces->size());
722 S2Polygon* result = (*pieces)[i];
723 pieces->erase(pieces->begin() + i);
724 return result;
725 }
726
SplitAndAssemble(S2Polygon const * polygon)727 static void SplitAndAssemble(S2Polygon const* polygon) {
728 S2PolygonBuilder builder(S2PolygonBuilderOptions::DIRECTED_XOR());
729 S2Polygon expected;
730 builder.AddPolygon(polygon);
731 ASSERT_TRUE(builder.AssemblePolygon(&expected, NULL));
732
733 for (int iter = 0; iter < 10; ++iter) {
734 S2RegionCoverer coverer;
735 // Compute the minimum level such that the polygon's bounding
736 // cap is guaranteed to be cut.
737 double diameter = 2 * polygon->GetCapBound().angle().radians();
738 int min_level = S2::kMaxWidth.GetMinLevel(diameter);
739
740 // TODO: Choose a level that will have up to 256 cells in the covering.
741 int level = min_level + S2Testing::rnd.Uniform(4);
742 coverer.set_min_level(min_level);
743 coverer.set_max_level(level);
744 coverer.set_max_cells(500);
745
746 vector<S2CellId> cells;
747 coverer.GetCovering(*polygon, &cells);
748 S2CellUnion covering;
749 covering.Init(cells);
750 S2Testing::CheckCovering(*polygon, covering, false);
751 VLOG(2) << cells.size() << " cells in covering";
752 vector<S2Polygon*> pieces;
753 for (int i = 0; i < cells.size(); ++i) {
754 S2Cell cell(cells[i]);
755 S2Polygon window(cell);
756 S2Polygon* piece = new S2Polygon;
757 piece->InitToIntersection(polygon, &window);
758 pieces.push_back(piece);
759 VLOG(4) << "\nPiece " << i << ":\n Window: "
760 << S2Testing::ToString(&window)
761 << "\n Piece: " << S2Testing::ToString(piece);
762 }
763
764 // Now we repeatedly remove two random pieces, compute their union, and
765 // insert the result as a new piece until only one piece is left.
766 //
767 // We don't use S2Polygon::DestructiveUnion() because it joins the pieces
768 // in a mostly deterministic order. We don't just call random_shuffle()
769 // on the pieces and repeatedly join the last two pieces in the vector
770 // because this always joins a single original piece to the current union
771 // rather than doing the unions according to a random tree structure.
772 while (pieces.size() > 1) {
773 scoped_ptr<S2Polygon> a(ChoosePiece(&pieces));
774 scoped_ptr<S2Polygon> b(ChoosePiece(&pieces));
775 S2Polygon* c = new S2Polygon;
776 c->InitToUnion(a.get(), b.get());
777 pieces.push_back(c);
778 VLOG(4) << "\nJoining piece a: " << S2Testing::ToString(a.get())
779 << "\n With piece b: " << S2Testing::ToString(b.get())
780 << "\n To get piece c: " << S2Testing::ToString(c);
781 }
782 scoped_ptr<S2Polygon> result(pieces[0]);
783 pieces.pop_back();
784
785 // The moment of truth!
786 EXPECT_TRUE(expected.BoundaryNear(result.get()))
787 << "\nActual:\n" << S2Testing::ToString(result.get())
788 << "\nExpected:\n" << S2Testing::ToString(&expected);
789 }
790 }
791
TEST_F(S2PolygonTestBase,Splitting)792 TEST_F(S2PolygonTestBase, Splitting) {
793 // It takes too long to test all the polygons in debug mode, so we just pick
794 // out some of the more interesting ones.
795
796 SplitAndAssemble(near_H3210);
797 SplitAndAssemble(far_H3210);
798 SplitAndAssemble(south_0ab);
799 SplitAndAssemble(south_210b);
800 SplitAndAssemble(south_H20abc);
801 SplitAndAssemble(nf1_n10_f2_s10abc);
802 SplitAndAssemble(nf2_n2_f210_s210ab);
803 SplitAndAssemble(far_H);
804 SplitAndAssemble(south_H);
805 SplitAndAssemble(far_H_south_H);
806 }
807
TEST(S2Polygon,InitToCellUnionBorder)808 TEST(S2Polygon, InitToCellUnionBorder) {
809 // Test S2Polygon::InitToCellUnionBorder().
810 // The main thing to check is that adjacent cells of different sizes get
811 // merged correctly. To do this we generate two random adjacent cells,
812 // convert to polygon, and make sure the polygon only has a single loop.
813 for (int iter = 0; iter < 500; ++iter) {
814 SCOPED_TRACE(StringPrintf("Iteration %d", iter));
815
816 // Choose a random non-leaf cell.
817 S2CellId big_cell =
818 S2Testing::GetRandomCellId(S2Testing::rnd.Uniform(S2CellId::kMaxLevel));
819 // Get all neighbors at some smaller level.
820 int small_level = big_cell.level() +
821 S2Testing::rnd.Uniform(min(16, S2CellId::kMaxLevel - big_cell.level()));
822 vector<S2CellId> neighbors;
823 big_cell.AppendAllNeighbors(small_level, &neighbors);
824 // Pick one at random.
825 S2CellId small_cell = neighbors[S2Testing::rnd.Uniform(neighbors.size())];
826 // If it's diagonally adjacent, bail out.
827 S2CellId edge_neighbors[4];
828 big_cell.GetEdgeNeighbors(edge_neighbors);
829 bool diagonal = true;
830 for (int i = 0; i < 4; ++i) {
831 if (edge_neighbors[i].contains(small_cell)) {
832 diagonal = false;
833 }
834 }
835 VLOG(3) << iter << ": big_cell " << big_cell <<
836 " small_cell " << small_cell;
837 if (diagonal) {
838 VLOG(3) << " diagonal - bailing out!";
839 continue;
840 }
841
842 vector<S2CellId> cells;
843 cells.push_back(big_cell);
844 cells.push_back(small_cell);
845 S2CellUnion cell_union;
846 cell_union.Init(cells);
847 EXPECT_EQ(2, cell_union.num_cells());
848 S2Polygon poly;
849 poly.InitToCellUnionBorder(cell_union);
850 EXPECT_EQ(1, poly.num_loops());
851 // If the conversion were perfect we could test containment, but due to
852 // rounding the polygon won't always exactly contain both cells. We can
853 // at least test intersection.
854 EXPECT_TRUE(poly.MayIntersect(S2Cell(big_cell)));
855 EXPECT_TRUE(poly.MayIntersect(S2Cell(small_cell)));
856 }
857 }
858
TEST_F(S2PolygonTestBase,TestEncodeDecode)859 TEST_F(S2PolygonTestBase, TestEncodeDecode) {
860 Encoder encoder;
861 cross1->Encode(&encoder);
862 Decoder decoder(encoder.base(), encoder.length());
863 S2Polygon decoded_polygon;
864 ASSERT_TRUE(decoded_polygon.Decode(&decoder));
865 EXPECT_TRUE(cross1->BoundaryEquals(&decoded_polygon));
866 EXPECT_EQ(cross1->GetRectBound(), decoded_polygon.GetRectBound());
867 }
868
869 // This test checks that S2Polygons created directly from S2Cells behave
870 // identically to S2Polygons created from the vertices of those cells; this
871 // previously was not the case, because S2Cells calculate their bounding
872 // rectangles slightly differently, and S2Polygons created from them just
873 // copied the S2Cell bounds.
TEST(S2Polygon,TestS2CellConstructorAndContains)874 TEST(S2Polygon, TestS2CellConstructorAndContains) {
875 S2LatLng latlng(S1Angle::E6(40565459), S1Angle::E6(-74645276));
876 S2Cell cell(S2CellId::FromLatLng(latlng));
877 S2Polygon cell_as_polygon(cell);
878 S2Polygon empty;
879 S2Polygon polygon_copy;
880 polygon_copy.InitToUnion(&cell_as_polygon, &empty);
881 EXPECT_TRUE(polygon_copy.Contains(&cell_as_polygon));
882 EXPECT_TRUE(polygon_copy.Contains(cell));
883 }
884
TEST(S2PolygonTest,Project)885 TEST(S2PolygonTest, Project) {
886 scoped_ptr<S2Polygon> polygon(MakePolygon(kNear0 + kNear2));
887 S2Point point;
888 S2Point projected;
889
890 // The point inside the polygon should be projected into itself.
891 point = S2Testing::MakePoint("1.1:0");
892 projected = polygon->Project(point);
893 EXPECT_TRUE(S2::ApproxEquals(point, projected));
894
895 // The point is on the outside of the polygon.
896 point = S2Testing::MakePoint("5.1:-2");
897 projected = polygon->Project(point);
898 EXPECT_TRUE(S2::ApproxEquals(S2Testing::MakePoint("5:-2"), projected));
899
900 // The point is inside the hole in the polygon.
901 point = S2Testing::MakePoint("-0.49:-0.49");
902 projected = polygon->Project(point);
903 EXPECT_TRUE(S2::ApproxEquals(S2Testing::MakePoint("-0.5:-0.5"),
904 projected, 1e-6));
905
906 point = S2Testing::MakePoint("0:-3");
907 projected = polygon->Project(point);
908 EXPECT_TRUE(S2::ApproxEquals(S2Testing::MakePoint("0:-2"), projected));
909 }
910
911 // Returns the distance of a point to a polygon (distance is 0 if the
912 // point is in the polygon).
DistanceToPolygonInDegrees(S2Point point,S2Polygon const & poly)913 double DistanceToPolygonInDegrees(S2Point point, S2Polygon const& poly) {
914 S1Angle distance = S1Angle(poly.Project(point), point);
915 return distance.degrees();
916 }
917
918 // Returns the diameter of a loop (maximum distance between any two
919 // points in the loop).
LoopDiameter(S2Loop const & loop)920 S1Angle LoopDiameter(S2Loop const& loop) {
921 S1Angle diameter = S1Angle();
922 for (int i = 0; i < loop.num_vertices(); ++i) {
923 S2Point test_point = loop.vertex(i);
924 for (int j = i + 1; j < loop.num_vertices(); ++j) {
925 diameter = max(diameter,
926 S2EdgeUtil::GetDistance(test_point, loop.vertex(j),
927 loop.vertex(j+1)));
928 }
929 }
930 return diameter;
931 }
932
933 // Returns the maximum distance from any vertex of poly_a to poly_b, that is,
934 // the directed Haussdorf distance of the set of vertices of poly_a to the
935 // boundary of poly_b.
936 //
937 // Doesn't consider loops from poly_a that have diameter less than min_diameter
938 // in degrees.
MaximumDistanceInDegrees(S2Polygon const & poly_a,S2Polygon const & poly_b,double min_diameter_in_degrees)939 double MaximumDistanceInDegrees(S2Polygon const& poly_a,
940 S2Polygon const& poly_b,
941 double min_diameter_in_degrees) {
942 double min_distance = 360;
943 bool has_big_loops = false;
944 for (int l = 0; l < poly_a.num_loops(); ++l) {
945 S2Loop* a_loop = poly_a.loop(l);
946 if (LoopDiameter(*a_loop).degrees() <= min_diameter_in_degrees) {
947 continue;
948 }
949 has_big_loops = true;
950 for (int v = 0; v < a_loop->num_vertices(); ++v) {
951 double distance =
952 DistanceToPolygonInDegrees(a_loop->vertex(v), poly_b);
953 if (distance < min_distance) {
954 min_distance = distance;
955 }
956 }
957 }
958 if (has_big_loops) {
959 return min_distance;
960 } else {
961 return 0.; // As if the first polygon were empty.
962 }
963 }
964
965 class S2PolygonSimplifierTest : public ::testing::Test {
966 protected:
S2PolygonSimplifierTest()967 S2PolygonSimplifierTest() {
968 simplified = NULL;
969 original = NULL;
970 }
971
~S2PolygonSimplifierTest()972 ~S2PolygonSimplifierTest() {
973 delete simplified;
974 delete original;
975 }
976
977 // Owns poly.
SetInput(S2Polygon * poly,double tolerance_in_degrees)978 void SetInput(S2Polygon* poly, double tolerance_in_degrees) {
979 delete original;
980 delete simplified;
981 original = poly;
982
983 simplified = new S2Polygon();
984 return simplified->InitToSimplified(original,
985 S1Angle::Degrees(tolerance_in_degrees));
986 }
987
SetInput(string const & poly,double tolerance_in_degrees)988 void SetInput(string const& poly, double tolerance_in_degrees) {
989 SetInput(S2Testing::MakePolygon(poly), tolerance_in_degrees);
990 }
991
992 S2Polygon* simplified;
993 S2Polygon* original;
994 };
995
TEST_F(S2PolygonSimplifierTest,NoSimplification)996 TEST_F(S2PolygonSimplifierTest, NoSimplification) {
997 SetInput("0:0, 0:20, 20:20, 20:0", 1.0);
998 EXPECT_EQ(4, simplified->num_vertices());
999
1000 EXPECT_EQ(0, MaximumDistanceInDegrees(*simplified, *original, 0));
1001 EXPECT_EQ(0, MaximumDistanceInDegrees(*original, *simplified, 0));
1002 }
1003
1004 // Here, 10:-2 will be removed and 0:0-20:0 will intersect two edges.
1005 // (The resulting polygon will in fact probably have more edges.)
TEST_F(S2PolygonSimplifierTest,SimplifiedLoopSelfIntersects)1006 TEST_F(S2PolygonSimplifierTest, SimplifiedLoopSelfIntersects) {
1007 SetInput("0:0, 0:20, 10:-0.1, 20:20, 20:0, 10:-0.2", 0.22);
1008
1009 // To make sure that this test does something, we check
1010 // that the vertex 10:-0.2 is not in the simplification anymore.
1011 S2Point test_point = S2Testing::MakePoint("10:-0.2");
1012 EXPECT_LT(0.05, DistanceToPolygonInDegrees(test_point, *simplified));
1013
1014 EXPECT_GE(0.22, MaximumDistanceInDegrees(*simplified, *original, 0));
1015 EXPECT_GE(0.22, MaximumDistanceInDegrees(*original, *simplified, 0.22));
1016 }
1017
TEST_F(S2PolygonSimplifierTest,NoSimplificationManyLoops)1018 TEST_F(S2PolygonSimplifierTest, NoSimplificationManyLoops) {
1019 SetInput("0:0, 0:1, 1:0; 0:20, 0:21, 1:20; "
1020 "20:20, 20:21, 21:20; 20:0, 20:1, 21:0", 0.01);
1021 EXPECT_EQ(0, MaximumDistanceInDegrees(*simplified, *original, 0));
1022 EXPECT_EQ(0, MaximumDistanceInDegrees(*original, *simplified, 0));
1023 }
1024
TEST_F(S2PolygonSimplifierTest,TinyLoopDisappears)1025 TEST_F(S2PolygonSimplifierTest, TinyLoopDisappears) {
1026 SetInput("0:0, 0:1, 1:1, 1:0", 1.1);
1027 EXPECT_EQ(0, simplified->num_vertices());
1028 }
1029
TEST_F(S2PolygonSimplifierTest,StraightLinesAreSimplified)1030 TEST_F(S2PolygonSimplifierTest, StraightLinesAreSimplified) {
1031 SetInput("0:0, 1:0, 2:0, 3:0, 4:0, 5:0, 6:0,"
1032 "6:1, 5:1, 4:1, 3:1, 2:1, 1:1, 0:1", 0.01);
1033 EXPECT_EQ(4, simplified->num_vertices());
1034 }
1035
TEST_F(S2PolygonSimplifierTest,EdgeSplitInManyPieces)1036 TEST_F(S2PolygonSimplifierTest, EdgeSplitInManyPieces) {
1037 // near_square's right four-point side will be simplified to a vertical
1038 // line at lng=7.9, that will cut the 9 teeth of the saw (the edge will
1039 // therefore be broken into 19 pieces).
1040 const string saw =
1041 "1:1, 1:8, 2:2, 2:8, 3:2, 3:8, 4:2, 4:8, 5:2, 5:8,"
1042 "6:2, 6:8, 7:2, 7:8, 8:2, 8:8, 9:2, 9:8, 10:1";
1043 const string near_square =
1044 "0:0, 0:7.9, 1:8.1, 10:8.1, 11:7.9, 11:0";
1045 SetInput(saw + ";" + near_square, 0.21);
1046
1047 EXPECT_TRUE(simplified->IsValid());
1048 EXPECT_GE(0.11, MaximumDistanceInDegrees(*simplified, *original, 0));
1049 EXPECT_GE(0.11, MaximumDistanceInDegrees(*original, *simplified, 0));
1050 // The resulting polygon's 9 little teeth are very small and disappear
1051 // due to the vertex_merge_radius of the polygon builder. There remains
1052 // nine loops.
1053 EXPECT_EQ(9, simplified->num_loops());
1054 }
1055
TEST_F(S2PolygonSimplifierTest,EdgesOverlap)1056 TEST_F(S2PolygonSimplifierTest, EdgesOverlap) {
1057 // Two loops, One edge of the second one ([0:1 - 0:2]) is part of an
1058 // edge of the first one..
1059 SetInput("0:0, 0:3, 1:0; 0:1, -1:1, 0:2", 0.01);
1060 scoped_ptr<S2Polygon> true_poly(
1061 S2Testing::MakePolygon("0:3, 1:0, 0:0, 0:1, -1:1, 0:2"));
1062 EXPECT_TRUE(simplified->BoundaryApproxEquals(true_poly.get()));
1063 }
1064
MakeRegularPolygon(const string & center,int num_points,double radius_in_degrees)1065 S2Polygon* MakeRegularPolygon(const string& center,
1066 int num_points, double radius_in_degrees) {
1067
1068 double radius_in_radians = S1Angle::Degrees(radius_in_degrees).radians();
1069 S2Loop* l = S2Testing::MakeRegularLoop(S2Testing::MakePoint(center),
1070 num_points,
1071 radius_in_radians);
1072 vector<S2Loop*> loops;
1073 loops.push_back(l);
1074 return new S2Polygon(&loops);
1075 }
1076
1077 // Tests that a regular polygon with many points gets simplified
1078 // enough.
TEST_F(S2PolygonSimplifierTest,LargeRegularPolygon)1079 TEST_F(S2PolygonSimplifierTest, LargeRegularPolygon) {
1080 const double kRadius = 2.; // in degrees
1081 const int num_initial_points = 1000;
1082 const int num_desired_points = 250;
1083 double tolerance = 1.05 * kRadius * (1 - cos(M_PI / num_desired_points));
1084
1085 S2Polygon* p = MakeRegularPolygon("0:0", num_initial_points, kRadius);
1086 SetInput(p, tolerance);
1087
1088 EXPECT_GE(tolerance, MaximumDistanceInDegrees(*simplified, *original, 0));
1089 EXPECT_GE(tolerance, MaximumDistanceInDegrees(*original, *simplified, 0));
1090 EXPECT_GE(250, simplified->num_vertices());
1091 EXPECT_LE(200, simplified->num_vertices());
1092 }
1093
GenerateInputForBenchmark(int num_vertices_per_loop_for_bm)1094 string GenerateInputForBenchmark(int num_vertices_per_loop_for_bm) {
1095 CHECK_LE(FLAGS_num_loops_per_polygon_for_bm, 90);
1096 vector<S2Loop*> loops;
1097 for (int li = 0; li < FLAGS_num_loops_per_polygon_for_bm; ++li) {
1098 vector<S2Point> vertices;
1099 double radius_degrees =
1100 1.0 + (50.0 * li) / FLAGS_num_loops_per_polygon_for_bm;
1101 for (int vi = 0; vi < num_vertices_per_loop_for_bm; ++vi) {
1102 double angle_radians = (2 * M_PI * vi) / num_vertices_per_loop_for_bm;
1103 double lat = radius_degrees * cos(angle_radians);
1104 double lng = radius_degrees * sin(angle_radians);
1105 vertices.push_back(S2LatLng::FromDegrees(lat, lng).ToPoint());
1106 }
1107 loops.push_back(new S2Loop(vertices));
1108 }
1109 S2Polygon polygon_to_encode(&loops);
1110
1111 Encoder encoder;
1112 polygon_to_encode.Encode(&encoder);
1113 string encoded(encoder.base(), encoder.length());
1114
1115 return encoded;
1116 }
1117
BM_S2Decoding(int iters,int num_vertices_per_loop_for_bm)1118 static void BM_S2Decoding(int iters, int num_vertices_per_loop_for_bm) {
1119 StopBenchmarkTiming();
1120 string encoded = GenerateInputForBenchmark(num_vertices_per_loop_for_bm);
1121 StartBenchmarkTiming();
1122 for (int i = 0; i < iters; ++i) {
1123 Decoder decoder(encoded.data(), encoded.size());
1124 S2Polygon decoded_polygon;
1125 decoded_polygon.Decode(&decoder);
1126 }
1127 }
1128 BENCHMARK_RANGE(BM_S2Decoding, 8, 131072);
1129
BM_S2DecodingWithinScope(int iters,int num_vertices_per_loop_for_bm)1130 static void BM_S2DecodingWithinScope(int iters,
1131 int num_vertices_per_loop_for_bm) {
1132 StopBenchmarkTiming();
1133 string encoded = GenerateInputForBenchmark(num_vertices_per_loop_for_bm);
1134 StartBenchmarkTiming();
1135 for (int i = 0; i < iters; ++i) {
1136 Decoder decoder(encoded.data(), encoded.size());
1137 S2Polygon decoded_polygon;
1138 decoded_polygon.DecodeWithinScope(&decoder);
1139 }
1140 }
1141 BENCHMARK_RANGE(BM_S2DecodingWithinScope, 8, 131072);
1142
1143
ConcentricLoops(S2Point center,int num_loops,int num_vertices_per_loop,S2Polygon * poly)1144 void ConcentricLoops(S2Point center,
1145 int num_loops,
1146 int num_vertices_per_loop,
1147 S2Polygon* poly) {
1148 Matrix3x3_d m;
1149 S2::GetFrame(center, &m);
1150 vector<S2Loop*> loops;
1151 for (int li = 0; li < num_loops; ++li) {
1152 vector<S2Point> vertices;
1153 double radius = 0.005 * (li + 1) / num_loops;
1154 double radian_step = 2 * M_PI / num_vertices_per_loop;
1155 for (int vi = 0; vi < num_vertices_per_loop; ++vi) {
1156 double angle = vi * radian_step;
1157 S2Point p(radius * cos(angle), radius * sin(angle), 1);
1158 vertices.push_back(S2::FromFrame(m, p.Normalize()));
1159 }
1160 loops.push_back(new S2Loop(vertices));
1161 }
1162 poly->Init(&loops);
1163 }
1164
UnionOfPolygons(int iters,int num_vertices_per_loop,double second_polygon_offset)1165 static void UnionOfPolygons(int iters,
1166 int num_vertices_per_loop,
1167 double second_polygon_offset) {
1168 for (int i = 0; i < iters; ++i) {
1169 StopBenchmarkTiming();
1170 S2Polygon p1, p2;
1171 S2Point center = S2Testing::RandomPoint();
1172 ConcentricLoops(center,
1173 FLAGS_num_loops_per_polygon_for_bm,
1174 num_vertices_per_loop,
1175 &p1);
1176 ConcentricLoops(
1177 (center + S2Point(second_polygon_offset,
1178 second_polygon_offset,
1179 second_polygon_offset)).Normalize(),
1180 FLAGS_num_loops_per_polygon_for_bm,
1181 num_vertices_per_loop,
1182 &p2);
1183 StartBenchmarkTiming();
1184 S2Polygon p_result;
1185 p_result.InitToUnion(&p1, &p2);
1186 }
1187 }
1188
BM_DeepPolygonUnion(int iters,int num_vertices_per_loop)1189 static void BM_DeepPolygonUnion(int iters, int num_vertices_per_loop) {
1190 UnionOfPolygons(iters, num_vertices_per_loop, 0.000001);
1191 }
1192 BENCHMARK(BM_DeepPolygonUnion)
1193 ->Arg(8)
1194 ->Arg(64)
1195 ->Arg(128)
1196 ->Arg(256)
1197 ->Arg(512)
1198 ->Arg(1024)
1199 ->Arg(4096)
1200 ->Arg(8192);
1201
BM_ShallowPolygonUnion(int iters,int num_vertices_per_loop)1202 static void BM_ShallowPolygonUnion(int iters, int num_vertices_per_loop) {
1203 UnionOfPolygons(iters, num_vertices_per_loop, 0.004);
1204 }
1205 BENCHMARK(BM_ShallowPolygonUnion)
1206 ->Arg(8)
1207 ->Arg(64)
1208 ->Arg(128)
1209 ->Arg(256)
1210 ->Arg(512)
1211 ->Arg(1024)
1212 ->Arg(4096)
1213 ->Arg(8192);
1214
BM_DisjointPolygonUnion(int iters,int num_vertices_per_loop)1215 static void BM_DisjointPolygonUnion(int iters, int num_vertices_per_loop) {
1216 UnionOfPolygons(iters, num_vertices_per_loop, 0.3);
1217 }
1218 BENCHMARK(BM_DisjointPolygonUnion)
1219 ->Arg(8)
1220 ->Arg(64)
1221 ->Arg(128)
1222 ->Arg(256)
1223 ->Arg(512)
1224 ->Arg(1024)
1225 ->Arg(4096)
1226 ->Arg(8192);
1227