socchastic somethingorother

Cargo.lock

@@ -129,6 +129,7 @@ version = "0.1.0"
 dependencies = [
  "minifb",
  "num",
+ "rand",
  "rayon",
 ]
 
@@ -230,6 +231,17 @@ dependencies = [
  "slab",
 ]
 
+[[package]]
+name = "getrandom"
+version = "0.2.10"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "be4136b2a15dd319360be1c07d9933517ccf0be8f16bf62a3bee4f0d618df427"
+dependencies = [
+ "cfg-if",
+ "libc",
+ "wasi",
+]
+
 [[package]]
 name = "hermit-abi"
 version = "0.2.6"
@@ -511,6 +523,12 @@ version = "0.3.27"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "26072860ba924cbfa98ea39c8c19b4dd6a4a25423dbdf219c1eca91aa0cf6964"
 
+[[package]]
+name = "ppv-lite86"
+version = "0.2.17"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "5b40af805b3121feab8a3c29f04d8ad262fa8e0561883e7653e024ae4479e6de"
+
 [[package]]
 name = "proc-macro2"
 version = "1.0.59"
@@ -529,6 +547,36 @@ dependencies = [
  "proc-macro2",
 ]
 
+[[package]]
+name = "rand"
+version = "0.8.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "34af8d1a0e25924bc5b7c43c079c942339d8f0a8b57c39049bef581b46327404"
+dependencies = [
+ "libc",
+ "rand_chacha",
+ "rand_core",
+]
+
+[[package]]
+name = "rand_chacha"
+version = "0.3.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "e6c10a63a0fa32252be49d21e7709d4d4baf8d231c2dbce1eaa8141b9b127d88"
+dependencies = [
+ "ppv-lite86",
+ "rand_core",
+]
+
+[[package]]
+name = "rand_core"
+version = "0.6.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ec0be4795e2f6a28069bec0b5ff3e2ac9bafc99e6a9a7dc3547996c5c816922c"
+dependencies = [
+ "getrandom",
+]
+
 [[package]]
 name = "raw-window-handle"
 version = "0.4.3"
@@ -689,6 +737,12 @@ version = "0.1.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "579a42fc0b8e0c63b76519a339be31bed574929511fa53c1a3acae26eb258f29"
 
+[[package]]
+name = "wasi"
+version = "0.11.0+wasi-snapshot-preview1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9c8d87e72b64a3b4db28d11ce29237c246188f4f51057d65a7eab63b7987e423"
+
 [[package]]
 name = "wasm-bindgen"
 version = "0.2.86"

Cargo.toml

@@ -6,4 +6,5 @@ edition = "2021"
 [dependencies]
 minifb = "0.24.0"
 num = "0.4.0"
+rand = "0.8.5"
 rayon = "1.7.0"

src/main.rs

@@ -1,10 +1,11 @@
 #![feature(core_intrinsics)]
-use std::f64::consts::PI;
+use rand::Rng;
+use std::f64::{consts::PI, INFINITY};
 
 use std::time::Duration;
 
 use minifb::{Key, Window, WindowOptions};
-use num::Complex;
+use num::{Complex, Float};
 use rayon::prelude::*;
 
 struct Timer {
@@ -55,20 +56,24 @@ struct Mandlebrot {
     pub view_x: f64,
     pub view_y: f64,
     pub view_zoom: f64,
+    pub supersampling_range: f64,
+    pub supersampling_min_dist: f64,
     pub max_iterations: u32,
 }
 
 impl Mandlebrot {
     fn new() -> Mandlebrot {
-        let height: usize = 800;
+        let height: usize = 1800;
-        let width: usize = 800;
+        let width: usize = 1800;
         Mandlebrot {
             view_width: width,
             view_height: height,
-            view_x: 0.0,
+            view_x: -0.6,
             view_y: 0.0,
-            view_zoom: 1.0,
+            view_zoom: 2.2,
-            max_iterations: 512,
+            supersampling_range: 4.0,
+            supersampling_min_dist: 0.5,
+            max_iterations: 1000,
         }
     }
 
@@ -83,19 +88,137 @@ impl Mandlebrot {
         );
     }
 
-    fn get_pixel(&self, x: u32, y: u32) -> Complex<f64> {
+    fn mandlebrot_sochastic_supersample(&self, x: u32, y: u32, limit: u32) -> Option<u32> {
+        // first we need to make the pisson disc of out sampling point and generate c for them
         let cx = self.view_zoom * (x as f64 / self.view_width as f64 - 0.5) + self.view_x;
         let cy = self.view_zoom * (y as f64 / self.view_height as f64 - 0.5) + self.view_y;
-        Complex { re: cx, im: cy }
+
+        // this is the fast bridson algorithm first outlined here:
+        // https://www.cs.ubc.ca/~rbridson/docs/bridson-siggraph07-poissondisk.pdf
+        let pisson_disc_sample = |cx: f64, cy: f64| -> Vec<Complex<f64>> {
+            let mut ret = Vec::<Complex<f64>>::new();
+            let mut rng = rand::thread_rng();
+
+            let cell_size = self.supersampling_min_dist / 1.41421356237; // sqrt of 2
+            let grid_size = (self.supersampling_range / cell_size).ceil();
+
+            let mut grid = Vec::<(f64, f64)>::with_capacity(
+                (self.supersampling_range * self.supersampling_range) as usize,
+            );
+            let mut proc = Vec::<(f64, f64)>::new(); // treat as a stack
+
+            let in_area = |p: (f64, f64)| -> bool {
+                p.0 > 0.0
+                    && p.0 < self.supersampling_range
+                    && p.1 > 0.0
+                    && p.1 < self.supersampling_range
+            };
+
+            let squared_distance = |a: (f64, f64), b: (f64, f64)| -> f64 {
+                let delta_x = a.0 - b.0;
+                let delta_y = a.1 - b.1;
+                delta_x * delta_x * delta_y * delta_y
+            };
+
+            let point_around = |p: (f64, f64)| -> (f64, f64) {
+                let radius =
+                    self.supersampling_min_dist * (rng.gen_range(0.0..3.0) + 1.0f64).sqrt();
+                let angle = rng.gen_range(0.0..PI * 2.0);
+                (p.0 + angle.cos() * radius, p.1 + angle.sin() * radius)
+            };
+
+            let mut set = |p: (f64, f64)| {
+                let cell = ((p.0 / cell_size) as usize, (p.1 / cell_size) as usize);
+                grid[cell.1 * self.supersampling_range as usize + cell.0] = p;
+            };
+
+            let mut add = |p: (f64, f64)| {
+                let c = Complex { re: p.0, im: p.1 };
+                ret.push(c);
+                proc.push(p);
+                set(p);
+            };
+
+            let mut p_too_close = |p: (f64, f64)| -> bool {
+                let xi = (p.0 / cell_size).floor() as usize;
+                let yi = (p.1 / cell_size).floor() as usize;
+
+                if grid[yi * grid_size as usize + xi].0 != INFINITY {
+                    return true;
+                }
+
+                let min_dist_sq = self.supersampling_min_dist * self.supersampling_min_dist;
+                let minx = (xi - 2).max(0);
+                let miny = (yi - 2).max(0);
+                let maxx = (xi + 2).min(grid_size as usize - 1);
+                let maxy = (yi + 2).min(grid_size as usize - 1);
+
+                for y in miny..=maxy {
+                    for x in minx..=maxx {
+                        let point = grid[y * grid_size as usize + x];
+                        let exists = point.0 != INFINITY;
+                        if exists && squared_distance(p, point) < min_dist_sq {
+                            return true;
+                        }
+                    }
+                }
+
+                false
+            };
+
+            let start_p = (
+                rng.gen_range(0.0..self.supersampling_range),
+                rng.gen_range(0.0..self.supersampling_range),
+            );
+            add(start_p);
+
+            while !proc.is_empty() {
+                let point = proc.pop().unwrap(); // we can safely unwrap here
+                for i in 0..30 {
+                    let p = point_around(point);
+                    if in_area(p) && !p_too_close(p) {
+                        add(p);
+                    }
+                }
+            }
+
+            ret
+        };
+
+        let samples = pisson_disc_sample(cx, cy);
+        let samples: Vec<Option<u32>> = samples
+            .iter()
+            .map(|c| -> Option<u32> {
+                let mut z = Complex { re: 0.0, im: 0.0 };
+
+                for i in 0..limit {
+                    z = z * z + c;
+                    if z.norm() > 2.0 {
+                        return Some(i);
+                    }
+                }
+                None
+            })
+            .collect();
+        let summed_iterations: u32 = samples
+            .iter()
+            .map(|x| if x.is_some() { x.unwrap() } else { 0 })
+            .sum();
+
+        Some(summed_iterations / samples.len() as u32)
     }
 
-    fn mandlebrot(&self, c: Complex<f64>, limit: u32) -> Option<(u32, Complex<f64>)> {
+    fn mandlebrot(&self, x: u32, y: u32, limit: u32) -> Option<u32> {
+        let cx = self.view_zoom * (x as f64 / self.view_width as f64 - 0.5) + self.view_x;
+        let cy = self.view_zoom * (y as f64 / self.view_height as f64 - 0.5) + self.view_y;
+
+        let c = Complex { re: cx, im: cy };
         let mut z = Complex { re: 0.0, im: 0.0 };
 
         for i in 0..limit {
             z = z * z + c;
             if z.norm() > 2.0 {
-                return Some((i, z));
+                return Some(i);
             }
         }
         None
@@ -107,10 +230,9 @@ impl Mandlebrot {
             (scaled_v, scaled_v, scaled_v)
         };
 
-        let c = self.get_pixel(x, y);
+        match self.mandlebrot(x, y, self.max_iterations) {
-        match self.mandlebrot(c, self.max_iterations) {
             None => 0,
-            Some((count, _z)) => {
+            Some(count) => {
                 let col = cyclic_shading(count, self.max_iterations);
                 (col.0 as u32) << 16 | (col.1 as u32) << 8 | col.2 as u32
             }