spirograph.rs
39.6 kB · rust · 983 lines
1use crate::factor::gcd;2use mrlycore::errors::{value_error, Result};3use mrlycore::Rng;4use std::collections::HashSet;5use std::f64::consts::{PI, TAU};67/// The largest radius of a ring or a wheel.8pub const RADIUS_CAP: usize = 96;9/// The fewest and the most sides of a polygon track.10pub const SIDES: (usize, usize) = (3, 12);11/// The most laps of a line or a polygon.12pub const LAPS_CAP: usize = 24;13/// The most pencils a wheel seats.14pub const PENCIL_CAP: usize = 4096;15/// The most points one trace returns.16pub const POINT_CAP: usize = 4_000_000;17const REACH: (f64, f64) = (0.05, 2.0);18const RING: usize = 360;1920// THE PENCILS2122/// What a pencil sits on: a filled cell, an empty cell, or a corner of a filled cell.23#[derive(Clone, Copy, Debug, PartialEq, Eq)]24pub enum Kind {25 /// The centre of a filled cell.26 Fill,27 /// The centre of an empty cell.28 Void,29 /// A corner of a filled cell.30 Corner,31}3233/// A pencil on the wheel: its seat in units of the wheel's radius with the tile centre at the origin, the exact seat it came from, and its kind.34#[derive(Clone, Copy, Debug, PartialEq)]35pub struct Pencil {36 /// The seat's abscissa, in wheel radii.37 pub x: f64,38 /// The seat's ordinate, up the page, in wheel radii.39 pub y: f64,40 /// The exact seat, twice the cell coordinates from the tile centre, before any jitter.41 pub seat: (i64, i64),42 /// What the pencil sits on.43 pub kind: Kind,44}4546/// The mass of a byte grid taken as a wheel: how many pencils of each kind it seats.47#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]48pub struct Seats {49 /// The pencils on filled cells.50 pub fills: usize,51 /// The pencils on empty cells.52 pub voids: usize,53 /// The pencils on corners.54 pub corners: usize,55}5657/// Seats one pencil per chosen site of a byte grid: `fill` the filled cells, `void` the empty ones, `both`, or `corners` the corners of the filled cells, each once. The tile is scaled so its circumradius is `reach` wheel radii, and `jitter` moves every seat by up to that fraction of a cell each way, seeded.58pub fn pencils(59 types: &[u8],60 width: usize,61 height: usize,62 mode: &str,63 reach: f64,64 jitter: f64,65 seed: u32,66) -> Result<Vec<Pencil>> {67 if width == 0 || height == 0 || types.len() != width * height {68 return value_error("the wheel must be a grid of width by height bytes.");69 }70 if !(REACH.0..=REACH.1).contains(&reach) {71 return value_error(format!(72 "the reach must be between {} and {} wheel radii.",73 REACH.0, REACH.174 ));75 }76 if !(0.0..=1.0).contains(&jitter) {77 return value_error("the jitter must be between 0 and 1 cells.");78 }79 if !["fill", "void", "both", "corners"].contains(&mode) {80 return value_error("the pencils sit on fill, void, both or corners.");81 }82 let (w, h) = (width as i64, height as i64);83 let mut seen = HashSet::new();84 let mut out = Vec::new();85 for i in 0..height {86 for j in 0..width {87 let on = types[i * width + j] != 0;88 let (i, j) = (i as i64, j as i64);89 if mode == "corners" {90 if !on {91 continue;92 }93 for (di, dj) in [(0, 0), (0, 1), (1, 0), (1, 1)] {94 let seat = (2 * (j + dj) - w, h - 2 * (i + di));95 if seen.insert(seat) {96 out.push(seat_pencil(seat, Kind::Corner));97 }98 }99 } else if (on && mode != "void") || (!on && mode != "fill") {100 let kind = if on { Kind::Fill } else { Kind::Void };101 out.push(seat_pencil((2 * j + 1 - w, h - 2 * i - 1), kind));102 }103 }104 }105 if out.len() > PENCIL_CAP {106 return value_error(format!(107 "{} pencils is past the cap of {PENCIL_CAP}: take a lower level or a smaller tile.",108 out.len()109 ));110 }111 let unit = reach / (width as f64 / 2.0).hypot(height as f64 / 2.0);112 let mut rng = Rng::new(u64::from(seed));113 for pencil in &mut out {114 let (dx, dy) = if jitter > 0.0 {115 ((rng.unit() - 0.5) * jitter, (rng.unit() - 0.5) * jitter)116 } else {117 (0.0, 0.0)118 };119 pencil.x = (pencil.seat.0 as f64 / 2.0 + dx) * unit;120 pencil.y = (pencil.seat.1 as f64 / 2.0 + dy) * unit;121 }122 Ok(out)123}124125/// The side of one cell in wheel radii at a reach, the number the page needs to draw the tile on the wheel.126pub fn cell(width: usize, height: usize, reach: f64) -> f64 {127 reach / (width as f64 / 2.0).hypot(height as f64 / 2.0)128}129130/// Counts the pencils by kind.131pub fn seats(pencils: &[Pencil]) -> Seats {132 let mut out = Seats::default();133 for pencil in pencils {134 match pencil.kind {135 Kind::Fill => out.fills += 1,136 Kind::Void => out.voids += 1,137 Kind::Corner => out.corners += 1,138 }139 }140 out141}142143fn seat_pencil(seat: (i64, i64), kind: Kind) -> Pencil {144 Pencil {145 x: 0.0,146 y: 0.0,147 seat,148 kind,149 }150}151152// THE TRACK153154/// One piece of the centre path: a straight run, or a turn about a point.155#[derive(Clone, Copy, Debug, PartialEq)]156pub enum Piece {157 /// A straight run from one point to another.158 Run {159 /// Where the run starts.160 from: (f64, f64),161 /// Where the run ends.162 to: (f64, f64),163 },164 /// A turn about a point at a radius, counterclockwise from one angle to another.165 Turn {166 /// The point turned about.167 about: (f64, f64),168 /// The radius of the turn.169 radius: f64,170 /// The angle the turn starts at.171 from: f64,172 /// The angle the turn ends at, past the start.173 to: f64,174 },175}176177impl Piece {178 fn len(&self) -> f64 {179 match *self {180 Piece::Run { from, to } => (to.0 - from.0).hypot(to.1 - from.1),181 Piece::Turn {182 radius, from, to, ..183 } => radius * (to - from),184 }185 }186187 fn at(&self, s: f64) -> (f64, f64) {188 match *self {189 Piece::Run { from, to } => {190 let len = self.len();191 let t = if len > 0.0 {192 (s / len).clamp(0.0, 1.0)193 } else {194 0.0195 };196 (from.0 + (to.0 - from.0) * t, from.1 + (to.1 - from.1) * t)197 }198 Piece::Turn {199 about,200 radius,201 from,202 to,203 } => {204 let angle = (from + s / radius).min(to);205 (206 about.0 + radius * angle.cos(),207 about.1 + radius * angle.sin(),208 )209 }210 }211 }212213 fn bounds(&self) -> [f64; 4] {214 match *self {215 Piece::Run { from, to } => [216 from.0.min(to.0),217 from.1.min(to.1),218 from.0.max(to.0),219 from.1.max(to.1),220 ],221 Piece::Turn { about, radius, .. } => [222 about.0 - radius,223 about.1 - radius,224 about.0 + radius,225 about.1 + radius,226 ],227 }228 }229}230231/// A track and the path the wheel's centre takes along it: the wheel of radius `wheel` rolls without slipping, on the left of the track when `side` is minus one and on the right when it is plus one, and turns by `side` times the centre's path length over the wheel's radius.232#[derive(Clone, Debug, PartialEq)]233pub struct Track {234 /// The kind: `line`, `in`, `out`, `polyin` or `polyout`.235 pub kind: String,236 /// The wheel's radius.237 pub wheel: f64,238 /// The pieces of the centre path, in order.239 pub pieces: Vec<Piece>,240 /// The length of the centre path.241 pub total: f64,242 /// Minus one inside the track, plus one outside it.243 pub side: f64,244 /// The track itself as a polyline, closed when `closed` says so.245 pub outline: Vec<(f64, f64)>,246 /// Whether the outline closes on itself.247 pub closed: bool,248 /// The ring's radius over the wheel's in lowest terms on a circle, zero over zero elsewhere.249 pub ratio: (usize, usize),250 /// How many times the centre goes round: the ratio's denominator on a circle, the laps on a polygon, one on a line.251 pub orbits: usize,252 /// The rotation order of the whole picture: the ratio's numerator on a circle, none elsewhere.253 pub fold: usize,254}255256/// Lays a track: `line` a straight line under the wheel for `laps` turns; `in` and `out` a circle of radius `ring` with the wheel inside or outside, closing after the reduced denominator of `ring` over `wheel` orbits; `polyin` and `polyout` a regular polygon of `sides` sides and circumradius `ring` for `laps` laps.257pub fn track(kind: &str, ring: usize, wheel: usize, sides: usize, laps: usize) -> Result<Track> {258 if !(1..=RADIUS_CAP).contains(&wheel) || !(1..=RADIUS_CAP).contains(&ring) {259 return value_error(format!(260 "the ring and the wheel must have radius between 1 and {RADIUS_CAP}."261 ));262 }263 if !(1..=LAPS_CAP).contains(&laps) {264 return value_error(format!("the laps must be between 1 and {LAPS_CAP}."));265 }266 let r = wheel as f64;267 let big = ring as f64;268 let mut out = Track {269 kind: kind.to_string(),270 wheel: r,271 pieces: Vec::new(),272 total: 0.0,273 side: -1.0,274 outline: Vec::new(),275 closed: false,276 ratio: (0, 0),277 orbits: laps,278 fold: 0,279 };280 match kind {281 "line" => {282 let run = laps as f64 * TAU * r;283 out.pieces.push(Piece::Run {284 from: (0.0, r),285 to: (run, r),286 });287 out.outline = vec![(-r, 0.0), (run + r, 0.0)];288 out.orbits = 1;289 }290 "in" | "out" => {291 let inside = kind == "in";292 if inside && ring <= wheel {293 return value_error("the wheel must be smaller than the ring it rolls inside.");294 }295 let g = gcd(ring, wheel);296 let (a, b) = (ring / g, wheel / g);297 let rho = if inside { big - r } else { big + r };298 out.side = if inside { -1.0 } else { 1.0 };299 out.pieces.push(Piece::Turn {300 about: (0.0, 0.0),301 radius: rho,302 from: 0.0,303 to: TAU * b as f64,304 });305 out.outline = (0..RING)306 .map(|k| {307 let angle = TAU * k as f64 / RING as f64;308 (big * angle.cos(), big * angle.sin())309 })310 .collect();311 out.closed = true;312 out.ratio = (a, b);313 out.orbits = b;314 out.fold = a;315 }316 "polyin" | "polyout" => {317 if !(SIDES.0..=SIDES.1).contains(&sides) {318 return value_error(format!(319 "the polygon must have between {} and {} sides.",320 SIDES.0, SIDES.1321 ));322 }323 let corners = polygon(big, sides);324 out.outline = corners.clone();325 out.closed = true;326 if kind == "polyin" {327 let inset = big * (PI / sides as f64).cos() - r;328 if inset <= 0.0 {329 return value_error("the wheel does not fit inside the polygon.");330 }331 let inner = polygon(inset / (PI / sides as f64).cos(), sides);332 for _ in 0..laps {333 for k in 0..sides {334 out.pieces.push(Piece::Run {335 from: inner[k],336 to: inner[(k + 1) % sides],337 });338 }339 }340 } else {341 out.side = 1.0;342 let normal = |k: usize| {343 let (p, q) = (corners[k], corners[(k + 1) % sides]);344 let (dx, dy) = (q.0 - p.0, q.1 - p.1);345 let len = dx.hypot(dy);346 (dy / len, -dx / len)347 };348 for _ in 0..laps {349 for k in 0..sides {350 let (p, q) = (corners[k], corners[(k + 1) % sides]);351 let (m, next) = (normal(k), normal((k + 1) % sides));352 out.pieces.push(Piece::Run {353 from: (p.0 + r * m.0, p.1 + r * m.1),354 to: (q.0 + r * m.0, q.1 + r * m.1),355 });356 let from = m.1.atan2(m.0);357 let mut to = next.1.atan2(next.0);358 while to < from {359 to += TAU;360 }361 out.pieces.push(Piece::Turn {362 about: q,363 radius: r,364 from,365 to,366 });367 }368 }369 }370 }371 _ => return value_error("the track is line, in, out, polyin or polyout."),372 }373 out.total = out.pieces.iter().map(Piece::len).sum();374 Ok(out)375}376377fn polygon(radius: f64, sides: usize) -> Vec<(f64, f64)> {378 (0..sides)379 .map(|k| {380 let angle = PI / sides as f64 + TAU * k as f64 / sides as f64;381 (radius * angle.cos(), radius * angle.sin())382 })383 .collect()384}385386/// The wheel's centre after `s` of path length.387pub fn pose(track: &Track, s: f64) -> (f64, f64) {388 let mut left = s.clamp(0.0, track.total);389 let last = track.pieces.len() - 1;390 for (k, piece) in track.pieces.iter().enumerate() {391 let len = piece.len();392 if left <= len || k == last {393 return piece.at(left);394 }395 left -= len;396 }397 (0.0, 0.0)398}399400/// The wheel's turn after `s` of path length, in radians: `side` times `s` over the wheel's radius.401pub fn turn(track: &Track, s: f64) -> f64 {402 let angle = track.side * s / track.wheel;403 if angle == 0.0 {404 0.0405 } else {406 angle407 }408}409410/// Where a pencil is after `s` of path length: the centre plus the seat turned with the wheel.411pub fn point(track: &Track, pencil: &Pencil, s: f64) -> (f64, f64) {412 let (cx, cy) = pose(track, s);413 let phi = turn(track, s);414 let (c, sn) = (phi.cos(), phi.sin());415 let r = track.wheel;416 (417 cx + r * (pencil.x * c - pencil.y * sn),418 cy + r * (pencil.x * sn + pencil.y * c),419 )420}421422/// Traces every pencil along the whole track at `samples` evenly spaced path lengths, first and last included: pencil by pencil, sample by sample, x then y.423pub fn trace(track: &Track, pencils: &[Pencil], samples: usize) -> Result<Vec<f32>> {424 if samples < 2 {425 return value_error("a trace needs at least two samples.");426 }427 if pencils.len() * samples > POINT_CAP {428 return value_error(format!(429 "{} points is past the cap of {POINT_CAP}: fewer samples or fewer pencils.",430 pencils.len() * samples431 ));432 }433 let mut out = Vec::with_capacity(pencils.len() * samples * 2);434 for pencil in pencils {435 for k in 0..samples {436 let s = track.total * k as f64 / (samples - 1) as f64;437 let (x, y) = point(track, pencil, s);438 out.push(x as f32);439 out.push(y as f32);440 }441 }442 Ok(out)443}444445/// How many classes `representatives` finds: the distinct curves on a circle track, the shapes up to a shift along a line track, one class per pencil on a polygon and under jitter.446pub fn distinct(track: &Track, pencils: &[Pencil], exact: bool) -> usize {447 representatives(track, pencils, exact).len()448}449450/// One pencil per class, the first index of every class in the order the pencils were seated. On a circle track the classes are the distinct curves, by the coincidence law on exact seats: two pencils draw one curve iff a rotation of a full turn over the ratio's denominator carries one seat to the other, which the square lattice allows only by half turns when the denominator is even and by quarter turns when four divides it. On a line track the classes are the seat radii, and those are shapes up to a shift, not curves: turning a seat by `gamma` slides its whole ribbon `gamma` wheel radii along the line while the ribbon's period is a full turn of the wheel, so two seats of one radius draw translates of one shape and share no point unless the seats are equal. On a polygon every pencil is its own class, and so is every pencil under jitter.451pub fn representatives(track: &Track, pencils: &[Pencil], exact: bool) -> Vec<usize> {452 if !exact {453 return (0..pencils.len()).collect();454 }455 let mut out = Vec::new();456 match track.kind.as_str() {457 "line" => {458 let mut seen = HashSet::new();459 for (k, p) in pencils.iter().enumerate() {460 if seen.insert(p.seat.0 * p.seat.0 + p.seat.1 * p.seat.1) {461 out.push(k);462 }463 }464 }465 "in" | "out" => {466 let order = gcd(track.ratio.1, 4);467 let quarter = |(u, v): (i64, i64)| (-v, u);468 let mut seen = HashSet::new();469 for (k, p) in pencils.iter().enumerate() {470 let mut least = p.seat;471 let mut seat = p.seat;472 for _ in 1..order {473 for _ in 0..4 / order {474 seat = quarter(seat);475 }476 least = least.min(seat);477 }478 if seen.insert(least) {479 out.push(k);480 }481 }482 }483 _ => out.extend(0..pencils.len()),484 }485 out486}487488/// The box the whole picture sits in: the centre path and the track, padded by the wheel's radius or the farthest seat, whichever reaches further.489pub fn frame(track: &Track, pencils: &[Pencil]) -> [f64; 4] {490 let reach = pencils491 .iter()492 .map(|p| p.x.hypot(p.y))493 .fold(1.0_f64, f64::max);494 let pad = track.wheel * reach;495 let mut box_ = [f64::MAX, f64::MAX, f64::MIN, f64::MIN];496 let mut take = |x: f64, y: f64| {497 box_[0] = box_[0].min(x);498 box_[1] = box_[1].min(y);499 box_[2] = box_[2].max(x);500 box_[3] = box_[3].max(y);501 };502 for piece in &track.pieces {503 let b = piece.bounds();504 take(b[0] - pad, b[1] - pad);505 take(b[2] + pad, b[3] + pad);506 }507 for &(x, y) in &track.outline {508 take(x, y);509 }510 box_511}512513// THE NODES514515/// The crossings of the whole roulette on a circle track, the generic count, with `R/r = a/b` in lowest terms. Write `|p|` for a seat's distance from the wheel's centre in wheel radii and `A` for the centre path's radius in the same units, `(a - b)/b` inside and `(a + b)/b` outside. Every seat must lie strictly inside the window `0 < |p| < min(1, A)`, which three hypotheses cut: `|p| > 0`, since a seat at the wheel's centre draws the centre circle `b` times over and never crosses; `|p| < 1`, the loop threshold, past which a curve loops; and `|p| < A`, the seat threshold, where the seat reaches the centre path, which comes before the loop threshold on every inside track with `a < 2b` and never bites outside. Inside that window two distinct curves cross exactly `2ab` times, one curve crosses itself `a(b - 1)` times, and `k` distinct curves cross `2ab k(k - 1) / 2 + k a (b - 1)` times, the design entering only through `k`. `exact` reads the coincidence law on the seats, as `distinct` does. `None` on a line or a polygon track, and `None` when any seat leaves the window, where neither count is the law's. At isolated reaches some crossings merge, so the count holds for the generic reach.516pub fn nodes(track: &Track, pencils: &[Pencil], exact: bool) -> Option<u64> {517 if track.kind != "in" && track.kind != "out" {518 return None;519 }520 let (a, b) = (track.ratio.0 as u64, track.ratio.1 as u64);521 let window = (a as f64 / b as f64 + track.side).abs().min(1.0);522 let out = pencils.iter().any(|p| {523 let seat = p.x.hypot(p.y);524 !(seat > 0.0 && seat < window)525 });526 if out {527 return None;528 }529 let k = distinct(track, pencils, exact) as u64;530 Some(2 * a * b * (k * k.saturating_sub(1) / 2) + k * a * (b - 1))531}532533// THE COVER534535/// The largest raster side a cover rasters.536pub const RASTER_CAP: usize = 4096;537538/// The disc a circle roulette sits in, in the track's units.539#[derive(Clone, Copy, Debug, PartialEq)]540pub struct Disc {541 /// The centre's abscissa in the track's units, the centre of the ring.542 pub x: f64,543 /// The centre's ordinate in the track's units.544 pub y: f64,545 /// The radius no curve leaves, in the track's units: `rho + d` with `rho` the centre circle's radius and `d = r abs(p)` the outermost seat, which is `(a - b)/b + abs(p)` wheel radii inside and `(a + b)/b + abs(p)` outside.546 pub radius: f64,547 /// The radius no curve enters, in the track's units: the least of `abs(rho - d)` over the seats, which is not `rho` less the outermost seat when the seats straddle `rho`, and `rho` itself when no pencil is seated.548 pub hole: f64,549}550551/// The shape between the walls on a raster, with its numbers.552#[derive(Clone, Debug, PartialEq)]553pub struct Cover {554 /// The raster row by row, four codes: zero outside the disc, one the flood from the raster's edge, two the flood from the centre, three the shape, the walls and their pockets included.555 pub mask: Vec<u8>,556 /// The raster's side in pixels.557 pub side: usize,558 /// The shape's share of the disc's pixels, the walls counted in.559 pub covered: f64,560 /// The centre flood's share of the disc's pixels.561 pub hole: f64,562 /// The share of the disc's pixels the polylines themselves mark, the boundary error `covered` carries and loses as the raster grows.563 pub wall: f64,564 /// The mean signed winding number of the disc's pixel centres, read by scanline off the polylines.565 pub winding: f64,566 /// The closed form `winding` converges to: the distinct curves' signed areas summed and divided by the disc's area.567 pub areas: f64,568}569570/// The signed area one pencil's closed trochoid sweeps over the whole track, counterclockwise positive and counted with multiplicity, so it is the winding number integrated over the plane: `pi b rho (rho - d^2/r)` inside and `pi b rho (rho + d^2/r)` outside, with `rho` the centre circle's radius `R -+ r`, `d = r |p|` the seat's distance from the wheel's centre and `R/r = a/b` in lowest terms. Green's theorem on `z(t) = rho e^(i t) + p r e^(-+ i (rho/r) t)` gives it, and the cross terms carry `e^(-+ i a t / b)` over `b` centre turns and integrate to zero. No hypothesis on the seat: loops are counted with their sign. `None` off a circle track, where the roulette need not close.571pub fn signed_area(track: &Track, pencil: &Pencil) -> Option<f64> {572 let (_, rho) = ring(track).ok()?;573 let r = track.wheel;574 let d = r * pencil.x.hypot(pencil.y);575 let b = track.ratio.1 as f64;576 Some(PI * b * rho * (rho + track.side * d * d / r))577}578579/// The disc a circle roulette sits in: the wheel's centre turns on a circle of radius `rho`, and a seat `d` from the wheel's centre puts the pencil at `|z|^2 = rho^2 + d^2 + 2 rho d cos(a t / b -+ arg p)`, whose phase runs over `a` full turns, so that curve lies in the closed annulus from `abs(rho - d)` to `rho + d` and attains both bounds. The whole roulette therefore never leaves the disc of radius `rho + max d` and enters no disc of radius under `min abs(rho - d)`, the least over the seats and not the outermost seat's own, since seats on both sides of `rho` each keep their own inner radius. Refuses a line or a polygon track, whose roulette need not close and has no wall.580pub fn disc(track: &Track, pencils: &[Pencil]) -> Result<Disc> {581 let ((x, y), rho) = ring(track)?;582 let reach = |p: &Pencil| track.wheel * p.x.hypot(p.y);583 let far = pencils.iter().map(reach).fold(0.0_f64, f64::max);584 let near = pencils585 .iter()586 .map(|p| (rho - reach(p)).abs())587 .fold(f64::MAX, f64::min);588 Ok(Disc {589 x,590 y,591 radius: rho + far,592 hole: if pencils.is_empty() { rho } else { near },593 })594}595596/// The shape between the walls of a circle roulette, on a raster of `side` by `side` pixels over the disc, row zero at the top and the ordinate falling down the rows. Every distinct curve under the coincidence law is drawn once as a polyline of at least `samples` points, and of enough points that consecutive points land in one pixel or in two of the eight that touch, so the polylines make a wall no four-connected flood crosses. One flood starts from every pixel of the raster's edge, the fluid poured from outside; one starts from the centre pixel, the fluid poured at the centre, and is empty when the centre is a wall or the outside already reached it; the shape is the rest of the disc, pockets included. `covered` is the shape's share of the disc's pixels, the wall's own pixels counted in and reported apart as `wall`, and `hole` is the centre flood's share. `winding` is the mean signed winding number of the disc's pixel centres, read off crossings of the same polylines by scanline and never off a flood, and `areas` is the closed form it converges to, the distinct curves' `signed_area` summed over the disc's area: the pair checks the polylines and the raster against Green's theorem and never the floods, which are guarded instead by the sample spacing of at most half a pixel, which makes the wall eight-connected and a four-connected flood unable to cross it. Every share carries a boundary error of the order of the polylines' length times the pixel side over the disc's area.597pub fn cover(598 track: &Track,599 pencils: &[Pencil],600 exact: bool,601 samples: usize,602 side: usize,603) -> Result<Cover> {604 if !(16..=RASTER_CAP).contains(&side) {605 return value_error(format!(606 "the raster side must be between 16 and {RASTER_CAP}."607 ));608 }609 if samples < 2 {610 return value_error("a cover needs at least two samples.");611 }612 let (_, rho) = ring(track)?;613 let bounds = disc(track, pencils)?;614 let (n, r) = (side, track.wheel);615 let h = 2.0 * bounds.radius / n as f64;616 let left = bounds.x - bounds.radius;617 let top = bounds.y + bounds.radius;618 let drawn = representatives(track, pencils, exact);619 let mut steps = Vec::with_capacity(drawn.len());620 for &k in &drawn {621 let d = r * pencils[k].x.hypot(pencils[k].y);622 let arc = TAU * track.ratio.1 as f64 * rho * (1.0 + d / r);623 steps.push(samples.max((2.0 * arc / h).ceil().max(2.0) as usize));624 }625 let work: usize = steps.iter().sum();626 if work > POINT_CAP {627 return value_error(format!(628 "{work} samples is past the cap of {POINT_CAP}: a coarser raster or fewer curves."629 ));630 }631 let mut wall = vec![false; n * n];632 let mut rows: Vec<Vec<(f64, i32)>> = vec![Vec::new(); n];633 for (&k, &count) in drawn.iter().zip(&steps) {634 let mut last = point(track, &pencils[k], 0.0);635 mark(&mut wall, n, left, top, h, last);636 for step in 1..=count {637 let s = track.total * step as f64 / count as f64;638 let next = point(track, &pencils[k], s);639 mark(&mut wall, n, left, top, h, next);640 scan(&mut rows, n, top, h, last, next);641 last = next;642 }643 }644 let mut code = vec![0u8; n * n];645 let mut stack = Vec::new();646 for j in 0..n {647 for at in [j, (n - 1) * n + j, j * n, j * n + n - 1] {648 if !wall[at] && code[at] == 0 {649 code[at] = 1;650 stack.push(at);651 }652 }653 }654 flood(&mut code, &wall, n, &mut stack);655 let middle = (n / 2) * n + n / 2;656 if !wall[middle] && code[middle] == 0 {657 code[middle] = 2;658 stack.push(middle);659 flood(&mut code, &wall, n, &mut stack);660 }661 let (mut total, mut shape, mut hollow, mut edge, mut turns) =662 (0usize, 0usize, 0usize, 0usize, 0i64);663 for (row, line) in rows.iter_mut().enumerate() {664 let y = top - (row as f64 + 0.5) * h;665 line.sort_by(|a, b| b.0.total_cmp(&a.0));666 let (mut seen, mut winding) = (0usize, 0i64);667 for column in (0..n).rev() {668 let x = left + (column as f64 + 0.5) * h;669 while seen < line.len() && line[seen].0 > x {670 winding += i64::from(line[seen].1);671 seen += 1;672 }673 let at = row * n + column;674 if (x - bounds.x).hypot(y - bounds.y) > bounds.radius {675 code[at] = 0;676 continue;677 }678 total += 1;679 turns += winding;680 if wall[at] {681 edge += 1;682 }683 match code[at] {684 1 => {}685 2 => hollow += 1,686 _ => {687 code[at] = 3;688 shape += 1;689 }690 }691 }692 }693 let count = total.max(1) as f64;694 let sum: f64 = drawn695 .iter()696 .filter_map(|&k| signed_area(track, &pencils[k]))697 .sum();698 Ok(Cover {699 mask: code,700 side: n,701 covered: shape as f64 / count,702 hole: hollow as f64 / count,703 wall: edge as f64 / count,704 winding: turns as f64 / count,705 areas: sum / (PI * bounds.radius * bounds.radius),706 })707}708709fn ring(track: &Track) -> Result<((f64, f64), f64)> {710 match track.pieces.first() {711 Some(&Piece::Turn { about, radius, .. }) if track.kind == "in" || track.kind == "out" => {712 Ok((about, radius))713 }714 _ => value_error("the cover needs a circle track: the wall of a line or a polygon roulette need not close."),715 }716}717718fn mark(wall: &mut [bool], n: usize, left: f64, top: f64, h: f64, p: (f64, f64)) {719 let column = (((p.0 - left) / h) as isize).clamp(0, n as isize - 1) as usize;720 let row = (((top - p.1) / h) as isize).clamp(0, n as isize - 1) as usize;721 wall[row * n + column] = true;722}723724fn scan(rows: &mut [Vec<(f64, i32)>], n: usize, top: f64, h: f64, p: (f64, f64), q: (f64, f64)) {725 if p.1 == q.1 {726 return;727 }728 let (sign, low, high) = if q.1 > p.1 {729 (1, p.1, q.1)730 } else {731 (-1, q.1, p.1)732 };733 let last = n as f64 - 1.0;734 let first = ((top - high) / h - 1.5).clamp(0.0, last) as usize;735 let stop = ((top - low) / h + 0.5).clamp(0.0, last) as usize;736 for (step, line) in rows[first..=stop].iter_mut().enumerate() {737 let y = top - ((first + step) as f64 + 0.5) * h;738 if y < low || y >= high {739 continue;740 }741 line.push((p.0 + (q.0 - p.0) * (y - p.1) / (q.1 - p.1), sign));742 }743}744745fn flood(code: &mut [u8], wall: &[bool], n: usize, stack: &mut Vec<usize>) {746 while let Some(at) = stack.pop() {747 let tint = code[at];748 let (row, column) = (at / n, at % n);749 let step = |to: usize, code: &mut [u8], stack: &mut Vec<usize>| {750 if !wall[to] && code[to] == 0 {751 code[to] = tint;752 stack.push(to);753 }754 };755 if row > 0 {756 step(at - n, code, stack);757 }758 if row + 1 < n {759 step(at + n, code, stack);760 }761 if column > 0 {762 step(at - 1, code, stack);763 }764 if column + 1 < n {765 step(at + 1, code, stack);766 }767 }768}769770#[cfg(test)]771mod tests {772 use super::*;773774 fn ring() -> Vec<u8> {775 vec![1, 1, 1, 1, 0, 1, 1, 1, 1]776 }777778 #[test]779 fn the_carpet_seats_eight_pencils_and_the_law_folds_them_into_curves_and_shifts() {780 let fills = pencils(&ring(), 3, 3, "fill", 0.9, 0.0, 1).unwrap();781 assert_eq!(fills.len(), 8);782 assert_eq!(seats(&fills).fills, 8);783 assert_eq!(784 pencils(&ring(), 3, 3, "void", 0.9, 0.0, 1).unwrap()[0].seat,785 (0, 0)786 );787 assert_eq!(788 pencils(&ring(), 3, 3, "corners", 0.9, 0.0, 1)789 .unwrap()790 .len(),791 16792 );793 let curves = |kind: &str, ring: usize, wheel: usize| {794 distinct(&track(kind, ring, wheel, 4, 1).unwrap(), &fills, true)795 };796 assert_eq!(curves("in", 7, 4), 2);797 assert_eq!(curves("in", 7, 2), 4);798 assert_eq!(curves("in", 7, 3), 8);799 assert_eq!(curves("out", 5, 8), 2);800 assert_eq!(curves("in", 7, 6), 4);801 let shapes_up_to_a_shift = curves("line", 7, 3);802 assert_eq!(shapes_up_to_a_shift, 2);803 let jittered = pencils(&ring(), 3, 3, "fill", 0.9, 0.3, 1).unwrap();804 assert_eq!(805 distinct(&track("in", 7, 4, 4, 1).unwrap(), &jittered, false),806 8807 );808 }809810 #[test]811 fn a_pencil_traces_the_textbook_hypotrochoid() {812 let path = track("in", 7, 3, 4, 1).unwrap();813 let pencil = Pencil {814 x: 0.5,815 y: 0.0,816 seat: (0, 0),817 kind: Kind::Fill,818 };819 let (big, r, d) = (7.0, 3.0, 1.5);820 for k in 0..12 {821 let theta = TAU * k as f64 / 7.0;822 let s = (big - r) * theta;823 let (x, y) = point(&path, &pencil, s);824 let want = (825 (big - r) * theta.cos() + d * ((big - r) / r * theta).cos(),826 (big - r) * theta.sin() - d * ((big - r) / r * theta).sin(),827 );828 assert!((x - want.0).abs() < 1e-9 && (y - want.1).abs() < 1e-9);829 }830 assert_eq!(path.ratio, (7, 3));831 assert_eq!((path.orbits, path.fold), (3, 7));832 assert!((path.total - TAU * 4.0 * 3.0).abs() < 1e-9);833 }834835 #[test]836 fn the_circle_closes_and_the_line_and_polygon_run_on() {837 let fills = pencils(&ring(), 3, 3, "fill", 0.9, 0.0, 1).unwrap();838 let closed = trace(&track("in", 7, 3, 4, 1).unwrap(), &fills, 500).unwrap();839 assert!((closed[0] - closed[998]).abs() < 1e-5 && (closed[1] - closed[999]).abs() < 1e-5);840 let open = trace(&track("line", 7, 3, 4, 2).unwrap(), &fills, 500).unwrap();841 assert!((open[0] - open[998]).abs() > 1.0);842 let square = track("polyout", 7, 3, 4, 1).unwrap();843 assert_eq!(square.pieces.len(), 8);844 assert!((square.total - (8.0 * 7.0 / 2f64.sqrt() + TAU * 3.0)).abs() < 1e-9);845 let inner = track("polyin", 7, 3, 4, 2).unwrap();846 let (x0, y0) = pose(&inner, 0.0);847 let (x1, y1) = pose(&inner, inner.total / 2.0);848 assert!((x0 - x1).abs() < 1e-9 && (y0 - y1).abs() < 1e-9);849 }850851 #[test]852 fn the_nodes_count_only_inside_the_seat_window() {853 let fills = pencils(&ring(), 3, 3, "fill", 0.9, 0.0, 1).unwrap();854 let seven = track("in", 7, 3, 4, 1).unwrap();855 assert_eq!(distinct(&seven, &fills, true), 8);856 assert_eq!(nodes(&seven, &fills, true), Some(1288));857 let four = track("in", 7, 4, 4, 1).unwrap();858 assert_eq!(nodes(&four, &fills, true), Some(98));859 let past = pencils(&ring(), 3, 3, "fill", 1.35, 0.0, 1).unwrap();860 assert_eq!(nodes(&four, &past, true), None);861 assert_eq!(nodes(&track("in", 7, 6, 4, 1).unwrap(), &fills, true), None);862 assert_eq!(863 nodes(&track("out", 5, 8, 4, 1).unwrap(), &fills, true),864 Some(150)865 );866 assert_eq!(867 nodes(&track("line", 7, 3, 4, 1).unwrap(), &fills, true),868 None869 );870 let far = pencils(&ring(), 3, 3, "corners", 1.2, 0.0, 1).unwrap();871 assert_eq!(nodes(&seven, &far, true), None);872 let centre = pencils(&ring(), 3, 3, "void", 0.9, 0.0, 1).unwrap();873 assert_eq!(centre[0].seat, (0, 0));874 assert_eq!(nodes(&seven, ¢re, true), None);875 let empty = pencils(&[1u8; 9], 3, 3, "void", 0.9, 0.0, 1).unwrap();876 assert_eq!(nodes(&seven, &empty, true), Some(0));877 }878879 #[test]880 fn the_simple_curve_holds_its_whole_inside_in_the_hole() {881 let d = 0.9;882 for (kind, rho) in [("in", 2.0), ("out", 4.0)] {883 let path = track(kind, 3, 1, 4, 1).unwrap();884 let pencil = Pencil {885 x: d,886 y: 0.0,887 seat: (1, 0),888 kind: Kind::Fill,889 };890 let bounds = disc(&path, &[pencil]).unwrap();891 assert!((bounds.radius - (rho + d)).abs() < 1e-12);892 assert!((bounds.hole - (rho - d)).abs() < 1e-12);893 let form = rho * (rho + path.side * d * d) / (rho + d).powi(2);894 let side = 512;895 let slack = 2.0 * TAU * rho * (1.0 + d) * (2.0 * (rho + d) / side as f64)896 / (PI * (rho + d).powi(2));897 let cover = cover(&path, &[pencil], true, 2, side).unwrap();898 assert_eq!(cover.mask.len(), side * side);899 assert!((cover.areas - form).abs() < 1e-12);900 assert!(cover.covered < slack);901 assert!((cover.hole - form).abs() < slack);902 assert!((cover.winding - form).abs() < slack);903 }904 }905906 #[test]907 fn the_disc_reads_the_innermost_curve_and_not_the_farthest_seat() {908 let corners = pencils(&ring(), 3, 3, "corners", 1.2, 0.0, 1).unwrap();909 let straddle = track("in", 3, 2, 4, 1).unwrap();910 let bounds = disc(&straddle, &corners).unwrap();911 assert_eq!(912 format!("{:.6} {:.6}", bounds.radius, bounds.hole),913 "3.400000 0.200000"914 );915 let inscribed = (bounds.hole / bounds.radius).powi(2);916 let cover = cover(&straddle, &corners, true, 2, 512).unwrap();917 assert!(cover.hole >= inscribed && cover.hole < 2.0 * inscribed);918 let out = track("in", 5, 4, 4, 1).unwrap();919 let far = Pencil {920 x: 0.9,921 y: 0.0,922 seat: (1, 0),923 kind: Kind::Fill,924 };925 let past = disc(&out, &[far]).unwrap();926 assert_eq!(927 format!("{:.6} {:.6}", past.radius, past.hole),928 "4.600000 2.600000"929 );930 let empty = disc(&straddle, &[]).unwrap();931 assert_eq!(932 format!("{:.6} {:.6}", empty.radius, empty.hole),933 "1.000000 1.000000"934 );935 }936937 #[test]938 fn the_mean_winding_is_the_signed_areas_over_the_disc() {939 let fills = pencils(&ring(), 3, 3, "fill", 0.9, 0.0, 1).unwrap();940 let seven = track("in", 7, 3, 4, 1).unwrap();941 let drawn = representatives(&seven, &fills, true);942 assert_eq!(drawn.len(), 8);943 let bounds = disc(&seven, &fills).unwrap();944 let area = PI * bounds.radius * bounds.radius;945 let want: f64 = drawn946 .iter()947 .map(|&k| signed_area(&seven, &fills[k]).unwrap())948 .sum::<f64>()949 / area;950 let side = 384;951 let arcs: f64 = drawn952 .iter()953 .map(|&k| TAU * 3.0 * 4.0 * (1.0 + fills[k].x.hypot(fills[k].y)))954 .sum();955 let slack = 2.0 * arcs * (2.0 * bounds.radius / side as f64) / area;956 let cover = cover(&seven, &fills, true, 2, side).unwrap();957 assert!((cover.areas - want).abs() < 1e-12);958 assert!((cover.winding - want).abs() < slack);959 assert!(cover.covered > 0.0 && cover.hole > 0.0);960 assert!((cover.covered + cover.hole) < 1.0);961 }962963 #[test]964 fn the_cover_refuses_what_has_no_wall() {965 let fills = pencils(&ring(), 3, 3, "fill", 0.9, 0.0, 1).unwrap();966 let straight = track("line", 7, 3, 4, 1).unwrap();967 assert!(disc(&straight, &fills).is_err());968 assert!(signed_area(&straight, &fills[0]).is_none());969 assert!(cover(&track("polyin", 7, 3, 4, 1).unwrap(), &fills, true, 2, 64).is_err());970 let seven = track("in", 7, 3, 4, 1).unwrap();971 assert!(cover(&seven, &fills, true, 2, 8).is_err());972 assert!(cover(&seven, &fills, true, 1, 64).is_err());973 }974975 #[test]976 fn a_wheel_too_big_for_its_track_is_refused() {977 assert!(track("in", 3, 3, 4, 1).is_err());978 assert!(track("polyin", 4, 3, 3, 1).is_err());979 assert!(track("orbit", 7, 3, 4, 1).is_err());980 assert!(pencils(&ring(), 3, 3, "edges", 0.9, 0.0, 1).is_err());981 assert!(pencils(&ring(), 3, 2, "fill", 0.9, 0.0, 1).is_err());982 }983}