Initial commit

.gitignore

@@ -0,0 +1 @@
+/target

Cargo.toml

@@ -0,0 +1,12 @@
+[package]
+name = "mandelia"
+version = "0.1.0"
+edition = "2024"
+
+[dependencies]
+image = { version = "0.25.6", features = ["png"] }
+rayon = "1.10.0"
+num = "0.4.3"
+nalgebra = "0.33.2"
+anyhow = "1.0.98"
+clap = { version = "4.5.40", features = ["derive"] }

src/main.rs

@@ -0,0 +1,111 @@
+//#![feature(portable_simd)]
+
+use std::path::PathBuf;
+use clap::{Arg, Parser};
+use image::{GrayImage, ImageBuffer, Luma};
+use image::buffer::ConvertBuffer;
+use rayon::iter::ParallelIterator;
+
+mod transform;
+
+#[derive(Parser)]
+pub struct Args {
+    #[clap(short, long, default_value_t = 3840)]
+    pub width: u32,
+    #[clap(short, long, default_value_t = 2160)]
+    pub height: u32,
+
+    #[clap(short, long, default_value_t = 100)]
+    pub maxiter: usize,
+
+    #[clap(short, long, default_value="out.png")]
+    pub output: PathBuf,
+
+    pub transform: Vec<String>
+}
+
+type Float = f32;
+type Point = Complex<Float>;
+use num::complex::Complex;
+pub type Transform = [(Point, Point); 3];
+
+trait ComplexExt<T> {
+    fn norm2(self) -> T;
+}
+
+impl ComplexExt<Float> for Complex<Float> {
+    fn norm2(self) -> Float {
+        self.im * self.im + self.re * self.re
+    }
+}
+
+fn transform_point(transform: &Transform, point: Point) -> (Point, Point) {
+    let start = transform[0].0 * point.re + transform[1].0 * point.im + transform[2].0;
+    let delta = transform[0].1 * point.re + transform[1].1 * point.im + transform[2].1;
+    (start, delta)
+}
+
+pub struct Config {
+    pub transform: Transform,
+    pub width: u32,
+    pub height: u32,
+    pub maxiter: usize,
+    // We represent escape as the square of its magnitude
+}
+
+
+
+impl Config {
+    fn render_point(&self, (base, delta): (Point, Point)) -> f32 {
+        const ESCAPE: Float = 2.;
+        const ESCAPE2: Float = ESCAPE * ESCAPE;
+
+        let mut base = base;
+        let result  =
+            std::iter::successors(Some(base), |i| Some(*i * *i + delta))
+                .take(self.maxiter)
+                .enumerate()
+                .skip_while(|(_, i)| i.norm2() <= ESCAPE2)
+                .next()
+                .map_or(0.,
+                        |(iter, val)|
+                            iter as Float + 1. - val.norm().ln().ln() / ESCAPE.ln()
+                            );
+        result as f32
+    }
+}
+
+fn main() -> anyhow::Result<()> {
+
+    // let transform = [
+    //     (Point::new(0., 0.), Point::new(1., 0.)),
+    //     (Point::new(0., 0.), Point::new(0., 1.)),
+    //     (Point::new(0., 0.), Point::new(0., 0.)),
+    // ];
+
+    let args = <Args as clap::Parser>::parse();
+
+    let config = Config {
+        transform: transform::parse_transforms(args.transform.iter())?,
+        width: args.width,
+        height: args.height,
+        maxiter: args.maxiter,
+    };
+
+    let scale = 2.0 / std::cmp::min(config.width, config.height) as Float;
+
+    let mut img = ImageBuffer::<image::Luma<f32>, _>::from_par_fn(config.width, config.height, |x, y| {
+        let x = (x as Float) * scale - 1.;
+        let y = (y as Float) * scale - 1.;
+        let point = transform_point(&config.transform, Complex::new(x, y));
+        image::Luma([config.render_point(point)])
+    });
+    let img_norm = img.pixels()
+        .map(|p: &Luma<f32>| p.0[0])
+        .reduce(f32::max)
+        .unwrap_or(1.);
+    img.par_pixels_mut().for_each(|Luma([p])| *p /= img_norm);
+    let img: GrayImage = img.convert();
+    img.save(&args.output)?;
+    Ok(())
+}

src/transform.rs

@@ -0,0 +1,124 @@
+use std::cmp::{max, min};
+use num::Num;
+use super::Float;
+type MTransform = nalgebra::Matrix5<Float>;
+
+#[derive(Debug)]
+enum Type {
+    Rotation(u8, u8),
+    Translation(u8),
+    CoordSwap(u8, u8),
+}
+
+fn parse_axis(c: char) -> Option<u8> {
+    Some(match c {
+        'x' => 0,
+        'y' => 1,
+        'X' => 2,
+        'Y' => 3,
+        _ => return None,
+    })
+}
+
+/// We take a type and distance.
+/// We have four axes, x,y,X,Y. The first two adjust the base point; the second two adjust the delta.
+///
+/// Given that, a plane can be defined with two axes, and a translation can be defined with one.
+pub fn parse_step(s: &str) -> Option<MTransform> {
+    if s.len() < 2 {
+        return None;
+    }
+    let mut chars = s.chars();
+    let axis0 = chars.next().and_then(parse_axis)?;
+    let chars0 = chars.clone();
+    let ttyp = if let Some(c) = chars.next().and_then(parse_axis) {
+        if c == axis0 {
+            return None;
+        }
+        Type::Rotation(min(axis0, c), max(axis0, c))
+    } else {
+        chars = chars0.clone();
+        if chars.next() == Some('=') {
+            let axis1 = chars.next().and_then(parse_axis)?;
+            eprintln!("Transform: coord swap {axis0} {axis1}");
+            if chars.next().is_some() {
+                return None;
+            }
+            let mut mat = nalgebra::Matrix5::identity();
+            mat[(axis0 as usize, axis0 as usize)] = 0.0;
+            mat[(axis1 as usize, axis1 as usize)] = 0.0;
+            mat[(axis0 as usize, axis1 as usize)] = 1.0;
+            mat[(axis1 as usize, axis0 as usize)] = 1.0;
+            return Some(mat);
+        }
+        chars = chars0;
+
+
+        Type::Translation(axis0)
+    };
+
+    // take what's left
+    let s = chars.as_str();
+    let delta = f32::from_str_radix(chars.as_str().trim(), 10).ok()?;
+
+    eprintln!("Transform: {:?} {}", ttyp, delta);
+
+    match ttyp {
+        Type::Rotation(a, b) => {
+            let (a0, a1) = match (a,b) {
+                (0,1) => (2,3),
+                (0,2) => (1,3),
+                (0,3) => (1,2),
+                (1,2) => (0,3),
+                (1,3) => (0,2),
+                (2,3) => (0,1),
+                _ => unreachable!("Invalid rotation axes"),
+            };
+            let delta = delta.to_radians();
+            let s = delta.sin();
+            let c = delta.cos();
+            let mut m = nalgebra::Matrix5::identity();
+            m[(a0,a0)] = c;
+            m[(a0,a1)] = -s;
+            m[(a1,a0)] = s;
+            m[(a1,a1)] = c;
+            Some(m)
+        }
+        Type::Translation(a) => {
+            let mut m = nalgebra::Matrix5::identity();
+            m[(4, a as usize)] = -delta;
+            Some(m)
+        }
+        Type::CoordSwap(a, b) => {
+            unreachable!()
+        }
+    }
+}
+
+pub fn parse_transforms<'a>(transforms: impl Iterator<Item=impl AsRef<str>>) ->anyhow::Result<super::Transform> {
+    use super::Point;
+    let mat = transforms
+        .map(|s| parse_step(s.as_ref())
+            .ok_or_else(|| anyhow::anyhow!("Invalid step: {:?}", s.as_ref())))
+            .try_fold(nalgebra::Matrix5::identity(),
+                      |acc, new| new.map(|new| acc * new))?;
+
+    let bx = Point::new(mat[(0,0)], mat[(0,1)]);
+    let by = Point::new(mat[(1,0)], mat[(1,1)]);
+    let bd = Point::new(mat[(4,0)], mat[(4,1)]);
+
+    let cx = Point::new(mat[(0,2)], mat[(0,3)]);
+    let cy = Point::new(mat[(1,2)], mat[(1,3)]);
+    let cd = Point::new(mat[(4,2)], mat[(4,3)]);
+
+    Ok([
+        (bx, cx),
+        (by, cy),
+        (bd, cd),
+    ])
+
+}
+
+pub fn parse_transform(s: &str) -> anyhow::Result<super::Transform> {
+    parse_transforms(s.split_whitespace())
+}