spirograph.rs
42.6 kB · rust · 1058 lines
1use crate::core::error::{value_error, Result};2use crate::core::Rng;3use crate::num::factor::gcd;4use serde::{Deserialize, Serialize};5use std::collections::HashSet;6use std::f64::consts::{PI, TAU};78/// The largest radius of a ring or a wheel.9pub const RADIUS_CAP: usize = 96;10/// The fewest and the most sides of a polygon track.11pub const SIDES: (usize, usize) = (3, 12);12/// The most laps of a line or a polygon.13pub const LAPS_CAP: usize = 24;14/// The most pencils a wheel seats.15pub const PENCIL_CAP: usize = 4096;16/// The most points one trace returns.17pub const POINT_CAP: usize = 4_000_000;18const REACH: (f64, f64) = (0.05, 2.0);19const RING: usize = 360;2021// THE PENCILS2223/// What a pencil sits on: a filled cell, an empty cell, or a corner of a filled cell.24#[derive(Clone, Copy, Debug, PartialEq, Eq, Serialize, Deserialize)]25pub enum Kind {26 /// The centre of a filled cell.27 Fill,28 /// The centre of an empty cell.29 Void,30 /// A corner of a filled cell.31 Corner,32}3334/// 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.35#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]36pub struct Pencil {37 /// The seat's abscissa, in wheel radii.38 pub x: f64,39 /// The seat's ordinate, up the page, in wheel radii.40 pub y: f64,41 /// The exact seat, twice the cell coordinates from the tile centre, before any jitter.42 pub seat: (i64, i64),43 /// What the pencil sits on.44 pub kind: Kind,45}4647/// The mass of a byte grid taken as a wheel: how many pencils of each kind it seats.48#[derive(Clone, Copy, Debug, Default, PartialEq, Eq, Serialize, Deserialize)]49pub struct Seats {50 /// The pencils on filled cells.51 pub fills: usize,52 /// The pencils on empty cells.53 pub voids: usize,54 /// The pencils on corners.55 pub corners: usize,56}5758/// 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.59///60/// # Errors61///62/// Errors on a grid that is not width by height, a reach or jitter out of range, or a mode that is not fill, void, both or corners.63pub fn pencils(64 types: &[u8],65 width: usize,66 height: usize,67 mode: &str,68 reach: f64,69 jitter: f64,70 seed: u32,71) -> Result<Vec<Pencil>> {72 if width == 0 || height == 0 || types.len() != width * height {73 return value_error("the wheel must be a grid of width by height bytes.");74 }75 if !(REACH.0..=REACH.1).contains(&reach) {76 return value_error(format!(77 "the reach must be between {} and {} wheel radii.",78 REACH.0, REACH.179 ));80 }81 if !(0.0..=1.0).contains(&jitter) {82 return value_error("the jitter must be between 0 and 1 cells.");83 }84 if !["fill", "void", "both", "corners"].contains(&mode) {85 return value_error("the pencils sit on fill, void, both or corners.");86 }87 let (w, h) = (width as i64, height as i64);88 let mut seen = HashSet::new();89 let mut out = Vec::new();90 for i in 0..height {91 for j in 0..width {92 let on = types[i * width + j] != 0;93 let (i, j) = (i as i64, j as i64);94 if mode == "corners" {95 if !on {96 continue;97 }98 for (di, dj) in [(0, 0), (0, 1), (1, 0), (1, 1)] {99 let seat = (2 * (j + dj) - w, h - 2 * (i + di));100 if seen.insert(seat) {101 out.push(seat_pencil(seat, Kind::Corner));102 }103 }104 } else if (on && mode != "void") || (!on && mode != "fill") {105 let kind = if on { Kind::Fill } else { Kind::Void };106 out.push(seat_pencil((2 * j + 1 - w, h - 2 * i - 1), kind));107 }108 }109 }110 if out.len() > PENCIL_CAP {111 return value_error(format!(112 "{} pencils is past the cap of {PENCIL_CAP}: take a lower level or a smaller tile.",113 out.len()114 ));115 }116 let unit = reach / (width as f64 / 2.0).hypot(height as f64 / 2.0);117 let mut rng = Rng::new(u64::from(seed));118 for pencil in &mut out {119 let (dx, dy) = if jitter > 0.0 {120 ((rng.unit() - 0.5) * jitter, (rng.unit() - 0.5) * jitter)121 } else {122 (0.0, 0.0)123 };124 pencil.x = (pencil.seat.0 as f64 / 2.0 + dx) * unit;125 pencil.y = (pencil.seat.1 as f64 / 2.0 + dy) * unit;126 }127 Ok(out)128}129130/// The side of one cell in wheel radii at a reach, the number the page needs to draw the tile on the wheel.131pub fn cell(width: usize, height: usize, reach: f64) -> f64 {132 reach / (width as f64 / 2.0).hypot(height as f64 / 2.0)133}134135/// Counts the pencils by kind.136pub fn seats(pencils: &[Pencil]) -> Seats {137 let mut out = Seats::default();138 for pencil in pencils {139 match pencil.kind {140 Kind::Fill => out.fills += 1,141 Kind::Void => out.voids += 1,142 Kind::Corner => out.corners += 1,143 }144 }145 out146}147148fn seat_pencil(seat: (i64, i64), kind: Kind) -> Pencil {149 Pencil {150 x: 0.0,151 y: 0.0,152 seat,153 kind,154 }155}156157// THE TRACK158159/// One piece of the centre path: a straight run, or a turn about a point.160#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]161pub enum Piece {162 /// A straight run from one point to another.163 Run {164 /// Where the run starts.165 from: (f64, f64),166 /// Where the run ends.167 to: (f64, f64),168 },169 /// A turn about a point at a radius, counterclockwise from one angle to another.170 Turn {171 /// The point turned about.172 about: (f64, f64),173 /// The radius of the turn.174 radius: f64,175 /// The angle the turn starts at.176 from: f64,177 /// The angle the turn ends at, past the start.178 to: f64,179 },180}181182impl Piece {183 fn len(&self) -> f64 {184 match *self {185 Piece::Run { from, to } => (to.0 - from.0).hypot(to.1 - from.1),186 Piece::Turn {187 radius, from, to, ..188 } => radius * (to - from),189 }190 }191192 fn at(&self, s: f64) -> (f64, f64) {193 match *self {194 Piece::Run { from, to } => {195 let len = self.len();196 let t = if len > 0.0 {197 (s / len).clamp(0.0, 1.0)198 } else {199 0.0200 };201 (from.0 + (to.0 - from.0) * t, from.1 + (to.1 - from.1) * t)202 }203 Piece::Turn {204 about,205 radius,206 from,207 to,208 } => {209 let angle = (from + s / radius).min(to);210 (211 about.0 + radius * angle.cos(),212 about.1 + radius * angle.sin(),213 )214 }215 }216 }217218 fn bounds(&self) -> [f64; 4] {219 match *self {220 Piece::Run { from, to } => [221 from.0.min(to.0),222 from.1.min(to.1),223 from.0.max(to.0),224 from.1.max(to.1),225 ],226 Piece::Turn { about, radius, .. } => [227 about.0 - radius,228 about.1 - radius,229 about.0 + radius,230 about.1 + radius,231 ],232 }233 }234}235236/// 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.237#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]238pub struct Track {239 /// The kind: `line`, `in`, `out`, `polyin` or `polyout`.240 pub kind: String,241 /// The wheel's radius.242 pub wheel: f64,243 /// The pieces of the centre path, in order.244 pub pieces: Vec<Piece>,245 /// The length of the centre path.246 pub total: f64,247 /// Minus one inside the track, plus one outside it.248 pub side: f64,249 /// The track itself as a polyline, closed when `closed` says so.250 pub outline: Vec<(f64, f64)>,251 /// Whether the outline closes on itself.252 pub closed: bool,253 /// The ring's radius over the wheel's in lowest terms on a circle, zero over zero elsewhere.254 pub ratio: (usize, usize),255 /// How many times the centre goes round: the ratio's denominator on a circle, the laps on a polygon, one on a line.256 pub orbits: usize,257 /// The rotation order of the whole picture: the ratio's numerator on a circle, none elsewhere.258 pub fold: usize,259}260261/// 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.262///263/// # Errors264///265/// Errors on a radius or lap count out of range, a wheel that does not fit its ring, or a kind that is not line, in, out, polyin or polyout.266pub fn track(kind: &str, ring: usize, wheel: usize, sides: usize, laps: usize) -> Result<Track> {267 if !(1..=RADIUS_CAP).contains(&wheel) || !(1..=RADIUS_CAP).contains(&ring) {268 return value_error(format!(269 "the ring and the wheel must have radius between 1 and {RADIUS_CAP}."270 ));271 }272 if !(1..=LAPS_CAP).contains(&laps) {273 return value_error(format!("the laps must be between 1 and {LAPS_CAP}."));274 }275 let r = wheel as f64;276 let big = ring as f64;277 let mut out = Track {278 kind: kind.to_string(),279 wheel: r,280 pieces: Vec::new(),281 total: 0.0,282 side: -1.0,283 outline: Vec::new(),284 closed: false,285 ratio: (0, 0),286 orbits: laps,287 fold: 0,288 };289 match kind {290 "line" => {291 let run = laps as f64 * TAU * r;292 out.pieces.push(Piece::Run {293 from: (0.0, r),294 to: (run, r),295 });296 out.outline = vec![(-r, 0.0), (run + r, 0.0)];297 out.orbits = 1;298 }299 "in" | "out" => {300 let inside = kind == "in";301 if inside && ring <= wheel {302 return value_error("the wheel must be smaller than the ring it rolls inside.");303 }304 let g = gcd(ring as u128, wheel as u128) as usize;305 let (a, b) = (ring / g, wheel / g);306 let rho = if inside { big - r } else { big + r };307 out.side = if inside { -1.0 } else { 1.0 };308 out.pieces.push(Piece::Turn {309 about: (0.0, 0.0),310 radius: rho,311 from: 0.0,312 to: TAU * b as f64,313 });314 out.outline = (0..RING)315 .map(|k| {316 let angle = TAU * k as f64 / RING as f64;317 (big * angle.cos(), big * angle.sin())318 })319 .collect();320 out.closed = true;321 out.ratio = (a, b);322 out.orbits = b;323 out.fold = a;324 }325 "polyin" | "polyout" => {326 if !(SIDES.0..=SIDES.1).contains(&sides) {327 return value_error(format!(328 "the polygon must have between {} and {} sides.",329 SIDES.0, SIDES.1330 ));331 }332 let corners = polygon(big, sides);333 out.outline = corners.clone();334 out.closed = true;335 if kind == "polyin" {336 let inset = big * (PI / sides as f64).cos() - r;337 if inset <= 0.0 {338 return value_error("the wheel does not fit inside the polygon.");339 }340 let inner = polygon(inset / (PI / sides as f64).cos(), sides);341 for _ in 0..laps {342 for k in 0..sides {343 out.pieces.push(Piece::Run {344 from: inner[k],345 to: inner[(k + 1) % sides],346 });347 }348 }349 } else {350 out.side = 1.0;351 let normal = |k: usize| {352 let (p, q) = (corners[k], corners[(k + 1) % sides]);353 let (dx, dy) = (q.0 - p.0, q.1 - p.1);354 let len = dx.hypot(dy);355 (dy / len, -dx / len)356 };357 for _ in 0..laps {358 for k in 0..sides {359 let (p, q) = (corners[k], corners[(k + 1) % sides]);360 let (m, next) = (normal(k), normal((k + 1) % sides));361 out.pieces.push(Piece::Run {362 from: (p.0 + r * m.0, p.1 + r * m.1),363 to: (q.0 + r * m.0, q.1 + r * m.1),364 });365 let from = m.1.atan2(m.0);366 let mut to = next.1.atan2(next.0);367 while to < from {368 to += TAU;369 }370 out.pieces.push(Piece::Turn {371 about: q,372 radius: r,373 from,374 to,375 });376 }377 }378 }379 }380 _ => return value_error("the track is line, in, out, polyin or polyout."),381 }382 out.total = out.pieces.iter().map(Piece::len).sum();383 Ok(out)384}385386fn polygon(radius: f64, sides: usize) -> Vec<(f64, f64)> {387 (0..sides)388 .map(|k| {389 let angle = PI / sides as f64 + TAU * k as f64 / sides as f64;390 (radius * angle.cos(), radius * angle.sin())391 })392 .collect()393}394395/// The wheel's centre after `s` of path length.396pub fn pose(track: &Track, s: f64) -> (f64, f64) {397 let mut left = s.clamp(0.0, track.total);398 let last = track.pieces.len() - 1;399 for (k, piece) in track.pieces.iter().enumerate() {400 let len = piece.len();401 if left <= len || k == last {402 return piece.at(left);403 }404 left -= len;405 }406 (0.0, 0.0)407}408409/// The wheel's turn after `s` of path length, in radians: `side` times `s` over the wheel's radius.410pub fn turn(track: &Track, s: f64) -> f64 {411 let angle = track.side * s / track.wheel;412 if angle == 0.0 {413 0.0414 } else {415 angle416 }417}418419/// Where a pencil is after `s` of path length: the centre plus the seat turned with the wheel.420pub fn point(track: &Track, pencil: &Pencil, s: f64) -> (f64, f64) {421 let (cx, cy) = pose(track, s);422 let phi = turn(track, s);423 let (c, sn) = (phi.cos(), phi.sin());424 let r = track.wheel;425 (426 cx + r * (pencil.x * c - pencil.y * sn),427 cy + r * (pencil.x * sn + pencil.y * c),428 )429}430431/// 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.432///433/// # Errors434///435/// Errors under two samples, or past the point cap.436pub fn trace(track: &Track, pencils: &[Pencil], samples: usize) -> Result<Vec<f32>> {437 if samples < 2 {438 return value_error("a trace needs at least two samples.");439 }440 if pencils.len() * samples > POINT_CAP {441 return value_error(format!(442 "{} points is past the cap of {POINT_CAP}: fewer samples or fewer pencils.",443 pencils.len() * samples444 ));445 }446 let mut out = Vec::with_capacity(pencils.len() * samples * 2);447 for pencil in pencils {448 for k in 0..samples {449 let s = track.total * k as f64 / (samples - 1) as f64;450 let (x, y) = point(track, pencil, s);451 out.push(x as f32);452 out.push(y as f32);453 }454 }455 Ok(out)456}457458/// 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.459pub fn distinct(track: &Track, pencils: &[Pencil], exact: bool) -> usize {460 representatives(track, pencils, exact).len()461}462463/// 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.464pub fn representatives(track: &Track, pencils: &[Pencil], exact: bool) -> Vec<usize> {465 if !exact {466 return (0..pencils.len()).collect();467 }468 let mut out = Vec::new();469 match track.kind.as_str() {470 "line" => {471 let mut seen = HashSet::new();472 for (k, p) in pencils.iter().enumerate() {473 if seen.insert(p.seat.0 * p.seat.0 + p.seat.1 * p.seat.1) {474 out.push(k);475 }476 }477 }478 "in" | "out" => {479 let order = gcd(track.ratio.1 as u128, 4) as usize;480 let quarter = |(u, v): (i64, i64)| (-v, u);481 let mut seen = HashSet::new();482 for (k, p) in pencils.iter().enumerate() {483 let mut least = p.seat;484 let mut seat = p.seat;485 for _ in 1..order {486 for _ in 0..4 / order {487 seat = quarter(seat);488 }489 least = least.min(seat);490 }491 if seen.insert(least) {492 out.push(k);493 }494 }495 }496 _ => out.extend(0..pencils.len()),497 }498 out499}500501/// 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.502pub fn frame(track: &Track, pencils: &[Pencil]) -> [f64; 4] {503 let reach = pencils504 .iter()505 .map(|p| p.x.hypot(p.y))506 .fold(1.0_f64, f64::max);507 let pad = track.wheel * reach;508 let mut box_ = [f64::MAX, f64::MAX, f64::MIN, f64::MIN];509 let mut take = |x: f64, y: f64| {510 box_[0] = box_[0].min(x);511 box_[1] = box_[1].min(y);512 box_[2] = box_[2].max(x);513 box_[3] = box_[3].max(y);514 };515 for piece in &track.pieces {516 let b = piece.bounds();517 take(b[0] - pad, b[1] - pad);518 take(b[2] + pad, b[3] + pad);519 }520 for &(x, y) in &track.outline {521 take(x, y);522 }523 box_524}525526// THE NODES527528/// 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.529pub fn nodes(track: &Track, pencils: &[Pencil], exact: bool) -> Option<u64> {530 if track.kind != "in" && track.kind != "out" {531 return None;532 }533 let (a, b) = (track.ratio.0 as u64, track.ratio.1 as u64);534 let window = (a as f64 / b as f64 + track.side).abs().min(1.0);535 let out = pencils.iter().any(|p| {536 let seat = p.x.hypot(p.y);537 !(seat > 0.0 && seat < window)538 });539 if out {540 return None;541 }542 let k = distinct(track, pencils, exact) as u64;543 Some(2 * a * b * (k * k.saturating_sub(1) / 2) + k * a * (b - 1))544}545546// THE COVER547548/// The largest raster side a cover rasters.549pub const RASTER_CAP: usize = 4096;550551/// The disc a circle roulette sits in, in the track's units.552#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]553pub struct Disc {554 /// The centre's abscissa in the track's units, the centre of the ring.555 pub x: f64,556 /// The centre's ordinate in the track's units.557 pub y: f64,558 /// 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.559 pub radius: f64,560 /// 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.561 pub hole: f64,562}563564/// The shape between the walls on a raster, with its numbers.565#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]566pub struct Cover {567 /// 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.568 pub mask: Vec<u8>,569 /// The raster's side in pixels.570 pub side: usize,571 /// The shape's share of the disc's pixels, the walls counted in.572 pub covered: f64,573 /// The centre flood's share of the disc's pixels.574 pub hole: f64,575 /// The share of the disc's pixels the polylines themselves mark, the boundary error `covered` carries and loses as the raster grows.576 pub wall: f64,577 /// The mean signed winding number of the disc's pixel centres, read by scanline off the polylines.578 pub winding: f64,579 /// The closed form `winding` converges to: the distinct curves' signed areas summed and divided by the disc's area.580 pub areas: f64,581}582583/// 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.584pub fn signed_area(track: &Track, pencil: &Pencil) -> Option<f64> {585 let (_, rho) = ring(track).ok()?;586 let r = track.wheel;587 let d = r * pencil.x.hypot(pencil.y);588 let b = track.ratio.1 as f64;589 Some(PI * b * rho * (rho + track.side * d * d / r))590}591592/// 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.593///594/// # Errors595///596/// Errors on a line or a polygon track.597pub fn disc(track: &Track, pencils: &[Pencil]) -> Result<Disc> {598 let ((x, y), rho) = ring(track)?;599 let reach = |p: &Pencil| track.wheel * p.x.hypot(p.y);600 let far = pencils.iter().map(reach).fold(0.0_f64, f64::max);601 let near = pencils602 .iter()603 .map(|p| (rho - reach(p)).abs())604 .fold(f64::MAX, f64::min);605 Ok(Disc {606 x,607 y,608 radius: rho + far,609 hole: if pencils.is_empty() { rho } else { near },610 })611}612613/// 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.614///615/// # Errors616///617/// Errors on a line or a polygon track, a raster side out of range, under two samples, or past the point cap.618pub fn cover(619 track: &Track,620 pencils: &[Pencil],621 exact: bool,622 samples: usize,623 side: usize,624) -> Result<Cover> {625 if !(16..=RASTER_CAP).contains(&side) {626 return value_error(format!(627 "the raster side must be between 16 and {RASTER_CAP}."628 ));629 }630 if samples < 2 {631 return value_error("a cover needs at least two samples.");632 }633 let (_, rho) = ring(track)?;634 let bounds = disc(track, pencils)?;635 let (n, r) = (side, track.wheel);636 let h = 2.0 * bounds.radius / n as f64;637 let left = bounds.x - bounds.radius;638 let top = bounds.y + bounds.radius;639 let drawn = representatives(track, pencils, exact);640 let mut steps = Vec::with_capacity(drawn.len());641 for &k in &drawn {642 let d = r * pencils[k].x.hypot(pencils[k].y);643 let arc = TAU * track.ratio.1 as f64 * rho * (1.0 + d / r);644 steps.push(samples.max((2.0 * arc / h).ceil().max(2.0) as usize));645 }646 let work: usize = steps.iter().sum();647 if work > POINT_CAP {648 return value_error(format!(649 "{work} samples is past the cap of {POINT_CAP}: a coarser raster or fewer curves."650 ));651 }652 let mut wall = vec![false; n * n];653 let mut rows: Vec<Vec<(f64, i32)>> = vec![Vec::new(); n];654 for (&k, &count) in drawn.iter().zip(&steps) {655 let mut last = point(track, &pencils[k], 0.0);656 mark(&mut wall, n, left, top, h, last);657 for step in 1..=count {658 let s = track.total * step as f64 / count as f64;659 let next = point(track, &pencils[k], s);660 mark(&mut wall, n, left, top, h, next);661 scan(&mut rows, n, top, h, last, next);662 last = next;663 }664 }665 let mut code = vec![0u8; n * n];666 let mut stack = Vec::new();667 for j in 0..n {668 for at in [j, (n - 1) * n + j, j * n, j * n + n - 1] {669 if !wall[at] && code[at] == 0 {670 code[at] = 1;671 stack.push(at);672 }673 }674 }675 flood(&mut code, &wall, n, &mut stack);676 let middle = (n / 2) * n + n / 2;677 if !wall[middle] && code[middle] == 0 {678 code[middle] = 2;679 stack.push(middle);680 flood(&mut code, &wall, n, &mut stack);681 }682 let (mut total, mut shape, mut hollow, mut edge, mut turns) =683 (0usize, 0usize, 0usize, 0usize, 0i64);684 for (row, line) in rows.iter_mut().enumerate() {685 let y = top - (row as f64 + 0.5) * h;686 line.sort_by(|a, b| b.0.total_cmp(&a.0));687 let (mut seen, mut winding) = (0usize, 0i64);688 for column in (0..n).rev() {689 let x = left + (column as f64 + 0.5) * h;690 while seen < line.len() && line[seen].0 > x {691 winding += i64::from(line[seen].1);692 seen += 1;693 }694 let at = row * n + column;695 if (x - bounds.x).hypot(y - bounds.y) > bounds.radius {696 code[at] = 0;697 continue;698 }699 total += 1;700 turns += winding;701 if wall[at] {702 edge += 1;703 }704 match code[at] {705 1 => {}706 2 => hollow += 1,707 _ => {708 code[at] = 3;709 shape += 1;710 }711 }712 }713 }714 let count = total.max(1) as f64;715 let sum: f64 = drawn716 .iter()717 .filter_map(|&k| signed_area(track, &pencils[k]))718 .sum();719 Ok(Cover {720 mask: code,721 side: n,722 covered: shape as f64 / count,723 hole: hollow as f64 / count,724 wall: edge as f64 / count,725 winding: turns as f64 / count,726 areas: sum / (PI * bounds.radius * bounds.radius),727 })728}729730fn ring(track: &Track) -> Result<((f64, f64), f64)> {731 match track.pieces.first() {732 Some(&Piece::Turn { about, radius, .. }) if track.kind == "in" || track.kind == "out" => {733 Ok((about, radius))734 }735 _ => value_error("the cover needs a circle track: the wall of a line or a polygon roulette need not close."),736 }737}738739fn mark(wall: &mut [bool], n: usize, left: f64, top: f64, h: f64, p: (f64, f64)) {740 let column = (((p.0 - left) / h) as isize).clamp(0, n as isize - 1) as usize;741 let row = (((top - p.1) / h) as isize).clamp(0, n as isize - 1) as usize;742 wall[row * n + column] = true;743}744745fn scan(rows: &mut [Vec<(f64, i32)>], n: usize, top: f64, h: f64, p: (f64, f64), q: (f64, f64)) {746 if p.1 == q.1 {747 return;748 }749 let (sign, low, high) = if q.1 > p.1 {750 (1, p.1, q.1)751 } else {752 (-1, q.1, p.1)753 };754 let last = n as f64 - 1.0;755 let first = ((top - high) / h - 1.5).clamp(0.0, last) as usize;756 let stop = ((top - low) / h + 0.5).clamp(0.0, last) as usize;757 for (step, line) in rows[first..=stop].iter_mut().enumerate() {758 let y = top - ((first + step) as f64 + 0.5) * h;759 if y < low || y >= high {760 continue;761 }762 line.push((p.0 + (q.0 - p.0) * (y - p.1) / (q.1 - p.1), sign));763 }764}765766fn flood(code: &mut [u8], wall: &[bool], n: usize, stack: &mut Vec<usize>) {767 while let Some(at) = stack.pop() {768 let tint = code[at];769 let (row, column) = (at / n, at % n);770 let step = |to: usize, code: &mut [u8], stack: &mut Vec<usize>| {771 if !wall[to] && code[to] == 0 {772 code[to] = tint;773 stack.push(to);774 }775 };776 if row > 0 {777 step(at - n, code, stack);778 }779 if row + 1 < n {780 step(at + n, code, stack);781 }782 if column > 0 {783 step(at - 1, code, stack);784 }785 if column + 1 < n {786 step(at + 1, code, stack);787 }788 }789}790791#[cfg(test)]792mod tests {793 use super::*;794795 fn ring() -> Vec<u8> {796 vec![1, 1, 1, 1, 0, 1, 1, 1, 1]797 }798799 #[test]800 fn the_carpet_seats_eight_pencils_and_the_law_folds_them_into_curves_and_shifts() {801 let fills = pencils(&ring(), 3, 3, "fill", 0.9, 0.0, 1).unwrap();802 assert_eq!(fills.len(), 8);803 assert_eq!(seats(&fills).fills, 8);804 let voids = pencils(&ring(), 3, 3, "void", 0.9, 0.0, 1).unwrap();805 assert_eq!(voids.len(), 1);806 assert_eq!(voids[0].seat, (0, 0));807 assert_eq!(808 seats(&pencils(&ring(), 3, 3, "corners", 0.9, 0.0, 1).unwrap()).corners,809 16810 );811 assert_eq!(812 pencils(&ring(), 3, 3, "corners", 0.9, 0.0, 1)813 .unwrap()814 .len(),815 16816 );817 let curves = |kind: &str, ring: usize, wheel: usize| {818 distinct(&track(kind, ring, wheel, 4, 1).unwrap(), &fills, true)819 };820 assert_eq!(curves("in", 7, 4), 2);821 assert_eq!(curves("in", 7, 2), 4);822 assert_eq!(curves("in", 7, 3), 8);823 assert_eq!(curves("out", 5, 8), 2);824 assert_eq!(curves("in", 7, 6), 4);825 let shapes_up_to_a_shift = curves("line", 7, 3);826 assert_eq!(shapes_up_to_a_shift, 2);827 let jittered = pencils(&ring(), 3, 3, "fill", 0.9, 0.3, 1).unwrap();828 assert_eq!(829 distinct(&track("in", 7, 4, 4, 1).unwrap(), &jittered, false),830 8831 );832 }833834 #[test]835 fn a_pencil_traces_the_textbook_hypotrochoid() {836 let path = track("in", 7, 3, 4, 1).unwrap();837 let pencil = Pencil {838 x: 0.5,839 y: 0.0,840 seat: (0, 0),841 kind: Kind::Fill,842 };843 let (big, r, d) = (7.0, 3.0, 1.5);844 for k in 0..12 {845 let theta = TAU * k as f64 / 7.0;846 let s = (big - r) * theta;847 let (x, y) = point(&path, &pencil, s);848 let want = (849 (big - r) * theta.cos() + d * ((big - r) / r * theta).cos(),850 (big - r) * theta.sin() - d * ((big - r) / r * theta).sin(),851 );852 assert!((x - want.0).abs() < 1e-9 && (y - want.1).abs() < 1e-9);853 }854 assert_eq!(path.ratio, (7, 3));855 assert_eq!((path.orbits, path.fold), (3, 7));856 assert!((path.total - TAU * 4.0 * 3.0).abs() < 1e-9);857 assert_eq!(858 format!("{:.6}", path.total / (TAU * path.wheel)),859 "4.000000"860 );861 let start = pose(&path, 0.0);862 assert_eq!(863 format!("{:.6} {:.6}", start.0, start.1),864 "4.000000 0.000000"865 );866 assert_eq!(format!("{:.6}", turn(&path, 0.0)), "0.000000");867 let end = pose(&path, path.total);868 assert_eq!(format!("{:.6}", end.0), "4.000000");869 assert_eq!(format!("{:.6}", turn(&path, path.total) / PI), "-8.000000");870 }871872 #[test]873 fn the_circle_closes_and_the_line_and_polygon_run_on() {874 let fills = pencils(&ring(), 3, 3, "fill", 0.9, 0.0, 1).unwrap();875 let closed = trace(&track("in", 7, 3, 4, 1).unwrap(), &fills, 500).unwrap();876 assert!((closed[0] - closed[998]).abs() < 1e-5 && (closed[1] - closed[999]).abs() < 1e-5);877 let open = trace(&track("line", 7, 3, 4, 2).unwrap(), &fills, 500).unwrap();878 assert!((open[0] - open[998]).abs() > 1.0);879 let square = track("polyout", 7, 3, 4, 1).unwrap();880 assert_eq!(square.pieces.len(), 8);881 assert!((square.total - (8.0 * 7.0 / 2f64.sqrt() + TAU * 3.0)).abs() < 1e-9);882 let inner = track("polyin", 7, 3, 4, 2).unwrap();883 let (x0, y0) = pose(&inner, 0.0);884 let (x1, y1) = pose(&inner, inner.total / 2.0);885 assert!((x0 - x1).abs() < 1e-9 && (y0 - y1).abs() < 1e-9);886 }887888 #[test]889 fn the_nodes_count_only_inside_the_seat_window() {890 let fills = pencils(&ring(), 3, 3, "fill", 0.9, 0.0, 1).unwrap();891 let seven = track("in", 7, 3, 4, 1).unwrap();892 assert_eq!(distinct(&seven, &fills, true), 8);893 assert_eq!(nodes(&seven, &fills, true), Some(1288));894 let four = track("in", 7, 4, 4, 1).unwrap();895 assert_eq!(nodes(&four, &fills, true), Some(98));896 let past = pencils(&ring(), 3, 3, "fill", 1.35, 0.0, 1).unwrap();897 assert_eq!(nodes(&four, &past, true), None);898 assert_eq!(nodes(&track("in", 7, 6, 4, 1).unwrap(), &fills, true), None);899 assert_eq!(900 nodes(&track("out", 5, 8, 4, 1).unwrap(), &fills, true),901 Some(150)902 );903 assert_eq!(904 nodes(&track("line", 7, 3, 4, 1).unwrap(), &fills, true),905 None906 );907 let far = pencils(&ring(), 3, 3, "corners", 1.2, 0.0, 1).unwrap();908 assert_eq!(nodes(&seven, &far, true), None);909 let centre = pencils(&ring(), 3, 3, "void", 0.9, 0.0, 1).unwrap();910 assert_eq!(centre[0].seat, (0, 0));911 assert_eq!(nodes(&seven, ¢re, true), None);912 let empty = pencils(&[1u8; 9], 3, 3, "void", 0.9, 0.0, 1).unwrap();913 assert_eq!(nodes(&seven, &empty, true), Some(0));914 }915916 #[test]917 fn the_simple_curve_holds_its_whole_inside_in_the_hole() {918 let d = 0.9;919 for (kind, rho) in [("in", 2.0), ("out", 4.0)] {920 let path = track(kind, 3, 1, 4, 1).unwrap();921 let pencil = Pencil {922 x: d,923 y: 0.0,924 seat: (1, 0),925 kind: Kind::Fill,926 };927 let bounds = disc(&path, &[pencil]).unwrap();928 assert!((bounds.radius - (rho + d)).abs() < 1e-12);929 assert!((bounds.hole - (rho - d)).abs() < 1e-12);930 let form = rho * (rho + path.side * d * d) / (rho + d).powi(2);931 let side = 512;932 let slack = 2.0 * TAU * rho * (1.0 + d) * (2.0 * (rho + d) / side as f64)933 / (PI * (rho + d).powi(2));934 let cover = cover(&path, &[pencil], true, 2, side).unwrap();935 assert_eq!(cover.mask.len(), side * side);936 assert!((cover.areas - form).abs() < 1e-12);937 assert!(cover.covered < slack);938 assert!((cover.hole - form).abs() < slack);939 assert!((cover.winding - form).abs() < slack);940 }941 }942943 #[test]944 fn the_disc_reads_the_innermost_curve_and_not_the_farthest_seat() {945 let corners = pencils(&ring(), 3, 3, "corners", 1.2, 0.0, 1).unwrap();946 let straddle = track("in", 3, 2, 4, 1).unwrap();947 let bounds = disc(&straddle, &corners).unwrap();948 assert_eq!(949 format!("{:.6} {:.6}", bounds.radius, bounds.hole),950 "3.400000 0.200000"951 );952 let inscribed = (bounds.hole / bounds.radius).powi(2);953 let cover = cover(&straddle, &corners, true, 2, 512).unwrap();954 assert!(cover.hole >= inscribed && cover.hole < 2.0 * inscribed);955 let out = track("in", 5, 4, 4, 1).unwrap();956 let far = Pencil {957 x: 0.9,958 y: 0.0,959 seat: (1, 0),960 kind: Kind::Fill,961 };962 let past = disc(&out, &[far]).unwrap();963 assert_eq!(964 format!("{:.6} {:.6}", past.radius, past.hole),965 "4.600000 2.600000"966 );967 let empty = disc(&straddle, &[]).unwrap();968 assert_eq!(969 format!("{:.6} {:.6}", empty.radius, empty.hole),970 "1.000000 1.000000"971 );972 }973974 #[test]975 fn the_mean_winding_is_the_signed_areas_over_the_disc() {976 let fills = pencils(&ring(), 3, 3, "fill", 0.9, 0.0, 1).unwrap();977 let seven = track("in", 7, 3, 4, 1).unwrap();978 let drawn = representatives(&seven, &fills, true);979 assert_eq!(drawn.len(), 8);980 let bounds = disc(&seven, &fills).unwrap();981 let area = PI * bounds.radius * bounds.radius;982 let want: f64 = drawn983 .iter()984 .map(|&k| signed_area(&seven, &fills[k]).unwrap())985 .sum::<f64>()986 / area;987 let side = 384;988 let arcs: f64 = drawn989 .iter()990 .map(|&k| TAU * 3.0 * 4.0 * (1.0 + fills[k].x.hypot(fills[k].y)))991 .sum();992 let slack = 2.0 * arcs * (2.0 * bounds.radius / side as f64) / area;993 let cover = cover(&seven, &fills, true, 2, side).unwrap();994 assert!((cover.areas - want).abs() < 1e-12);995 assert!((cover.winding - want).abs() < slack);996 assert!(cover.covered > 0.0 && cover.hole > 0.0);997 assert!((cover.covered + cover.hole) < 1.0);998 assert_eq!(999 format!("{:.6} {:.6}", bounds.radius, bounds.hole),1000 "5.800000 2.200000"1001 );1002 let walled = super::cover(&seven, &fills, true, 2, 256).unwrap();1003 assert_eq!(walled.mask.len(), 256 * 256);1004 assert_eq!(1005 format!(1006 "{:.6} {:.6} {:.6} {:.6} {:.6}",1007 walled.covered, walled.hole, walled.wall, walled.winding, walled.areas1008 ),1009 "0.814487 0.140825 0.269255 9.104609 9.103448"1010 );1011 let inside = walled.mask.iter().filter(|&&code| code > 0).count();1012 let shape = walled.mask.iter().filter(|&&code| code == 3).count();1013 assert_eq!(format!("{:.6}", shape as f64 / inside as f64), "0.814487");1014 }10151016 #[test]1017 fn one_curve_of_a_wheel_covers_the_area_its_closed_form_names() {1018 let one = pencils(&[1, 0, 0], 3, 1, "fill", 0.9, 0.0, 1).unwrap();1019 let d = one[0].x.hypot(one[0].y);1020 for (kind, rho, want) in [1021 ("in", 2.0, "0.014922 0.500622 0.508005 0.507814"),1022 ("out", 4.0, "0.016671 0.819733 0.828554 0.828445"),1023 ] {1024 let path = track(kind, 3, 1, 4, 1).unwrap();1025 let form = rho * (rho + path.side * d * d) / (rho + d).powi(2);1026 let drawn = cover(&path, &one, true, 2, 256).unwrap();1027 assert!((drawn.areas - form).abs() < 1e-12);1028 assert_eq!(1029 format!(1030 "{:.6} {:.6} {:.6} {:.6}",1031 drawn.covered, drawn.hole, drawn.winding, form1032 ),1033 want1034 );1035 }1036 }10371038 #[test]1039 fn the_cover_refuses_what_has_no_wall() {1040 let fills = pencils(&ring(), 3, 3, "fill", 0.9, 0.0, 1).unwrap();1041 let straight = track("line", 7, 3, 4, 1).unwrap();1042 assert!(disc(&straight, &fills).is_err());1043 assert!(signed_area(&straight, &fills[0]).is_none());1044 assert!(cover(&track("polyin", 7, 3, 4, 1).unwrap(), &fills, true, 2, 64).is_err());1045 let seven = track("in", 7, 3, 4, 1).unwrap();1046 assert!(cover(&seven, &fills, true, 2, 8).is_err());1047 assert!(cover(&seven, &fills, true, 1, 64).is_err());1048 }10491050 #[test]1051 fn a_wheel_too_big_for_its_track_is_refused() {1052 assert!(track("in", 3, 3, 4, 1).is_err());1053 assert!(track("polyin", 4, 3, 3, 1).is_err());1054 assert!(track("orbit", 7, 3, 4, 1).is_err());1055 assert!(pencils(&ring(), 3, 3, "edges", 0.9, 0.0, 1).is_err());1056 assert!(pencils(&ring(), 3, 2, "fill", 0.9, 0.0, 1).is_err());1057 }1058}