Oregami
Repositories/oxedyne/fe2o3

oxedyne/fe2o3/fe2o3_geom/src/cell.rs

24.5 KiB, 8 runs

created by r1870400018:44200, which is this file's identity for as long as the history lasts, whatever it is later renamed to

download · who wrote it · its history

1//! A global cell-index grid: a cube-sphere quadtree with great-circle cell edges.
2//!
3//! The sphere is wrapped in a cube. A direction picks one of the cube's six faces by its
4//! dominant axis; the other two coordinates, divided by the dominant one, give a gnomonic
5//! `(u, v)` on that face. A tangent warp `s = (2/pi)*atan(u) + 1/2` (and the matching one
6//! for `v`) evens the cell areas out, so that every cell at a given level has an area within
7//! a factor of sqrt(2) of every other. Each face is then a `2^level x 2^level` grid indexed
8//! by `(i, j)`, and a cell is a spherical quadrilateral whose four edges are great circles.
9//!
10//! Because the warp and the grid are applied per face in closed form, the scheme has two
11//! properties Uber's H3 cannot offer: a cell's parent is an exact prefix of its index (a
12//! bit mask, never a lookup), and containment is a closed-form four-half-space test. The
13//! poles are ordinary interior points of the +z and -z faces, so there is no polar special
14//! case.
15//!
16//! An id packs into one `u64`:
17//!
18//! ```text
19//! [ 4 bits scheme | 3 bits face | 5 bits level | 52 bits interleaved Morton(i, j) ]
20//! ```
21//!
22//! The Morton field is left-justified, so masking its low bits yields an ancestor's field
23//! directly. The scheme nibble is `0x1` ("cube-tan-v1"); other values are reserved for
24//! future warps or face layouts. The string form is sixteen lowercase hexadecimal digits.
25//!
26//! # Provenance
27//!
28//! The face and `(u, v)` conventions are those of Google's S2 geometry library, chosen
29//! because adjacent faces share edges continuously under them; the tangent warp is S2's as
30//! well. The quadtree indexing, the `u64` layout, the exact-parent property and the
31//! re-indexing neighbour walk are this library's own.
32
33use crate::proj::EARTH_RADIUS_M;
34
35use oxedyne_fe2o3_core::prelude::*;
36
37use std::{
38 f64::consts::{FRAC_2_PI, FRAC_PI_2, PI},
39 fmt,
40 str::FromStr,
41};
42
43/// The deepest level; a level-26 cell is roughly 14 cm across at the equator.
44pub const MAX_LEVEL: u8 = 26;
45
46// Bit layout of the u64 id.
47const SCHEME_SHIFT: u64 = 60;
48const FACE_SHIFT: u64 = 57;
49const LEVEL_SHIFT: u64 = 52;
50const FACE_MASK: u64 = 0x7 << FACE_SHIFT;
51const LEVEL_MASK: u64 = 0x1F << LEVEL_SHIFT;
52const FIELD_MASK: u64 = (1u64 << LEVEL_SHIFT) - 1; // low 52 bits
53const SCHEME_CUBE_TAN_V1: u64 = 0x1;
54
55// ---------------------------------------------------------------------------------------------
56// Cube faces
57// ---------------------------------------------------------------------------------------------
58
59/// One of the six faces of the cube enclosing the sphere.
60///
61/// The discriminants are the on-wire face indices, so `PosX` is 0 and `NegZ` is 5.
62#[derive(Clone, Copy, Debug, PartialEq, Eq)]
63enum Face {
64 PosX, // 0, +x dominant
65 PosY, // 1, +y dominant
66 PosZ, // 2, +z dominant (north pole is its interior)
67 NegX, // 3
68 NegY, // 4
69 NegZ, // 5, -z dominant (south pole is its interior)
70}
71
72impl Face {
73 fn index(self) -> u8 {
74 match self {
75 Self::PosX => 0,
76 Self::PosY => 1,
77 Self::PosZ => 2,
78 Self::NegX => 3,
79 Self::NegY => 4,
80 Self::NegZ => 5,
81 }
82 }
83
84 fn from_index(n: u8) -> Outcome<Face> {
85 match n {
86 0 => Ok(Self::PosX),
87 1 => Ok(Self::PosY),
88 2 => Ok(Self::PosZ),
89 3 => Ok(Self::NegX),
90 4 => Ok(Self::NegY),
91 5 => Ok(Self::NegZ),
92 _ => Err(err!("Cell face index {} is out of range 0..=5.", n; Invalid, Input)),
93 }
94 }
95
96 /// The face whose dominant axis best matches the direction `v`.
97 ///
98 /// On a tie between axes -- a point exactly on a cube edge or vertex -- the lowest axis
99 /// index wins (x before y before z). This is the single tie rule the whole grid shares,
100 /// so that boundary points have exactly one owning cell.
101 fn of_vec(v: &[f64; 3]) -> Face {
102 let (ax, ay, az) = (v[0].abs(), v[1].abs(), v[2].abs());
103 if ax >= ay && ax >= az {
104 if v[0] >= 0.0 { Self::PosX } else { Self::NegX }
105 } else if ay >= az {
106 if v[1] >= 0.0 { Self::PosY } else { Self::NegY }
107 } else if v[2] >= 0.0 {
108 Self::PosZ
109 } else {
110 Self::NegZ
111 }
112 }
113
114 /// Turns a face-local `(u, v)` into an (unnormalised) direction.
115 ///
116 /// The inverse of [`Face::vec_to_uv`] up to length. These formulae are S2's, so that
117 /// the shared edge between two faces is parameterised identically from both sides.
118 fn uv_to_vec(self, u: f64, v: f64) -> [f64; 3] {
119 match self {
120 Self::PosX => [ 1.0, u, v],
121 Self::PosY => [ -u, 1.0, v],
122 Self::PosZ => [ -u, -v, 1.0],
123 Self::NegX => [-1.0, -v, -u],
124 Self::NegY => [ v, -1.0, -u],
125 Self::NegZ => [ v, u, -1.0],
126 }
127 }
128
129 /// Turns a direction known to lie on this face into its `(u, v)`.
130 ///
131 /// The dominant component is the divisor, so the result lies in `[-1, 1]^2`.
132 fn vec_to_uv(self, p: &[f64; 3]) -> (f64, f64) {
133 let (x, y, z) = (p[0], p[1], p[2]);
134 match self {
135 Self::PosX => ( y / x, z / x),
136 Self::PosY => (-x / y, z / y),
137 Self::PosZ => (-x / z, -y / z),
138 Self::NegX => ( z / x, y / x),
139 Self::NegY => ( z / y, -x / y),
140 Self::NegZ => (-y / z, -x / z),
141 }
142 }
143}
144
145// ---------------------------------------------------------------------------------------------
146// Scalar geometry
147// ---------------------------------------------------------------------------------------------
148
149/// The tangent warp: face coordinate `u` in `[-1, 1]` to grid coordinate `s` in `[0, 1]`.
150fn uv_to_st(u: f64) -> f64 { FRAC_2_PI * u.atan() + 0.5 }
151
152/// The inverse tangent warp: grid coordinate `s` to face coordinate `u`.
153fn st_to_uv(s: f64) -> f64 { ((s - 0.5) * FRAC_PI_2).tan() }
154
155fn latlon_to_vec(lat_deg: f64, lon_deg: f64) -> [f64; 3] {
156 let (la, lo) = (lat_deg.to_radians(), lon_deg.to_radians());
157 let cl = la.cos();
158 [cl * lo.cos(), cl * lo.sin(), la.sin()]
159}
160
161fn vec_to_latlon(v: &[f64; 3]) -> (f64, f64) {
162 let len = (v[0]*v[0] + v[1]*v[1] + v[2]*v[2]).sqrt();
163 let z = (v[2] / len).clamp(-1.0, 1.0);
164 (z.asin().to_degrees(), v[1].atan2(v[0]).to_degrees())
165}
166
167fn normalise(v: &[f64; 3]) -> [f64; 3] {
168 let len = (v[0]*v[0] + v[1]*v[1] + v[2]*v[2]).sqrt();
169 [v[0]/len, v[1]/len, v[2]/len]
170}
171
172fn dot(a: &[f64; 3], b: &[f64; 3]) -> f64 { a[0]*b[0] + a[1]*b[1] + a[2]*b[2] }
173
174fn cross(a: &[f64; 3], b: &[f64; 3]) -> [f64; 3] {
175 [a[1]*b[2] - a[2]*b[1], a[2]*b[0] - a[0]*b[2], a[0]*b[1] - a[1]*b[0]]
176}
177
178// ---------------------------------------------------------------------------------------------
179// Morton interleave
180// ---------------------------------------------------------------------------------------------
181
182/// Interleaves two indices, `i` in the even bit positions and `j` in the odd.
183///
184/// Each input carries at most [`MAX_LEVEL`] bits, so the result carries at most 52.
185fn interleave(i: u32, j: u32) -> u64 {
186 let mut m = 0u64;
187 let mut b = 0u32;
188 while b < MAX_LEVEL as u32 {
189 m |= (((i >> b) & 1) as u64) << (2 * b);
190 m |= (((j >> b) & 1) as u64) << (2 * b + 1);
191 b += 1;
192 }
193 m
194}
195
196/// The inverse of [`interleave`].
197fn deinterleave(m: u64) -> (u32, u32) {
198 let mut i = 0u32;
199 let mut j = 0u32;
200 let mut b = 0u32;
201 while b < MAX_LEVEL as u32 {
202 i |= (((m >> (2 * b)) & 1) as u32) << b;
203 j |= (((m >> (2 * b + 1)) & 1) as u32) << b;
204 b += 1;
205 }
206 (i, j)
207}
208
209// ---------------------------------------------------------------------------------------------
210// Corner
211// ---------------------------------------------------------------------------------------------
212
213/// A cell corner as both a unit direction and its geodetic coordinates.
214#[derive(Clone, Copy, Debug)]
215pub struct Corner {
216 pub vec: [f64; 3], // unit direction from the sphere's centre
217 pub lat: f64, // degrees, positive north
218 pub lon: f64, // degrees, positive east
219}
220
221// ---------------------------------------------------------------------------------------------
222// Cell
223// ---------------------------------------------------------------------------------------------
224
225/// A single cell of the cube-sphere quadtree grid, addressed by one packed `u64`.
226#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash, PartialOrd, Ord)]
227pub struct Cell(u64);
228
229impl Cell {
230 /// The cell of the given level containing the point at `lat_deg`, `lon_deg`.
231 ///
232 /// Latitude is degrees north, longitude degrees east; neither need be pre-wrapped. A
233 /// point on a cell boundary is resolved by the shared tie rule (see [`Face::of_vec`] and
234 /// the half-open grid convention), so every point has exactly one owning cell.
235 pub fn at(lat_deg: f64, lon_deg: f64, level: u8) -> Outcome<Cell> {
236 if level > MAX_LEVEL {
237 return Err(err!("Cell level {} exceeds the maximum {}.", level, MAX_LEVEL;
238 Invalid, Input, Range));
239 }
240 let v = latlon_to_vec(lat_deg, lon_deg);
241 let face = Face::of_vec(&v);
242 let (u, w) = face.vec_to_uv(&v);
243 let (s, t) = (uv_to_st(u), uv_to_st(w));
244 let n = 1u32 << level;
245 let (i, j) = (clamp_index(s, n), clamp_index(t, n));
246 Cell::from_face_ij(face.index(), level, i, j)
247 }
248
249 /// Builds a cell directly from its face, level and grid indices.
250 pub fn from_face_ij(face: u8, level: u8, i: u32, j: u32) -> Outcome<Cell> {
251 let f = res!(Face::from_index(face));
252 if level > MAX_LEVEL {
253 return Err(err!("Cell level {} exceeds the maximum {}.", level, MAX_LEVEL;
254 Invalid, Input, Range));
255 }
256 let n = 1u32 << level;
257 if i >= n || j >= n {
258 return Err(err!("Cell index ({}, {}) is out of range for level {} (0..{}).",
259 i, j, level, n; Invalid, Input, Range));
260 }
261 let field = interleave(i, j) << (LEVEL_SHIFT - 2 * level as u64);
262 let bits = (SCHEME_CUBE_TAN_V1 << SCHEME_SHIFT)
263 | ((f.index() as u64) << FACE_SHIFT)
264 | ((level as u64) << LEVEL_SHIFT)
265 | (field & FIELD_MASK);
266 Ok(Cell(bits))
267 }
268
269 /// Validates and wraps a raw `u64`, rejecting every malformed field.
270 pub fn from_bits(bits: u64) -> Outcome<Cell> {
271 let scheme = bits >> SCHEME_SHIFT;
272 if scheme != SCHEME_CUBE_TAN_V1 {
273 return Err(err!("Cell scheme nibble {:#x} is not the cube-tan-v1 scheme {:#x}.",
274 scheme, SCHEME_CUBE_TAN_V1; Invalid, Input));
275 }
276 let face = ((bits & FACE_MASK) >> FACE_SHIFT) as u8;
277 res!(Face::from_index(face));
278 let level = ((bits & LEVEL_MASK) >> LEVEL_SHIFT) as u8;
279 if level > MAX_LEVEL {
280 return Err(err!("Cell level {} exceeds the maximum {}.", level, MAX_LEVEL;
281 Invalid, Input, Range));
282 }
283 let unused = (1u64 << (LEVEL_SHIFT - 2 * level as u64)) - 1; // low bits that must be zero
284 if bits & unused != 0 {
285 return Err(err!("Cell id has non-zero unused low bits below level {}.", level;
286 Invalid, Input));
287 }
288 Ok(Cell(bits))
289 }
290
291 /// The raw packed id.
292 pub fn bits(&self) -> u64 { self.0 }
293
294 /// The level, `0..=MAX_LEVEL`.
295 pub fn level(&self) -> u8 { ((self.0 & LEVEL_MASK) >> LEVEL_SHIFT) as u8 }
296
297 /// The face index, `0..=5`.
298 pub fn face(&self) -> u8 { ((self.0 & FACE_MASK) >> FACE_SHIFT) as u8 }
299
300 fn face_enum(&self) -> Face {
301 // The stored face is validated on every construction path, so it is always in range.
302 match Face::from_index(self.face()) {
303 Ok(f) => f,
304 Err(_) => Face::PosX,
305 }
306 }
307
308 /// The grid indices `(i, j)` at the cell's own level.
309 pub fn coords(&self) -> (u32, u32) {
310 let level = self.level();
311 let morton = (self.0 & FIELD_MASK) >> (LEVEL_SHIFT - 2 * level as u64);
312 deinterleave(morton)
313 }
314
315 /// The ancestor at level `level`, obtained by masking the Morton field.
316 ///
317 /// Because the field is left-justified, this is an exact prefix operation -- no
318 /// re-projection and no rounding. The target level must not exceed this cell's.
319 pub fn parent(&self, level: u8) -> Outcome<Cell> {
320 let own = self.level();
321 if level > own {
322 return Err(err!("Cannot take a level-{} parent of a level-{} cell.", level, own;
323 Invalid, Input, Range));
324 }
325 // Clear the Morton bits below the ancestor level, then restamp the level field.
326 let clear = if level == 0 { FIELD_MASK } else { (1u64 << (LEVEL_SHIFT - 2 * level as u64)) - 1 };
327 let bits = (self.0 & !LEVEL_MASK & !clear) | ((level as u64) << LEVEL_SHIFT);
328 Ok(Cell(bits))
329 }
330
331 /// The four child cells at the next level, in Morton order.
332 pub fn children(&self) -> Outcome<[Cell; 4]> {
333 let level = self.level();
334 if level >= MAX_LEVEL {
335 return Err(err!("A level-{} cell has no children.", level; Invalid, Input, Range));
336 }
337 let (i, j) = self.coords();
338 let face = self.face();
339 Ok([
340 res!(Cell::from_face_ij(face, level + 1, 2*i, 2*j)),
341 res!(Cell::from_face_ij(face, level + 1, 2*i + 1, 2*j)),
342 res!(Cell::from_face_ij(face, level + 1, 2*i, 2*j + 1)),
343 res!(Cell::from_face_ij(face, level + 1, 2*i + 1, 2*j + 1)),
344 ])
345 }
346
347 /// The centre of the cell as a unit direction.
348 pub fn centre_vec(&self) -> [f64; 3] {
349 let level = self.level();
350 let (i, j) = self.coords();
351 let n = (1u32 << level) as f64;
352 let s = (i as f64 + 0.5) / n;
353 let t = (j as f64 + 0.5) / n;
354 normalise(&self.face_enum().uv_to_vec(st_to_uv(s), st_to_uv(t)))
355 }
356
357 /// The centre of the cell as `(lat, lon)` in degrees.
358 pub fn centre(&self) -> (f64, f64) { vec_to_latlon(&self.centre_vec()) }
359
360 /// The four corners, counter-clockwise in the face's `(s, t)` frame.
361 pub fn corners(&self) -> [Corner; 4] {
362 let level = self.level();
363 let (i, j) = self.coords();
364 let n = (1u32 << level) as f64;
365 let face = self.face_enum();
366 let st = [
367 (i as f64 / n, j as f64 / n),
368 ((i as f64 + 1.0) / n, j as f64 / n),
369 ((i as f64 + 1.0) / n, (j as f64 + 1.0) / n),
370 (i as f64 / n, (j as f64 + 1.0) / n),
371 ];
372 let mut out = [Corner { vec: [0.0; 3], lat: 0.0, lon: 0.0 }; 4];
373 for (k, (s, t)) in st.iter().enumerate() {
374 let vec = normalise(&face.uv_to_vec(st_to_uv(*s), st_to_uv(*t)));
375 let (lat, lon) = vec_to_latlon(&vec);
376 out[k] = Corner { vec, lat, lon };
377 }
378 out
379 }
380
381 /// Does the cell contain the point at `lat_deg`, `lon_deg`?
382 ///
383 /// A cell's four edges are great circles, so containment is exact: the point belongs to
384 /// exactly the cell that [`Cell::at`] returns for it at this level. Both share the one
385 /// tie rule, so a point on a shared boundary is contained by exactly one cell.
386 pub fn contains(&self, lat_deg: f64, lon_deg: f64) -> Outcome<bool> {
387 let owner = res!(Cell::at(lat_deg, lon_deg, self.level()));
388 Ok(owner == *self)
389 }
390
391 /// The cell's area in steradians, by spherical excess of its two triangles.
392 pub fn area(&self) -> f64 {
393 let c = self.corners();
394 tri_area(&c[0].vec, &c[1].vec, &c[2].vec) + tri_area(&c[0].vec, &c[2].vec, &c[3].vec)
395 }
396
397 /// The eight (or, at a cube vertex, seven) cells sharing an edge or corner with this one.
398 ///
399 /// A step within the face is exact index arithmetic; a step that leaves the face is
400 /// re-indexed through the sphere, which resolves the face change and its orientation flip
401 /// automatically. At the 24 cells per level that touch a cube vertex, the diagonal step
402 /// across the vertex lands on an already-listed cell, so those cells have seven neighbours.
403 pub fn neighbours(&self) -> Outcome<Vec<Cell>> {
404 const DIRS: [(i32, i32); 8] = [
405 (1, 0), (-1, 0), (0, 1), (0, -1), // edge-adjacent
406 (1, 1), (1, -1), (-1, 1), (-1, -1), // corner-adjacent
407 ];
408 let level = self.level();
409 let (i, j) = self.coords();
410 let face = self.face_enum();
411 let n = 1i64 << level;
412 let mut out: Vec<Cell> = Vec::with_capacity(8);
413 for (di, dj) in DIRS.iter() {
414 let ni = i as i64 + *di as i64;
415 let nj = j as i64 + *dj as i64;
416 let cell = if ni >= 0 && ni < n && nj >= 0 && nj < n {
417 res!(Cell::from_face_ij(face.index(), level, ni as u32, nj as u32))
418 } else {
419 // Leave the face: a representative point a quarter-cell past the crossed edge,
420 // centred on the in-range axis. A quarter cell never reaches the far edge
421 // (which would be a singular tangent), so no infinity ever arises.
422 let nf = n as f64;
423 let s = axis_rep(ni, n, nf);
424 let t = axis_rep(nj, n, nf);
425 let vec = face.uv_to_vec(st_to_uv(s), st_to_uv(t));
426 let (lat, lon) = vec_to_latlon(&vec);
427 res!(Cell::at(lat, lon, level))
428 };
429 if cell != *self && !out.contains(&cell) {
430 out.push(cell);
431 }
432 }
433 Ok(out)
434 }
435
436 /// The cells at exactly graph-distance `k` from this one, over the neighbour relation.
437 ///
438 /// `ring(0)` is the cell itself; `ring(1)` is [`Cell::neighbours`].
439 pub fn ring(&self, k: u32) -> Outcome<Vec<Cell>> {
440 let mut seen: std::collections::HashSet<Cell> = std::collections::HashSet::new();
441 let mut frontier = vec![*self];
442 seen.insert(*self);
443 for _ in 0..k {
444 let mut next = Vec::new();
445 for cell in &frontier {
446 for nb in res!(cell.neighbours()) {
447 if seen.insert(nb) {
448 next.push(nb);
449 }
450 }
451 }
452 frontier = next;
453 if frontier.is_empty() {
454 break;
455 }
456 }
457 Ok(frontier)
458 }
459
460 /// The cell's boundary as unit vectors, each edge divided into `segs` pieces along its
461 /// great circle, counter-clockwise from the first corner and not closed.
462 ///
463 /// The edges are great circles, so on the globe a straight line between corners is
464 /// already the edge's chord; on a flat map a large cell's edge is visibly curved. About
465 /// eight pieces serve below level 8 and one above it, where an edge is a few kilometres
466 /// and no projection bends it by a pixel.
467 pub fn outline(&self, segs: u32) -> Vec<[f64; 3]> {
468 let segs = segs.max(1) as usize;
469 let c = self.corners();
470 let mut out = Vec::with_capacity(4 * segs);
471 for k in 0..4 {
472 let a = c[k].vec;
473 let b = c[(k + 1) % 4].vec;
474 for s in 0..segs {
475 // A normalised straight blend of two points stays on their great circle.
476 let t = s as f64 / segs as f64;
477 out.push(normalise(&[
478 a[0] + t * (b[0] - a[0]),
479 a[1] + t * (b[1] - a[1]),
480 a[2] + t * (b[2] - a[2]),
481 ]));
482 }
483 }
484 out
485 }
486
487 /// The smallest cap about the cell's centre that holds the whole cell, as the centre and
488 /// the cap's angular radius in radians.
489 ///
490 /// The furthest point of a quadrilateral with great-circle edges from a point inside it is
491 /// a corner, so the radius is the furthest corner's.
492 pub fn bounding_cap(&self) -> ([f64; 3], f64) {
493 let centre = self.centre_vec();
494 let mut far: f64 = 0.0;
495 for corner in self.corners().iter() {
496 far = far.max(dot(&centre, &corner.vec).clamp(-1.0, 1.0).acos());
497 }
498 (centre, far)
499 }
500}
501
502/// The cells of a level that a spherical cap may touch: every cell whose bounding cap meets
503/// the query cap.
504///
505/// A breadth-first walk over [`Cell::neighbours`] from the cell holding the cap's centre, so
506/// the cost is the size of the answer and not of the level, and the answer comes nearest
507/// first. `max` bounds that answer: a cap that would cover more cells is refused rather than
508/// walked, because a caller that asked for the cells on a screen and got millions has asked
509/// at the wrong level.
510///
511/// # Arguments
512/// * `centre` - The cap's centre as a vector; its length does not matter.
513/// * `radius_rad` - The cap's angular radius, in radians of arc.
514pub fn cover_cap(
515 centre: [f64; 3],
516 radius_rad: f64,
517 level: u8,
518 max: usize,
519)
520 -> Outcome<Vec<Cell>>
521{
522 let len = dot(&centre, &centre).sqrt();
523 if !(len > 0.0 && len.is_finite()) {
524 return Err(err!("A cap centred on {:?} has no direction.", centre; Invalid, Input));
525 }
526 if !(radius_rad >= 0.0 && radius_rad.is_finite()) {
527 return Err(err!("A cap of radius {} radians is not a cap.", radius_rad;
528 Invalid, Input, Range));
529 }
530 if max == 0 {
531 return Err(err!("A cover of at most no cells cannot hold the cap's own centre.";
532 Invalid, Input, Range));
533 }
534 let c = normalise(&centre);
535 let (lat, lon) = vec_to_latlon(&c);
536 let first = res!(Cell::at(lat, lon, level));
537 let meets = |cell: &Cell| -> bool {
538 let (v, r) = cell.bounding_cap();
539 dot(&v, &c).clamp(-1.0, 1.0).acos() <= r + radius_rad + 1.0e-12
540 };
541 let mut seen: std::collections::HashSet<Cell> = std::collections::HashSet::new();
542 let mut out: Vec<Cell> = Vec::new();
543 let mut queue: std::collections::VecDeque<Cell> = std::collections::VecDeque::new();
544 seen.insert(first);
545 out.push(first);
546 queue.push_back(first);
547 while let Some(cell) = queue.pop_front() {
548 for nb in res!(cell.neighbours()) {
549 if seen.insert(nb) && meets(&nb) {
550 if out.len() >= max {
551 return Err(err!(
552 "A cap of {} radians covers more than {} level-{} cells.",
553 radius_rad, max, level;
554 Excessive, Size));
555 }
556 out.push(nb);
557 queue.push_back(nb);
558 }
559 }
560 }
561 Ok(out)
562}
563
564/// The mean side of a level's cells in metres, on a sphere of [`EARTH_RADIUS_M`]: the square
565/// root of the mean cell area, about 9,220 km at level 0 and halving with each level.
566pub fn mean_side_m(level: u8) -> f64 {
567 EARTH_RADIUS_M * (2.0 * PI / 3.0).sqrt() / (1u64 << level.min(MAX_LEVEL)) as f64
568}
569
570/// The finest level whose mean cell side spans at least `min_px` pixels at a scale of
571/// `m_per_px` ground metres per pixel, which is the level a map draws its grid at.
572///
573/// Level 0 when even that is smaller, and [`MAX_LEVEL`] when every level is larger.
574pub fn level_for_scale(m_per_px: f64, min_px: f64) -> u8 {
575 if !(m_per_px > 0.0 && min_px > 0.0) || !(m_per_px * min_px).is_finite() {
576 return 0;
577 }
578 let want = m_per_px * min_px;
579 let mut level = 0u8;
580 while level < MAX_LEVEL && mean_side_m(level + 1) >= want {
581 level += 1;
582 }
583 level
584}
585
586/// The `(s, t)` coordinate representing a stepped neighbour along one axis.
587///
588/// In range, the cell centre; below the face, a quarter-cell before the near edge; above it,
589/// a quarter-cell past the far edge.
590fn axis_rep(idx: i64, n: i64, nf: f64) -> f64 {
591 if idx < 0 {
592 -0.25 / nf
593 } else if idx >= n {
594 1.0 + 0.25 / nf
595 } else {
596 (idx as f64 + 0.5) / nf
597 }
598}
599
600/// Floors a warped coordinate `s` in `[0, 1]` to a grid index, clamped to `0..n`.
601fn clamp_index(s: f64, n: u32) -> u32 {
602 let raw = (s * n as f64).floor();
603 if raw < 0.0 {
604 0
605 } else if raw >= n as f64 {
606 n - 1
607 } else {
608 raw as u32
609 }
610}
611
612/// Spherical-excess area of a unit-vector triangle (Van Oosterom & Strackee).
613fn tri_area(a: &[f64; 3], b: &[f64; 3], c: &[f64; 3]) -> f64 {
614 let triple = dot(a, &cross(b, c)).abs();
615 let den = 1.0 + dot(a, b) + dot(b, c) + dot(c, a);
616 2.0 * triple.atan2(den)
617}
618
619impl fmt::Display for Cell {
620 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
621 write!(f, "{:016x}", self.0)
622 }
623}
624
625impl FromStr for Cell {
626 type Err = Error<ErrTag>;
627
628 fn from_str(s: &str) -> Outcome<Cell> {
629 if s.len() != 16 {
630 return Err(err!("A cell id must be 16 hex characters, got {}.", s.len();
631 Invalid, Input));
632 }
633 let bits = res!(u64::from_str_radix(s, 16), Invalid, Input);
634 Cell::from_bits(bits)
635 }
636}