1// Copyright 2009 Google Inc. All Rights Reserved.
2// julienbasch@google.com (Julien Basch)
3
4// Implementation of class S2EdgeIndex, a fast lookup structure for edges in S2.
5//
6// An object of this class contains a set S of edges called the test edges.
7// For a query edge q, you want to compute a superset of all test edges that
8// intersect q.
9//
10// The idea is roughly that of
11// Each edge is covered by one or several S2 cells, stored in a multimap
12// cell -> edge*.
13// To perform a query, you cover the query edge with a set of cells. For
14// each such cell c, you find all test edges that are in c,in an ancestor of c
15// or in a child of c.
16//
17// This is simple, but there are two complications:
18//
19// 1. For containment queries, the query edge is very long (from S2::Origin()
20// to the query point). A standard cell covering of q is either useless or
21// too large. The covering needs to be adapted to S: if a cell contains too
22// many edges from S, you subdivide it and keep only the subcells that
23// intersect q. See comments for FindCandidateCrossings().
24//
25// 2. To decide if edge q could possibly cross edge e, we end up comparing
26// both with edges that bound s2 cells. Numerical inaccuracies
27// can lead to inconcistencies, e.g.: there may be an edge b at the
28// boundary of two cells such that q and e are on opposite sides of b,
29// yet they cross each other. This special case happens a lot if your
30// test and query edges are cell boundaries themselves, and this in turn
31// is a common case when regions are approximated by cell unions.
32//
33// We expand here on the solution to the second problem. Two components:
34//
35// 1. Each test edge is thickened to a rectangle before it is S2-covered.
36// See the comment for GetThickenedEdgeCovering().
37//
38// 2. When recursing through the children of a cell c for a query edge q,
39// we test q against the boundaries of c's children in a 'lenient'
40// way. That is, instead of testing e.g. area(abc)*area(abd) < 0,
41// we check if it is 'approximately negative'.
42//
43// To see how the second point is necessary, imagine that your query
44// edge q is the North boundary of cell x. We recurse into the four
45// children a,b,c,d of x. To do so, we check if q crosses or touches any
46// of a,b,c or d boundaries. As all the situations are degenerate, it is
47// possible that all crossing tests return false, thus making q suddenly
48// 'disappear'. Using the lenient crossing test, we are guaranteed that q
49// will intersect one of the four edges of the cross that bounds a,b,c,d.
50// The same holds true if q passes through the cell center of x.
51
52
53
54#include "s2edgeindex.h"
55
56#include <algorithm>
57using std::min;
58using std::max;
59using std::swap;
60using std::reverse;
61
62#include <set>
63using std::set;
64using std::multiset;
65
66#include <utility>
67using std::pair;
68using std::make_pair;
69
70
71// #include "base/commandlineflags.h"
72#include "base/logging.h"
73#include "s2cell.h"
74#include "s2edgeutil.h"
75#include "s2polyline.h"
76#include "s2regioncoverer.h"
77
78
79// DEFINE_bool(always_recurse_on_children, false,
80// "When we test a query edge against a cell, we don't "
81// "recurse if there are only a few test edges in it. "
82// "For testing, it is useful to always recurse to the end. "
83// "You don't want to use this flag anywhere but in tests.");
84static bool FLAGS_always_recurse_on_children = false;
85
86void S2EdgeIndex::Reset() {
87 minimum_s2_level_used_ = S2CellId::kMaxLevel;
88 index_computed_ = false;
89 query_count_ = 0;
90 mapping_.clear();
91}
92
93void S2EdgeIndex::ComputeIndex() {
94 DCHECK(!index_computed_);
95
96 for (int i = 0; i < num_edges(); ++i) {
97 S2Point from, to;
98 vector<S2CellId> cover;
99 int level = GetCovering(*edge_from(i), *edge_to(i),
100 true, &cover);
101 minimum_s2_level_used_ = min(minimum_s2_level_used_, level);
102
103 for (vector<S2CellId>::const_iterator it = cover.begin(); it != cover.end();
104 ++it) {
105 mapping_.insert(make_pair(*it, i));
106 }
107 }
108 index_computed_ = true;
109}
110
111bool S2EdgeIndex::IsIndexComputed() const {
112 return index_computed_;
113}
114
115void S2EdgeIndex::IncrementQueryCount() {
116 query_count_++;
117}
118
119
120// If we have m data edges and n query edges, then the brute force cost is
121// m * n * test_cost
122// where test_cost is taken to be the cost of EdgeCrosser::RobustCrossing,
123// measured to be about 30ns at the time of this writing.
124//
125// If we compute the index, the cost becomes:
126// m * cost_insert + n * cost_find(m)
127//
128// - cost_insert can be expected to be reasonably stable, and was measured
129// at 1200ns with the BM_QuadEdgeInsertionCost benchmark.
130//
131// - cost_find depends on the length of the edge . For m=1000 edges,
132// we got timings ranging from 1ms (edge the length of the polygon) to
133// 40ms. The latter is for very long query edges, and needs to be
134// optimized. We will assume for the rest of the discussion that
135// cost_find is roughly 3ms.
136//
137// When doing one additional query, the differential cost is
138// m * test_cost - cost_find(m)
139// With the numbers above, it is better to use the quad tree (if we have it)
140// if m >= 100.
141//
142// If m = 100, 30 queries will give m*n*test_cost = m*cost_insert = 100ms,
143// while the marginal cost to find is 3ms. Thus, this is a reasonable
144// thing to do.
145void S2EdgeIndex::PredictAdditionalCalls(int n) {
146 if (index_computed_) return;
147 if (num_edges() > 100 && (query_count_ + n) > 30) {
148 ComputeIndex();
149 }
150}
151
152void S2EdgeIndex::GetEdgesInParentCells(
153 const vector<S2CellId>& cover,
154 const CellEdgeMultimap& mapping,
155 int minimum_s2_level_used,
156 vector<int>* candidate_crossings) {
157 // Find all parent cells of covering cells.
158 set<S2CellId> parent_cells;
159 for (vector<S2CellId>::const_iterator it = cover.begin(); it != cover.end();
160 ++it) {
161 for (int parent_level = it->level() - 1;
162 parent_level >= minimum_s2_level_used;
163 --parent_level) {
164 if (!parent_cells.insert(it->parent(parent_level)).second) {
165 break; // cell is already in => parents are too.
166 }
167 }
168 }
169
170 // Put parent cell edge references into result.
171 for (set<S2CellId>::const_iterator it = parent_cells.begin(); it
172 != parent_cells.end(); ++it) {
173 pair<CellEdgeMultimap::const_iterator,
174 CellEdgeMultimap::const_iterator> range =
175 mapping.equal_range(*it);
176 for (CellEdgeMultimap::const_iterator it2 = range.first;
177 it2 != range.second; ++it2) {
178 candidate_crossings->push_back(it2->second);
179 }
180 }
181}
182
183// Returns true if ab possibly crosses cd, by clipping tiny angles to
184// zero.
185static bool LenientCrossing(S2Point const& a, S2Point const& b,
186 S2Point const& c, S2Point const& d) {
187 DCHECK(S2::IsUnitLength(a));
188 DCHECK(S2::IsUnitLength(b));
189 DCHECK(S2::IsUnitLength(c));
190 // See comment for RobustCCW() in s2.h
191 const double kMaxDetError = 1.e-14;
192 double acb = a.CrossProd(c).DotProd(b);
193 double bda = b.CrossProd(d).DotProd(a);
194 if (fabs(acb) < kMaxDetError || fabs(bda) < kMaxDetError) {
195 return true;
196 }
197 if (acb * bda < 0) return false;
198 double cbd = c.CrossProd(b).DotProd(d);
199 double dac = d.CrossProd(a).DotProd(c);
200 if (fabs(cbd) < kMaxDetError || fabs(dac) < kMaxDetError) {
201 return true;
202 }
203 return (acb * cbd >= 0) && (acb * dac >= 0);
204}
205
206bool S2EdgeIndex::EdgeIntersectsCellBoundary(
207 S2Point const& a, S2Point const& b, const S2Cell& cell) {
208 S2Point start_vertex = cell.GetVertex(0);
209
210 S2Point vertices[4];
211 for (int i = 0; i < 4; ++i) {
212 vertices[i] = cell.GetVertex(i);
213 }
214 for (int i = 0; i < 4; ++i) {
215 S2Point from_point = vertices[i];
216 S2Point to_point = vertices[(i+1) % 4];
217 if (LenientCrossing(a, b, from_point, to_point)) {
218 return true;
219 }
220 }
221 return false;
222}
223
224void S2EdgeIndex::GetEdgesInChildrenCells(
225 S2Point const& a, S2Point const& b,
226 vector<S2CellId>* cover,
227 const CellEdgeMultimap& mapping,
228 vector<int>* candidate_crossings) {
229 CellEdgeMultimap::const_iterator it, start, end;
230
231 int num_cells = 0;
232
233 // Put all edge references of (covering cells + descendant cells) into result.
234 // This relies on the natural ordering of S2CellIds.
235 while (!cover->empty()) {
236 S2CellId cell = cover->back();
237 cover->pop_back();
238 num_cells++;
239 start = mapping.lower_bound(cell.range_min());
240 end = mapping.upper_bound(cell.range_max());
241 int num_edges = 0;
242 bool rewind = FLAGS_always_recurse_on_children;
243 // TODO(user): Maybe distinguish between edges in current cell, that
244 // are going to be added anyhow, and edges in subcells, and rewind only
245 // those.
246 if (!rewind) {
247 for (it = start; it != end; ++it) {
248 candidate_crossings->push_back(it->second);
249 ++num_edges;
250 if (num_edges == 16 && !cell.is_leaf()) {
251 rewind = true;
252 break;
253 }
254 }
255 }
256 // If there are too many to insert, uninsert and recurse.
257 if (rewind) {
258 for (int i = 0; i < num_edges; ++i) {
259 candidate_crossings->pop_back();
260 }
261 // Add cells at this level
262 pair<CellEdgeMultimap::const_iterator,
263 CellEdgeMultimap::const_iterator> eq =
264 mapping.equal_range(cell);
265 for (it = eq.first; it != eq.second; ++it) {
266 candidate_crossings->push_back(it->second);
267 }
268 // Recurse on the children -- hopefully some will be empty.
269 if (eq.first != start || eq.second != end) {
270 S2Cell children[4];
271 S2Cell c(cell);
272 c.Subdivide(children);
273 for (int i = 0; i < 4; ++i) {
274 // TODO(user): Do the check for the four cells at once,
275 // as it is enough to check the four edges between the cells. At
276 // this time, we are checking 16 edges, 4 times too many.
277 //
278 // Note that given the guarantee of AppendCovering, it is enough
279 // to check that the edge intersect with the cell boundary as it
280 // cannot be fully contained in a cell.
281 if (EdgeIntersectsCellBoundary(a, b, children[i])) {
282 cover->push_back(children[i].id());
283 }
284 }
285 }
286 }
287 }
288 VLOG(1) << "Num cells traversed: " << num_cells;
289}
290
291// Appends to "candidate_crossings" all edge references which may cross the
292// given edge. This is done by covering the edge and then finding all
293// references of edges whose coverings overlap this covering. Parent cells
294// are checked level by level. Child cells are checked all at once by taking
295// advantage of the natural ordering of S2CellIds.
296void S2EdgeIndex::FindCandidateCrossings(
297 S2Point const& a, S2Point const& b,
298 vector<int>* candidate_crossings) const {
299 DCHECK(index_computed_);
300 vector<S2CellId> cover;
301 GetCovering(a, b, false, &cover);
302 GetEdgesInParentCells(cover, mapping_, minimum_s2_level_used_,
303 candidate_crossings);
304
305 // TODO(user): An important optimization for long query
306 // edges (Contains queries): keep a bounding cap and clip the query
307 // edge to the cap before starting the descent.
308 GetEdgesInChildrenCells(a, b, &cover, mapping_, candidate_crossings);
309
310 // Remove duplicates: This is necessary because edge references are
311 // inserted into the map once for each covering cell. (Testing shows
312 // this to be at least as fast as using a set.)
313 sort(candidate_crossings->begin(), candidate_crossings->end());
314 candidate_crossings->erase(
315 unique(candidate_crossings->begin(), candidate_crossings->end()),
316 candidate_crossings->end());
317}
318
319
320// Returns the smallest cell containing all four points, or Sentinel
321// if they are not all on the same face.
322// The points don't need to be normalized.
323static S2CellId ContainingCell(S2Point const& pa, S2Point const& pb,
324 S2Point const& pc, S2Point const& pd) {
325 S2CellId a = S2CellId::FromPoint(pa);
326 S2CellId b = S2CellId::FromPoint(pb);
327 S2CellId c = S2CellId::FromPoint(pc);
328 S2CellId d = S2CellId::FromPoint(pd);
329
330 if (a.face() != b.face() || a.face() != c.face() || a.face() != d.face()) {
331 return S2CellId::Sentinel();
332 }
333
334 while (a != b || a != c || a != d) {
335 a = a.parent();
336 b = b.parent();
337 c = c.parent();
338 d = d.parent();
339 }
340 return a;
341}
342
343// Returns the smallest cell containing both points, or Sentinel
344// if they are not all on the same face.
345// The points don't need to be normalized.
346static S2CellId ContainingCell(S2Point const& pa, S2Point const& pb) {
347 S2CellId a = S2CellId::FromPoint(pa);
348 S2CellId b = S2CellId::FromPoint(pb);
349
350 if (a.face() != b.face()) return S2CellId::Sentinel();
351
352 while (a != b) {
353 a = a.parent();
354 b = b.parent();
355 }
356 return a;
357}
358
359int S2EdgeIndex::GetCovering(
360 S2Point const& a, S2Point const& b,
361 bool thicken_edge,
362 vector<S2CellId>* edge_covering) const {
363 edge_covering->clear();
364
365 // Thicken the edge in all directions by roughly 1% of the edge length when
366 // thicken_edge is true.
367 static const double kThickening = 0.01;
368
369 // Selects the ideal s2 level at which to cover the edge, this will be the
370 // level whose S2 cells have a width roughly commensurate to the length of
371 // the edge. We multiply the edge length by 2*kThickening to guarantee the
372 // thickening is honored (it's not a big deal if we honor it when we don't
373 // request it) when doing the covering-by-cap trick.
374 const double edge_length = a.Angle(b);
375 const int ideal_level = S2::kMinWidth.GetMaxLevel(
376 edge_length * (1 + 2 * kThickening));
377
378 S2CellId containing_cell;
379 if (!thicken_edge) {
380 containing_cell = ContainingCell(a, b);
381 } else {
382 if (ideal_level == S2CellId::kMaxLevel) {
383 // If the edge is tiny, instabilities are more likely, so we
384 // want to limit the number of operations.
385 // We pretend we are in a cell much larger so as to trigger the
386 // 'needs covering' case, so we won't try to thicken the edge.
387 containing_cell = S2CellId(0xFFF0).parent(3);
388 } else {
389 S2Point pq = (b - a) * kThickening;
390 S2Point ortho = pq.CrossProd(a).Normalize() *
391 edge_length * kThickening;
392 S2Point p = a - pq;
393 S2Point q = b + pq;
394 // If p and q were antipodal, the edge wouldn't be lengthened,
395 // and it could even flip! This is not a problem because
396 // ideal_level != 0 here. The farther p and q can be is roughly
397 // a quarter Earth away from each other, so we remain
398 // Theta(kThickening).
399 containing_cell = ContainingCell(p - ortho, p + ortho,
400 q - ortho, q + ortho);
401 }
402 }
403
404 // Best case: edge is fully contained in a cell that's not too big.
405 if (containing_cell != S2CellId::Sentinel() &&
406 containing_cell.level() >= ideal_level - 2) {
407 edge_covering->push_back(containing_cell);
408 return containing_cell.level();
409 }
410
411 if (ideal_level == 0) {
412 // Edge is very long, maybe even longer than a face width, so the
413 // trick below doesn't work. For now, we will add the whole S2 sphere.
414 // TODO(user): Do something a tad smarter (and beware of the
415 // antipodal case).
416 for (S2CellId cellid = S2CellId::Begin(0); cellid != S2CellId::End(0);
417 cellid = cellid.next()) {
418 edge_covering->push_back(cellid);
419 }
420 return 0;
421 }
422 // TODO(user): Check trick below works even when vertex is at interface
423 // between three faces.
424
425 // Use trick as in S2PolygonBuilder::PointIndex::FindNearbyPoint:
426 // Cover the edge by a cap centered at the edge midpoint, then cover
427 // the cap by four big-enough cells around the cell vertex closest to the
428 // cap center.
429 S2Point middle = ((a + b) / 2).Normalize();
430 int actual_level = min(ideal_level, S2CellId::kMaxLevel-1);
431 S2CellId::FromPoint(middle).AppendVertexNeighbors(
432 actual_level, edge_covering);
433 return actual_level;
434}
435
436void S2EdgeIndex::Iterator::GetCandidates(S2Point const& a, S2Point const& b) {
437 edge_index_->PredictAdditionalCalls(1);
438 is_brute_force_ = !edge_index_->IsIndexComputed();
439 if (is_brute_force_) {
440 edge_index_->IncrementQueryCount();
441 current_index_ = 0;
442 num_edges_ = edge_index_->num_edges();
443 } else {
444 candidates_.clear();
445 edge_index_->FindCandidateCrossings(a, b, &candidates_);
446 current_index_in_candidates_ = 0;
447 if (!candidates_.empty()) {
448 current_index_ = candidates_[0];
449 }
450 }
451}
452