Implement a basic Sierpinski Gasket IFS

Reviewed-on: #1

Cargo.toml

@@ -21,6 +21,10 @@ spirv-std = { git = "https://github.com/Rust-GPU/rust-gpu.git", rev = "67f1ff2"
 
 anyhow = "1.0.102"
 bytemuck = { version = "1.25.0", features = ["derive"] }
-glam = { version = "0.33.1", default-features = false, features = ["libm"] }
+glam = { version = "0.33.1", default-features = false, features = ["bytemuck", "scalar-math"] }
+image = { version = "0.25.10", default-features = false, features = ["default-formats"]}
+libm = "0.2.16"
+rand = { version = "0.10.1", default-features = false }
 rspirv = "0.13.0"
-
+rand_xoshiro = "0.8.1"
+tempfile = "3.27.0"

enkou-shaders/Cargo.toml

@@ -10,6 +10,14 @@ repository.workspace = true
 workspace = true
 
 [dependencies]
-spirv-std.workspace = true
-glam.workspace = true
 bytemuck.workspace = true
+glam.workspace = true
+libm.workspace = true
+rand.workspace = true
+spirv-std.workspace = true
+
+[dev-dependencies]
+anyhow.workspace = true
+image.workspace = true
+rand_xoshiro.workspace = true
+tempfile.workspace = true

enkou-shaders/examples/gasket.rs

@@ -0,0 +1,75 @@
+use anyhow::{Context, Result};
+use enkou_shaders::Coefficients2;
+use enkou_shaders::camera::Camera;
+use enkou_shaders::chaos_game::ChaosGame;
+use enkou_shaders::transform::Transform;
+use glam::{Affine2, Vec2, uvec2, UVec2};
+use image::{GrayImage, Luma};
+use rand::SeedableRng;
+use rand_xoshiro::Xoshiro256StarStar;
+use std::mem;
+use std::process::Command;
+use tempfile::{NamedTempFile};
+
+const ITERATIONS: u32 = 50_000;
+const ITERATIONS_DISCARD: u32 = 20;
+const IMAGE_DIMENSION: UVec2 = uvec2(600, 600);
+
+pub fn main() -> Result<()> {
+    let seed: u64 = 4; // chosen by fair dice roll
+    let mut rng = Xoshiro256StarStar::seed_from_u64(seed);
+
+    let transforms = [
+        // F_0: (x / 2, y / 2)
+        Transform::new(Affine2::from_coefficients(0.5, 0.0, 0.0, 0.0, 0.5, 0.0)),
+        // F_1: ((x + 1) / 2, y / 2)
+        Transform::new(Affine2::from_coefficients(0.5, 0.0, 0.5, 0.0, 0.5, 0.0)),
+        // F_2: (x / 2, (y + 1) / 2)
+        Transform::new(Affine2::from_coefficients(0.5, 0.0, 0.0, 0.0, 0.5, 0.5)),
+    ];
+
+    let weights = [1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0];
+
+    let mut image = GrayImage::new(IMAGE_DIMENSION.x, IMAGE_DIMENSION.y);
+
+    // The gasket is defined on the range [0, 1] for both X and Y
+    let camera = Camera::new(
+        IMAGE_DIMENSION,
+        Vec2::ONE * 0.5,
+        0.0,
+        Vec2::ZERO,
+        IMAGE_DIMENSION.as_vec2(),
+    );
+
+    let mut chaos_game = ChaosGame::new(&mut rng, &transforms, &weights);
+    for i in 0..ITERATIONS {
+        let next_point = chaos_game.next().unwrap();
+
+        if i < ITERATIONS_DISCARD {
+            continue;
+        }
+
+        if let Some(next_point) = camera.transform_point_to_image(next_point) {
+            image.put_pixel(next_point.x, next_point.y, Luma([255u8]))
+        }
+    }
+
+    let temp = NamedTempFile::with_suffix(".png").context("Unable to create file for image")?;
+    image.save(temp.path()).context("Unable to save image")?;
+
+    let open_program = cfg_select! {
+        unix => "xdg-open",
+        _ => panic!("Unknown system"),
+    };
+
+    Command::new(open_program)
+        .arg(temp.path())
+        .spawn()?
+        .wait()?;
+
+    // In case the image viewer forks and gives control back prior to reading the file,
+    // drop it and don't run the constructor
+    mem::forget(temp);
+
+    Ok(())
+}

enkou-shaders/src/camera.rs

@@ -0,0 +1,175 @@
+use bytemuck::{Pod, Zeroable};
+use glam::{Affine2, IVec2, UVec2, Vec2, vec2};
+use libm::powf;
+
+#[derive(Copy, Clone, Pod, Zeroable)]
+#[repr(C)]
+pub struct Camera {
+    dimensions: UVec2,
+    transform: Affine2,
+}
+
+impl Camera {
+    /// Construct a new camera that maps IFS coordinates to pixel coordinates.
+    ///
+    /// The camera object is itself an affine transformation, but it's helpful to express
+    /// the parameters in individual steps, and compose them internally.
+    ///
+    /// # Arguments
+    ///
+    /// * `dimensions` - Width and height of the output image (in pixels).
+    /// * `center` - Location of the origin in IFS coordinates. Positive `x` shifts the image
+    ///     left, and positive `y` position shifts the image up.
+    /// * `rotate` - Rotation angle (in radians) of IFS coordinates. Rotation is applied after the
+    ///     `center` translation, so it is about the new origin.
+    /// * `zoom` - Zoom factor applied to IFS coordinates. IFS coordinates are scaled by
+    ///     `pow(2, zoom)`, so a zoom factor of 0 is the identity.
+    /// * `scale` - Pixels per unit of IFS coordinates. By default, this parameter is chosen such
+    ///     that the largest dimension will cover the range `[-2, 2]`, but values higher or lower
+    ///     can be used as a secondary zoom.
+    pub fn new(dimensions: UVec2, center: Vec2, rotate: f32, zoom: Vec2, scale: Vec2) -> Camera {
+        let ifs_center_transform = Affine2::from_translation(-center);
+        let zoom_transform = Affine2::from_scale(vec2(powf(2.0, zoom.x), powf(2.0, zoom.y)));
+        let scale_transform = Affine2::from_scale(scale);
+        let rotate_transform = Affine2::from_angle(rotate);
+        let image_center_transform = Affine2::from_translation((dimensions / 2).as_vec2());
+
+        let transform = image_center_transform
+            * rotate_transform
+            * scale_transform
+            * zoom_transform
+            * ifs_center_transform;
+
+        Camera {
+            dimensions,
+            transform,
+        }
+    }
+
+    /// Map a point from IFS coordinates to pixel coordinates.
+    ///
+    /// ```
+    /// # use glam::{vec2, ivec2, uvec2, Vec2};
+    /// # use crate::enkou_shaders::camera::Camera;
+    /// // Output image is 600x600 pixels, centered at the origin, no rotation, no zoom,
+    /// // and scaled such that it covers the range [-2, 2].
+    /// // Use the origin as the IFS coordinate, so the pixel coordinate is the center of the image
+    /// let camera = Camera::new(
+    ///     uvec2(600, 600),
+    ///     Vec2::ZERO,
+    ///     0.0,
+    ///     Vec2::ZERO,
+    ///     vec2(150.0, 150.0)
+    /// );
+    /// assert_eq!(camera.transform_point(vec2(0.0, 0.0)), ivec2(300, 300));
+    /// ```
+    pub fn transform_point(&self, point: Vec2) -> IVec2 {
+        self.transform.transform_point2(point).as_ivec2()
+    }
+
+    /// Map a point from IFS coordinates to pixel coordinates (like [`transform_point`]),
+    /// and check that the result is within the provided image dimensions.
+    pub fn transform_point_to_image(&self, point: Vec2) -> Option<UVec2> {
+        let pixel_coordinates = self.transform_point(point);
+        if pixel_coordinates.x < 0
+            || pixel_coordinates.y < 0
+            || (pixel_coordinates.x as u32) >= self.dimensions.x
+            || (pixel_coordinates.y as u32) >= self.dimensions.y
+        {
+            None
+        } else {
+            Some(pixel_coordinates.as_uvec2())
+        }
+    }
+}
+
+#[cfg(test)]
+mod test {
+    use crate::camera::Camera;
+    use glam::{Affine2, Vec2, ivec2, uvec2, vec2};
+    use libm::powf;
+
+    #[test]
+    pub fn manual_camera() {
+        let starting_point = vec2(1.0, 1.0);
+
+        // Move the origin; points move right and up by one unit, giving us (2.0, 2.0)
+        let center = vec2(-1.0, -1.0);
+        let point = starting_point - center;
+
+        // Rotate about the new origin; points move counter-clockwise, giving us (-2.0, 2.0)
+        let rotate = 90.0f32.to_radians();
+        let point = Affine2::from_angle(rotate).transform_point2(point);
+
+        // Zoom in by a factor of 1; points will be twice as far from the origin,
+        // giving us (-4.0, 4.0)
+        let zoom = vec2(1.0, 1.0);
+        let point = point * vec2(powf(2.0, zoom.x), powf(2.0, zoom.y));
+
+        // Apply scaling; scale 100 in a 1000 x 1000 image is an effective range
+        // of [-5, 5] in IFS coordinates.
+        // After scaling, the point is (-400.0, 400.0)
+        let scale = vec2(100.0, 100.0);
+        let point = point * scale;
+
+        // Move the origin from (0, 0) to image center,
+        // giving us (100.0, 900.0)
+        let dimensions = uvec2(1000, 1000);
+        let point = point.as_ivec2() + dimensions.as_ivec2() / 2;
+
+        // Check that the camera implementation ends up at the same point
+        let camera = Camera::new(dimensions, center, rotate, zoom, scale);
+
+        // The camera is implemented by composing affine transforms,
+        // which ends up with a slightly different result because of rounding.
+        let error = camera.transform_point(starting_point) - point;
+        assert!(error.x.abs() <= 1);
+        assert!(error.y.abs() <= 1);
+    }
+
+    #[test]
+    pub fn point_outside_camera() {
+        // Scale 250 for an image 1000 x 1000 gives an effective range of [-2, 2]
+        let camera = Camera::new(
+            uvec2(1000, 1000),
+            Vec2::ZERO,
+            0.0,
+            Vec2::ZERO,
+            vec2(250.0, 250.0),
+        );
+
+        // Converting a point outside the effective range is legal, but outside the image bounds
+        assert_eq!(camera.transform_point(vec2(3.0, 3.0)), ivec2(1250, 1250));
+    }
+
+    #[test]
+    pub fn point_outside_camera_negative() {
+        // Scale 250 for an image 1000 x 1000 gives an effective range of [-2, 2]
+        let camera = Camera::new(
+            uvec2(1000, 1000),
+            Vec2::ZERO,
+            0.0,
+            Vec2::ZERO,
+            vec2(250.0, 250.0),
+        );
+
+        // Converting a point outside the effective range is legal, but outside the image bounds
+        assert_eq!(camera.transform_point(vec2(-3.0, -3.0)), ivec2(-250, -250));
+    }
+
+    #[test]
+    pub fn aspect_ratio() {
+        // Scale 100 for an image 1600 x 900 gives an effective X range of [-8, 8],
+        // and effective Y range of [-4.5, 4.5]
+        let camera = Camera::new(
+            uvec2(1600, 900),
+            Vec2::ZERO,
+            0.0,
+            Vec2::ZERO,
+            vec2(100.0, 100.0),
+        );
+
+        // This point is inside the image width, but outside its height
+        assert_eq!(camera.transform_point(vec2(6.0, 6.0)), ivec2(1400, 1050));
+    }
+}

enkou-shaders/src/chaos_game.rs

@@ -0,0 +1,68 @@
+use glam::{vec2, Vec2};
+use rand::distr::{Distribution, StandardUniform};
+use rand::{Rng, RngExt};
+use crate::transform::Transform;
+
+struct BiUnit;
+impl Distribution<f32> for BiUnit {
+    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f32 {
+        rng.sample::<f32, _>(StandardUniform) * 2.0 - 1.0
+    }
+}
+
+/// Iterate one step in the chaos game; choose the next transform, apply it,
+/// and return the resulting point. Also returns the transform index so that
+/// path-dependent weights (the "Xaos" table in Apophysis) can be chosen
+/// for the next iteration step.
+///
+/// # Arguments
+///
+/// * `weights` - Weights are assumed to be normalized; adding all elements together should return the value 1
+pub fn step_chaos_game<R: Rng>(
+    point: Vec2,
+    rng: &mut R,
+    transforms: &[Transform],
+    weights: &[f32],
+) -> (Vec2, u32) {
+    let mut choice_weight = rng.sample::<f32, _>(StandardUniform);
+    let mut transform_index: u32 = 0;
+
+    for weight in weights {
+        choice_weight -= weight;
+        if choice_weight <= 0.0 {
+            break;
+        }
+
+        transform_index += 1;
+    }
+
+    (
+        transforms[transform_index as usize].transform_point(point),
+        transform_index,
+    )
+}
+
+pub struct ChaosGame<'a, R: Rng> {
+    current_point: Vec2,
+    rng: &'a mut R,
+    transforms: &'a [Transform],
+    weights: &'a [f32],
+}
+
+impl<'a, R: Rng> ChaosGame<'a, R> {
+    pub fn new(rng: &'a mut R, transforms: &'a [Transform], weights: &'a [f32]) -> Self {
+        let current_point = vec2(rng.sample(BiUnit), rng.sample(BiUnit));
+        ChaosGame { current_point, rng, transforms, weights }
+    }
+}
+
+impl<'a, R: Rng> Iterator for ChaosGame<'a, R> {
+    type Item = Vec2;
+
+    fn next(&mut self) -> Option<Self::Item> {
+        let (next_point, _) = step_chaos_game(self.current_point, self.rng, self.transforms, self.weights);
+        self.current_point = next_point;
+
+        Some(next_point)
+    }
+}

enkou-shaders/src/lib.rs

@@ -1,12 +1,80 @@
+//! # Enkou
 #![no_std]
+#![warn(missing_docs)]
+
+pub mod camera;
+pub mod chaos_game;
+pub mod transform;
 
 use bytemuck::{Pod, Zeroable};
 use core::f32::consts::PI;
-use glam::{Vec3, Vec4, vec2, vec3};
+use glam::{Affine2, Vec3, Vec4, vec2, vec3};
 #[cfg(target_arch = "spirv")]
 use spirv_std::num_traits::Float;
 use spirv_std::spirv;
 
+/// Utility trait for [`Affine2`] to convert between `flam3` notation and [`glam`].
+pub trait Coefficients2 {
+    /// Convert affine transformation coefficients to the [`Affine2`] representation.
+    /// Parameters use the following form:
+    ///
+    /// ```text
+    /// (a * x + b * y + c, d * x + e * y + f)
+    /// ```
+    ///
+    /// ```
+    /// # use glam::{Affine2, vec2};
+    /// # use crate::enkou_shaders::Coefficients2;
+    /// let coefs = Affine2::from_coefficients(1.0, 2.0, 3.0, 4.0, 5.0, 6.0);
+    /// let (x, y) = (7.0, 8.0);
+    /// assert_eq!(
+    ///     coefs.transform_point2(vec2(x, y)),
+    ///     vec2(
+    ///         coefs.a() * x + coefs.b() * y + coefs.c(),
+    ///         coefs.d() * x + coefs.e() * y + coefs.f()
+    ///     )
+    /// );
+    /// ```
+    fn from_coefficients(a: f32, b: f32, c: f32, d: f32, e: f32, f: f32) -> Affine2;
+
+    fn a(&self) -> f32;
+    fn b(&self) -> f32;
+    fn c(&self) -> f32;
+    fn d(&self) -> f32;
+    fn e(&self) -> f32;
+    fn f(&self) -> f32;
+}
+
+impl Coefficients2 for Affine2 {
+    fn from_coefficients(a: f32, b: f32, c: f32, d: f32, e: f32, f: f32) -> Affine2 {
+        Affine2::from_cols_array(&[a, d, b, e, c, f])
+    }
+
+    fn a(&self) -> f32 {
+        self.matrix2.x_axis.x
+    }
+
+    fn b(&self) -> f32 {
+        self.matrix2.y_axis.x
+    }
+
+    fn c(&self) -> f32 {
+        self.translation.x
+    }
+
+    fn d(&self) -> f32 {
+        self.matrix2.x_axis.y
+    }
+
+    fn e(&self) -> f32 {
+        self.matrix2.y_axis.y
+    }
+
+    fn f(&self) -> f32 {
+        self.translation.y
+    }
+}
+
 #[derive(Copy, Clone, Pod, Zeroable)]
 #[repr(C)]
 pub struct ShaderConstants {

enkou-shaders/src/transform.rs

@@ -0,0 +1,18 @@
+use bytemuck::{Pod, Zeroable};
+use glam::{Affine2, Vec2};
+
+#[derive(Copy, Clone, Pod, Zeroable)]
+#[repr(C)]
+pub struct Transform {
+    pub coefficients: Affine2,
+}
+
+impl Transform {
+    pub fn new(coefficients: Affine2) -> Self {
+        Transform { coefficients }
+    }
+
+    pub fn transform_point(&self, point: Vec2) -> Vec2 {
+        self.coefficients.transform_point2(point)
+    }
+}