tourbillon.rs
21.1 kB · rust · 602 lines
1use crate::classics::primes;2use crate::factor::{factorize, gcd, mobius, squarefree};3use crate::prime::{is_prime, squares};4use crate::spin::Blend;5use mrlycore::errors::{value_error, Result};6use mrlycore::Rng;78const SCALE_CAP: usize = 199;9const SIZE_FLOOR: usize = 16;10const SIZE_CAP: usize = 1024;11const SHOWN: usize = 8;12const PEAKS: usize = 3;13const LADDER: usize = 1000;14const PLANE_CAP: usize = 64_000_000;15const PERIOD_CAP: usize = 360;16const TIGHT: f64 = 1e-9;1718// THE LAYERS1920/// One carpet of a stack: the odd scale it is drawn at, the weight the linear blends carry it at and the turn it takes about the centre, in degrees.21#[derive(Clone, Copy, Debug, PartialEq)]22pub struct Layer {23 /// The odd scale, the number of cells across the carpet.24 pub scale: usize,25 /// The weight, scaled so the magnitudes average one.26 pub weight: f64,27 /// The turn about the centre, in degrees.28 pub degrees: f64,29}3031/// The readings of a spun stack: its layers, the first eight scales and angles, the mean and RMS contrast over the disc, that contrast times the root of the layer count, the exact centre value, whether the blend carries the weights, the span the raster covers, the sites inside the disc and the brightest three.32#[derive(Clone, Debug, PartialEq)]33pub struct Stats {34 /// The count of layers in the stack.35 pub layers: usize,36 /// The first eight scales.37 pub scales: Vec<usize>,38 /// The first eight angles, in degrees.39 pub angles: Vec<f64>,40 /// The mean over the disc.41 pub mean: f64,42 /// The RMS contrast over the disc.43 pub rms: f64,44 /// The RMS contrast times the root of the layer count.45 pub faded: f64,46 /// The exact value at the centre, computed and never sampled.47 pub centre: f64,48 /// Whether the blend carries the weights.49 pub weighted: bool,50 /// The smallest value over the disc.51 pub low: f64,52 /// The largest value over the disc.53 pub high: f64,54 /// The count of sites inside the disc.55 pub inside: usize,56 /// The brightest three sites, each as its unit coordinates and its value.57 pub peaks: Vec<[f64; 3]>,58 /// The least whole number of increments that closes a quarter turn, none past the cap.59 pub period: Option<usize>,60 /// The count of distinct angle classes the layers fall in, read a quarter turn apart.61 pub classes: usize,62 /// The count of layer pairs sharing an angle class.63 pub pairs: usize,64}6566/// One angle of the quarter-turn lattice: the turn in degrees and the ninety a over q that names it.67#[derive(Clone, Copy, Debug, PartialEq)]68pub struct Eye {69 /// The turn in degrees.70 pub angle: f64,71 /// The numerator, ninety times a.72 pub numer: usize,73 /// The denominator q.74 pub denom: usize,75}7677fn golden() -> f64 {78 180.0 * (3.0 - 5f64.sqrt())79}8081// THE QUARTER TURN8283/// The angles a quarter turn shares with itself: ninety a over q for every q up to the cap and every a from zero to four q coprime to it, sorted by angle.84///85/// The odd carpet is unchanged by a quarter turn, so two layers whose angles agree modulo ninety stand on one exact lattice of nodes and the stack stops moving there.86///87/// ```88/// let list = mrlynum::tourbillon::eyes(2);89/// assert_eq!(list.len(), 9);90/// assert_eq!(list[1].angle, 45.0);91/// ```92pub fn eyes(qmax: usize) -> Vec<Eye> {93 let mut out = Vec::new();94 for denom in 1..=qmax {95 for step in 0..=4 * denom {96 if gcd(step, denom) != 1 {97 continue;98 }99 out.push(Eye {100 angle: 90.0 * step as f64 / denom as f64,101 numer: 90 * step,102 denom,103 });104 }105 }106 out.sort_by(|a, b| a.angle.total_cmp(&b.angle));107 out108}109110/// The least whole number of increments that closes a quarter turn, none once the count passes the cap.111///112/// ```113/// assert_eq!(mrlynum::tourbillon::period(18.0), Some(5));114/// assert_eq!(mrlynum::tourbillon::period(1.0), Some(90));115/// ```116pub fn period(increment: f64) -> Option<usize> {117 (1..=PERIOD_CAP).find(|count| {118 let turns = *count as f64 * increment / 90.0;119 (turns - turns.round()).abs() <= TIGHT120 })121}122123fn quarter(degrees: f64) -> f64 {124 let angle = degrees.rem_euclid(90.0);125 if 90.0 - angle <= TIGHT {126 0.0127 } else {128 angle129 }130}131132/// The angle classes of a stack read a quarter turn apart: how many the layers fall in, and how many layer pairs share one.133pub fn sharing(list: &[Layer]) -> (usize, usize) {134 let mut angles: Vec<f64> = list.iter().map(|layer| quarter(layer.degrees)).collect();135 angles.sort_by(f64::total_cmp);136 let mut counts: Vec<usize> = Vec::new();137 let mut held = f64::NEG_INFINITY;138 for angle in angles {139 match counts.last_mut() {140 Some(count) if angle - held <= TIGHT => *count += 1,141 _ => {142 counts.push(1);143 held = angle;144 }145 }146 }147 let pairs = counts.iter().map(|count| count * (count - 1) / 2).sum();148 (counts.len(), pairs)149}150151fn odd_top(top: usize) -> Result<usize> {152 if !(3..=SCALE_CAP).contains(&top) || top.is_multiple_of(2) {153 return value_error(format!(154 "the top scale must be odd between 3 and {SCALE_CAP}."155 ));156 }157 Ok(top)158}159160fn side(size: usize) -> Result<usize> {161 if !(SIZE_FLOOR..=SIZE_CAP).contains(&size) {162 return value_error(format!(163 "the raster must be between {SIZE_FLOOR} and {SIZE_CAP} pixels."164 ));165 }166 Ok(size)167}168169fn blend_of(name: &str) -> Result<Blend> {170 match Blend::named(name) {171 Some(blend) => Ok(blend),172 None => value_error(format!(173 "blend {name:?} is not mean, sum, union, meet, parity or difference."174 )),175 }176}177178fn linear(blend: Blend) -> bool {179 matches!(blend, Blend::Mean | Blend::Sum)180}181182fn chosen(scale: usize, set: &str) -> Result<bool> {183 match set {184 "odd" => Ok(true),185 "primes" => Ok(is_prime(scale)),186 "squarefree" => Ok(squarefree(scale)),187 "prime powers" => Ok(factorize(scale).len() == 1),188 _ => value_error(format!(189 "layer set {set:?} is not odd, primes, squarefree or prime powers."190 )),191 }192}193194fn weight_of(scale: usize, weights: &str) -> Result<f64> {195 match weights {196 "plain" => Ok(1.0),197 "mobius" => Ok(f64::from(mobius(scale))),198 "harmonic" => Ok(1.0 / scale as f64),199 _ => value_error(format!(200 "weights {weights:?} are not plain, mobius or harmonic."201 )),202 }203}204205/// The layers of a stack: every scale one, three, five up to the top the set keeps, each with its weight and its angle.206///207/// The schedule is unspun, degrees, golden, primes, random or gaussian, the layer set is odd, primes, squarefree or prime powers and the weights are plain, mobius or harmonic. The weights are rescaled so their magnitudes sum to the layer count, so the mean blend stays paper coverage and the sum blend stays the inked layer count.208pub fn layers(209 top: usize,210 schedule: &str,211 increment: f64,212 set: &str,213 weights: &str,214 seed: u32,215) -> Result<Vec<Layer>> {216 let top = odd_top(top)?;217 let ladder = primes(LADDER);218 let mut rng = Rng::new(u64::from(seed));219 let mut out: Vec<Layer> = Vec::new();220 let mut index = 0usize;221 for scale in (1..=top).step_by(2) {222 if !chosen(scale, set)? {223 continue;224 }225 index += 1;226 let degrees = match schedule {227 "unspun" => 0.0,228 "degrees" => index as f64 * increment,229 "golden" => index as f64 * golden(),230 "primes" => match index {231 1 => 0.0,232 _ => ladder.get(index - 2).copied().unwrap_or(0) as f64,233 },234 "random" => rng.unit() * 360.0,235 "gaussian" => {236 squares(scale).map_or(0.0, |(a, b)| (b as f64).atan2(a as f64).to_degrees())237 }238 _ => {239 return value_error(format!(240 "schedule {schedule:?} is not unspun, degrees, golden, primes, random or gaussian."241 ))242 }243 };244 out.push(Layer {245 scale,246 weight: weight_of(scale, weights)?,247 degrees,248 });249 }250 let mass: f64 = out.iter().map(|layer| layer.weight.abs()).sum();251 if mass > 0.0 {252 let factor = out.len() as f64 / mass;253 for layer in &mut out {254 layer.weight *= factor;255 }256 }257 Ok(out)258}259260fn lit(scale: usize) -> f32 {261 f32::from(u8::from(((scale as f64 * 0.5).floor() as i64) & 1 == 1))262}263264fn sample(list: &[Layer], blend: Blend) -> Vec<f32> {265 list.iter()266 .map(|layer| {267 let ink = lit(layer.scale);268 if linear(blend) {269 layer.weight as f32 * ink270 } else {271 ink272 }273 })274 .collect()275}276277// THE RASTER278279fn disc(size: usize) -> Vec<bool> {280 let mut out = vec![false; size * size];281 for row in 0..size {282 let v = (row as f64 + 0.5) / size as f64 - 0.5;283 for column in 0..size {284 let u = (column as f64 + 0.5) / size as f64 - 0.5;285 out[row * size + column] = u * u + v * v <= 0.25;286 }287 }288 out289}290291fn cells(layer: &Layer, size: usize, inside: &[bool], ink: &mut [u8]) {292 let (sin, cos) = layer.degrees.to_radians().sin_cos();293 let scale = layer.scale as f64;294 let step = 1.0 / size as f64;295 let (dx, dy) = (cos * step, -sin * step);296 let edge = 0.5 * step - 0.5;297 for row in 0..size {298 let v = (row as f64 + 0.5) * step - 0.5;299 let mut x = cos * edge + sin * v + 0.5;300 let mut y = cos * v - sin * edge + 0.5;301 for column in 0..size {302 let at = row * size + column;303 let held = inside[at] && (0.0..1.0).contains(&x) && (0.0..1.0).contains(&y);304 ink[at] = u8::from(305 held && (scale * x).floor() as i64 & 1 == 1 && (scale * y).floor() as i64 & 1 == 1,306 );307 x += dx;308 y += dy;309 }310 }311}312313fn edges(size: usize, inside: &[bool], ink: &[u8], mark: &mut [u8]) {314 for row in 0..size {315 for column in 0..size {316 let at = row * size + column;317 if !inside[at] {318 mark[at] = 0;319 continue;320 }321 let right = column + 1 < size && inside[at + 1] && ink[at + 1] != ink[at];322 let below = row + 1 < size && inside[at + size] && ink[at + size] != ink[at];323 mark[at] = u8::from(right || below);324 }325 }326}327328fn corners(layer: &Layer, size: usize, inside: &[bool], mark: &mut [u8]) {329 mark.fill(0);330 let (sin, cos) = layer.degrees.to_radians().sin_cos();331 let scale = layer.scale;332 for down in (1..scale).step_by(2) {333 for across in (1..scale).step_by(2) {334 for (vx, vy) in [335 (across, down),336 (across + 1, down),337 (across, down + 1),338 (across + 1, down + 1),339 ] {340 let x = vx as f64 / scale as f64 - 0.5;341 let y = vy as f64 / scale as f64 - 0.5;342 let u = cos * x - sin * y;343 let v = sin * x + cos * y;344 if u * u + v * v > 0.25 {345 continue;346 }347 let column = ((u + 0.5) * size as f64 - 0.5).round() as i64;348 let row = ((v + 0.5) * size as f64 - 0.5).round() as i64;349 for dy in -1..=1i64 {350 for dx in -1..=1i64 {351 let (a, b) = (column + dx, row + dy);352 if a < 0 || b < 0 || a >= size as i64 || b >= size as i64 {353 continue;354 }355 let at = b as usize * size + a as usize;356 if inside[at] {357 mark[at] = 1;358 }359 }360 }361 }362 }363 }364}365366/// Rasters the layers onto a square of the size, every one turned about the centre by its own angle and masked to the inscribed disc, then merged site by site.367///368/// The render mode is cells, the lit squares of the carpet, edges, the boundaries between its cells, or corners, the vertices of its lit squares. The weights ride on the layers under the linear blends only. Sites outside the disc are NaN.369pub fn stack(list: &[Layer], size: usize, mode: &str, blend: Blend) -> Result<Vec<f32>> {370 if !["cells", "edges", "corners"].contains(&mode) {371 return value_error(format!(372 "render mode {mode:?} is not cells, edges or corners."373 ));374 }375 if list.len() * size * size > PLANE_CAP {376 return value_error(format!(377 "{} layers on a {size} raster passes the {PLANE_CAP} site budget.",378 list.len()379 ));380 }381 let inside = disc(size);382 let mut out = vec![f32::NAN; size * size];383 if list.is_empty() {384 return Ok(out);385 }386 let mut ink = vec![0u8; size * size];387 let mut planes: Vec<Vec<u8>> = Vec::with_capacity(list.len());388 for layer in list {389 let mut mark = vec![0u8; size * size];390 match mode {391 "corners" => corners(layer, size, &inside, &mut mark),392 "edges" => {393 cells(layer, size, &inside, &mut ink);394 edges(size, &inside, &ink, &mut mark);395 }396 _ => cells(layer, size, &inside, &mut mark),397 }398 planes.push(mark);399 }400 let scaled: Vec<f32> = list401 .iter()402 .map(|layer| {403 if linear(blend) {404 layer.weight as f32405 } else {406 1.0407 }408 })409 .collect();410 let mut buffer = vec![0f32; list.len()];411 for at in 0..size * size {412 if !inside[at] {413 continue;414 }415 for (slot, plane) in planes.iter().enumerate() {416 buffer[slot] = scaled[slot] * f32::from(plane[at]);417 }418 out[at] = blend.fold(&buffer);419 }420 Ok(out)421}422423fn peaks(field: &[f32], size: usize, top: usize) -> Vec<[f64; 3]> {424 let block = (size / top.max(1)).max(1);425 let wide = size.div_ceil(block);426 let mut best: Vec<(f64, usize, usize)> = Vec::with_capacity(wide * wide);427 for tile in 0..wide * wide {428 let (bx, by) = (tile % wide * block, tile / wide * block);429 let mut hit: Option<(f64, usize, usize)> = None;430 for row in by..(by + block).min(size) {431 for column in bx..(bx + block).min(size) {432 let value = field[row * size + column];433 if value.is_nan() {434 continue;435 }436 let value = f64::from(value);437 if hit.is_none_or(|(seen, _, _)| value > seen) {438 hit = Some((value, column, row));439 }440 }441 }442 if let Some(found) = hit {443 best.push(found);444 }445 }446 best.sort_by(|a, b| b.0.total_cmp(&a.0));447 let apart = block as f64 / size as f64;448 let mut out: Vec<[f64; 3]> = Vec::new();449 for (value, column, row) in best {450 let x = (column as f64 + 0.5) / size as f64;451 let y = (row as f64 + 0.5) / size as f64;452 if out453 .iter()454 .all(|site| (site[0] - x).hypot(site[1] - y) >= apart)455 {456 out.push([x, y, value]);457 }458 if out.len() == PEAKS {459 break;460 }461 }462 out463}464465// THE STACK466467/// Spins the odd parity carpets at the scales one, three, five up to the top into one stack on a square of the size, every layer turned about the centre by its own angle and masked to the inscribed disc, so every pixel sees every layer.468///469/// The schedule is unspun, degrees, golden, primes, random or gaussian, the layer set is odd, primes, squarefree or prime powers, the weights are plain, mobius or harmonic, the render mode is cells, edges or corners and the blend is mean, sum, union, meet, parity or difference. The weights are scaled to average magnitude one and apply only under the linear blends, so mean is paper coverage and sum is the inked layer count. Sites outside the disc are NaN.470#[allow(clippy::too_many_arguments)]471pub fn field(472 top: usize,473 size: usize,474 schedule: &str,475 increment: f64,476 set: &str,477 weights: &str,478 mode: &str,479 blend: &str,480 seed: u32,481) -> Result<Vec<f32>> {482 let size = side(size)?;483 let blend = blend_of(blend)?;484 let list = layers(top, schedule, increment, set, weights, seed)?;485 stack(&list, size, mode, blend)486}487488/// Reads a spun stack against the schedule that made it: the layer count, the first eight scales and angles, the mean and RMS contrast over the disc, that contrast times the root of the layer count, the exact centre value, whether the blend carries the weights, the span the raster covers and the brightest three sites.489///490/// The centre is computed from the layers at the exact half-half point, never sampled off the raster, so no turn of the layers moves it.491#[allow(clippy::too_many_arguments)]492pub fn stats(493 field: &[f32],494 size: usize,495 top: usize,496 schedule: &str,497 increment: f64,498 set: &str,499 weights: &str,500 blend: &str,501 seed: u32,502) -> Result<Stats> {503 let size = side(size)?;504 if field.len() != size * size {505 return value_error("the field must be size by size with size at least 1.");506 }507 let blend = blend_of(blend)?;508 let list = layers(top, schedule, increment, set, weights, seed)?;509 let seen: Vec<f64> = field510 .iter()511 .filter(|value| !value.is_nan())512 .map(|value| f64::from(*value))513 .collect();514 let count = seen.len().max(1) as f64;515 let mean = seen.iter().sum::<f64>() / count;516 let rms = (seen.iter().map(|value| (value - mean).powi(2)).sum::<f64>() / count).sqrt();517 let low = seen.iter().cloned().fold(f64::INFINITY, f64::min);518 let high = seen.iter().cloned().fold(f64::NEG_INFINITY, f64::max);519 let centre = if list.is_empty() {520 0.0521 } else {522 f64::from(blend.fold(&sample(&list, blend)))523 };524 let (classes, pairs) = sharing(&list);525 Ok(Stats {526 layers: list.len(),527 scales: list.iter().take(SHOWN).map(|layer| layer.scale).collect(),528 angles: list.iter().take(SHOWN).map(|layer| layer.degrees).collect(),529 mean,530 rms,531 faded: rms * (list.len() as f64).sqrt(),532 centre,533 weighted: linear(blend),534 low: if low.is_finite() { low } else { 0.0 },535 high: if high.is_finite() { high } else { 0.0 },536 inside: seen.len(),537 peaks: peaks(field, size, top),538 period: period(increment),539 classes,540 pairs,541 })542}543544#[cfg(test)]545mod tests {546 use super::*;547548 #[test]549 fn the_schedule_keeps_the_odd_scales_and_turns_them_by_the_primes() {550 let odd = layers(55, "primes", 0.0, "odd", "plain", 1).unwrap();551 assert_eq!(odd.len(), 28);552 let scales: Vec<usize> = odd.iter().take(SHOWN).map(|layer| layer.scale).collect();553 assert_eq!(scales, vec![1, 3, 5, 7, 9, 11, 13, 15]);554 let angles: Vec<f64> = odd.iter().take(SHOWN).map(|layer| layer.degrees).collect();555 assert_eq!(angles, vec![0.0, 2.0, 3.0, 5.0, 7.0, 11.0, 13.0, 17.0]);556 assert_eq!(557 layers(199, "unspun", 0.0, "primes", "plain", 1)558 .unwrap()559 .len(),560 45561 );562 let counts: Vec<usize> = ["odd", "primes", "squarefree", "prime powers"]563 .iter()564 .map(|set| layers(55, "unspun", 0.0, set, "plain", 1).unwrap().len())565 .collect();566 assert_eq!(counts, vec![28, 15, 23, 19]);567 assert!(layers(54, "unspun", 0.0, "odd", "plain", 1).is_err());568 assert!(layers(55, "spiral", 0.0, "odd", "plain", 1).is_err());569 }570571 #[test]572 fn the_unspun_stack_counts_eighteen_layers_at_its_brightest() {573 let list = layers(55, "unspun", 0.0, "odd", "plain", 1).unwrap();574 let raster = stack(&list, 512, "cells", Blend::Sum).unwrap();575 assert_eq!(raster.len(), 512 * 512);576 let high = raster577 .iter()578 .filter(|value| !value.is_nan())579 .cloned()580 .fold(f32::NEG_INFINITY, f32::max);581 assert_eq!(high, 18.0);582 assert!(raster[0].is_nan());583 assert!(stack(&list, 512, "dots", Blend::Sum).is_err());584 }585586 #[test]587 fn the_field_and_the_stats_hold_the_centre_at_half_the_layers() {588 let raster = field(55, 512, "unspun", 0.0, "odd", "plain", "cells", "mean", 1).unwrap();589 let read = stats(&raster, 512, 55, "unspun", 0.0, "odd", "plain", "mean", 1).unwrap();590 assert_eq!(read.layers, 28);591 assert_eq!(read.centre, 14.0 / 28.0);592 assert_eq!(format!("{:.6}", read.high), format!("{:.6}", 18.0 / 28.0));593 assert!(read.weighted);594 let counted = stats(&raster, 512, 55, "unspun", 0.0, "odd", "plain", "sum", 1).unwrap();595 assert_eq!(counted.centre, 14.0);596 let folded = stats(&raster, 512, 55, "unspun", 0.0, "odd", "plain", "parity", 1).unwrap();597 assert_eq!(folded.centre, 0.0);598 assert!(!folded.weighted);599 assert!(field(55, 4096, "unspun", 0.0, "odd", "plain", "cells", "mean", 1).is_err());600 assert!(stats(&raster, 256, 55, "unspun", 0.0, "odd", "plain", "mean", 1).is_err());601 }602}