spin.rs

1use std::f64::consts::{PI, SQRT_2};23fn centre(size: usize) -> f64 {4    size as f64 / 2.05}67/// The radius of the corner circle of a square raster of the side, the last radius a profile reads.8pub fn reach(size: usize) -> f64 {9    size as f64 / SQRT_210}1112/// The arcs of the circle of the radius about the raster's centre: each as its start angle, end angle and the value of the one cell it lies in, zero outside.13pub fn arcs(data: &[f32], size: usize, radius: f64) -> Vec<(f64, f64, f32)> {14    assert_eq!(data.len(), size * size, "data must be size*size");15    if size == 0 {16        return Vec::new();17    }18    let c = centre(size);19    let r = radius.max(1e-9);20    let mut cuts = vec![0.0, 2.0 * PI];21    let low = (c - r).floor().max(0.0) as usize;22    let high = ((c + r).ceil() as usize).min(size);23    for line in low..=high {24        let u = (line as f64 - c) / r;25        if u.abs() < 1.0 {26            let a = u.acos();27            cuts.push(a);28            cuts.push(2.0 * PI - a);29            let b = u.asin();30            cuts.push(b.rem_euclid(2.0 * PI));31            cuts.push((PI - b).rem_euclid(2.0 * PI));32        }33    }34    cuts.sort_by(|a, b| a.partial_cmp(b).unwrap());35    let mut out = Vec::with_capacity(cuts.len());36    for pair in cuts.windows(2) {37        let (a, b) = (pair[0], pair[1]);38        if b <= a {39            continue;40        }41        let mid = (a + b) / 2.0;42        let x = c + r * mid.cos();43        let y = c + r * mid.sin();44        let mut value = 0.0;45        if x >= 0.0 && y >= 0.0 {46            let (col, row) = (x as usize, y as usize);47            if col < size && row < size {48                value = data[row * size + col];49            }50        }51        out.push((a, b, value));52    }53    out54}5556/// The exact mean of a square raster over the circle of the radius about its centre, each cell read as a constant and the outside as zero.57///58/// ```59/// let solid = vec![1.0; 16];60/// assert!((mrlynum::spin::ring(&solid, 4, 1.0) - 1.0).abs() < 1e-12);61/// assert!(mrlynum::spin::ring(&solid, 4, 3.0).abs() < 1e-12);62/// ```63pub fn ring(data: &[f32], size: usize, radius: f64) -> f64 {64    arcs(data, size, radius)65        .iter()66        .map(|&(a, b, v)| v as f64 * (b - a))67        .sum::<f64>()68        / (2.0 * PI)69}7071/// The circular-harmonic power of a raster: for every order `m` up to the last, the energy `sum |c_m(r)|^2 2 pi r dr` of its `m`-th harmonic over rings radii, each ring's coefficient exact from its arcs.72pub fn harmonics(data: &[f32], size: usize, rings: usize, orders: usize) -> Vec<f64> {73    let rings = rings.max(2);74    let far = reach(size);75    let step = far / (rings - 1) as f64;76    let mut power = vec![0.0; orders + 1];77    for k in 0..rings {78        let r = k as f64 * step;79        let pieces = arcs(data, size, r);80        let mut re = vec![0.0; orders + 1];81        let mut im = vec![0.0; orders + 1];82        for &(a, b, v) in &pieces {83            if v == 0.0 {84                continue;85            }86            let v = v as f64;87            re[0] += v * (b - a);88            for m in 1..=orders {89                let f = m as f64;90                re[m] += v * ((f * b).sin() - (f * a).sin()) / f;91                im[m] += v * ((f * b).cos() - (f * a).cos()) / f;92            }93        }94        for m in 0..=orders {95            let (x, y) = (re[m] / (2.0 * PI), im[m] / (2.0 * PI));96            power[m] += (x * x + y * y) * 2.0 * PI * r * step;97        }98    }99    power100}101102fn gcd(a: usize, b: usize) -> usize {103    if b == 0 {104        a105    } else {106        gcd(b, a % b)107    }108}109110/// The rotation order a harmonic power spectrum reveals: the gcd of the orders carrying more than a ten-thousandth of the power, the share pixel aliasing stays under, or zero when none does.111pub fn turns(power: &[f64]) -> usize {112    let total: f64 = power.iter().sum();113    power114        .iter()115        .enumerate()116        .skip(1)117        .filter(|&(_, &p)| p > 1e-4 * total)118        .fold(0, |g, (m, _)| gcd(g, m))119}120121/// The petals a full radial stack of the copies shows on a design of the rotation order: their least common multiple.122pub fn petals(copies: usize, order: usize) -> usize {123    if copies == 0 || order == 0 {124        return 0;125    }126    copies / gcd(copies, order) * order127}128129/// The way radial copies merge: their mean, their sum, their union, their meet, their parity or what the first keeps that no other has.130#[derive(Clone, Copy, Debug, PartialEq, Eq)]131pub enum Blend {132    /// The mean of the copies.133    Mean,134    /// The sum of the copies.135    Sum,136    /// The largest copy.137    Union,138    /// The smallest copy.139    Meet,140    /// The sum folded to its parity.141    Parity,142    /// The first copy less the largest of the rest, floored at zero.143    Difference,144}145146impl Blend {147    /// Reads a blend by name: mean, sum, union, meet, parity or difference.148    pub fn named(name: &str) -> Option<Blend> {149        match name {150            "mean" => Some(Blend::Mean),151            "sum" => Some(Blend::Sum),152            "union" => Some(Blend::Union),153            "meet" => Some(Blend::Meet),154            "parity" => Some(Blend::Parity),155            "difference" => Some(Blend::Difference),156            _ => None,157        }158    }159160    /// Merges one site's copies into the blended value.161    pub fn fold(self, values: &[f32]) -> f32 {162        let sum: f32 = values.iter().sum();163        let max = values.iter().cloned().fold(f32::NEG_INFINITY, f32::max);164        match self {165            Blend::Mean => sum / values.len() as f32,166            Blend::Sum => sum,167            Blend::Union => max,168            Blend::Meet => values.iter().cloned().fold(f32::INFINITY, f32::min),169            Blend::Parity => 1.0 - (sum.rem_euclid(2.0) - 1.0).abs(),170            Blend::Difference => {171                let rest = values[1..].iter().cloned().fold(0.0, f32::max);172                (values[0] - rest).max(0.0)173            }174        }175    }176}177178/// Stacks a raster radially: copies turned by multiples of the step, in turns, about the centre and merged by the blend, on an output raster of the side whose inscribed circle is the source's corner circle, every pixel the mean of samples by samples points.179pub fn radial(180    data: &[f32],181    size: usize,182    out: usize,183    copies: usize,184    step: f64,185    blend: Blend,186    samples: usize,187) -> Vec<f32> {188    assert_eq!(data.len(), size * size, "data must be size*size");189    let copies = copies.max(1);190    let samples = samples.max(1);191    let c = centre(size);192    let scale = 2.0 * reach(size) / out as f64;193    let turns: Vec<(f64, f64)> = (0..copies)194        .map(|k| {195            let angle = 2.0 * PI * step * k as f64;196            (angle.cos(), angle.sin())197        })198        .collect();199    let mut values = vec![0.0f32; copies];200    let mut field = Vec::with_capacity(out * out);201    for i in 0..out {202        for j in 0..out {203            let mut total = 0.0;204            for a in 0..samples {205                for b in 0..samples {206                    let px =207                        (j as f64 + (b as f64 + 0.5) / samples as f64 - out as f64 / 2.0) * scale;208                    let py =209                        (i as f64 + (a as f64 + 0.5) / samples as f64 - out as f64 / 2.0) * scale;210                    for (k, &(cos, sin)) in turns.iter().enumerate() {211                        let x = c + px * cos + py * sin;212                        let y = c - px * sin + py * cos;213                        values[k] = if x >= 0.0 && y >= 0.0 && x < size as f64 && y < size as f64 {214                            data[y as usize * size + x as usize]215                        } else {216                            0.0217                        };218                    }219                    total += blend.fold(&values);220                }221            }222            field.push(total / (samples * samples) as f32);223        }224    }225    field226}227228/// The ring profile: the circle means at steps radii spaced evenly from the centre to the corner circle.229pub fn profile(data: &[f32], size: usize, steps: usize) -> Vec<f32> {230    let steps = steps.max(2);231    let far = reach(size);232    (0..steps)233        .map(|k| ring(data, size, far * k as f64 / (steps - 1) as f64) as f32)234        .collect()235}236237/// The wheel: a profile spread over a square raster of the side, the corner circle it ends on drawn as the inscribed circle, every pixel reading the profile at its own radius.238pub fn wheel(profile: &[f32], size: usize) -> Vec<f32> {239    let last = profile.len().saturating_sub(1);240    if last == 0 {241        return vec![profile.first().copied().unwrap_or(0.0); size * size];242    }243    let c = centre(size);244    let scale = 2.0 * last as f64 / size as f64;245    let mut out = Vec::with_capacity(size * size);246    for row in 0..size {247        for col in 0..size {248            let (dx, dy) = (col as f64 + 0.5 - c, row as f64 + 0.5 - c);249            let t = (dx * dx + dy * dy).sqrt() * scale;250            let i = (t.floor() as usize).min(last);251            let f = (t - i as f64) as f32;252            let value = if i == last {253                profile[last]254            } else {255                profile[i] * (1.0 - f) + profile[i + 1] * f256            };257            out.push(value);258        }259    }260    out261}262263/// The mass a profile carries, the trapezoid integral of `2 pi r F(r)` in cells of the raster it came from.264pub fn mass(profile: &[f32], size: usize) -> f64 {265    let last = profile.len().saturating_sub(1);266    if last == 0 {267        return 0.0;268    }269    let step = reach(size) / last as f64;270    let weight = |k: usize| 2.0 * PI * (k as f64 * step) * profile[k] as f64;271    let inner: f64 = (1..last).map(weight).sum();272    step * (inner + (weight(0) + weight(last)) / 2.0)273}274275/// The mass a profile carries inside the radius, the trapezoid integral of `2 pi r F(r)` from the centre out, in cells of the raster it came from.276///277/// ```278/// let solid = vec![1.0; 64];279/// let rings = mrlynum::spin::profile(&solid, 8, 4000);280/// let inner = mrlynum::spin::mass_within(&rings, 8, 3.0);281/// assert!((inner / (9.0 * std::f64::consts::PI) - 1.0).abs() < 1e-3);282/// ```283pub fn mass_within(profile: &[f32], size: usize, radius: f64) -> f64 {284    let last = profile.len().saturating_sub(1);285    if last == 0 || radius <= 0.0 {286        return 0.0;287    }288    let step = reach(size) / last as f64;289    let weight = |k: usize| 2.0 * PI * (k as f64 * step) * profile[k] as f64;290    let full = (radius / step).floor().min(last as f64) as usize;291    let mut total = 0.0;292    for k in 1..=full {293        total += (weight(k - 1) + weight(k)) / 2.0 * step;294    }295    if full < last {296        let rest = radius - full as f64 * step;297        let share = rest / step;298        let edge = weight(full) + (weight(full + 1) - weight(full)) * share;299        total += (weight(full) + edge) / 2.0 * rest;300    }301    total302}303304#[cfg(test)]305mod tests {306    use super::*;307    use mrlycore::atoms;308309    fn floats(grid: &mrlycore::tensor::Tensor) -> Vec<f32> {310        grid.bytes().iter().map(|&b| b as f32).collect()311    }312313    #[test]314    fn a_solid_square_leaves_through_four_arcs() {315        let side = 8usize;316        let solid = floats(&atoms::ones_2d(side));317        for r in [0.5, 2.0, 3.99] {318            assert!((ring(&solid, side, r) - 1.0).abs() < 1e-12);319        }320        for r in [4.5, 5.0, 5.5] {321            let expect = 1.0 - 4.0 / PI * (side as f64 / 2.0 / r).acos();322            assert!((ring(&solid, side, r) - expect).abs() < 1e-12);323        }324        assert!(ring(&solid, side, reach(side) + 0.01).abs() < 1e-12);325    }326327    #[test]328    fn the_carpet_opens_on_a_black_disc() {329        let carpet = atoms::carpet_nd(3, 2).kron(&atoms::carpet_nd(3, 2));330        let data = floats(&carpet);331        assert!(ring(&data, 9, 1.4).abs() < 1e-12);332        assert!(ring(&data, 9, 2.0) > 0.0);333        assert!((ring(&data, 9, 1.0) - 0.0).abs() < 1e-12);334    }335336    #[test]337    fn the_mass_of_the_profile_is_the_fill() {338        let carpet = atoms::carpet_nd(3, 2)339            .kron(&atoms::carpet_nd(3, 2))340            .kron(&atoms::carpet_nd(3, 2));341        let data = floats(&carpet);342        let fills = data.iter().sum::<f32>() as f64;343        assert_eq!(fills, 512.0);344        let rings = profile(&data, 27, 4000);345        assert!((mass(&rings, 27) - fills).abs() / fills < 0.002);346    }347348    #[test]349    fn the_spin_mass_scales_by_the_fill_about_a_filled_corner() {350        let tile = atoms::carpet_nd(3, 2);351        let carpet = tile.kron(&tile).kron(&tile).kron(&tile);352        let side = 81usize;353        let wide = 2 * side;354        let mut data = vec![0.0f32; wide * wide];355        let bytes = carpet.bytes();356        for row in 0..side {357            for col in 0..side {358                data[(row + side) * wide + col + side] = bytes[row * side + col] as f32;359            }360        }361        let rings = profile(&data, wide, 4000);362        assert!(mass_within(&rings, wide, 1e-6).abs() < 1e-9);363        for radius in [12.0f64, 18.0, 27.0] {364            let near = mass_within(&rings, wide, radius);365            let far = mass_within(&rings, wide, 3.0 * radius);366            assert!(near > 0.0);367            assert!((far / (8.0 * near) - 1.0).abs() < 0.08, "radius {radius}");368        }369        let whole = mass_within(&rings, wide, reach(wide));370        assert!((whole - mass(&rings, wide)).abs() / whole < 1e-9);371    }372373    #[test]374    fn the_wheel_reads_the_profile_by_radius() {375        let rings = profile(&floats(&atoms::ones_2d(8)), 8, 64);376        let spun = wheel(&rings, 16);377        assert_eq!(spun.len(), 256);378        assert!((spun[8 * 16 + 8] - 1.0).abs() < 1e-6);379        assert!((spun[8 * 16 + 12] - 1.0).abs() < 1e-6);380        assert!(spun[8 * 16 + 15] < 0.2);381        assert_eq!(spun[0], 0.0);382        assert_eq!(profile(&[], 0, 0).len(), 2);383    }384385    fn at(field: &[f32], out: usize, size: usize, px: f64, py: f64) -> f32 {386        let scale = 2.0 * reach(size) / out as f64;387        let col = (px / scale + out as f64 / 2.0) as usize;388        let row = (py / scale + out as f64 / 2.0) as usize;389        field[row * out + col]390    }391392    #[test]393    fn two_squares_at_an_eighth_turn_make_a_star() {394        let solid = floats(&atoms::ones_2d(8));395        let stack = |blend| radial(&solid, 8, 64, 2, 0.125, blend, 1);396        let probe = |blend, px, py| at(&stack(blend), 64, 8, px, py);397        assert_eq!(probe(Blend::Union, 0.0, 0.0), 1.0);398        assert_eq!(probe(Blend::Union, 0.0, -5.0), 1.0);399        assert_eq!(probe(Blend::Meet, 0.0, -5.0), 0.0);400        assert_eq!(probe(Blend::Mean, 0.0, -5.0), 0.5);401        assert_eq!(probe(Blend::Sum, 0.0, 0.0), 2.0);402        assert_eq!(probe(Blend::Parity, 0.0, -5.0), 1.0);403        assert_eq!(probe(Blend::Parity, 0.0, 0.0), 0.0);404        assert_eq!(probe(Blend::Difference, 0.0, -5.0), 0.0);405        assert_eq!(probe(Blend::Difference, 3.8, -3.8), 1.0);406        assert_eq!(Blend::named("soup"), None);407    }408409    #[test]410    fn the_harmonics_read_the_rotation_order() {411        let square = harmonics(&floats(&atoms::ones_2d(8)), 8, 128, 12);412        assert!(square[0] > 0.0);413        assert!(square[4] > 1e-3 * square[0]);414        assert!(square[8] > 1e-3 * square[0]);415        for m in [1, 2, 3, 5, 6, 7, 9, 10, 11] {416            assert!(square[m] < 1e-9 * square[0], "m {m}");417        }418        assert_eq!(turns(&square), 4);419        let mut bar = vec![0.0f32; 16];420        for i in [5, 6, 9, 10, 4, 11] {421            bar[i] = 1.0;422        }423        assert_eq!(turns(&harmonics(&bar, 4, 64, 8)), 2);424        let mut blob = vec![0.0f32; 16];425        blob[0] = 1.0;426        assert_eq!(turns(&harmonics(&blob, 4, 64, 8)), 1);427        assert_eq!(turns(&[1.0, 0.0, 0.0]), 0);428        assert_eq!(429            (petals(6, 4), petals(8, 4), petals(5, 1), petals(3, 0)),430            (12, 8, 5, 0)431        );432    }433}