Oregami
Repositories/oxedyne/fe2o3

oxedyne/fe2o3/fe2o3_geom/src/world.rs

22.8 KiB, 5 runs

created by r1870400018:59718, 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//! An offline world map: coastlines, lakes, borders and place names as unprojected positions.
2//!
3//! # Provenance
4//!
5//! The file layout, the delta coding and the spherical simplifier are Ochre's
6//! (`web/apps/ochre/dev/gen_world.py` and `src/web/globe.rs`, August 2026), moved here when
7//! Oxegen became the second program to draw the world. Version 1 is Ochre's layout byte for
8//! byte with the magic generalised from `OCHRWRLD` to `FE2O3WLD`; version 2 adds label layers.
9//!
10//! # The layout
11//!
12//! ```text
13//! magic "FE2O3WLD" | u16 version | u16 layer count
14//! one 20-byte directory entry a layer:
15//! u8 kind (0 fill, 1 stroke, 2 label) | u8 level of detail | u8 name length
16//! | 7 bytes of ASCII name, zero padded | u16 tolerance in hundreds of metres
17//! | u32 first ring (or label) | u32 ring (or label) count
18//! one u32 vertex count a ring, for every ring of every fill and stroke layer
19//! the rings: i32 lat, i32 lng in whole hundred-thousandths of a degree, then one step a
20//! vertex as i16 dlat, i16 dlng; a step of dlat = -32768 is a jump, and the eight bytes
21//! after it are the next position in full
22//! version 2 only, the labels: i32 lat, i32 lng, u8 rank, u8 name length, UTF-8 name
23//! ```
24//!
25//! All numbers are little-endian. A hundred-thousandth of a degree is 1.11 m of latitude. A
26//! step is taken between rounded positions rather than rounded itself, so error cannot build
27//! up along a coastline, and a step too long for sixteen bits is written as a jump rather than
28//! broken up, so every position in the file is a source position. Longitude steps take the
29//! short way round, so a ring's running longitude may walk past 180 degrees; to a sine that is
30//! no trouble.
31//!
32//! # Simplification
33//!
34//! [`simplify_sphere`] is Douglas and Peucker's, measured in three dimensions against the
35//! chord. The orthographic projection is parallel, hence linear, so the straight line a
36//! canvas draws between two kept vertices is exactly the chord the tolerance was measured
37//! against: nothing needs densifying and nothing is cut at the antimeridian.
38
39use crate::proj::{
40 EARTH_RADIUS_M,
41 Viewport,
42 unit_vec,
43};
44
45use oxedyne_fe2o3_core::prelude::*;
46
47pub const MAGIC: &[u8; 8] = b"FE2O3WLD";
48pub const VERSION_RINGS: u16 = 1; // fill and stroke layers only: Ochre's layout
49pub const VERSION_LABELS: u16 = 2; // adds label layers
50pub const Q: f64 = 100_000.0; // positions are whole hundred-thousandths of a degree
51
52const HEADER: usize = 12;
53const DIRENT: usize = 20;
54const NAME_MAX: usize = 7;
55const STEP_LIMIT: i64 = 32_767;
56const JUMP: i16 = -32_768; // a step that is not a step: a position in full follows
57const HALF_TURN: i64 = 180 * 100_000;
58
59/// How a layer is drawn.
60#[derive(Clone, Copy, Debug, PartialEq, Eq)]
61pub enum LayerKind {
62 Fill, // closed rings, painted
63 Stroke, // open lines, stroked
64 Label, // named points
65}
66
67impl LayerKind {
68 fn code(self) -> u8 {
69 match self {
70 Self::Fill => 0,
71 Self::Stroke => 1,
72 Self::Label => 2,
73 }
74 }
75
76 fn from_code(c: u8) -> Outcome<Self> {
77 match c {
78 0 => Ok(Self::Fill),
79 1 => Ok(Self::Stroke),
80 2 => Ok(Self::Label),
81 _ => Err(err!("A world layer of kind {} is none this reads.", c; Invalid, Input)),
82 }
83 }
84}
85
86/// A named point: a town, a city, a capital.
87#[derive(Clone, Debug, PartialEq)]
88pub struct Label {
89 pub lat: f64, // degrees
90 pub lng: f64, // degrees
91 pub rank: u8, // importance, 0 the most: Natural Earth's `scalerank`
92 pub name: String, // UTF-8, at most 255 bytes
93}
94
95/// One layer of the world at one level of detail.
96#[derive(Clone, Debug, PartialEq)]
97pub struct Layer {
98 pub name: String, // at most 7 ASCII bytes: `land`, `lakes`, `borders`
99 pub kind: LayerKind,
100 pub detail: u8, // level of detail, 0 the coarsest
101 pub tol_m: f64, // simplification tolerance, kept in hundreds of metres
102 pub rings: Vec<Vec<(f64, f64)>>, // (lat, lng) in degrees, for Fill and Stroke
103 pub labels: Vec<Label>, // for Label
104}
105
106impl Layer {
107 /// The rings as unit vectors, which is what [`crate::proj::Viewport::project_rings`]
108 /// draws. A caller drawing every frame converts once and keeps the result.
109 pub fn unit_rings(&self) -> Vec<Vec<[f64; 3]>> {
110 self.rings.iter()
111 .map(|r| r.iter().map(|(lat, lng)| unit_vec(*lat, *lng)).collect())
112 .collect()
113 }
114}
115
116/// A world: its layers, in the order the file holds them.
117#[derive(Clone, Debug, Default, PartialEq)]
118pub struct World {
119 pub layers: Vec<Layer>,
120}
121
122impl World {
123 pub fn layer(&self, name: &str, detail: u8) -> Option<&Layer> {
124 self.layers.iter().find(|l| l.name == name && l.detail == detail)
125 }
126
127 /// The number of levels of detail, one more than the finest.
128 pub fn details(&self) -> u8 {
129 self.layers.iter().map(|l| l.detail.saturating_add(1)).max().unwrap_or(0)
130 }
131
132 /// The level of detail to draw at `m_per_px` ground metres per pixel: the coarsest whose
133 /// rings stray from the truth by at most `max_px` pixels, or the finest held when none is
134 /// that fine. `None` for a world with no rings.
135 ///
136 /// Past the finest level's tolerance the coastline is still drawn, but the eye begins to
137 /// see the simplification rather than the coast, which is the caller's cue to fade it.
138 pub fn detail_for(&self, m_per_px: f64, max_px: f64) -> Option<u8> {
139 let mut levels: Vec<(u8, f64)> = self.layers.iter()
140 .filter(|l| l.kind != LayerKind::Label)
141 .map(|l| (l.detail, l.tol_m))
142 .collect();
143 levels.sort_by(|a, b| a.0.cmp(&b.0));
144 levels.dedup_by_key(|l| l.0);
145 let finest = match levels.last() {
146 Some(l) => l.0,
147 None => return None,
148 };
149 let reach = m_per_px * max_px;
150 for (detail, tol_m) in levels.iter() {
151 if *tol_m <= reach {
152 return Some(*detail);
153 }
154 }
155 Some(finest)
156 }
157
158 /// Adds another world's layers to this one, a layer of the same name and level of detail
159 /// replacing the one held. A coarse world loaded first and a finer one fetched later
160 /// become one world this way, and a finer file's labels supersede a coarse file's.
161 pub fn merge(&mut self, other: World) {
162 for layer in other.layers {
163 match self.layers.iter_mut().find(|l| l.name == layer.name && l.detail == layer.detail) {
164 Some(held) => *held = layer,
165 None => self.layers.push(layer),
166 }
167 }
168 }
169}
170
171/// A label the screen shows, placed by [`Viewport::place_labels`].
172#[derive(Clone, Copy, Debug, PartialEq)]
173pub struct Placed {
174 pub x: f32, // screen pixels
175 pub y: f32,
176 pub index: usize, // into the labels given
177}
178
179impl Viewport {
180 /// The labels to draw, in the order given, each kept only if it lands on the screen and
181 /// at least `gap_px` from every label already kept, until `max` are kept.
182 ///
183 /// Given most important first, as a world file holds them, this is the greedy
184 /// declutter that keeps a capital over the suburbs around it. Text widths are the
185 /// painter's, so the gap is a radius rather than a box.
186 pub fn place_labels(&self, labels: &[Label], gap_px: f64, max: usize) -> Vec<Placed> {
187 let mut kept: Vec<Placed> = Vec::new();
188 if self.check().is_err() || max == 0 {
189 return kept;
190 }
191 let f = self.frame();
192 let g2 = (gap_px.max(0.0) * gap_px.max(0.0)) as f32;
193 for (index, label) in labels.iter().enumerate() {
194 let p = match self.forward_in(&f, label.lat, label.lng) {
195 Some(p) => p,
196 None => continue,
197 };
198 if !(p.x >= 0.0 && p.y >= 0.0 && p.x <= self.w && p.y <= self.h) {
199 continue;
200 }
201 let (x, y) = (p.x as f32, p.y as f32);
202 if kept.iter().any(|k| (k.x - x) * (k.x - x) + (k.y - y) * (k.y - y) < g2) {
203 continue;
204 }
205 kept.push(Placed { x, y, index });
206 if kept.len() >= max {
207 break;
208 }
209 }
210 kept
211 }
212}
213
214/// Reads every level of a world file.
215pub fn read(bytes: &[u8]) -> Outcome<World> {
216 read_detail(bytes, None)
217}
218
219/// Reads a world file, building only the layers of one level of detail when `detail` names
220/// one.
221///
222/// Every ring is walked whatever level it belongs to, because a jump makes a ring's length in
223/// bytes unknowable from its length in vertices, but only the wanted rings are built.
224pub fn read_detail(bytes: &[u8], detail: Option<u8>) -> Outcome<World> {
225 if bytes.len() < HEADER || &bytes[..8] != MAGIC {
226 return Err(err!("The world file does not begin as one."; Invalid, Input));
227 }
228 let version = u16::from_le_bytes([bytes[8], bytes[9]]);
229 if version != VERSION_RINGS && version != VERSION_LABELS {
230 return Err(err!("The world file is version {}, and this reads versions {} and {}.",
231 version, VERSION_RINGS, VERSION_LABELS; Invalid, Input, Version));
232 }
233 let count = u16::from_le_bytes([bytes[10], bytes[11]]) as usize;
234 let table_at = HEADER + count * DIRENT;
235 if bytes.len() < table_at {
236 return Err(err!("The world file's directory is cut short at {} of {} bytes.",
237 bytes.len(), table_at; Invalid, Input));
238 }
239
240 struct Dir {
241 kind: LayerKind,
242 detail: u8,
243 name: String,
244 tol_m: f64,
245 first: usize,
246 n: usize,
247 }
248 let mut dirs: Vec<Dir> = Vec::with_capacity(count);
249 let mut rings = 0usize;
250 let mut labels = 0usize;
251 for i in 0..count {
252 let at = HEADER + i * DIRENT;
253 let kind = res!(LayerKind::from_code(bytes[at]));
254 if kind == LayerKind::Label && version == VERSION_RINGS {
255 return Err(err!("Layer {} of a version {} world file holds labels.", i, version;
256 Invalid, Input));
257 }
258 let len = bytes[at + 2] as usize;
259 if len > NAME_MAX {
260 return Err(err!("Layer {} claims a name of {} bytes, more than {}.", i, len, NAME_MAX;
261 Invalid, Input));
262 }
263 let name = match std::str::from_utf8(&bytes[at + 3..at + 3 + len]) {
264 Ok(s) => s.to_string(),
265 Err(_) => return Err(err!("Layer {} has a name that is not text.", i; Invalid, Input)),
266 };
267 let tol_m = u16::from_le_bytes([bytes[at + 10], bytes[at + 11]]) as f64 * 100.0;
268 let first = res!(u32_at(bytes, at + 12)) as usize;
269 let n = res!(u32_at(bytes, at + 16)) as usize;
270 match kind {
271 LayerKind::Label => labels = labels.max(first + n),
272 _ => rings = rings.max(first + n),
273 }
274 dirs.push(Dir { kind, detail: bytes[at + 1], name, tol_m, first, n });
275 }
276 let mut at = table_at + rings * 4;
277 if bytes.len() < at {
278 return Err(err!("The world file's ring table is cut short at {} of {} bytes.",
279 bytes.len(), at; Invalid, Input));
280 }
281
282 let mut wanted = vec![false; rings];
283 for d in &dirs {
284 if d.kind != LayerKind::Label && detail.map_or(true, |w| w == d.detail) {
285 for r in d.first..(d.first + d.n).min(rings) {
286 wanted[r] = true;
287 }
288 }
289 }
290 let mut built: Vec<Vec<(f64, f64)>> = vec![Vec::new(); rings];
291 for r in 0..rings {
292 let n = res!(u32_at(bytes, table_at + r * 4)) as usize;
293 let (drawn, next) = res!(read_ring(bytes, at, n, wanted[r], r));
294 built[r] = drawn;
295 at = next;
296 }
297
298 let mut read_labels: Vec<Label> = Vec::with_capacity(labels);
299 for i in 0..labels {
300 if at + 10 > bytes.len() {
301 return Err(err!("The world file ends inside label {} of {}.", i, labels;
302 Invalid, Input));
303 }
304 let lat = res!(u32_at(bytes, at)) as i32;
305 let lng = res!(u32_at(bytes, at + 4)) as i32;
306 let rank = bytes[at + 8];
307 let len = bytes[at + 9] as usize;
308 let end = at + 10 + len;
309 if end > bytes.len() {
310 return Err(err!("The world file ends inside the name of label {}.", i; Invalid, Input));
311 }
312 let name = match std::str::from_utf8(&bytes[at + 10..end]) {
313 Ok(s) => s.to_string(),
314 Err(_) => return Err(err!("Label {} has a name that is not UTF-8.", i; Invalid, Input)),
315 };
316 read_labels.push(Label { lat: lat as f64 / Q, lng: lng as f64 / Q, rank, name });
317 at = end;
318 }
319 if at != bytes.len() {
320 return Err(err!("The world file has {} bytes after its last position.",
321 bytes.len() as i64 - at as i64; Invalid, Input));
322 }
323
324 let mut world = World::default();
325 for d in dirs {
326 if !detail.map_or(true, |w| w == d.detail) {
327 continue;
328 }
329 let mut layer = Layer {
330 name: d.name, kind: d.kind, detail: d.detail, tol_m: d.tol_m,
331 rings: Vec::new(), labels: Vec::new(),
332 };
333 match d.kind {
334 LayerKind::Label => {
335 layer.labels = read_labels[d.first..d.first + d.n].to_vec();
336 },
337 _ => {
338 layer.rings = built[d.first..d.first + d.n].to_vec();
339 },
340 }
341 world.layers.push(layer);
342 }
343 Ok(world)
344}
345
346/// One ring as degrees, and where the next one begins.
347fn read_ring(bytes: &[u8], at: usize, n: usize, keep: bool, r: usize)
348 -> Outcome<(Vec<(f64, f64)>, usize)>
349{
350 if at + 8 > bytes.len() || n == 0 {
351 return Err(err!("The world file ends inside ring {} at byte {}.", r, at; Invalid, Input));
352 }
353 let mut lat = res!(u32_at(bytes, at)) as i32;
354 let mut lng = res!(u32_at(bytes, at + 4)) as i32;
355 let mut out = Vec::with_capacity(if keep { n } else { 0 });
356 if keep {
357 out.push((lat as f64 / Q, lng as f64 / Q));
358 }
359 let mut step = at + 8;
360 for _ in 1..n {
361 if step + 4 > bytes.len() {
362 return Err(err!("The world file ends inside ring {} at byte {}.", r, step;
363 Invalid, Input));
364 }
365 let dlat = i16::from_le_bytes([bytes[step], bytes[step + 1]]);
366 if dlat == JUMP {
367 if step + 12 > bytes.len() {
368 return Err(err!("The world file ends inside a jump in ring {}.", r; Invalid, Input));
369 }
370 lat = res!(u32_at(bytes, step + 4)) as i32;
371 lng = res!(u32_at(bytes, step + 8)) as i32;
372 step += 12;
373 } else {
374 lat = lat.wrapping_add(dlat as i32);
375 lng = lng.wrapping_add(i16::from_le_bytes([bytes[step + 2], bytes[step + 3]]) as i32);
376 step += 4;
377 }
378 if keep {
379 out.push((lat as f64 / Q, lng as f64 / Q));
380 }
381 }
382 Ok((out, step))
383}
384
385fn u32_at(bytes: &[u8], at: usize) -> Outcome<u32> {
386 match bytes.get(at..at + 4) {
387 Some(b) => Ok(u32::from_le_bytes([b[0], b[1], b[2], b[3]])),
388 None => Err(err!("The world file ends inside a number at byte {}.", at; Invalid, Input)),
389 }
390}
391
392/// Writes a world file: version 1, Ochre's layout, when it holds no labels, and version 2
393/// otherwise.
394///
395/// Positions are rounded to the nearest hundred-thousandth of a degree, ties to even, and a
396/// vertex that rounds onto the one before it is dropped; a ring left with fewer than two
397/// vertices is dropped with it. Reading the result back and writing it again gives the same
398/// bytes.
399pub fn write(world: &World) -> Outcome<Vec<u8>> {
400 let labelled = world.layers.iter().any(|l| l.kind == LayerKind::Label);
401 let version = if labelled { VERSION_LABELS } else { VERSION_RINGS };
402 if world.layers.len() > u16::MAX as usize {
403 return Err(err!("A world of {} layers is more than a file can list.", world.layers.len();
404 Invalid, Input, Excessive));
405 }
406 let mut head: Vec<u8> = Vec::with_capacity(HEADER + DIRENT * world.layers.len());
407 head.extend_from_slice(MAGIC);
408 head.extend_from_slice(&version.to_le_bytes());
409 head.extend_from_slice(&(world.layers.len() as u16).to_le_bytes());
410 let mut table: Vec<u8> = Vec::new();
411 let mut coords: Vec<u8> = Vec::new();
412 let mut names: Vec<u8> = Vec::new();
413 let mut ring_at = 0u32;
414 let mut label_at = 0u32;
415 for (i, layer) in world.layers.iter().enumerate() {
416 let name = layer.name.as_bytes();
417 if name.len() > NAME_MAX || !layer.name.is_ascii() {
418 return Err(err!("Layer {} is named {:?}; a name is at most {} ASCII bytes.",
419 i, layer.name, NAME_MAX; Invalid, Input));
420 }
421 let tol = (layer.tol_m / 100.0).round();
422 if !(tol >= 0.0 && tol <= u16::MAX as f64) {
423 return Err(err!("Layer {} has a tolerance of {} m, outside 0 to {} m.",
424 i, layer.tol_m, u16::MAX as f64 * 100.0; Invalid, Input, Range));
425 }
426 let (first, n) = match layer.kind {
427 LayerKind::Label => {
428 for (k, label) in layer.labels.iter().enumerate() {
429 let text = label.name.as_bytes();
430 if text.len() > u8::MAX as usize {
431 return Err(err!("Label {} of layer {} has a name of {} bytes, more than {}.",
432 k, i, text.len(), u8::MAX; Invalid, Input, Excessive));
433 }
434 names.extend_from_slice(&res!(quantise(label.lat, i)).to_le_bytes());
435 names.extend_from_slice(&res!(quantise(label.lng, i)).to_le_bytes());
436 names.push(label.rank);
437 names.push(text.len() as u8);
438 names.extend_from_slice(text);
439 }
440 let first = label_at;
441 label_at += layer.labels.len() as u32;
442 (first, layer.labels.len() as u32)
443 },
444 _ => {
445 let mut kept = 0u32;
446 for ring in &layer.rings {
447 if let Some(n) = res!(write_ring(ring, i, &mut coords)) {
448 table.extend_from_slice(&n.to_le_bytes());
449 kept += 1;
450 }
451 }
452 let first = ring_at;
453 ring_at += kept;
454 (first, kept)
455 },
456 };
457 head.push(layer.kind.code());
458 head.push(layer.detail);
459 head.push(name.len() as u8);
460 let mut padded = [0u8; NAME_MAX];
461 padded[..name.len()].copy_from_slice(name);
462 head.extend_from_slice(&padded);
463 head.extend_from_slice(&(tol as u16).to_le_bytes());
464 head.extend_from_slice(&first.to_le_bytes());
465 head.extend_from_slice(&n.to_le_bytes());
466 }
467 head.extend_from_slice(&table);
468 head.extend_from_slice(&coords);
469 head.extend_from_slice(&names);
470 Ok(head)
471}
472
473/// A coordinate in whole hundred-thousandths of a degree, rounded half to even.
474fn quantise(deg: f64, layer: usize) -> Outcome<i32> {
475 let q = (deg * Q).round_ties_even();
476 if !(q >= i32::MIN as f64 && q <= i32::MAX as f64) {
477 return Err(err!("A coordinate of {} degrees in layer {} cannot be written.", deg, layer;
478 Invalid, Input, Range));
479 }
480 Ok(q as i32)
481}
482
483/// Appends one ring, returning how many vertices it came to, or `None` if it came to fewer
484/// than two and was not written.
485fn write_ring(pts: &[(f64, f64)], layer: usize, out: &mut Vec<u8>) -> Outcome<Option<u32>> {
486 let mut grid: Vec<(i64, i64)> = Vec::with_capacity(pts.len());
487 for (lat, lng) in pts {
488 let q = (res!(quantise(*lat, layer)) as i64, res!(quantise(*lng, layer)) as i64);
489 if grid.last() != Some(&q) {
490 grid.push(q);
491 }
492 }
493 if grid.len() < 2 {
494 return Ok(None);
495 }
496 out.extend_from_slice(&(grid[0].0 as i32).to_le_bytes());
497 out.extend_from_slice(&(grid[0].1 as i32).to_le_bytes());
498 let (mut py, mut px) = grid[0];
499 for (qy, qx) in grid.iter().skip(1) {
500 let dy = qy - py;
501 // The short way round: a ring crossing the antimeridian steps a whole turn in the data
502 // and a few hundred metres on the ground.
503 let dx = (qx - px + HALF_TURN).rem_euclid(2 * HALF_TURN) - HALF_TURN;
504 let (ny, nx) = (py + dy, px + dx);
505 if dy.abs() > STEP_LIMIT || dx.abs() > STEP_LIMIT {
506 if ny < i32::MIN as i64 || ny > i32::MAX as i64 || nx < i32::MIN as i64 || nx > i32::MAX as i64 {
507 return Err(err!("A ring in layer {} wanders past what a jump can hold.", layer;
508 Invalid, Input, Range));
509 }
510 out.extend_from_slice(&JUMP.to_le_bytes());
511 out.extend_from_slice(&0i16.to_le_bytes());
512 out.extend_from_slice(&(ny as i32).to_le_bytes());
513 out.extend_from_slice(&(nx as i32).to_le_bytes());
514 } else {
515 out.extend_from_slice(&(dy as i16).to_le_bytes());
516 out.extend_from_slice(&(dx as i16).to_le_bytes());
517 }
518 py = ny;
519 px = nx;
520 }
521 Ok(Some(grid.len() as u32))
522}
523
524/// Douglas and Peucker's simplification of a run of unit vectors, against the chord, returning
525/// the indices kept in order.
526///
527/// The distance measured is the perpendicular in three dimensions from a vertex to the straight
528/// segment through the kept pair either side of it, on a sphere of [`EARTH_RADIUS_M`], which is
529/// the distance from the edge a globe actually draws. The first and last vertices are always
530/// kept. Iterative, so a coastline of a million vertices cannot overflow the stack.
531pub fn simplify_sphere(pts: &[[f64; 3]], eps_m: f64) -> Vec<usize> {
532 let n = pts.len();
533 if n < 3 || !(eps_m > 0.0) {
534 return (0..n).collect();
535 }
536 let eps = eps_m / EARTH_RADIUS_M;
537 let e2 = eps * eps;
538 let mut keep = vec![false; n];
539 keep[0] = true;
540 keep[n - 1] = true;
541 let mut stack: Vec<(usize, usize)> = vec![(0, n - 1)];
542 while let Some((i, j)) = stack.pop() {
543 if j <= i + 1 {
544 continue;
545 }
546 let a = pts[i];
547 let b = pts[j];
548 let (dx, dy, dz) = (b[0] - a[0], b[1] - a[1], b[2] - a[2]);
549 let dd = dx * dx + dy * dy + dz * dz;
550 let mut best = -1.0;
551 let mut at = i;
552 for k in (i + 1)..j {
553 let p = pts[k];
554 let d2 = if dd == 0.0 {
555 sq(p[0] - a[0]) + sq(p[1] - a[1]) + sq(p[2] - a[2])
556 } else {
557 let t = (((p[0] - a[0]) * dx + (p[1] - a[1]) * dy + (p[2] - a[2]) * dz) / dd)
558 .clamp(0.0, 1.0);
559 sq(p[0] - a[0] - t * dx) + sq(p[1] - a[1] - t * dy) + sq(p[2] - a[2] - t * dz)
560 };
561 if d2 > best {
562 best = d2;
563 at = k;
564 }
565 }
566 if best > e2 {
567 keep[at] = true;
568 stack.push((i, at));
569 stack.push((at, j));
570 }
571 }
572 (0..n).filter(|k| keep[*k]).collect()
573}
574
575fn sq(x: f64) -> f64 { x * x }