Initial commit. WIP blur implementation. Grid struct is tentatively ready. Model struct is in its nascency.

src/blur.rs

@@ -0,0 +1,124 @@
+/// Approximate 1D Gaussian filter of standard deviation sigma with N box filter passes. Each
+/// element in the output array contains the radius of the box filter for the corresponding pass.
+pub fn boxes_for_gaussian<const N: usize>(sigma: f32) -> ([usize; N]) {
+    let w_ideal = (12.0 * sigma * sigma / N as f32 + 1.0).sqrt();
+    let mut w = w_ideal as usize;
+    w -= 1 - w & 1;
+
+    let mut m = ((w * w + 4 * w + 3) * N) as f32;
+    m -= 12.0 * sigma * sigma;
+    m *= 0.25;
+    m /= (w + 1) as f32;
+    let m = m.round() as usize;
+
+    let mut result = [0; N];
+    for (i, value) in result.iter_mut().enumerate() {
+        *value = (if i < m { w - 1 } else { w + 1 }) / 2;
+    }
+    result
+}
+
+/// Blur an image with 3 box filter passes. The result will be written to the src slice, while the
+/// buf slice is used as a scratch space.
+pub fn approximate_gauss_blur(
+    src: &mut [f32],
+    buf: &mut [f32],
+    width: usize,
+    height: usize,
+    sigma: f32,
+    decay: f32,
+) {
+    let boxes = boxes_for_gaussian::<3>(sigma);
+    box_blur(src, buf, width, height, boxes[0], 1.0);
+    box_blur(src, buf, width, height, boxes[1], 1.0);
+    box_blur(src, buf, width, height, boxes[2], decay);
+}
+
+/// Perform one pass of the 2D box filter of the given radius. The result will be written to the src
+/// slice, while the buf slice is used as a scratch space.
+fn box_blur(
+    src: &mut [f32],
+    buf: &mut [f32],
+    width: usize,
+    height: usize,
+    radius: usize,
+    decay: f32,
+) {
+    box_blur_h(src, buf, width, height, radius, 1.0);
+    box_blur_v(buf, src, width, height, radius, decay);
+}
+
+/// Perform one pass of the 1D box filter of the given radius along x axis. Applies the decay factor
+/// to the destination buffer.
+fn box_blur_h(
+    src: &[f32],
+    dst: &mut [f32],
+    width: usize,
+    height: usize,
+    radius: usize,
+    decay: f32,
+) {
+    let weight = decay / (2 * radius + 1) as f32;
+
+    // TODO: Parallelize with rayon
+    for i in 0..height {
+        // First we build a value for the beginning of each row. We assume periodic boundary
+        // conditions, so we need to push the left index to the opposite side of the row.
+        let mut value = src[(i + 1) * width - radius - 1];
+        for j in 0..radius {
+            value += src[(i + 1) * width - radius + j] + src[i * width + j];
+        }
+        // At this point "value" contains the unweighted sum for the right-most row element.
+
+        for current_id in i * width..(i + 1) * width {
+            let left_id = ((current_id + width - radius - 1) & (width - 1)) + i * width;
+            let right_id = ((current_id + radius) & (width - 1)) + i * width;
+            value += src[right_id] - src[left_id];
+            dst[current_id] = value * weight;
+        }
+    }
+}
+
+/// Perform one pass of the 1D box filter of the given radius along y axis. Applies the decay factor
+/// to the destination buffer.
+fn box_blur_v(
+    src: &[f32],
+    dst: &mut [f32],
+    width: usize,
+    height: usize,
+    radius: usize,
+    decay: f32,
+) {
+    let weight = decay / (2 * radius + 1) as f32;
+
+    // TODO: Parallelize with rayon
+    for i in 0..width {
+        // First we build a value for the beginning of each column. We assume periodic boundary
+        // conditions, so we need to push the bottom index to the opposite side of the column.
+        let mut value = src[i + (height - radius - 1) * width];
+        for j in 0..radius {
+            value += src[i + (height - radius + j) * width] + src[i + j * width];
+        }
+        // At this point "value" contains the unweighted sum for the top-most column element.
+
+        for current_id in (i..i + height * width).step_by(width) {
+            let bottom_id = (current_id + (height - radius - 1) * width) & (width * height - 1);
+            let top_id = (current_id + radius * width) & (width * height - 1);
+            value += src[top_id] - src[bottom_id];
+            dst[current_id] = value * weight;
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_blur() {
+        let src = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0];
+        let mut dst = vec![0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
+        box_blur_v(&src, &mut dst, 2, 4, 1, 1.0);
+        println!("Out: {:?}", dst);
+    }
+}

src/grid.rs

@@ -0,0 +1,95 @@
+use crate::blur::approximate_gauss_blur;
+use rand::{distributions::Uniform, Rng};
+
+/// A 2D grid with a scalar value per each grid block.
+#[derive(Debug)]
+pub struct Grid {
+    width: usize,
+    height: usize,
+    data: Vec<f32>,
+    buf: Vec<f32>,
+}
+
+#[inline(always)]
+fn is_power_of_two(x: usize) -> bool {
+    (x & (x - 1)) == 0
+}
+
+impl Grid {
+    /// Create a new grid filled with random floats in the [0.0..1.0) range.
+    pub fn new(width: usize, height: usize) -> Self {
+        if !is_power_of_two(width) || !is_power_of_two(height) {
+            panic!("Grid dimensitions must be a power of two.");
+        }
+        let rng = rand::thread_rng();
+        let range = Uniform::from(0.0..1.0);
+        let data = rng.sample_iter(range).take(width * height).collect();
+
+        Grid {
+            width,
+            height,
+            data,
+            buf: vec![0.0; width * height],
+        }
+    }
+
+    /// Truncate x and y and return a corresponding index into the data slice.
+    fn index(&self, x: f32, y: f32) -> usize {
+        let i = (x as usize + self.width) & (self.width - 1);
+        let j = (y as usize + self.height) & (self.height - 1);
+        j * self.width + i
+    }
+
+    /// Get the data value at a given position. The implementation effectively treats data as
+    /// periodic, hence any finite position will produce a value.
+    pub fn get(&self, x: f32, y: f32) -> f32 {
+        self.data[self.index(x, y)]
+    }
+
+    /// Get the buffer value at a given position. The implementation effectively treats data as
+    /// periodic, hence any finite position will produce a value.
+    pub fn get_buf(&self, x: f32, y: f32) -> f32 {
+        self.buf[self.index(x, y)]
+    }
+
+    /// Add a value to the grid data at a given position.
+    pub fn add(&mut self, x: f32, y: f32, value: f32) {
+        let idx = self.index(x, y);
+        self.data[idx] += value
+    }
+
+    /// Diffuse grid data and apply a decay multiplier.
+    pub fn diffuse(&mut self, radius: usize, decay_factor: f32) {
+        approximate_gauss_blur(
+            &mut self.data,
+            &mut self.buf,
+            self.width,
+            self.height,
+            radius as f32,
+            decay_factor,
+        );
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    #[should_panic]
+    fn test_grid_new_panics() {
+        let _ = Grid::new(5, 5);
+    }
+
+    #[test]
+    fn test_grid_new() {
+        let grid = Grid::new(8, 8);
+        assert_eq!(grid.index(0.5, 0.6), 0);
+        assert_eq!(grid.index(1.5, 0.6), 1);
+        assert_eq!(grid.index(0.5, 1.6), 8);
+        assert_eq!(grid.index(2.5, 0.6), 2);
+        assert_eq!(grid.index(2.5, 1.6), 10);
+        assert_eq!(grid.index(7.9, 7.9), 63);
+        assert_eq!(grid.index(-0.5, -0.6), 0);
+    }
+}

src/main.rs

@@ -0,0 +1,8 @@
+mod blur;
+mod grid;
+mod model;
+
+fn main() {
+    let boxes = blur::boxes_for_gaussian::<3>(2.5);
+    println!("boxes: {:?}", boxes);
+}

src/model.rs

@@ -0,0 +1,32 @@
+use crate::grid::Grid;
+
+/// A single Physarum agent. The x and y positions are continuous, hence we use floating point
+/// numbers instead of integers.
+#[derive(Debug)]
+pub struct Agent {
+    pub x: f32,
+    pub y: f32,
+    pub angle: f32,
+}
+
+/// A model configuration. We make it into a separate type, because we will eventually have multiple
+/// configurations in one model.
+#[derive(Debug)]
+pub struct PopulationConfig {
+    pub sensor_angle: f32,
+    pub sensor_distance: f32,
+    pub rotation_angle: f32,
+    pub step_distance: f32,
+    pub decay_factor: f32,
+}
+
+impl PopulationConfig {}
+
+/// Top-level simulation class.
+#[derive(Debug)]
+pub struct Model {
+    grid: Grid,
+    agents: Vec<Agent>,
+
+    config: PopulationConfig,
+}