Oregami
Repositories/oxedyne/fe2o3

oxedyne/fe2o3/fe2o3_geom/tests/view.rs

32.2 KiB, 2 runs

created by r1870400018:59722, 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//! The viewport, clipped ring projection, cell outlines and cap covering.
2//!
3//! The screen transforms are held to PROJ 9.7.1's numbers, the covering to a brute-force
4//! enumeration, and the ring clip to places that are sea or land as a matter of record: the
5//! "sea stays sea at every turn" check Ochre's globe was held to, moved here with the clip and
6//! extended to the flat map and to a globe turned and zoomed in.
7
8use oxedyne_fe2o3_geom::{
9 cell::{
10 self,
11 Cell,
12 },
13 planar::Pt,
14 proj::{
15 EARTH_RADIUS_M,
16 Projection,
17 RingMode,
18 ScreenPaths,
19 Viewport,
20 unit_vec,
21 vec_lat_lng,
22 },
23 world,
24};
25
26use oxedyne_fe2o3_core::prelude::*;
27
28use std::{
29 collections::HashSet,
30 f64::consts::{
31 PI,
32 TAU,
33 },
34};
35
36const WORLD: &[u8] = include_bytes!("data/world_coarse.bin");
37
38struct Rng(u64);
39
40impl Rng {
41 fn next_u64(&mut self) -> u64 {
42 self.0 = self.0.wrapping_add(0x9E37_79B9_7F4A_7C15);
43 let mut z = self.0;
44 z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
45 z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
46 z ^ (z >> 31)
47 }
48
49 fn unit(&mut self) -> f64 { (self.next_u64() >> 11) as f64 / ((1u64 << 53) as f64) }
50}
51
52fn dot(a: &[f64; 3], b: &[f64; 3]) -> f64 { a[0] * b[0] + a[1] * b[1] + a[2] * b[2] }
53
54fn angle(a: &[f64; 3], b: &[f64; 3]) -> f64 { dot(a, b).clamp(-1.0, 1.0).acos() }
55
56/// How many times the closed paths wind round a point, which is what a nonzero fill asks.
57fn winding(paths: &ScreenPaths, x: f64, y: f64) -> f64 {
58 let mut sum = 0.0;
59 for i in 0..paths.len() {
60 let (pts, closed) = match paths.path(i) {
61 Some(p) => p,
62 None => continue,
63 };
64 if !closed {
65 continue;
66 }
67 let n = pts.len() / 2;
68 for k in 0..n {
69 let (ax, ay) = (pts[2 * k] as f64 - x, pts[2 * k + 1] as f64 - y);
70 let m = (k + 1) % n;
71 let (bx, by) = (pts[2 * m] as f64 - x, pts[2 * m + 1] as f64 - y);
72 sum += (ax * by - ay * bx).atan2(ax * bx + ay * by);
73 }
74 }
75 sum / TAU
76}
77
78fn view(kind: Projection, lat_0: f64, lon_0: f64, heading: f64, m_per_px: f64, w: f64, h: f64)
79 -> Outcome<Viewport>
80{
81 Viewport::new(kind, lat_0, lon_0, heading, m_per_px, w, h)
82}
83
84fn ring(pts: &[(f64, f64)]) -> Vec<[f64; 3]> {
85 pts.iter().map(|(lat, lng)| unit_vec(*lat, *lng)).collect()
86}
87
88// ---------------------------------------------------------------------------------------------
89// The camera
90// ---------------------------------------------------------------------------------------------
91
92#[test]
93fn test_the_viewport_puts_cities_where_proj_does_00() -> Outcome<()> {
94 // Orthographic about Perth at 300 px to the Earth's radius; the unit-sphere numbers are
95 // PROJ's (`+proj=ortho +R=1 +lat_0=-31.9535 +lon_0=115.8571`), with y flipped for a screen.
96 let v = res!(view(Projection::Orthographic, -31.9535, 115.8571, 0.0,
97 EARTH_RADIUS_M / 300.0, 800.0, 700.0));
98 for (lat, lng, x, y) in [
99 (-33.8688, 151.2093, 0.480421544470, -0.114447989448),
100 ( 35.6895, 139.6917, 0.328204336296, 0.888173527162),
101 ( 3.1390, 101.6869, -0.244435840733, 0.558819302390),
102 (-90.0, 0.0, 0.0, -0.848477887693),
103 ] {
104 let p = res!(v.forward(lat, lng).ok_or_else(|| err!("{}, {} was hidden.", lat, lng; Test)));
105 let near = (p.x - (400.0 + 300.0 * x)).abs() < 1.0e-6
106 && (p.y - (350.0 - 300.0 * y)).abs() < 1.0e-6;
107 req!(near, true, "{}, {} landed at {}, {}.", lat, lng, p.x, p.y);
108 }
109 // London is behind Perth.
110 req!(v.forward(51.5074, -0.1278).is_none(), true, "London showed through the globe.");
111
112 // Web Mercator about Perth at 100 px to the Earth's radius at the centre's parallel, so the
113 // equator's radius is 100 / cos(lat_0). PROJ: `+proj=webmerc +R=1 +lon_0=115.8571`.
114 let phi = (-31.9535f64).to_radians();
115 let v = res!(view(Projection::WebMercator, -31.9535, 115.8571, 0.0,
116 EARTH_RADIUS_M * phi.cos() / 100.0, 800.0, 700.0));
117 let perth_y = -0.589076140273;
118 for (lat, lng, x, y) in [
119 ( 35.6895, 139.6917, 0.415992245896, 0.667590039565),
120 ( 51.5074, -0.1278, -2.024318387596, 1.052273175629),
121 (-33.9189, 18.4233, -1.700540612730, -0.629951578195),
122 (-31.9535, 115.8571, 0.0, perth_y),
123 ] {
124 let p = res!(v.forward(lat, lng).ok_or_else(|| err!("{}, {} was lost.", lat, lng; Test)));
125 let near = (p.x - (400.0 + 100.0 * x)).abs() < 1.0e-6
126 && (p.y - (350.0 - 100.0 * (y - perth_y))).abs() < 1.0e-6;
127 req!(near, true, "{}, {} landed at {}, {}.", lat, lng, p.x, p.y);
128 }
129 Ok(())
130}
131
132#[test]
133fn test_a_heading_turns_the_picture_01() -> Outcome<()> {
134 for kind in [Projection::Orthographic, Projection::WebMercator] {
135 // East up: a place due east is straight above the centre, and north is to the left.
136 let v = res!(view(kind, 0.0, 0.0, 90.0, 1000.0, 400.0, 400.0));
137 let e = res!(v.forward(0.0, 0.5).ok_or_else(|| err!("East was hidden."; Test)));
138 let up = e.y < 200.0 - 10.0 && (e.x - 200.0).abs() < 1.0e-6;
139 req!(up, true, "{:?}: east landed at {}, {}.", kind, e.x, e.y);
140 let n = res!(v.forward(0.5, 0.0).ok_or_else(|| err!("North was hidden."; Test)));
141 let left = n.x < 200.0 - 10.0 && (n.y - 200.0).abs() < 1.0e-3;
142 req!(left, true, "{:?}: north landed at {}, {}.", kind, n.x, n.y);
143 }
144 Ok(())
145}
146
147#[test]
148fn test_forward_and_inverse_undo_each_other_02() -> Outcome<()> {
149 let mut rng = Rng(7);
150 for kind in [Projection::Orthographic, Projection::WebMercator] {
151 for (lat_0, lon_0, heading, m) in [
152 (-31.9535, 115.8571, 37.0, 50.0),
153 (64.0, -150.0, -120.0, 2_000.0),
154 (0.0, 179.9, 0.0, 0.3),
155 (-75.0, 10.0, 200.0, 20_000.0),
156 ] {
157 let v = res!(view(kind, lat_0, lon_0, heading, m, 390.0, 700.0));
158 let wraps = kind == Projection::WebMercator && v.bounding_cap().1 >= PI;
159 for _ in 0..200 {
160 let pt = Pt::new(rng.unit() * 390.0, rng.unit() * 700.0);
161 let (lat, lng) = match v.inverse(pt) {
162 Some(p) => p,
163 None => continue,
164 };
165 let back = res!(v.forward(lat, lng).ok_or_else(|| err!(
166 "{:?}: {}, {} came from the screen and then hid.", kind, lat, lng; Test)));
167 let same = (back.x - pt.x).abs() < 1.0e-5 && (back.y - pt.y).abs() < 1.0e-5;
168 if wraps && !same {
169 // A flat map wider than the world shows a place more than once, and
170 // forward answers with the copy nearest the centre.
171 let (a, b) = res!(v.inverse(back).ok_or_else(|| err!("Lost."; Test)));
172 let twin = (a - lat).abs() < 1.0e-9 && ((b - lng + 540.0) % 360.0 - 180.0).abs() < 1.0e-9;
173 req!(twin, true, "{:?} about {}, {}: ({}, {}) came back at ({}, {}), \
174 which is {}, {}.", kind, lat_0, lon_0, pt.x, pt.y, back.x, back.y, a, b);
175 continue;
176 }
177 req!(same, true, "{:?} about {}, {}: ({}, {}) came back at ({}, {}).",
178 kind, lat_0, lon_0, pt.x, pt.y, back.x, back.y);
179 }
180 }
181 }
182 Ok(())
183}
184
185#[test]
186fn test_the_bounding_cap_holds_the_screen_03() -> Outcome<()> {
187 for kind in [Projection::Orthographic, Projection::WebMercator] {
188 for (lat_0, lon_0, heading, m) in [
189 (-33.87, 151.21, 0.0, 1.0),
190 (60.0, 10.0, 30.0, 300.0),
191 (0.0, -179.0, 0.0, 30_000.0),
192 (-80.0, 0.0, 0.0, 5_000.0),
193 (20.0, 40.0, 0.0, 200_000.0),
194 ] {
195 let v = res!(view(kind, lat_0, lon_0, heading, m, 390.0, 700.0));
196 let (c, r) = v.bounding_cap();
197 for i in 0..=20 {
198 for j in 0..=20 {
199 let pt = Pt::new(390.0 * i as f64 / 20.0, 700.0 * j as f64 / 20.0);
200 if let Some((lat, lng)) = v.inverse(pt) {
201 let a = angle(&c, &unit_vec(lat, lng));
202 req!((a <= r + 1.0e-12), true, "{:?} about {}, {} at {} m/px: a screen \
203 point is {} rad out, the cap {}.", kind, lat_0, lon_0, m, a, r);
204 }
205 }
206 }
207 }
208 }
209 Ok(())
210}
211
212// ---------------------------------------------------------------------------------------------
213// Rings
214// ---------------------------------------------------------------------------------------------
215
216#[test]
217fn test_the_flat_map_joins_a_ring_across_the_antimeridian_04() -> Outcome<()> {
218 let fiji = ring(&[(-17.0, 179.5), (-17.0, -179.5), (-16.0, -179.5), (-16.0, 179.5)]);
219 let v = res!(view(Projection::WebMercator, -16.5, 180.0, 0.0, 400.0, 800.0, 800.0));
220 let mut out = ScreenPaths::new();
221 res!(v.project_rings(&[fiji], RingMode::Fill, 0.0, 0.25, &mut out));
222 req!(out.len(), 1, "The island came out in {} pieces.", out.len());
223 let (pts, closed) = res!(out.path(0).ok_or_else(|| err!("No path."; Test)));
224 req!(closed, true);
225 let xs: Vec<f64> = pts.chunks(2).map(|p| p[0] as f64).collect();
226 let span = xs.iter().cloned().fold(f64::MIN, f64::max) - xs.iter().cloned().fold(f64::MAX, f64::min);
227 let narrow = span < 400.0;
228 req!(narrow, true, "An island two degrees wide spans {} px.", span);
229 let inside = winding(&out, 400.0, 400.0).abs() > 0.5;
230 req!(inside, true, "The island's middle is not filled.");
231 Ok(())
232}
233
234#[test]
235fn test_the_flat_map_closes_a_ring_round_the_pole_05() -> Outcome<()> {
236 let east: Vec<(f64, f64)> = (0..36).map(|k| (-70.0, -180.0 + 10.0 * k as f64)).collect();
237 let west: Vec<(f64, f64)> = east.iter().rev().cloned().collect();
238 // The whole map, square, 800 px on a side.
239 let v = res!(view(Projection::WebMercator, 0.0, 0.0, 0.0, TAU * EARTH_RADIUS_M / 800.0,
240 800.0, 800.0));
241 for pts in [east, west] {
242 let mut out = ScreenPaths::new();
243 res!(v.project_rings(&[ring(&pts)], RingMode::Fill, 0.0, 0.25, &mut out));
244 for (lat, lng, filled) in [(-80.0, 45.0, true), (-75.0, -170.0, true), (-60.0, 0.0, false),
245 (0.0, 90.0, false), (70.0, 0.0, false)]
246 {
247 let p = res!(v.forward(lat, lng).ok_or_else(|| err!("Lost."; Test)));
248 let got = winding(&out, p.x, p.y).abs() > 0.5;
249 req!(got, filled, "{}, {} came out {}.", lat, lng, if got { "filled" } else { "empty" });
250 }
251 }
252 // A cell with the north pole in it, and one with the pole as a corner.
253 let face = res!(Cell::from_face_ij(2, 0, 0, 0)).outline(8);
254 let corner = res!(Cell::from_face_ij(2, 1, 0, 0)).outline(8);
255 let mut out = ScreenPaths::new();
256 res!(v.project_rings(&[face], RingMode::Fill, 0.0, 0.25, &mut out));
257 for (lat, lng, filled) in [(80.0, 45.0, true), (40.0, 45.0, true), (40.0, 0.0, false),
258 (0.0, 0.0, false)]
259 {
260 let p = res!(v.forward(lat, lng).ok_or_else(|| err!("Lost."; Test)));
261 let got = winding(&out, p.x, p.y).abs() > 0.5;
262 req!(got, filled, "Face 2: {}, {} came out {}.", lat, lng, if got { "filled" } else { "empty" });
263 }
264 let mut out = ScreenPaths::new();
265 res!(v.project_rings(&[corner], RingMode::Fill, 0.0, 0.25, &mut out));
266 for (lat, lng, filled) in [(80.0, 45.0, true), (80.0, 135.0, false), (80.0, -45.0, false)] {
267 let p = res!(v.forward(lat, lng).ok_or_else(|| err!("Lost."; Test)));
268 let got = winding(&out, p.x, p.y).abs() > 0.5;
269 req!(got, filled, "Pole corner: {}, {} came out {}.", lat, lng,
270 if got { "filled" } else { "empty" });
271 }
272 Ok(())
273}
274
275#[test]
276fn test_the_flat_map_repeats_the_world_06() -> Outcome<()> {
277 // 1200 px of a world 800 px round: a ring at 170 degrees east shows twice.
278 let v = res!(view(Projection::WebMercator, 0.0, 0.0, 0.0, TAU * EARTH_RADIUS_M / 800.0,
279 1200.0, 400.0));
280 let isle = ring(&[(-1.0, 169.0), (-1.0, 171.0), (1.0, 171.0), (1.0, 169.0)]);
281 let mut out = ScreenPaths::new();
282 res!(v.project_rings(&[isle], RingMode::Fill, 0.0, 0.25, &mut out));
283 req!(out.len(), 2, "The island showed {} times.", out.len());
284 let (_, r) = v.bounding_cap();
285 req!(r, PI, "A screen wider than the world has a cap of {} rad.", r);
286 Ok(())
287}
288
289#[test]
290fn test_the_globe_fills_the_view_from_inside_a_ring_07() -> Outcome<()> {
291 let mut big = Vec::new();
292 for k in 0..40 {
293 big.push((-40.0, -40.0 + 2.0 * k as f64));
294 }
295 for k in 0..40 {
296 big.push((-40.0 + 2.0 * k as f64, 40.0));
297 }
298 for k in 0..40 {
299 big.push((40.0, 40.0 - 2.0 * k as f64));
300 }
301 for k in 0..40 {
302 big.push((40.0 - 2.0 * k as f64, -40.0));
303 }
304 let rev: Vec<(f64, f64)> = big.iter().rev().cloned().collect();
305 for pts in [big, rev] {
306 let r = ring(&pts);
307 // From right in the middle, close in: no vertex shows, and the screen is filled.
308 let v = res!(view(Projection::Orthographic, 0.0, 0.0, 0.0, 100.0, 390.0, 700.0));
309 let mut out = ScreenPaths::new();
310 res!(v.project_rings(&[r.clone()], RingMode::Fill, 0.0, 0.25, &mut out));
311 let filled = winding(&out, 195.0, 350.0).abs() > 0.5 && winding(&out, 1.0, 1.0).abs() > 0.5;
312 req!(filled, true, "The view from inside the ring is not filled.");
313 // From outside it, close in, nothing is drawn.
314 let v = res!(view(Projection::Orthographic, 0.0, 90.0, 0.0, 100.0, 390.0, 700.0));
315 let mut out = ScreenPaths::new();
316 res!(v.project_rings(&[r.clone()], RingMode::Fill, 0.0, 0.25, &mut out));
317 req!(out.len(), 0, "The view from outside drew {} paths.", out.len());
318 // From the antipode, whole globe: the ring is behind it, and nothing is drawn.
319 let v = res!(view(Projection::Orthographic, 0.0, 180.0, 0.0, EARTH_RADIUS_M / 300.0,
320 700.0, 700.0));
321 let mut out = ScreenPaths::new();
322 res!(v.project_rings(&[r], RingMode::Fill, 0.0, 0.25, &mut out));
323 req!(out.len(), 0, "The far side drew {} paths.", out.len());
324 }
325 Ok(())
326}
327
328#[test]
329fn test_clipped_paths_stay_by_the_screen_08() -> Outcome<()> {
330 let w = res!(world::read(WORLD));
331 let land = res!(w.layer("land", 0).ok_or_else(|| err!("No land."; Test))).unit_rings();
332 let borders = res!(w.layer("borders", 0).ok_or_else(|| err!("No borders."; Test))).unit_rings();
333 let mut rng = Rng(11);
334 for _ in 0..40 {
335 let lat = rng.unit() * 170.0 - 85.0;
336 let lng = rng.unit() * 360.0 - 180.0;
337 let m = 10f64.powf(1.0 + rng.unit() * 4.5);
338 let hd = rng.unit() * 360.0;
339 for kind in [Projection::Orthographic, Projection::WebMercator] {
340 let v = res!(view(kind, lat, lng, hd, m, 390.0, 700.0));
341 let mut out = ScreenPaths::new();
342 res!(v.project_rings(&land, RingMode::Fill, 30_000.0, 0.25, &mut out));
343 res!(v.project_rings(&borders, RingMode::Line, 30_000.0, 0.25, &mut out));
344 // The clip lies eight pixels outside the screen, or round it on the globe.
345 let reach = 195.0f64.hypot(350.0) + 8.0 + 0.01;
346 for p in out.xy.chunks(2) {
347 let (x, y) = (p[0] as f64, p[1] as f64);
348 let ok = match kind {
349 Projection::Orthographic => (x - 195.0).hypot(y - 350.0) <= reach
350 .max(EARTH_RADIUS_M / m + 0.01),
351 Projection::WebMercator => x >= -8.01 && x <= 398.01 && y >= -8.01 && y <= 708.01,
352 };
353 req!(ok, true, "{:?} about {}, {} at {} m/px drew a point at {}, {}.",
354 kind, lat, lng, m, x, y);
355 }
356 }
357 }
358 Ok(())
359}
360
361#[test]
362fn test_a_stroked_ring_is_cut_at_the_horizon_09() -> Outcome<()> {
363 // A parallel at 10 degrees north, seen from the equator: the half facing us is drawn and
364 // its two ends sit on the rim.
365 let pts: Vec<(f64, f64)> = (0..72).map(|k| (10.0, -180.0 + 5.0 * k as f64)).collect();
366 let v = res!(view(Projection::Orthographic, 0.0, 0.0, 0.0, EARTH_RADIUS_M / 300.0,
367 700.0, 700.0));
368 let mut out = ScreenPaths::new();
369 res!(v.project_rings(&[ring(&pts)], RingMode::Outline, 0.0, 0.0, &mut out));
370 req!(out.len(), 1, "The parallel came out in {} pieces.", out.len());
371 let (xy, closed) = res!(out.path(0).ok_or_else(|| err!("No path."; Test)));
372 req!(closed, false);
373 let n = xy.len() / 2;
374 for k in [0, n - 1] {
375 let r = (xy[2 * k] as f64 - 350.0).hypot(xy[2 * k + 1] as f64 - 350.0);
376 let on = (r - 300.0).abs() < 0.01;
377 req!(on, true, "An end of the cut parallel is {} px from the centre.", r);
378 }
379 Ok(())
380}
381
382#[test]
383fn test_the_sea_stays_sea_at_every_turn_10() -> Outcome<()> {
384 // Open water and dry ground, both a matter of record and neither from this repository;
385 // the ground is far enough from any coast that a thirty-kilometre outline cannot reach it.
386 // The lists are Ochre's (`globe.rs`, `dev/verify_world.py`).
387 let sea = [
388 ("the middle of the Indian Ocean", -30.0, 80.0),
389 ("the middle of the Pacific", 0.0, -150.0),
390 ("the middle of the Atlantic", -20.0, -25.0),
391 ("Point Nemo", -48.8767, -123.3933),
392 ("the north Pacific", 40.0, -170.0),
393 ("the Sargasso Sea", 35.0, -45.0),
394 ("the south Atlantic", -40.0, -20.0),
395 ("the south Indian Ocean", -45.0, 90.0),
396 ("the Southern Ocean", -60.0, 0.0),
397 ("the south Pacific", -30.0, -140.0),
398 ("the Arabian Sea", 15.0, 65.0),
399 ("the Arctic Ocean", 85.0, 0.0),
400 ];
401 let land = [
402 ("Moscow", 55.7558, 37.6173),
403 ("Ulaanbaatar", 47.8864, 106.9057),
404 ("Cairo", 30.0444, 31.2357),
405 ("Nairobi", -1.2921, 36.8219),
406 ("Fairbanks", 64.8378, -147.7164),
407 ("Denver", 39.7392, -104.9903),
408 ("Brasilia", -15.7939, -47.8828),
409 ("Novosibirsk", 55.0084, 82.9357),
410 ("Alice Springs", -23.6980, 133.8807),
411 ];
412 let places: Vec<(&str, f64, f64, bool)> = sea.iter().map(|(n, a, b)| (*n, *a, *b, false))
413 .chain(land.iter().map(|(n, a, b)| (*n, *a, *b, true)))
414 .collect();
415 let w = res!(world::read(WORLD));
416 let layer = res!(w.layer("land", 0).ok_or_else(|| err!("No land."; Test)));
417 let rings = layer.unit_rings();
418 let mut asked = 0usize;
419
420 // The whole globe at 36 turns, and again with the picture turned.
421 for heading in [0.0, 30.0] {
422 for lat_0 in [-75.0, -45.0, -15.0, 15.0, 45.0, 75.0] {
423 for lon_0 in [-165.0, -105.0, -45.0, 15.0, 75.0, 135.0] {
424 let v = res!(view(Projection::Orthographic, lat_0, lon_0, heading,
425 EARTH_RADIUS_M / 300.0, 700.0, 700.0));
426 let mut out = ScreenPaths::new();
427 res!(v.project_rings(&rings, RingMode::Fill, layer.tol_m, 0.25, &mut out));
428 for (name, lat, lng, dry) in &places {
429 let p = match v.forward(*lat, *lng) {
430 Some(p) => p,
431 None => continue,
432 };
433 // A band at the rim where a coastline is edge-on means nothing either way.
434 if (p.x - 350.0).hypot(p.y - 350.0) > 0.92 * 300.0 {
435 continue;
436 }
437 asked += 1;
438 let filled = winding(&out, p.x, p.y).abs() > 0.5;
439 req!(filled, *dry, "Globe turned to {}, {} heading {}: {} came out {}.",
440 lat_0, lon_0, heading, name, if filled { "land" } else { "sea" });
441 }
442 }
443 }
444 }
445 // Close in over each place, where the clip is the circle round the screen and a place
446 // inland has no coastline in view at all.
447 for (name, lat, lng, dry) in &places {
448 for m in [200.0, 5_000.0] {
449 let v = res!(view(Projection::Orthographic, *lat, *lng, 15.0, m, 390.0, 700.0));
450 let mut out = ScreenPaths::new();
451 res!(v.project_rings(&rings, RingMode::Fill, layer.tol_m, 0.25, &mut out));
452 asked += 1;
453 let filled = winding(&out, 195.0, 350.0).abs() > 0.5;
454 req!(filled, *dry, "Close over {} at {} m/px it came out {}.", name, m,
455 if filled { "land" } else { "sea" });
456 }
457 }
458 // The flat map, whole and close in.
459 for lon_0 in [-165.0, -105.0, -45.0, 15.0, 75.0, 135.0] {
460 let v = res!(view(Projection::WebMercator, 0.0, lon_0, 0.0,
461 TAU * EARTH_RADIUS_M / 1000.0, 1000.0, 1000.0));
462 let mut out = ScreenPaths::new();
463 res!(v.project_rings(&rings, RingMode::Fill, layer.tol_m, 0.25, &mut out));
464 for (name, lat, lng, dry) in &places {
465 if lat.abs() > 80.0 {
466 continue;
467 }
468 let p = res!(v.forward(*lat, *lng).ok_or_else(|| err!("{} was lost.", name; Test)));
469 asked += 1;
470 let filled = winding(&out, p.x, p.y).abs() > 0.5;
471 req!(filled, *dry, "Map about {}: {} came out {}.", lon_0, name,
472 if filled { "land" } else { "sea" });
473 }
474 }
475 for (name, lat, lng, dry) in &places {
476 let v = res!(view(Projection::WebMercator, *lat, *lng, 0.0, 2_000.0, 390.0, 700.0));
477 let mut out = ScreenPaths::new();
478 res!(v.project_rings(&rings, RingMode::Fill, layer.tol_m, 0.25, &mut out));
479 asked += 1;
480 let filled = winding(&out, 195.0, 350.0).abs() > 0.5;
481 req!(filled, *dry, "Map close over {} came out {}.", name, if filled { "land" } else { "sea" });
482 }
483 let plenty = asked > 400;
484 req!(plenty, true, "Only {} places were asked about.", asked);
485 Ok(())
486}
487
488// ---------------------------------------------------------------------------------------------
489// Cells
490// ---------------------------------------------------------------------------------------------
491
492#[test]
493fn test_a_cell_outline_runs_along_its_edges_11() -> Outcome<()> {
494 let mut rng = Rng(3);
495 for level in [0u8, 2, 7, 15, 22] {
496 for _ in 0..20 {
497 let n = 1u32 << level;
498 let face = (rng.next_u64() % 6) as u8;
499 let cell = res!(Cell::from_face_ij(face, level,
500 (rng.next_u64() % n as u64) as u32, (rng.next_u64() % n as u64) as u32));
501 let segs = 1 + (rng.next_u64() % 9) as u32;
502 let out = cell.outline(segs);
503 req!(out.len(), 4 * segs as usize);
504 let c = cell.corners();
505 let centre = cell.centre_vec();
506 for (k, p) in out.iter().enumerate() {
507 let (a, b) = (c[k / segs as usize].vec, c[(k / segs as usize + 1) % 4].vec);
508 // On the great circle through the edge's two corners.
509 let nrm = [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]];
510 let len = dot(&nrm, &nrm).sqrt();
511 let off = (dot(p, &nrm) / len).abs();
512 // The normal of a short edge is itself only good to a unit in the last
513 // place over the edge's length, so the tolerance grows as cells shrink.
514 let tol = 1.0e-12 + 1.0e-15 / len;
515 req!((off < tol), true, "Level {} outline point {} is {} off its edge.", level, k, off);
516 // And a hair inside it is the cell itself.
517 let t = 1.0e-4; // a ten-thousandth of the way to the centre
518 let q = [p[0] + t * (centre[0] - p[0]), p[1] + t * (centre[1] - p[1]),
519 p[2] + t * (centre[2] - p[2])];
520 let (lat, lng) = vec_lat_lng(&q);
521 let owner = res!(Cell::at(lat, lng, level));
522 req!(owner, cell, "Level {} outline point {} belongs to {}.", level, k, owner);
523 }
524 }
525 }
526 Ok(())
527}
528
529/// Every cell of a level, by enumeration, whose bounding cap meets the cap.
530fn brute(centre: &[f64; 3], radius: f64, level: u8) -> Outcome<HashSet<Cell>> {
531 let n = 1u32 << level;
532 let mut out = HashSet::new();
533 for face in 0..6u8 {
534 for i in 0..n {
535 for j in 0..n {
536 let cell = res!(Cell::from_face_ij(face, level, i, j));
537 let (v, r) = cell.bounding_cap();
538 if angle(&v, centre) <= r + radius + 1.0e-12 {
539 out.insert(cell);
540 }
541 }
542 }
543 }
544 Ok(out)
545}
546
547/// A random point within a cap, uniform over its area.
548fn in_cap(rng: &mut Rng, c: &[f64; 3], radius: f64) -> [f64; 3] {
549 let z = 1.0 - rng.unit() * (1.0 - radius.cos());
550 let t = rng.unit() * TAU;
551 let s = (1.0 - z * z).max(0.0).sqrt();
552 // A frame about the centre.
553 let up = if c[2].abs() < 0.9 { [0.0, 0.0, 1.0] } else { [1.0, 0.0, 0.0] };
554 let mut e = [up[1] * c[2] - up[2] * c[1], up[2] * c[0] - up[0] * c[2], up[0] * c[1] - up[1] * c[0]];
555 let el = dot(&e, &e).sqrt();
556 e = [e[0] / el, e[1] / el, e[2] / el];
557 let nn = [c[1] * e[2] - c[2] * e[1], c[2] * e[0] - c[0] * e[2], c[0] * e[1] - c[1] * e[0]];
558 let (a, b) = (s * t.cos(), s * t.sin());
559 [z * c[0] + a * e[0] + b * nn[0], z * c[1] + a * e[1] + b * nn[1], z * c[2] + a * e[2] + b * nn[2]]
560}
561
562fn check_samples(rng: &mut Rng, c: &[f64; 3], radius: f64, level: u8, got: &[Cell]) -> Outcome<()> {
563 let set: HashSet<Cell> = got.iter().cloned().collect();
564 for _ in 0..500 {
565 let p = in_cap(rng, c, radius);
566 let (lat, lng) = vec_lat_lng(&p);
567 let owner = res!(Cell::at(lat, lng, level));
568 req!(set.contains(&owner), true, "A point in the cap lies in {}, not covered.", owner);
569 let mut holders = 0;
570 for cell in got {
571 if res!(cell.contains(lat, lng)) {
572 holders += 1;
573 }
574 }
575 req!(holders, 1, "A point in the cap lies in {} covered cells.", holders);
576 }
577 Ok(())
578}
579
580#[test]
581fn test_cover_cap_agrees_with_enumeration_12() -> Outcome<()> {
582 let mut rng = Rng(5);
583 // Level 3, every cell, caps from tiny to a third of the sphere, some on cube edges and
584 // corners.
585 let mut caps: Vec<([f64; 3], f64)> = vec![
586 (unit_vec(35.26439, 45.0), 0.2), // a cube vertex
587 (unit_vec(0.0, 45.0), 0.05), // a cube edge
588 (unit_vec(90.0, 0.0), 0.7), // the north pole
589 (unit_vec(-31.9535, 115.8571), 1.1),
590 ];
591 for _ in 0..8 {
592 let v = [rng.unit() - 0.5, rng.unit() - 0.5, rng.unit() - 0.5];
593 caps.push((v, rng.unit()));
594 }
595 for (c, r) in &caps {
596 let len = dot(c, c).sqrt();
597 let cn = [c[0] / len, c[1] / len, c[2] / len];
598 let got = res!(cell::cover_cap(*c, *r, 3, 10_000));
599 let set: HashSet<Cell> = got.iter().cloned().collect();
600 req!(set.len(), got.len(), "The cover listed a cell twice.");
601 let want = res!(brute(&cn, *r, 3));
602 req!(set, want, "Level 3 cover of {:?} radius {} differs from enumeration.", c, r);
603 res!(check_samples(&mut rng, &cn, *r, 3, &got));
604 }
605 // Level 9, every one of its 1.6 million cells, for a cap across a cube vertex.
606 let c = unit_vec(35.0, 44.0);
607 let got = res!(cell::cover_cap(c, 0.02, 9, 100_000));
608 let set: HashSet<Cell> = got.iter().cloned().collect();
609 let want = res!(brute(&c, 0.02, 9));
610 req!(set, want, "Level 9 cover across a cube vertex differs from enumeration.");
611 let faces: HashSet<u8> = got.iter().map(|c| c.face()).collect();
612 req!(faces.len(), 3, "The level 9 cap was meant to touch three faces, touched {}.", faces.len());
613 res!(check_samples(&mut rng, &c, 0.02, 9, &got));
614 Ok(())
615}
616
617#[test]
618fn test_cover_cap_at_street_level_13() -> Outcome<()> {
619 // Level 15 about Sydney, 1.5 km: enumerate a window of the face round the centre wide
620 // enough that nothing on its border qualifies.
621 let mut rng = Rng(9);
622 let c = unit_vec(-33.8688, 151.2093);
623 let r = 1_500.0 / EARTH_RADIUS_M;
624 let got = res!(cell::cover_cap(c, r, 15, 10_000));
625 let home = res!(Cell::at(-33.8688, 151.2093, 15));
626 let (hi, hj) = home.coords();
627 let k = 20i64;
628 let mut want = HashSet::new();
629 for di in -k..=k {
630 for dj in -k..=k {
631 let cell = res!(Cell::from_face_ij(home.face(), 15, (hi as i64 + di) as u32,
632 (hj as i64 + dj) as u32));
633 let (v, cr) = cell.bounding_cap();
634 if angle(&v, &c) <= cr + r + 1.0e-12 {
635 let edge = di.abs() == k || dj.abs() == k;
636 req!(edge, false, "The enumeration window was too small.");
637 want.insert(cell);
638 }
639 }
640 }
641 let set: HashSet<Cell> = got.iter().cloned().collect();
642 req!(set, want, "Level 15 cover about Sydney differs from enumeration.");
643 let ok = got.len() > 50 && got.len() < 400;
644 req!(ok, true, "A 1.5 km cap took {} level 15 cells.", got.len());
645 res!(check_samples(&mut rng, &c, r, 15, &got));
646 // Nearest first: the first cell holds the centre.
647 req!(got[0], home);
648 // A cap too big for the bound is refused rather than walked.
649 let refused = cell::cover_cap(c, 1.0, 9, 100).is_err();
650 req!(refused, true, "A cap of a radian at level 9 was walked into a cover of 100.");
651 Ok(())
652}
653
654#[test]
655fn test_a_level_side_is_the_root_of_its_mean_area_14() -> Outcome<()> {
656 // The oracle is the cells' own spherical areas: every one of the 384 level-3 cells.
657 let mut sum = 0.0;
658 let mut n = 0.0;
659 for face in 0..6u8 {
660 for i in 0..8u32 {
661 for j in 0..8u32 {
662 sum += res!(Cell::from_face_ij(face, 3, i, j)).area();
663 n += 1.0;
664 }
665 }
666 }
667 req!((((sum / n).sqrt() * EARTH_RADIUS_M - cell::mean_side_m(3)).abs() < 1.0), true);
668 req!(((cell::mean_side_m(0) - 9_220_000.0).abs() < 1_000.0), true);
669 // At about 9 m a pixel, a 28-pixel grid is level 15 (281 m cells), as the map plan has it.
670 req!(cell::level_for_scale(9.0, 28.0), 15);
671 for m_per_px in [0.01, 0.3, 9.0, 1_000.0, 100_000.0] {
672 let l = cell::level_for_scale(m_per_px, 28.0);
673 req!((cell::mean_side_m(l) >= m_per_px * 28.0 || l == 0), true);
674 req!((l == cell::MAX_LEVEL || cell::mean_side_m(l + 1) < m_per_px * 28.0), true);
675 }
676 req!(cell::level_for_scale(1.0e9, 28.0), 0);
677 req!(cell::level_for_scale(1.0e-9, 28.0), cell::MAX_LEVEL);
678 Ok(())
679}
680
681#[test]
682fn test_placed_labels_are_on_screen_apart_and_in_rank_order_15() -> Outcome<()> {
683 let mut rng = Rng(15);
684 let mut labels = Vec::new();
685 for k in 0..3000 {
686 let lat = (2.0 * rng.unit() - 1.0).asin().to_degrees();
687 let lng = 360.0 * rng.unit() - 180.0;
688 labels.push(world::Label { lat, lng, rank: (k / 300) as u8, name: fmt!("p{}", k) });
689 }
690 for kind in [Projection::WebMercator, Projection::Orthographic] {
691 for m_per_px in [20_000.0, 2_000.0] {
692 let vp = res!(Viewport::new(kind, -30.0, 140.0, 20.0, m_per_px, 390.0, 700.0));
693 let kept = vp.place_labels(&labels, 40.0, 60);
694 req!(kept.is_empty(), false);
695 req!((kept.len() <= 60), true);
696 let mut last = 0usize;
697 for (a, p) in kept.iter().enumerate() {
698 req!((p.x >= 0.0 && p.y >= 0.0 && p.x <= 390.0 && p.y <= 700.0), true);
699 // Where it says the label is, the projection agrees.
700 let l = &labels[p.index];
701 let q = res!(vp.forward(l.lat, l.lng).ok_or_else(|| err!("behind"; Missing)));
702 req!(((q.x as f32 - p.x).abs() < 0.01 && (q.y as f32 - p.y).abs() < 0.01), true);
703 req!((a == 0 || p.index > last), true);
704 last = p.index;
705 for b in kept.iter().take(a) {
706 req!(((b.x - p.x).hypot(b.y - p.y) >= 40.0), true);
707 }
708 }
709 // Nothing dropped that could have been kept: every unkept visible label is near a
710 // kept one of lower index, unless the budget ran out.
711 if kept.len() < 60 {
712 for (i, l) in labels.iter().enumerate() {
713 if let Some(q) = vp.forward(l.lat, l.lng) {
714 if q.x >= 0.0 && q.y >= 0.0 && q.x <= 390.0 && q.y <= 700.0 {
715 let near = kept.iter().any(|k| k.index <= i
716 && (k.x - q.x as f32).hypot(k.y - q.y as f32) < 40.0);
717 req!((near), true);
718 }
719 }
720 }
721 }
722 }
723 }
724 Ok(())
725}