spin.rs

1use crate::core::error::{shape_error, Result};2use crate::num::factor::gcd;3use serde::{Deserialize, Serialize};4use std::f64::consts::{PI, SQRT_2};56fn centre(size: usize) -> f64 {7    size as f64 / 2.08}910/// The radius of the corner circle of a square raster of the side, the last radius a profile reads.11pub fn reach(size: usize) -> f64 {12    size as f64 / SQRT_213}1415/// 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.16///17/// ```18/// let arcs = mrlyrs::math::spin::arcs(&[1.0, 1.0, 1.0, 1.0], 2, 0.5).unwrap();19/// assert!(arcs.iter().all(|arc| arc.2 == 1.0));20/// ```21///22/// # Errors23///24/// Errors when the sample count is not the side squared.25pub fn arcs(data: &[f32], size: usize, radius: f64) -> Result<Vec<(f64, f64, f32)>> {26    if data.len() != size * size {27        return shape_error(format!(28            "a raster of side {size} needs {} samples, got {}.",29            size * size,30            data.len()31        ));32    }33    if size == 0 {34        return Ok(Vec::new());35    }36    let c = centre(size);37    let r = radius.max(1e-9);38    let mut cuts = vec![0.0, 2.0 * PI];39    let low = (c - r).floor().max(0.0) as usize;40    let high = ((c + r).ceil() as usize).min(size);41    for line in low..=high {42        let u = (line as f64 - c) / r;43        if u.abs() < 1.0 {44            let a = u.acos();45            cuts.push(a);46            cuts.push(2.0 * PI - a);47            let b = u.asin();48            cuts.push(b.rem_euclid(2.0 * PI));49            cuts.push((PI - b).rem_euclid(2.0 * PI));50        }51    }52    cuts.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));53    let mut out = Vec::with_capacity(cuts.len());54    for pair in cuts.windows(2) {55        let (a, b) = (pair[0], pair[1]);56        if b <= a {57            continue;58        }59        let mid = (a + b) / 2.0;60        let x = c + r * mid.cos();61        let y = c + r * mid.sin();62        let mut value = 0.0;63        if x >= 0.0 && y >= 0.0 {64            let (col, row) = (x as usize, y as usize);65            if col < size && row < size {66                value = data[row * size + col];67            }68        }69        out.push((a, b, value));70    }71    Ok(out)72}7374/// 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.75///76/// ```77/// let solid = vec![1.0; 16];78/// assert!((mrlyrs::math::spin::ring(&solid, 4, 1.0).unwrap() - 1.0).abs() < 1e-12);79/// assert!(mrlyrs::math::spin::ring(&solid, 4, 3.0).unwrap().abs() < 1e-12);80/// ```81///82/// # Errors83///84/// Errors when the sample count is not the side squared.85pub fn ring(data: &[f32], size: usize, radius: f64) -> Result<f64> {86    Ok(arcs(data, size, radius)?87        .iter()88        .map(|&(a, b, v)| v as f64 * (b - a))89        .sum::<f64>()90        / (2.0 * PI))91}9293/// 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.94///95/// # Errors96///97/// Errors when the sample count is not the side squared.98pub fn harmonics(data: &[f32], size: usize, rings: usize, orders: usize) -> Result<Vec<f64>> {99    let rings = rings.max(2);100    let far = reach(size);101    let step = far / (rings - 1) as f64;102    let mut power = vec![0.0; orders + 1];103    for k in 0..rings {104        let r = k as f64 * step;105        let pieces = arcs(data, size, r)?;106        let mut re = vec![0.0; orders + 1];107        let mut im = vec![0.0; orders + 1];108        for &(a, b, v) in &pieces {109            if v == 0.0 {110                continue;111            }112            let v = v as f64;113            re[0] += v * (b - a);114            for m in 1..=orders {115                let f = m as f64;116                re[m] += v * ((f * b).sin() - (f * a).sin()) / f;117                im[m] += v * ((f * b).cos() - (f * a).cos()) / f;118            }119        }120        for m in 0..=orders {121            let (x, y) = (re[m] / (2.0 * PI), im[m] / (2.0 * PI));122            power[m] += (x * x + y * y) * 2.0 * PI * r * step;123        }124    }125    Ok(power)126}127128/// 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.129pub fn turns(power: &[f64]) -> usize {130    let total: f64 = power.iter().sum();131    power132        .iter()133        .enumerate()134        .skip(1)135        .filter(|&(_, &p)| p > 1e-4 * total)136        .fold(0, |g, (m, _)| gcd(g as u128, m as u128) as usize)137}138139/// The petals a full radial stack of the copies shows on a design of the rotation order: their least common multiple.140pub fn petals(copies: usize, order: usize) -> usize {141    if copies == 0 || order == 0 {142        return 0;143    }144    copies / gcd(copies as u128, order as u128) as usize * order145}146147/// The way radial copies merge: their mean, their sum, their union, their meet, their parity or what the first keeps that no other has.148#[derive(Clone, Copy, Debug, PartialEq, Eq, Serialize, Deserialize)]149pub enum Blend {150    /// The mean of the copies.151    Mean,152    /// The sum of the copies.153    Sum,154    /// The largest copy.155    Union,156    /// The smallest copy.157    Meet,158    /// The sum folded to its parity.159    Parity,160    /// The first copy less the largest of the rest, floored at zero.161    Difference,162}163164impl Blend {165    /// Reads a blend by name: mean, sum, union, meet, parity or difference.166    pub fn named(name: &str) -> Option<Blend> {167        match name {168            "mean" => Some(Blend::Mean),169            "sum" => Some(Blend::Sum),170            "union" => Some(Blend::Union),171            "meet" => Some(Blend::Meet),172            "parity" => Some(Blend::Parity),173            "difference" => Some(Blend::Difference),174            _ => None,175        }176    }177178    /// Merges one site's copies into the blended value.179    pub fn fold(self, values: &[f32]) -> f32 {180        let sum: f32 = values.iter().sum();181        let max = values.iter().cloned().fold(f32::NEG_INFINITY, f32::max);182        match self {183            Blend::Mean => sum / values.len() as f32,184            Blend::Sum => sum,185            Blend::Union => max,186            Blend::Meet => values.iter().cloned().fold(f32::INFINITY, f32::min),187            Blend::Parity => 1.0 - (sum.rem_euclid(2.0) - 1.0).abs(),188            Blend::Difference => {189                let rest = values[1..].iter().cloned().fold(0.0, f32::max);190                (values[0] - rest).max(0.0)191            }192        }193    }194}195196/// 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.197///198/// # Errors199///200/// Errors when the sample count is not the side squared.201pub fn radial(202    data: &[f32],203    size: usize,204    out: usize,205    copies: usize,206    step: f64,207    blend: Blend,208    samples: usize,209) -> Result<Vec<f32>> {210    if data.len() != size * size {211        return shape_error(format!(212            "a raster of side {size} needs {} samples, got {}.",213            size * size,214            data.len()215        ));216    }217    let copies = copies.max(1);218    let samples = samples.max(1);219    let c = centre(size);220    let scale = 2.0 * reach(size) / out as f64;221    let turns: Vec<(f64, f64)> = (0..copies)222        .map(|k| {223            let angle = 2.0 * PI * step * k as f64;224            (angle.cos(), angle.sin())225        })226        .collect();227    let mut values = vec![0.0f32; copies];228    let mut field = Vec::with_capacity(out * out);229    for i in 0..out {230        for j in 0..out {231            let mut total = 0.0;232            for a in 0..samples {233                for b in 0..samples {234                    let px =235                        (j as f64 + (b as f64 + 0.5) / samples as f64 - out as f64 / 2.0) * scale;236                    let py =237                        (i as f64 + (a as f64 + 0.5) / samples as f64 - out as f64 / 2.0) * scale;238                    for (k, &(cos, sin)) in turns.iter().enumerate() {239                        let x = c + px * cos + py * sin;240                        let y = c - px * sin + py * cos;241                        values[k] = if x >= 0.0 && y >= 0.0 && x < size as f64 && y < size as f64 {242                            data[y as usize * size + x as usize]243                        } else {244                            0.0245                        };246                    }247                    total += blend.fold(&values);248                }249            }250            field.push(total / (samples * samples) as f32);251        }252    }253    Ok(field)254}255256/// The ring profile: the circle means at steps radii spaced evenly from the centre to the corner circle.257///258/// # Errors259///260/// Errors when the sample count is not the side squared.261pub fn profile(data: &[f32], size: usize, steps: usize) -> Result<Vec<f32>> {262    let steps = steps.max(2);263    let far = reach(size);264    (0..steps)265        .map(|k| Ok(ring(data, size, far * k as f64 / (steps - 1) as f64)? as f32))266        .collect()267}268269/// 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.270pub fn wheel(profile: &[f32], size: usize) -> Vec<f32> {271    let last = profile.len().saturating_sub(1);272    if last == 0 {273        return vec![profile.first().copied().unwrap_or(0.0); size * size];274    }275    let c = centre(size);276    let scale = 2.0 * last as f64 / size as f64;277    let mut out = Vec::with_capacity(size * size);278    for row in 0..size {279        for col in 0..size {280            let (dx, dy) = (col as f64 + 0.5 - c, row as f64 + 0.5 - c);281            let t = (dx * dx + dy * dy).sqrt() * scale;282            let i = (t.floor() as usize).min(last);283            let f = (t - i as f64) as f32;284            let value = if i == last {285                profile[last]286            } else {287                profile[i] * (1.0 - f) + profile[i + 1] * f288            };289            out.push(value);290        }291    }292    out293}294295/// The mass a profile carries, the trapezoid integral of `2 pi r F(r)` in cells of the raster it came from.296pub fn mass(profile: &[f32], size: usize) -> f64 {297    let last = profile.len().saturating_sub(1);298    if last == 0 {299        return 0.0;300    }301    let step = reach(size) / last as f64;302    let weight = |k: usize| 2.0 * PI * (k as f64 * step) * profile[k] as f64;303    let inner: f64 = (1..last).map(weight).sum();304    step * (inner + (weight(0) + weight(last)) / 2.0)305}306307/// 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.308///309/// ```310/// let solid = vec![1.0; 64];311/// let rings = mrlyrs::math::spin::profile(&solid, 8, 4000).unwrap();312/// let inner = mrlyrs::math::spin::mass_within(&rings, 8, 3.0);313/// assert!((inner / (9.0 * std::f64::consts::PI) - 1.0).abs() < 1e-3);314/// ```315pub fn mass_within(profile: &[f32], size: usize, radius: f64) -> f64 {316    let last = profile.len().saturating_sub(1);317    if last == 0 || radius <= 0.0 {318        return 0.0;319    }320    let step = reach(size) / last as f64;321    let weight = |k: usize| 2.0 * PI * (k as f64 * step) * profile[k] as f64;322    let full = (radius / step).floor().min(last as f64) as usize;323    let mut total = 0.0;324    for k in 1..=full {325        total += (weight(k - 1) + weight(k)) / 2.0 * step;326    }327    if full < last {328        let rest = radius - full as f64 * step;329        let share = rest / step;330        let edge = weight(full) + (weight(full + 1) - weight(full)) * share;331        total += (weight(full) + edge) / 2.0 * rest;332    }333    total334}335336#[cfg(test)]337mod tests {338    use super::*;339    use crate::math::atoms;340341    fn floats(grid: &crate::core::tensor::Tensor) -> Vec<f32> {342        grid.bytes().unwrap().iter().map(|&b| b as f32).collect()343    }344345    #[test]346    fn a_solid_square_leaves_through_four_arcs() {347        let side = 8usize;348        let solid = floats(&atoms::ones_2d(side));349        for r in [0.5, 2.0, 3.99] {350            assert!((ring(&solid, side, r).unwrap() - 1.0).abs() < 1e-12);351        }352        for r in [4.5, 5.0, 5.5] {353            let expect = 1.0 - 4.0 / PI * (side as f64 / 2.0 / r).acos();354            assert!((ring(&solid, side, r).unwrap() - expect).abs() < 1e-12);355        }356        assert!(ring(&solid, side, reach(side) + 0.01).unwrap().abs() < 1e-12);357    }358359    #[test]360    fn the_carpet_opens_on_a_black_disc() {361        let carpet = atoms::carpet_nd(3, 2).kron(&atoms::carpet_nd(3, 2));362        let data = floats(&carpet);363        assert!(ring(&data, 9, 1.4).unwrap().abs() < 1e-12);364        assert!(ring(&data, 9, 2.0).unwrap() > 0.0);365        assert!((ring(&data, 9, 1.0).unwrap() - 0.0).abs() < 1e-12);366    }367368    #[test]369    fn the_mass_of_the_profile_is_the_fill() {370        let carpet = atoms::carpet_nd(3, 2)371            .kron(&atoms::carpet_nd(3, 2))372            .kron(&atoms::carpet_nd(3, 2));373        let data = floats(&carpet);374        let fills = data.iter().sum::<f32>() as f64;375        assert_eq!(fills, 512.0);376        let rings = profile(&data, 27, 4000).unwrap();377        assert!((mass(&rings, 27) - fills).abs() / fills < 0.002);378    }379380    #[test]381    fn the_spin_mass_scales_by_the_fill_about_a_filled_corner() {382        let tile = atoms::carpet_nd(3, 2);383        let carpet = tile.kron(&tile).kron(&tile).kron(&tile);384        let side = 81usize;385        let wide = 2 * side;386        let mut data = vec![0.0f32; wide * wide];387        let bytes = carpet.bytes().unwrap();388        for row in 0..side {389            for col in 0..side {390                data[(row + side) * wide + col + side] = bytes[row * side + col] as f32;391            }392        }393        let rings = profile(&data, wide, 4000).unwrap();394        assert!(mass_within(&rings, wide, 1e-6).abs() < 1e-9);395        for radius in [12.0f64, 18.0, 27.0] {396            let near = mass_within(&rings, wide, radius);397            let far = mass_within(&rings, wide, 3.0 * radius);398            assert!(near > 0.0);399            assert!((far / (8.0 * near) - 1.0).abs() < 0.08, "radius {radius}");400        }401        let whole = mass_within(&rings, wide, reach(wide));402        assert!((whole - mass(&rings, wide)).abs() / whole < 1e-9);403    }404405    #[test]406    fn the_wheel_reads_the_profile_by_radius() {407        let rings = profile(&floats(&atoms::ones_2d(8)), 8, 64).unwrap();408        let spun = wheel(&rings, 16);409        assert_eq!(spun.len(), 256);410        assert!((spun[8 * 16 + 8] - 1.0).abs() < 1e-6);411        assert!((spun[8 * 16 + 12] - 1.0).abs() < 1e-6);412        assert!(spun[8 * 16 + 15] < 0.2);413        assert_eq!(spun[0], 0.0);414        assert_eq!(profile(&[], 0, 0).unwrap().len(), 2);415    }416417    fn at(field: &[f32], out: usize, size: usize, px: f64, py: f64) -> f32 {418        let scale = 2.0 * reach(size) / out as f64;419        let col = (px / scale + out as f64 / 2.0) as usize;420        let row = (py / scale + out as f64 / 2.0) as usize;421        field[row * out + col]422    }423424    #[test]425    fn two_squares_at_an_eighth_turn_make_a_star() {426        let solid = floats(&atoms::ones_2d(8));427        let stack = |blend| radial(&solid, 8, 64, 2, 0.125, blend, 1).unwrap();428        let probe = |blend, px, py| at(&stack(blend), 64, 8, px, py);429        assert_eq!(probe(Blend::Union, 0.0, 0.0), 1.0);430        assert_eq!(probe(Blend::Union, 0.0, -5.0), 1.0);431        assert_eq!(probe(Blend::Meet, 0.0, -5.0), 0.0);432        assert_eq!(probe(Blend::Mean, 0.0, -5.0), 0.5);433        assert_eq!(probe(Blend::Sum, 0.0, 0.0), 2.0);434        assert_eq!(probe(Blend::Parity, 0.0, -5.0), 1.0);435        assert_eq!(probe(Blend::Parity, 0.0, 0.0), 0.0);436        assert_eq!(probe(Blend::Difference, 0.0, -5.0), 0.0);437        assert_eq!(probe(Blend::Difference, 3.8, -3.8), 1.0);438        assert_eq!(Blend::named("soup"), None);439    }440441    #[test]442    fn the_harmonics_read_the_rotation_order() {443        let square = harmonics(&floats(&atoms::ones_2d(8)), 8, 128, 12).unwrap();444        assert!(square[0] > 0.0);445        assert!(square[4] > 1e-3 * square[0]);446        assert!(square[8] > 1e-3 * square[0]);447        for m in [1, 2, 3, 5, 6, 7, 9, 10, 11] {448            assert!(square[m] < 1e-9 * square[0], "m {m}");449        }450        assert_eq!(turns(&square), 4);451        let mut bar = vec![0.0f32; 16];452        for i in [5, 6, 9, 10, 4, 11] {453            bar[i] = 1.0;454        }455        assert_eq!(turns(&harmonics(&bar, 4, 64, 8).unwrap()), 2);456        let mut blob = vec![0.0f32; 16];457        blob[0] = 1.0;458        assert_eq!(turns(&harmonics(&blob, 4, 64, 8).unwrap()), 1);459        assert_eq!(turns(&[1.0, 0.0, 0.0]), 0);460        assert_eq!(461            (petals(6, 4), petals(8, 4), petals(5, 1), petals(3, 0)),462            (12, 8, 5, 0)463        );464    }465466    #[test]467    fn refuses_a_raster_that_is_not_the_side_squared() {468        let short = vec![1.0f32; 15];469        assert!(arcs(&short, 4, 1.0).is_err());470        assert!(ring(&short, 4, 1.0).is_err());471        assert!(harmonics(&short, 4, 8, 4).is_err());472        assert!(profile(&short, 4, 8).is_err());473        assert!(radial(&short, 4, 4, 2, 0.5, Blend::Mean, 1).is_err());474    }475}