1// Copyright 2005 Google Inc. All Rights Reserved.
2
3#include "s2regioncoverer.h"
4
5#include <pthread.h>
6
7#include <algorithm>
8using std::min;
9using std::max;
10using std::swap;
11using std::reverse;
12
13#include <functional>
14using std::less;
15
16#include <hash_set>
17using __gnu_cxx::hash_set;
18
19#include <queue>
20using std::priority_queue;
21
22#include <vector>
23using std::vector;
24
25
26#include "base/logging.h"
27#include "s2.h"
28#include "s2cap.h"
29#include "s2cellunion.h"
30
31// Define storage for header file constants (the values are not needed here).
32int const S2RegionCoverer::kDefaultMaxCells;
33
34// We define our own own comparison function on QueueEntries in order to
35// make the results deterministic. Using the default less<QueueEntry>,
36// entries of equal priority would be sorted according to the memory address
37// of the candidate.
38
39struct S2RegionCoverer::CompareQueueEntries : public less<QueueEntry> {
40 bool operator()(QueueEntry const& x, QueueEntry const& y) {
41 return x.first < y.first;
42 }
43};
44
45static S2Cell face_cells[6];
46
47void Init() {
48 for (int face = 0; face < 6; ++face) {
49 face_cells[face] = S2Cell::FromFacePosLevel(face, 0, 0);
50 }
51}
52
53static pthread_once_t init_once = PTHREAD_ONCE_INIT;
54inline static void MaybeInit() {
55 pthread_once(&init_once, Init);
56}
57
58S2RegionCoverer::S2RegionCoverer() :
59 min_level_(0),
60 max_level_(S2CellId::kMaxLevel),
61 level_mod_(1),
62 max_cells_(kDefaultMaxCells),
63 region_(NULL),
64 result_(new vector<S2CellId>),
65 pq_(new CandidateQueue) {
66 // Initialize the constants
67 MaybeInit();
68}
69
70S2RegionCoverer::~S2RegionCoverer() {
71 // Need to declare explicitly because of scoped pointers.
72}
73
74void S2RegionCoverer::set_min_level(int min_level) {
75 DCHECK_GE(min_level, 0);
76 DCHECK_LE(min_level, S2CellId::kMaxLevel);
77 min_level_ = max(0, min(S2CellId::kMaxLevel, min_level));
78}
79
80void S2RegionCoverer::set_max_level(int max_level) {
81 DCHECK_GE(max_level, 0);
82 DCHECK_LE(max_level, S2CellId::kMaxLevel);
83 max_level_ = max(0, min(S2CellId::kMaxLevel, max_level));
84}
85
86void S2RegionCoverer::set_level_mod(int level_mod) {
87 DCHECK_GE(level_mod, 1);
88 DCHECK_LE(level_mod, 3);
89 level_mod_ = max(1, min(3, level_mod));
90}
91
92void S2RegionCoverer::set_max_cells(int max_cells) {
93 max_cells_ = max_cells;
94}
95
96S2RegionCoverer::Candidate* S2RegionCoverer::NewCandidate(S2Cell const& cell) {
97 if (!region_->MayIntersect(cell)) return NULL;
98
99 bool is_terminal = false;
100 size_t size = sizeof(Candidate);
101 if (cell.level() >= min_level_) {
102 if (interior_covering_) {
103 if (region_->Contains(cell)) {
104 is_terminal = true;
105 } else if (cell.level() + level_mod_ > max_level_) {
106 return NULL;
107 }
108 } else {
109 if (cell.level() + level_mod_ > max_level_ || region_->Contains(cell)) {
110 is_terminal = true;
111 }
112 }
113 }
114 if (!is_terminal) {
115 size += sizeof(Candidate*) << max_children_shift();
116 }
117 Candidate* candidate = static_cast<Candidate*>(malloc(size));
118 memset(candidate, 0, size);
119 candidate->cell = cell;
120 candidate->is_terminal = is_terminal;
121 ++candidates_created_counter_;
122 return candidate;
123}
124
125void S2RegionCoverer::DeleteCandidate(Candidate* candidate,
126 bool delete_children) {
127 if (delete_children) {
128 for (int i = 0; i < candidate->num_children; ++i)
129 DeleteCandidate(candidate->children[i], true);
130 }
131 free(candidate);
132}
133
134int S2RegionCoverer::ExpandChildren(Candidate* candidate,
135 S2Cell const& cell, int num_levels) {
136 num_levels--;
137 S2Cell child_cells[4];
138 cell.Subdivide(child_cells);
139 int num_terminals = 0;
140 for (int i = 0; i < 4; ++i) {
141 if (num_levels > 0) {
142 if (region_->MayIntersect(child_cells[i])) {
143 num_terminals += ExpandChildren(candidate, child_cells[i], num_levels);
144 }
145 continue;
146 }
147 Candidate* child = NewCandidate(child_cells[i]);
148 if (child) {
149 candidate->children[candidate->num_children++] = child;
150 if (child->is_terminal) ++num_terminals;
151 }
152 }
153 return num_terminals;
154}
155
156void S2RegionCoverer::AddCandidate(Candidate* candidate) {
157 if (candidate == NULL) return;
158
159 if (candidate->is_terminal) {
160 result_->push_back(candidate->cell.id());
161 DeleteCandidate(candidate, true);
162 return;
163 }
164 DCHECK_EQ(0, candidate->num_children);
165
166 // Expand one level at a time until we hit min_level_ to ensure that
167 // we don't skip over it.
168 int num_levels = (candidate->cell.level() < min_level_) ? 1 : level_mod_;
169 int num_terminals = ExpandChildren(candidate, candidate->cell, num_levels);
170
171 if (candidate->num_children == 0) {
172 DeleteCandidate(candidate, false);
173
174 } else if (!interior_covering_ &&
175 num_terminals == 1 << max_children_shift() &&
176 candidate->cell.level() >= min_level_) {
177 // Optimization: add the parent cell rather than all of its children.
178 // We can't do this for interior coverings, since the children just
179 // intersect the region, but may not be contained by it - we need to
180 // subdivide them further.
181 candidate->is_terminal = true;
182 AddCandidate(candidate);
183
184 } else {
185 // We negate the priority so that smaller absolute priorities are returned
186 // first. The heuristic is designed to refine the largest cells first,
187 // since those are where we have the largest potential gain. Among cells
188 // at the same level, we prefer the cells with the smallest number of
189 // intersecting children. Finally, we prefer cells that have the smallest
190 // number of children that cannot be refined any further.
191 int priority = -((((candidate->cell.level() << max_children_shift())
192 + candidate->num_children) << max_children_shift())
193 + num_terminals);
194 pq_->push(make_pair(priority, candidate));
195 VLOG(2) << "Push: " << candidate->cell.id() << " (" << priority << ") ";
196 }
197}
198
199void S2RegionCoverer::GetInitialCandidates() {
200 // Optimization: if at least 4 cells are desired (the normal case),
201 // start with a 4-cell covering of the region's bounding cap. This
202 // lets us skip quite a few levels of refinement when the region to
203 // be covered is relatively small.
204 if (max_cells() >= 4) {
205 // Find the maximum level such that the bounding cap contains at most one
206 // cell vertex at that level.
207 S2Cap cap = region_->GetCapBound();
208 int level = min(S2::kMinWidth.GetMaxLevel(2 * cap.angle().radians()),
209 min(max_level(), S2CellId::kMaxLevel - 1));
210 if (level_mod() > 1 && level > min_level()) {
211 level -= (level - min_level()) % level_mod();
212 }
213 // We don't bother trying to optimize the level == 0 case, since more than
214 // four face cells may be required.
215 if (level > 0) {
216 // Find the leaf cell containing the cap axis, and determine which
217 // subcell of the parent cell contains it.
218 vector<S2CellId> base;
219 base.reserve(4);
220 S2CellId id = S2CellId::FromPoint(cap.axis());
221 id.AppendVertexNeighbors(level, &base);
222 for (int i = 0; i < base.size(); ++i) {
223 AddCandidate(NewCandidate(S2Cell(base[i])));
224 }
225 return;
226 }
227 }
228 // Default: start with all six cube faces.
229 for (int face = 0; face < 6; ++face) {
230 AddCandidate(NewCandidate(face_cells[face]));
231 }
232}
233
234void S2RegionCoverer::GetCoveringInternal(S2Region const& region) {
235 // Strategy: Start with the 6 faces of the cube. Discard any
236 // that do not intersect the shape. Then repeatedly choose the
237 // largest cell that intersects the shape and subdivide it.
238 //
239 // result_ contains the cells that will be part of the output, while pq_
240 // contains cells that we may still subdivide further. Cells that are
241 // entirely contained within the region are immediately added to the output,
242 // while cells that do not intersect the region are immediately discarded.
243 // Therefore pq_ only contains cells that partially intersect the region.
244 // Candidates are prioritized first according to cell size (larger cells
245 // first), then by the number of intersecting children they have (fewest
246 // children first), and then by the number of fully contained children
247 // (fewest children first).
248
249 DCHECK(pq_->empty());
250 DCHECK(result_->empty());
251 region_ = &region;
252 candidates_created_counter_ = 0;
253
254 GetInitialCandidates();
255 while (!pq_->empty() &&
256 (!interior_covering_ || result_->size() < max_cells_)) {
257 Candidate* candidate = pq_->top().second;
258 pq_->pop();
259 VLOG(2) << "Pop: " << candidate->cell.id();
260 if (candidate->cell.level() < min_level_ ||
261 candidate->num_children == 1 ||
262 result_->size() + (interior_covering_ ? 0 : pq_->size()) +
263 candidate->num_children <= max_cells_) {
264 // Expand this candidate into its children.
265 for (int i = 0; i < candidate->num_children; ++i) {
266 AddCandidate(candidate->children[i]);
267 }
268 DeleteCandidate(candidate, false);
269 } else if (interior_covering_) {
270 DeleteCandidate(candidate, true);
271 } else {
272 candidate->is_terminal = true;
273 AddCandidate(candidate);
274 }
275 }
276 VLOG(2) << "Created " << result_->size() << " cells, " <<
277 candidates_created_counter_ << " candidates created, " <<
278 pq_->size() << " left";
279 while (!pq_->empty()) {
280 DeleteCandidate(pq_->top().second, true);
281 pq_->pop();
282 }
283 region_ = NULL;
284}
285
286void S2RegionCoverer::GetCovering(S2Region const& region,
287 vector<S2CellId>* covering) {
288
289 // Rather than just returning the raw list of cell ids generated by
290 // GetCoveringInternal(), we construct a cell union and then denormalize it.
291 // This has the effect of replacing four child cells with their parent
292 // whenever this does not violate the covering parameters specified
293 // (min_level, level_mod, etc). This strategy significantly reduces the
294 // number of cells returned in many cases, and it is cheap compared to
295 // computing the covering in the first place.
296
297 S2CellUnion tmp;
298 GetCellUnion(region, &tmp);
299 tmp.Denormalize(min_level(), level_mod(), covering);
300}
301
302void S2RegionCoverer::GetInteriorCovering(S2Region const& region,
303 vector<S2CellId>* interior) {
304 S2CellUnion tmp;
305 GetInteriorCellUnion(region, &tmp);
306 tmp.Denormalize(min_level(), level_mod(), interior);
307}
308
309void S2RegionCoverer::GetCellUnion(S2Region const& region,
310 S2CellUnion* covering) {
311 interior_covering_ = false;
312 GetCoveringInternal(region);
313 covering->InitSwap(result_.get());
314}
315
316void S2RegionCoverer::GetInteriorCellUnion(S2Region const& region,
317 S2CellUnion* interior) {
318 interior_covering_ = true;
319 GetCoveringInternal(region);
320 interior->InitSwap(result_.get());
321}
322
323void S2RegionCoverer::FloodFill(
324 S2Region const& region, S2CellId const& start, vector<S2CellId>* output) {
325 hash_set<S2CellId> all;
326 vector<S2CellId> frontier;
327 output->clear();
328 all.insert(start);
329 frontier.push_back(start);
330 while (!frontier.empty()) {
331 S2CellId id = frontier.back();
332 frontier.pop_back();
333 if (!region.MayIntersect(S2Cell(id))) continue;
334 output->push_back(id);
335
336 S2CellId neighbors[4];
337 id.GetEdgeNeighbors(neighbors);
338 for (int edge = 0; edge < 4; ++edge) {
339 S2CellId nbr = neighbors[edge];
340 if (all.insert(nbr).second) {
341 frontier.push_back(nbr);
342 }
343 }
344 }
345}
346
347void S2RegionCoverer::GetSimpleCovering(
348 S2Region const& region, S2Point const& start,
349 int level, vector<S2CellId>* output) {
350 return FloodFill(region, S2CellId::FromPoint(start).parent(level), output);
351}
352