volume.rs

1use super::sample::{membership, pack};2use super::stack::merge;3use super::{Combine, Spec};4use mrlycore::errors::{value_error, Result};5use mrlycore::tensor::Tensor;67/// A cubic grid of f32 samples, x-major.8#[derive(Clone, Debug, PartialEq)]9pub struct Volume {10    /// The samples, x-major, then y, then z.11    pub data: Vec<f32>,12    /// The side in samples.13    pub size: usize,14}1516/// 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`.17#[derive(Clone, Debug, PartialEq)]18pub struct Frame {19    /// The point of the plane the window is centred on.20    pub centre: [f64; 3],21    /// The unit axis the window's columns run along.22    pub u: [f64; 3],23    /// The unit axis the window's rows run along.24    pub v: [f64; 3],25    /// The unit normal.26    pub normal: [f64; 3],27    /// The side of the square window.28    pub width: f64,29}3031fn dot(a: [f64; 3], b: [f64; 3]) -> f64 {32    a[0] * b[0] + a[1] * b[1] + a[2] * b[2]33}3435fn cross(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {36    [37        a[1] * b[2] - a[2] * b[1],38        a[2] * b[0] - a[0] * b[2],39        a[0] * b[1] - a[1] * b[0],40    ]41}4243fn unit(a: [f64; 3]) -> Option<[f64; 3]> {44    let n = dot(a, a).sqrt();45    (n > 1e-12).then(|| [a[0] / n, a[1] / n, a[2] / n])46}4748impl Volume {49    /// Builds a zeroed volume of the side.50    pub fn new(size: usize) -> Volume {51        Volume {52            data: vec![0.0; size * size * size],53            size,54        }55    }56    /// Wraps x-major samples of the side.57    pub fn from_data(data: Vec<f32>, size: usize) -> Result<Volume> {58        if data.len() != size * size * size {59            return value_error("data must be size*size*size.");60        }61        Ok(Volume { data, size })62    }63    /// Returns the smallest sample.64    pub fn min(&self) -> f32 {65        self.data.iter().cloned().fold(f32::INFINITY, f32::min)66    }67    /// Returns the largest sample.68    pub fn max(&self) -> f32 {69        self.data.iter().cloned().fold(f32::NEG_INFINITY, f32::max)70    }71    /// Reads the sample at a voxel.72    pub fn at(&self, x: usize, y: usize, z: usize) -> f32 {73        self.data[(x * self.size + y) * self.size + z]74    }75    /// Reads the voxel a point of the unit cube falls in, or zero outside it.76    pub fn sample(&self, p: [f64; 3]) -> Option<f32> {77        if p.iter().any(|&c| !(0.0..1.0).contains(&c)) {78            return None;79        }80        let s = self.size as f64;81        Some(self.at(82            (p[0] * s) as usize,83            (p[1] * s) as usize,84            (p[2] * s) as usize,85        ))86    }87    /// Thresholds into a byte tensor: one where a sample reaches the level, zero below.88    pub fn solid(&self, level: f32) -> Tensor {89        let mut grid = Tensor::new(vec![self.size; 3]);90        for (site, &v) in grid.bytes_mut().iter_mut().zip(self.data.iter()) {91            *site = (v >= level) as u8;92        }93        grid94    }95    /// Counts the samples at or above the level.96    pub fn count(&self, level: f32) -> usize {97        self.data.iter().filter(|&&v| v >= level).count()98    }99    /// 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.100    pub fn plane(&self, frame: &Frame, out: usize) -> (Vec<f32>, Vec<u8>) {101        let mut values = Vec::with_capacity(out * out);102        let mut inside = Vec::with_capacity(out * out);103        for i in 0..out {104            for j in 0..out {105                let a = ((j as f64 + 0.5) / out as f64 - 0.5) * frame.width;106                let b = ((i as f64 + 0.5) / out as f64 - 0.5) * frame.width;107                let p: Vec<f64> = (0..3)108                    .map(|k| (frame.centre[k] + a * frame.u[k] + b * frame.v[k] + 1.0) / 2.0)109                    .collect();110                match self.sample([p[0], p[1], p[2]]) {111                    Some(v) => {112                        values.push(v);113                        inside.push(1);114                    }115                    None => {116                        values.push(0.0);117                        inside.push(0);118                    }119                }120            }121        }122        (values, inside)123    }124}125126/// 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.127pub fn frame(normal: [f64; 3], offset: f64) -> Result<Frame> {128    let Some(n) = unit(normal) else {129        return value_error("the normal must not be zero.");130    };131    let seed = if n[0].abs() < 0.9 {132        [1.0, 0.0, 0.0]133    } else {134        [0.0, 1.0, 0.0]135    };136    let u = unit(cross(seed, n)).unwrap();137    let v = cross(n, u);138    let corners = (0..8).map(|c| {139        [140            if c & 1 == 0 { -1.0 } else { 1.0 },141            if c & 2 == 0 { -1.0 } else { 1.0 },142            if c & 4 == 0 { -1.0 } else { 1.0 },143        ]144    });145    let mut span = [[f64::INFINITY, f64::NEG_INFINITY]; 3];146    for c in corners {147        for (k, axis) in [u, v, n].iter().enumerate() {148            let t = dot(c, *axis);149            span[k][0] = span[k][0].min(t);150            span[k][1] = span[k][1].max(t);151        }152    }153    let mid = |k: usize| (span[k][0] + span[k][1]) / 2.0;154    let depth = span[2][0] + offset.clamp(0.0, 1.0) * (span[2][1] - span[2][0]);155    let centre = [156        mid(0) * u[0] + mid(1) * v[0] + depth * n[0],157        mid(0) * u[1] + mid(1) * v[1] + depth * n[1],158        mid(0) * u[2] + mid(1) * v[2] + depth * n[2],159    ];160    Ok(Frame {161        centre,162        u,163        v,164        normal: n,165        width: (span[0][1] - span[0][0]).max(span[1][1] - span[1][0]),166    })167}168169fn layer(spec: Spec, number: usize, level: usize, size: usize) -> Result<Vec<bool>> {170    let Spec {171        code,172        base: q,173        dimension: d,174    } = spec;175    if d != 3 {176        return value_error("a volume needs dimension 3.");177    }178    let table = membership(code, q, 3)?;179    let mut mask = vec![true; size * size * size];180    let inv = 1.0 / size as f64;181    for k in 0..level.max(1) {182        let s = (number * q.pow(k as u32)) as f64;183        let residue =184            |i: usize| ((s * (i as f64 + 0.5) * inv).floor() as i64).rem_euclid(q as i64) as usize;185        let digits: Vec<usize> = (0..size).map(residue).collect();186        for x in 0..size {187            for y in 0..size {188                for z in 0..size {189                    let cell = (x * size + y) * size + z;190                    if mask[cell] && !table[pack(&[digits[x], digits[y], digits[z]], q)] {191                        mask[cell] = false;192                    }193                }194            }195        }196    }197    Ok(mask)198}199200/// Layers one cube design at several side numbers into a volume under the chosen combine.201pub fn volume(202    spec: Spec,203    numbers: &[usize],204    combine: Combine,205    level: usize,206    size: usize,207) -> Result<Volume> {208    if size == 0 {209        return value_error("size must be at least 1.");210    }211    let mut acc = vec![0.0f32; size * size * size];212    let mut first = true;213    for &n in numbers {214        let mask = layer(spec, n, level, size)?;215        merge(&mut acc, &mask, combine, first);216        first = false;217    }218    Volume::from_data(acc, size)219}220221#[cfg(test)]222mod tests {223    use super::*;224    use mrlymath::bang::corners_to_code;225226    fn low() -> Spec {227        Spec::new(corners_to_code(&[vec![0, 0, 0]], 3, 2), 2, 3)228    }229230    #[test]231    fn the_low_corner_at_scale_one_fills_the_first_octant() {232        let v = volume(low(), &[1], Combine::Sum, 1, 4).unwrap();233        assert_eq!(v.count(1.0), 64);234        let v = volume(low(), &[2], Combine::Sum, 1, 4).unwrap();235        assert_eq!(v.count(1.0), 8);236        assert_eq!(v.at(0, 0, 0), 1.0);237        assert_eq!(v.at(3, 0, 0), 0.0);238    }239240    #[test]241    fn the_sponge_stack_counts_its_layers() {242        let sponge = Spec::new(23, 2, 3);243        let v = volume(sponge, &[1, 3], Combine::Sum, 1, 9).unwrap();244        assert_eq!((v.min(), v.max()), (1.0, 2.0));245        assert_eq!(v.count(2.0), 20 * 27);246        let x = volume(sponge, &[1, 3], Combine::Xor, 1, 9).unwrap();247        assert_eq!(x.count(1.0), 7 * 27);248        let a = volume(sponge, &[3, 9], Combine::And, 1, 9).unwrap();249        assert_eq!(a.count(1.0), 64 + 240 + 48);250        assert!(volume(Spec::new(1, 2, 2), &[1], Combine::Sum, 1, 4).is_err());251        assert!(volume(sponge, &[1], Combine::Sum, 1, 0).is_err());252    }253254    #[test]255    fn the_diagonal_frame_holds_the_hexagon() {256        let f = frame([1.0, 1.0, 1.0], 0.5).unwrap();257        assert!((dot(f.centre, f.normal)).abs() < 1e-12);258        assert!((dot(f.u, f.v)).abs() < 1e-12);259        assert!((f.width - 2.0 * (8.0f64 / 3.0).sqrt()).abs() < 1e-9);260        let solid = Volume::from_data(vec![1.0; 27], 3).unwrap();261        let (_, inside) = solid.plane(&f, 200);262        let share = inside.iter().map(|&b| b as f64).sum::<f64>() / 40000.0;263        let hexagon = 3.0 * 3f64.sqrt() / (f.width * f.width);264        assert!((share - hexagon).abs() < 0.01, "{share} {hexagon}");265        let x = frame([1.0, 0.0, 0.0], 0.25).unwrap();266        assert!((x.centre[0] + 0.5).abs() < 1e-12);267        assert_eq!(x.width, 2.0);268        assert!(frame([0.0, 0.0, 0.0], 0.5).is_err());269    }270}