volume.rs

1use super::sample::{membership, pack};2use super::stack::merge;3use super::{Combine, Spec};4use crate::core::error::{value_error, Result};5use crate::core::tensor::Tensor;6use serde::{Deserialize, Serialize};78/// A cubic grid of f32 samples, x-major.9#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]10pub struct Volume {11    /// The samples, x-major, then y, then z.12    pub data: Vec<f32>,13    /// The side in samples.14    pub size: usize,15}1617/// A plane through the unit box, framed for sampling: its centre, its two in-plane axes and the width of the square window that holds the whole section, all in the box `[-1, 1]^3`.18#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]19pub struct Frame {20    /// The point of the plane the window is centred on.21    pub centre: [f64; 3],22    /// The unit axis the window's columns run along.23    pub u: [f64; 3],24    /// The unit axis the window's rows run along.25    pub v: [f64; 3],26    /// The unit normal.27    pub normal: [f64; 3],28    /// The side of the square window.29    pub width: f64,30}3132fn dot(a: [f64; 3], b: [f64; 3]) -> f64 {33    a[0] * b[0] + a[1] * b[1] + a[2] * b[2]34}3536fn cross(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {37    [38        a[1] * b[2] - a[2] * b[1],39        a[2] * b[0] - a[0] * b[2],40        a[0] * b[1] - a[1] * b[0],41    ]42}4344fn unit(a: [f64; 3]) -> Option<[f64; 3]> {45    let n = dot(a, a).sqrt();46    (n > 1e-12).then(|| [a[0] / n, a[1] / n, a[2] / n])47}4849impl Volume {50    /// Builds a zeroed volume of the side.51    pub fn new(size: usize) -> Volume {52        Volume {53            data: vec![0.0; size * size * size],54            size,55        }56    }57    /// Wraps x-major samples of the side.58    ///59    /// # Errors60    ///61    /// Errors when the count is not the side cubed.62    pub fn from_data(data: Vec<f32>, size: usize) -> Result<Volume> {63        if data.len() != size * size * size {64            return value_error("data must be size*size*size.");65        }66        Ok(Volume { data, size })67    }68    /// Returns the smallest sample.69    pub fn min(&self) -> f32 {70        self.data.iter().cloned().fold(f32::INFINITY, f32::min)71    }72    /// Returns the largest sample.73    pub fn max(&self) -> f32 {74        self.data.iter().cloned().fold(f32::NEG_INFINITY, f32::max)75    }76    /// Reads the sample at a voxel.77    pub fn at(&self, x: usize, y: usize, z: usize) -> f32 {78        self.data[(x * self.size + y) * self.size + z]79    }80    /// Reads the voxel a point of the unit cube falls in, or zero outside it.81    pub fn sample(&self, p: [f64; 3]) -> Option<f32> {82        if p.iter().any(|&c| !(0.0..1.0).contains(&c)) {83            return None;84        }85        let s = self.size as f64;86        Some(self.at(87            (p[0] * s) as usize,88            (p[1] * s) as usize,89            (p[2] * s) as usize,90        ))91    }92    /// Thresholds into a byte tensor: one where a sample reaches the level, zero below.93    pub fn solid(&self, level: f32) -> Tensor {94        let mut grid = Tensor::new(vec![self.size; 3]);95        for (flat, &v) in self.data.iter().take(grid.size()).enumerate() {96            grid.put(flat, i64::from(v >= level));97        }98        grid99    }100    /// Counts the samples at or above the level.101    pub fn count(&self, level: f32) -> usize {102        self.data.iter().filter(|&&v| v >= level).count()103    }104    /// Samples the plane of the frame on an out by out window: the values row by row, and one byte per pixel saying whether it lies inside the cube.105    pub fn plane(&self, frame: &Frame, out: usize) -> (Vec<f32>, Vec<u8>) {106        let mut values = Vec::with_capacity(out * out);107        let mut inside = Vec::with_capacity(out * out);108        for i in 0..out {109            for j in 0..out {110                let a = ((j as f64 + 0.5) / out as f64 - 0.5) * frame.width;111                let b = ((i as f64 + 0.5) / out as f64 - 0.5) * frame.width;112                let p: Vec<f64> = (0..3)113                    .map(|k| (frame.centre[k] + a * frame.u[k] + b * frame.v[k] + 1.0) / 2.0)114                    .collect();115                match self.sample([p[0], p[1], p[2]]) {116                    Some(v) => {117                        values.push(v);118                        inside.push(1);119                    }120                    None => {121                        values.push(0.0);122                        inside.push(0);123                    }124                }125            }126        }127        (values, inside)128    }129}130131/// Frames the plane normal to the direction, at the offset from zero to one across the box along it; the window is the smallest square holding every section on that normal.132///133/// # Errors134///135/// Errors on a zero normal.136pub fn frame(normal: [f64; 3], offset: f64) -> Result<Frame> {137    let Some(n) = unit(normal) else {138        return value_error("the normal must not be zero.");139    };140    let seed = if n[0].abs() < 0.9 {141        [1.0, 0.0, 0.0]142    } else {143        [0.0, 1.0, 0.0]144    };145    let u = unit(cross(seed, n)).unwrap();146    let v = cross(n, u);147    let corners = (0..8).map(|c| {148        [149            if c & 1 == 0 { -1.0 } else { 1.0 },150            if c & 2 == 0 { -1.0 } else { 1.0 },151            if c & 4 == 0 { -1.0 } else { 1.0 },152        ]153    });154    let mut span = [[f64::INFINITY, f64::NEG_INFINITY]; 3];155    for c in corners {156        for (k, axis) in [u, v, n].iter().enumerate() {157            let t = dot(c, *axis);158            span[k][0] = span[k][0].min(t);159            span[k][1] = span[k][1].max(t);160        }161    }162    let mid = |k: usize| (span[k][0] + span[k][1]) / 2.0;163    let depth = span[2][0] + offset.clamp(0.0, 1.0) * (span[2][1] - span[2][0]);164    let centre = [165        mid(0) * u[0] + mid(1) * v[0] + depth * n[0],166        mid(0) * u[1] + mid(1) * v[1] + depth * n[1],167        mid(0) * u[2] + mid(1) * v[2] + depth * n[2],168    ];169    Ok(Frame {170        centre,171        u,172        v,173        normal: n,174        width: (span[0][1] - span[0][0]).max(span[1][1] - span[1][0]),175    })176}177178fn layer(spec: Spec, number: usize, level: usize, size: usize) -> Result<Vec<bool>> {179    let Spec {180        code,181        base: q,182        dimension: d,183    } = spec;184    if d != 3 {185        return value_error("a volume needs dimension 3.");186    }187    let table = membership(code, q, 3)?;188    let mut mask = vec![true; size * size * size];189    let inv = 1.0 / size as f64;190    for k in 0..level.max(1) {191        let s = (number * q.pow(k as u32)) as f64;192        let residue =193            |i: usize| ((s * (i as f64 + 0.5) * inv).floor() as i64).rem_euclid(q as i64) as usize;194        let digits: Vec<usize> = (0..size).map(residue).collect();195        for x in 0..size {196            for y in 0..size {197                for z in 0..size {198                    let cell = (x * size + y) * size + z;199                    if mask[cell] && !table[pack(&[digits[x], digits[y], digits[z]], q)] {200                        mask[cell] = false;201                    }202                }203            }204        }205    }206    Ok(mask)207}208209/// Layers one cube design at several side numbers into a volume under the chosen combine.210///211/// # Errors212///213/// Errors at size zero, or on a code out of range.214pub fn volume(215    spec: Spec,216    numbers: &[usize],217    combine: Combine,218    level: usize,219    size: usize,220) -> Result<Volume> {221    if size == 0 {222        return value_error("size must be at least 1.");223    }224    let mut acc = vec![0.0f32; size * size * size];225    let mut first = true;226    for &n in numbers {227        let mask = layer(spec, n, level, size)?;228        merge(&mut acc, &mask, combine, first);229        first = false;230    }231    Volume::from_data(acc, size)232}233234#[cfg(test)]235mod tests {236    use super::*;237    use crate::math::bang::corners_to_code;238239    fn low() -> Spec {240        Spec::new(corners_to_code(&[vec![0, 0, 0]], 3, 2).get(), 2, 3)241    }242243    #[test]244    fn the_low_corner_at_scale_one_fills_the_first_octant() {245        let v = volume(low(), &[1], Combine::Sum, 1, 4).unwrap();246        assert_eq!(v.count(1.0), 64);247        let v = volume(low(), &[2], Combine::Sum, 1, 4).unwrap();248        assert_eq!(v.count(1.0), 8);249        assert_eq!(v.at(0, 0, 0), 1.0);250        assert_eq!(v.at(3, 0, 0), 0.0);251    }252253    #[test]254    fn the_sponge_stack_counts_its_layers() {255        let sponge = Spec::new(23, 2, 3);256        let v = volume(sponge, &[1, 3], Combine::Sum, 1, 9).unwrap();257        assert_eq!((v.min(), v.max()), (1.0, 2.0));258        assert_eq!(v.count(2.0), 20 * 27);259        let x = volume(sponge, &[1, 3], Combine::Xor, 1, 9).unwrap();260        assert_eq!(x.count(1.0), 7 * 27);261        let a = volume(sponge, &[3, 9], Combine::And, 1, 9).unwrap();262        assert_eq!(a.count(1.0), 64 + 240 + 48);263        assert!(volume(Spec::new(1, 2, 2), &[1], Combine::Sum, 1, 4).is_err());264        assert!(volume(sponge, &[1], Combine::Sum, 1, 0).is_err());265    }266267    #[test]268    fn the_diagonal_frame_holds_the_hexagon() {269        let f = frame([1.0, 1.0, 1.0], 0.5).unwrap();270        assert!((dot(f.centre, f.normal)).abs() < 1e-12);271        assert!((dot(f.u, f.v)).abs() < 1e-12);272        assert!((f.width - 2.0 * (8.0f64 / 3.0).sqrt()).abs() < 1e-9);273        let solid = Volume::from_data(vec![1.0; 27], 3).unwrap();274        let (_, inside) = solid.plane(&f, 200);275        let share = inside.iter().map(|&b| b as f64).sum::<f64>() / 40000.0;276        let hexagon = 3.0 * 3f64.sqrt() / (f.width * f.width);277        assert!((share - hexagon).abs() < 0.01, "{share} {hexagon}");278        let x = frame([1.0, 0.0, 0.0], 0.25).unwrap();279        assert!((x.centre[0] + 0.5).abs() < 1e-12);280        assert_eq!(x.width, 2.0);281        assert!(frame([0.0, 0.0, 0.0], 0.5).is_err());282    }283}