research-spirograph.rs

1use mrlycore::errors::Result;2use mrlyfig::{ink, plot, save, Board};3use mrlylab::roulette;4use mrlymath::two;5use mrlynum::spirograph::{6    cover, disc, distinct, nodes as law, pencils, trace, track, Cover, Pencil, Track,7};8use std::collections::HashSet;910const CODE: u128 = 495;11const SIDE: usize = 3;12const LEVEL: usize = 1;13const BASE: usize = 3;14const RING: usize = 7;15const WHEEL: usize = 3;16const REACH: f64 = 0.9;17const SAMPLES: usize = 8000;18const CURVES: usize = 8;19const CROSSINGS: usize = 1288;20const RASTER: usize = 1024;21const TOL: f64 = 4e-4;22const MARGIN: f64 = 0.08;23const GRID: usize = 256;24const SUB: usize = 4;25const WASH: f64 = 0.3;26const THIN: f64 = 1.4;27const HAIR: f64 = 1.1;28const DOT: f64 = 2.6;2930fn main() -> Result<()> {31    let tile = two::create(CODE, SIDE, LEVEL, 0, BASE)?;32    let types = tile.types().bytes().to_vec();33    let path = track("in", RING, WHEEL, 4, 1)?;34    let pens = pencils(&types, tile.width(), tile.height(), "fill", REACH, 0.0, 1)?;35    assert_eq!(distinct(&path, &pens, true), CURVES);36    assert_eq!(law(&path, &pens, true), Some(CROSSINGS as u64));37    assert_eq!(38        roulette::nodes(&path, &pens, SAMPLES, TOL)?.total(),39        CROSSINGS40    );41    let marks = crossings(&path, &pens, SAMPLES)?;42    assert_eq!(marks.len(), CROSSINGS);43    let bounds = disc(&path, &pens)?;44    let shape = cover(&path, &pens, true, 2, RASTER)?;45    let points = trace(&path, &pens, SAMPLES)?;4647    let mut board = Board::square();48    let area = board.frame(MARGIN);49    let (cx, cy) = area.center();50    let scale = area.radius() / RING as f64;51    let at = |x: f64, y: f64| (cx + (x - bounds.x) * scale, cy - (y - bounds.y) * scale);52    let span = 2.0 * bounds.radius * scale;53    wash(&mut board, &shape, (cx - span / 2.0, cy - span / 2.0), span);54    board.ring(cx, cy, RING as f64 * scale, HAIR, ink::line());55    for k in 0..CURVES {56        let curve: Vec<(f64, f64)> = (0..SAMPLES)57            .map(|i| {58                let j = 2 * (k * SAMPLES + i);59                at(f64::from(points[j]), f64::from(points[j + 1]))60            })61            .collect();62        board.polyline(&curve, THIN, ink::blue());63    }64    let dots: Vec<(f64, f64)> = marks.iter().map(|&(x, y)| at(x, y)).collect();65    plot::dots(&mut board, &dots, DOT, ink::orange());66    save("research-spirograph", &board)?;67    Ok(())68}6970// THE SHAPE7172fn wash(board: &mut Board, shape: &Cover, corner: (f64, f64), span: f64) {73    let n = shape.side;74    let cell = span / n as f64;75    let x0 = corner.0.floor().max(0.0) as usize;76    let y0 = corner.1.floor().max(0.0) as usize;77    let x1 = (corner.0 + span).ceil().max(0.0) as usize;78    let y1 = (corner.1 + span).ceil().max(0.0) as usize;79    for py in y0..y1 {80        for px in x0..x1 {81            let mut hit = 0;82            for sy in 0..SUB {83                for sx in 0..SUB {84                    let u = (px as f64 + (sx as f64 + 0.5) / SUB as f64 - corner.0) / cell;85                    let v = (py as f64 + (sy as f64 + 0.5) / SUB as f64 - corner.1) / cell;86                    if u < 0.0 || v < 0.0 {87                        continue;88                    }89                    let (column, row) = (u as usize, v as usize);90                    if column < n && row < n && shape.mask[row * n + column] == 3 {91                        hit += 1;92                    }93                }94            }95            if hit > 0 {96                board.blend(97                    px,98                    py,99                    ink::dim(),100                    WASH * f64::from(hit) / (SUB * SUB) as f64,101                );102            }103        }104    }105}106107// THE NODES108109fn crossings(path: &Track, pens: &[Pencil], samples: usize) -> Result<Vec<(f64, f64)>> {110    let points = trace(path, pens, samples)?;111    let steps = samples - 1;112    let at = |k: usize, i: usize| {113        let base = 2 * (k * samples + i);114        (f64::from(points[base]), f64::from(points[base + 1]))115    };116    let (mut low, mut high) = ((f64::MAX, f64::MAX), (f64::MIN, f64::MIN));117    for pair in points.chunks_exact(2) {118        let (x, y) = (f64::from(pair[0]), f64::from(pair[1]));119        low = (low.0.min(x), low.1.min(y));120        high = (high.0.max(x), high.1.max(y));121    }122    let size = (high.0 - low.0, high.1 - low.1);123    let cell = |p: (f64, f64)| {124        let column = ((p.0 - low.0) / size.0 * GRID as f64) as usize;125        let row = ((p.1 - low.1) / size.1 * GRID as f64) as usize;126        (column.min(GRID - 1), row.min(GRID - 1))127    };128    let mut buckets: Vec<Vec<u32>> = vec![Vec::new(); GRID * GRID];129    for k in 0..pens.len() {130        for i in 0..steps {131            let (a, b) = (at(k, i), at(k, i + 1));132            let (x0, y0) = cell((a.0.min(b.0), a.1.min(b.1)));133            let (x1, y1) = cell((a.0.max(b.0), a.1.max(b.1)));134            for x in x0..=x1 {135                for y in y0..=y1 {136                    buckets[x * GRID + y].push((k * steps + i) as u32);137                }138            }139        }140    }141    let mut seen: HashSet<(u32, u32)> = HashSet::new();142    let mut out = Vec::new();143    for bucket in &buckets {144        for (u, &left) in bucket.iter().enumerate() {145            for &right in &bucket[u + 1..] {146                let (p, q) = (left as usize, right as usize);147                let ((kp, ip), (kq, iq)) = ((p / steps, p % steps), (q / steps, q % steps));148                let gap = ip.abs_diff(iq);149                if kp == kq && (gap == 1 || gap == steps - 1) {150                    continue;151                }152                let (a, b) = (at(kp, ip), at(kp, ip + 1));153                let (c, d) = (at(kq, iq), at(kq, iq + 1));154                if side(a, b, c) * side(a, b, d) < 0155                    && side(c, d, a) * side(c, d, b) < 0156                    && seen.insert((left.min(right), left.max(right)))157                {158                    out.push(meet(a, b, c, d));159                }160            }161        }162    }163    Ok(out)164}165166fn side(a: (f64, f64), b: (f64, f64), c: (f64, f64)) -> i32 {167    let turn = (b.0 - a.0) * (c.1 - a.1) - (b.1 - a.1) * (c.0 - a.0);168    match turn.partial_cmp(&0.0) {169        Some(std::cmp::Ordering::Greater) => 1,170        Some(std::cmp::Ordering::Less) => -1,171        _ => 0,172    }173}174175fn meet(a: (f64, f64), b: (f64, f64), c: (f64, f64), d: (f64, f64)) -> (f64, f64) {176    let run = (b.0 - a.0, b.1 - a.1);177    let other = (d.0 - c.0, d.1 - c.1);178    let denominator = run.0 * other.1 - run.1 * other.0;179    let step = ((c.0 - a.0) * other.1 - (c.1 - a.1) * other.0) / denominator;180    (a.0 + step * run.0, a.1 + step * run.1)181}