Compare commits

..
1 Commits
Author SHA1 Message Date
surv ea2a23999d feat: BLA iteration-skipping for deep-zoom Mandelbrot 2026-09-15 15:50:57 +02:00
26 changed files with 1060 additions and 5748 deletions
-3
View File
@@ -1,6 +1,3 @@
/target
Cargo.lock
dist
frames*
out.mp4
-216
View File
@@ -1,216 +0,0 @@
# CLAUDE.md
This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository.
## What this is
A deep-zoom fractal explorer (Rust + wgpu + egui + WGSL). It zooms past the
~10¹³× limit of plain `f64` using **perturbation theory**: one high-precision
reference orbit is computed on the CPU (arbitrary precision via `dashu-float`),
and every pixel is rendered on the GPU as a cheap `f32` delta from it, with
rebasing to avoid glitches. The `f32` GPU tier reaches roughly 10³⁰×. Runs
natively (Vulkan/Metal/DX12) and in the browser (WebGPU only — WebGL2 can't do
storage buffers, which the fragment shader needs for the reference orbit).
## Commands
```sh
cargo run --release # native, run (release matters: fractal math is hot)
cargo test # reference-orbit math, share-link round-trip, WGSL validation
cargo test --test shader_valid # just the WGSL parse/validate tests (naga, no GPU needed)
cargo clippy
cargo fmt # rustfmt.toml just pins edition = "2024"
```
Web build (WebGPU):
```sh
rustup target add wasm32-unknown-unknown
cargo install wasm-bindgen-cli --version 0.2.128 # must match the wasm-bindgen crate version
./build-web.sh # -> ./dist
python3 -m http.server -d dist 8080
```
Native CLI flags (`src/cli.rs`, applied in `FractalApp::apply_cli`): `--kind`,
`--power`, `--julia re,im`, `--phoenix-p re,im`, `--lambda-l re,im`,
`--palette`, `--share <fragment>`,
`--view re,im,half_height[,iterations]`, `--de`, `--buddhabrot`,
`--buddha-palette`. `--headless` (`src/headless.rs`) skips the window
entirely: it builds the same view from the other flags, creates its own
offscreen wgpu device, and renders straight to a PNG (`--width`/`--height`,
default 1920×1080, `--export-path out.png`) without needing a GPU-backed
window/event loop. Not yet supported with `--buddhabrot`. Run
`mandelbrot --help` for the full list.
`--headless` also has an animation mode, for feeding into `ffmpeg`: add
`--to-view re,im,half_height[,iterations]` (or `--to-share <fragment>`, which
only pulls position/zoom/iterations out of the link) alongside a start view
(`--view`/`--share`/`--kind`/`--julia`), plus `--frames N` or
`--fps`/`--duration`. `--export-path` then names an output *directory* of
`frame-00001.png`, `frame-00002.png`, ... instead of a single file. Only the
camera (center + half-height) is animated — kind, colors, and per-kind
constants stay fixed at whatever the start flags set. `view::interpolate_view`
does the interpolation: half-height geometrically (log-linear, since zoom
spans many decades), center linearly through the complex plane at full
`Big` precision; `--linear` swaps the default smoothstep easing for constant
pacing. Iteration count auto-scales with zoom depth per frame (same
`auto_iteration_count` the interactive app uses while zooming), overriding
any iteration count from `--view`/`--share`/`--to-view`/`--to-share`.
There's no GPU in most sandboxes: `cargo check`/`cargo test --test shader_valid`
are the fast, headless way to validate a change. `cargo test` also runs but
doesn't touch the GPU — the reference-orbit tests are pure CPU math (see
below), and `shader_valid` parses/validates WGSL with `naga` statically instead
of creating a pipeline.
## Architecture
### The perturbation pipeline (the core mechanism, spans several files)
For a pixel at parameter `c = C_ref + dc`, its orbit is written as
`y_n = X_n + e_n`, where `X_n` is the (shared, high-precision) reference orbit
and `e_n` is a small `f32` delta. Whenever `|y_n| < |e_n|` (or the reference
runs out), rebase: `e ← y_n − X_0`, restart the reference index at 0. This is
what makes deep zoom cheap — one expensive high-precision orbit, then every
pixel is a handful of `f32` complex multiplies.
- `src/view.rs` — `ViewState`; center is arbitrary-precision `FBig` (`Big`
type alias), pixel scale stays `f64` (still in-range at 10³⁰×). Precision
(bits) scales with zoom depth (`precision_for`).
- `src/fractal/kind.rs` — the `FractalKind` enum (Mandelbrot, Burning Ship,
Tricorn, Multibrot, Celtic, Perpendicular, Buffalo, Phoenix, Lambda,
Complex Multibrot) plus everything that only needs to switch on it:
`label`/`description`/`formula` (UI text), `share_tag`/`from_share_tag`
(share-link encoding), `default_set_view` (per-kind starting view), and the
`ALL` array used to enumerate every kind.
- `src/fractal/reference.rs` — `compute_reference`/`compute_set_reference`:
iterate the chosen formula at high precision on the CPU, emitting `Z_n` as
`f32` pairs — that's the reference orbit the GPU perturbs from. At
precision ≤ `F64_MAX_PRECISION` (80 bits, i.e. shallow views) it takes a
plain-`f64` fast path (`compute_reference_f64`), so each kind's formula
exists twice in this file (f64 + `FBig`) and both must stay in sync;
`f64_fast_path_matches_big` checks they agree. Requests are made with 1.5×
iteration headroom (`reference_iterations` in `app.rs`), so auto-iterations
creeping up during a zoom doesn't recompute the orbit every frame.
- `src/shaders/*.wgsl` — none of these are standalone WGSL modules; WGSL has
no `#include`, so each is compiled by concatenating plain-text fragments
with `concat!`/`include_str!` at the `create_shader_module` call site (see
`renderer.rs`, `buddhabrot.rs`, and `tests/shader_valid.rs`, which must
concatenate the same pieces to validate what actually gets built).
`common.wgsl` (fullscreen-triangle vertex helper, `cmul`/`cpow`, `KIND_*`
constants) is prepended to every shader. `iterate_uniforms.wgsl` (the
perturbation-pipeline `Uniforms` struct + `palette()`) is additionally
prepended to `mandelbrot.wgsl` and `colorize.wgsl`, which share that layout.
Because there's no namespacing, a definition must live in exactly one file
among those concatenated together for a given shader — don't redefine a
`common.wgsl`/`iterate_uniforms.wgsl` symbol locally.
- `src/shaders/mandelbrot.wgsl` — the perturbation fragment shader. It is
**specialized per pipeline** through WGSL `override` constants (`KIND`,
`IS_JULIA`, `DE`), so the per-iteration kind/Julia/DE branches fold away at
pipeline creation. Read those constants in the shader, never `u.kind` /
`u.is_julia` / `u.de_coloring` (they're still uploaded for layout reasons).
`renderer.rs` builds one pipeline set per `PipelineKey` lazily on first
use, and `tests/shader_valid.rs` compiles every kind × Julia × DE variant to
SPIR-V. So a new kind needs no pipeline-list change, only its `KIND_*`
constant. `buddhabrot.wgsl` does the same with its own `override KIND`.
Interior pixels exit early through **periodicity detection**. It uses
Brent-style checkpoints plus two guards: the cycle's multiplier must be
clearly attracting (`PERIOD_MAX_MULT2`), and the contracting return must
repeat in `PERIOD_CONFIRMATIONS` consecutive windows. Both guards are
needed: without them, exterior pixels at cusps and minibrot edges turned
black. Retune them only against f64 ground truth on such views. Phoenix is
excluded (two-term map).
`advance_delta(z, e)` is the per-kind delta step (`z` = reference point,
`e` = current delta); the caller adds `step_add` (= `dc`) afterward — this
relies on `c` being additive in every current kind's formula (a kind where
it isn't, e.g. a rational map with `c` in a denominator, would need its own
step function that consumes `dc` internally instead, plus extra per-step
reference data since the orbit point alone wouldn't be enough to recover an
exact delta). `fprime(z)` is the derivative used for distance-estimation
(DE) shading; exact for holomorphic kinds, an approximation (`~2Z`) for the
abs-based ones. A `KIND_*` constant (from `common.wgsl`) must match the
matching `FractalKind` variant's discriminant exactly.
- `src/fractal/renderer.rs` — `FractalRenderer` (wgpu pipelines, uniform +
storage buffers, bind groups), `Uniforms` (repr(C) layout that must match
the WGSL `Uniforms` struct field-for-field, including padding; it includes
CPU-precomputed data: `cm_coef`, the Complex Multibrot binomial
coefficients from `app.rs::complex_binomials`, and `light_count` for the
packed `GpuLight` buffer from `lights.rs::gpu_lights`), and
`FractalCallback` (the `egui_wgpu::CallbackTrait` impl: `prepare()` uploads
changed buffers and decides whether to re-run the iterate pass, the cheap
colourise pass, or just blit the cached texture). Also `ExportRender`, a
self-contained tiled renderer used for PNG export off the UI thread.
- `src/worker.rs` — native background thread for reference-orbit computation
(coalesces bursts of requests so a fast drag doesn't compute every
intermediate view). The wasm32 build computes inline instead (see the
`#[cfg(target_arch = "wasm32")]` branch in `app.rs::ensure_reference`) —
**any signature change to `compute_reference`/`compute_set_reference` or
`RefRequest`/`RefResult` must be applied to both call sites.**
- `src/app.rs` — `FractalApp` (the egui app + all UI). Key methods:
`should_request`/`ensure_reference` (decide when the reference is stale and
dispatch/collect it), `make_uniforms` (assemble the per-frame `Uniforms`),
`tick_animations` (drives the "morph c/p/λ" and auto-zoom animations),
`default_view_for` (wraps `FractalKind::default_set_view`, adding the
kind-independent Julia case). `JULIA_PRESETS` and `SET_PRESETS` are sized as
`[T; FractalKind::<last variant> as usize + 1]` — adding a new `FractalKind`
means bumping both (and adding an empty `&[]` slot to each if the kind has
none), plus adding it to `FractalKind::ALL` in `kind.rs`.
- `src/fractal/share.rs` — `ShareState`: encodes the full view (mode, kind,
full-precision decimal center, zoom, iterations, per-kind constants,
coloring) as a `#`-fragment URL for bookmarking/sharing deep-zoom locations.
### Adding a new `FractalKind`
Touches, in order: `kind.rs` (enum variant + `ALL` slot + `label`/
`description`/`formula`/`share_tag`/`from_share_tag`/`default_set_view`
arms), `reference.rs` (CPU iteration formula arm, and a test comparing
against a naive `f64` iteration), `common.wgsl` (matching `KIND_*` const),
`mandelbrot.wgsl` (matching `advance_delta`/`fprime` arms), `buddhabrot.wgsl`
(matching arm in `advance()`, if the kind makes sense as a Buddhabrot),
`renderer.rs` `Uniforms` (only if the kind needs a new per-kind constant,
e.g. Phoenix's `phoenix_p`), `app.rs` (`JULIA_PRESETS`/`SET_PRESETS` slot,
and optionally a UI control for its constant + an animation toggle,
following the Phoenix/Lambda pattern). If `c` doesn't enter the formula
additively (e.g. a rational map with `c` in a denominator), the
`advance_delta`/`step_add` split doesn't work — that needs its own step
function plus extra per-step reference data uploaded in a second GPU buffer
alongside the orbit.
### Buddhabrot is a separate pipeline
`src/fractal/buddhabrot.rs` + `src/shaders/buddhabrot.wgsl` implement the
Monte-Carlo orbit-density histogram. It does **not** use the perturbation/
reference-orbit machinery: a Buddhabrot sample's orbit scatters across the
whole image rather than staying in one pixel, so it's plain `f32` iteration
from the live view (no deep zoom) via a compute pass that accumulates into a
histogram buffer, tone-mapped by a fragment pass every frame. Its own
`KIND_*` iteration formulas in `advance()` must be kept in sync with
`reference.rs` by hand (there's no shared code path).
### Two-pass render + caching (`renderer.rs`)
The interactive path splits iteration (expensive, perturbation) from
colourising (cheap, palette remap) into separate offscreen textures, so
palette/color-scale/offset tweaks skip re-iteration entirely (`geom_differs`
vs `color_differs` in `renderer.rs` decide which pass reruns). A frame where
neither differs uploads and renders nothing and only blits. So any new
uniform field must go into one of those two functions (or the lights
comparison), or changing it won't redraw.
AA is **adaptive** on the interactive path. `fs_data` always iterates 1
sample per pixel. When AA is on, `fs_refine` reads that texture and runs the
2×2 grid only on pixels whose 4-neighbours differ (interior/exterior edge, or
`ci`/DE beyond `AA_CI_EPS`/`AA_DE_EPS`), copying the rest. Colourise then
reads the refined texture. PNG export (`fs_color`) still supersamples every
pixel.
The 3D view (`colorize.wgsl::ray_marching`) sphere-traces the DE height
field straight from the data texture. It's cheap: rays start on the z = 0
plane, and most hit within a few steps (about 4 on average). A min-height
mip pyramid (quadtree height-field tracing) was tried and measured about 3×
slower, because it needs about 12 costlier steps per ray. Don't reintroduce it.
Orbiting the camera only re-runs the colourise pass, never iteration.
While the user is actively panning/zooming, the app renders downscaled with
AA off (`INTERACT_DOWNSCALE`) and snaps back to full resolution once input
settles (`INTERACT_SETTLE`).
+2 -7
View File
@@ -8,17 +8,14 @@ bytemuck = { version = "1.25.2", features = ["derive"] }
dashu-float = "0.6.0"
eframe = { version = "0.36.2", default-features = false, features = ["wgpu", "default_fonts", "x11", "wayland", "accesskit"] }
egui = "0.36.2"
glam = "0.33.8"
futures-channel = { version = "0.3.34", default-features = false, features = ["alloc", "std"] }
log = "0.4.34"
png = "0.18.1"
[target.'cfg(not(target_arch = "wasm32"))'.dependencies]
env_logger = "0.11.11"
clap = { version = "4.5.51", features = ["derive"] }
pollster = "1.0.1"
[target.'cfg(target_arch = "wasm32")'.dependencies]
futures-channel = { version = "0.3.34", default-features = false, features = ["alloc", "std"] }
console_error_panic_hook = "0.1.7"
console_log = "1.1.0"
js-sys = "0.3.105"
@@ -29,8 +26,6 @@ web-sys = { version = "0.3.105", features = ["Window", "Location", "Url", "UrlSe
# Release: optimize hard (fractal math is hot).
[profile.release]
opt-level = 3
# codegen-units = 1
debug = true
# Dev: keep our own crate debuggable, but optimize dependencies (dashu, wgpu,
# egui) so the explorer is actually interactive during development.
@@ -41,4 +36,4 @@ opt-level = 1
opt-level = 3
[dev-dependencies]
naga = { version = "30", features = ["wgsl-in", "spv-out"] }
naga = { version = "30", features = ["wgsl-in"] }
-1
View File
@@ -20,7 +20,6 @@ wasm-bindgen \
target/wasm32-unknown-unknown/release/mandelbrot.wasm
cp index.html "$OUT/index.html"
cp favicon.ico "$OUT/favicon.ico"
echo "==> done: $OUT/ (index.html, mandelbrot.js, mandelbrot_bg.wasm)"
echo " serve: python3 -m http.server -d $OUT 8080"
BIN
View File
Binary file not shown.

Before

Width:  |  Height:  |  Size: 422 KiB

+256 -1555
View File
File diff suppressed because it is too large Load Diff
-74
View File
@@ -1,74 +0,0 @@
use std::f32::consts::{PI, TAU};
use glam::Vec3;
#[derive(Default, Clone)]
pub struct Camera {
pub position: glam::Vec3,
pub yaw: f32,
pub pitch: f32,
pub aspect_ratio: f32,
z_near: f32,
z_far: f32,
}
impl Camera {
pub fn new() -> Self {
Self {
position: Vec3::new(0., 0., -1.),
yaw: 0. * PI / 180.,
pitch: -30. * PI / 180.,
aspect_ratio: 1.,
z_near: 0.1,
z_far: 100.,
}
}
pub fn set_aspect_ratio(&mut self, aspect_ratio: f32) {
self.aspect_ratio = aspect_ratio;
}
/// Adjust yaw/pitch by the given deltas (radians). Pitch is clamped just
/// short of straight up/down to avoid the view flipping past the pole.
/// Yaw is wrapped to [-π, π) so the 2D <-> 3D transition (which scales
/// yaw by `t`) always unwinds the short way instead of every past turn.
pub fn rotate(&mut self, dyaw: f32, dpitch: f32) {
const PITCH_LIMIT: f32 = PI / 2.0 - 0.01;
self.yaw = (self.yaw + dyaw + PI).rem_euclid(TAU) - PI;
self.pitch = (self.pitch + dpitch).clamp(-PITCH_LIMIT, PITCH_LIMIT);
}
pub fn orthographic(&self, t: f32) -> glam::Mat4 {
let yaw = self.yaw * t;
let pitch = self.pitch * t;
let zoom = 0.5 * (1. + t);
// Orbit pivot: the center of the fractal texture, which the raymarcher's
// `sdf` lays out over world x ∈ [0, aspect], y ∈ [0, 1] on the z = 0 plane.
let view = glam::Mat4::from_translation(Vec3::new(0.5 * self.aspect_ratio, 0.5, 0.))
* glam::Mat4::from_rotation_z(-yaw)
* glam::Mat4::from_rotation_x(-pitch)
* glam::Mat4::from_translation(self.position);
glam::camera::lh::proj::directx::orthographic(
-self.aspect_ratio / 4. / zoom,
self.aspect_ratio / 4. / zoom,
-0.25 / zoom,
0.25 / zoom,
self.z_near,
self.z_far,
) * view.inverse()
}
pub fn direction(&self, t: f32) -> glam::Vec3 {
let yaw = self.yaw * t;
let pitch = self.pitch * t;
let forward =
glam::Mat3::from_rotation_z(-yaw) * glam::Mat3::from_rotation_x(-pitch) * glam::Vec3::Z;
forward.normalize()
}
}
-174
View File
@@ -1,174 +0,0 @@
// Native command-line arguments. Currently mirrors the old `MANDEL_*` debug
// env vars one-for-one; this is the foundation a future headless (no-window,
// render-to-file) mode will build on.
use clap::{Parser, ValueEnum};
use crate::fractal::FractalKind;
#[derive(Parser, Debug, Default)]
#[command(name = "mandelbrot", about = "Deep-zoom fractal explorer", version)]
pub struct Cli {
/// Fractal formula to render.
#[arg(long, value_enum)]
pub kind: Option<KindArg>,
/// Start in Julia mode with this seed constant.
#[arg(long, value_name = "RE,IM")]
pub julia: Option<String>,
/// Switch to the Buddhabrot renderer.
#[arg(long)]
pub buddhabrot: bool,
/// Rendering mode to use.
#[arg(long)]
pub rendering_kind: Option<RenderingKindArg>,
/// Exponent for the Multibrot kind (z -> z^power + c), clamped to [2, 8].
#[arg(long)]
pub power: Option<u32>,
/// Complex exponent for the Complex Multibrot kind (z -> z^power + c).
#[arg(long, value_name = "RE,IM")]
pub complex_power: Option<String>,
/// Phoenix constant p for the Phoenix kind (z -> z^2 + c + p*z_prev).
#[arg(long, value_name = "RE,IM")]
pub phoenix_p: Option<String>,
/// Lambda constant λ for the Lambda kind (z -> λ*z*(1 - z)).
#[arg(long, value_name = "RE,IM")]
pub lambda_l: Option<String>,
/// Restore a view from a share-link fragment (the part after '#').
#[arg(long, value_name = "FRAGMENT")]
pub share: Option<String>,
/// Jump to a view on startup.
#[arg(long, value_name = "RE,IM,HALF_HEIGHT[,ITERATIONS]")]
pub view: Option<String>,
/// Jump to a specific position on startup.
#[arg(long, short('p'), value_name = "RE,IM")]
pub position: Option<String>,
/// Set a maximum iterations count on startup.
#[arg(long, short('i'))]
pub iterations: Option<u32>,
/// Set the zoom level on startup.
#[arg(long("zoom"), short('z'))]
pub half_height: Option<String>,
/// Enable distance-estimation shading.
#[arg(long)]
pub de: bool,
/// Coloring palette index.
#[arg(long, value_name = "INDEX")]
pub palette: Option<u32>,
/// Output path for --headless (default: fractal-<timestamp>.png). When
/// animating (--to-view/--to-share), this is a directory of
/// frame-00001.png, frame-00002.png, ... instead (default:
/// frames-<timestamp>/).
#[arg(long, value_name = "PATH")]
pub export_path: Option<String>,
/// End view for an animation: "re,im,half_height[,iterations]", the same
/// syntax as --view. Combine with --view (or --share, --kind, --julia...)
/// for the start view; headless then renders a sequence of frames
/// interpolating the camera from start to end instead of a single PNG.
#[arg(long, value_name = "RE,IM,HALF_HEIGHT[,ITERATIONS]")]
pub to_view: Option<String>,
/// End view for an animation, as a share-link fragment (only the
/// position/zoom/iterations are used; alternative to --to-view for
/// pasting a location copied from the app's "Copy share link").
#[arg(long, value_name = "FRAGMENT")]
pub to_share: Option<String>,
/// Set a maximum iterations count at animation end.
#[arg(long)]
pub to_iterations: Option<u32>,
/// Number of frames to render for an animation. Alternative to --fps +
/// --duration.
#[arg(long, value_name = "N")]
pub frames: Option<u32>,
/// Frames per second, used with --duration to compute the frame count
/// (ignored if --frames is given). Also used in the ffmpeg command
/// hint printed after rendering.
#[arg(long, value_name = "N", default_value_t = 30.0)]
pub fps: f64,
/// Animation duration in seconds, used with --fps to compute the frame
/// count (ignored if --frames is given).
#[arg(long, value_name = "SECONDS")]
pub duration: Option<f64>,
/// Pace animation frames linearly instead of easing in/out (smoothstep).
#[arg(long)]
pub linear: bool,
/// Run without opening a window: render the current view to a PNG and
/// exit. Combine with --kind/--julia/--share/--view etc. to pick what to
/// render, or --to-view/--to-share to render an animation instead of a
/// single frame. Not yet supported with --buddhabrot.
#[arg(long)]
pub headless: bool,
/// Output image width in pixels (--headless only).
#[arg(long, value_name = "PX", default_value_t = 1920)]
pub width: u32,
/// Output image height in pixels (--headless only).
#[arg(long, value_name = "PX", default_value_t = 1080)]
pub height: u32,
}
#[derive(Copy, Clone, Debug, ValueEnum)]
pub enum KindArg {
Mandelbrot,
#[value(alias = "ship")]
BurningShip,
#[value(alias = "mandelbar")]
Tricorn,
#[value(alias = "multi")]
Multibrot,
Celtic,
#[value(alias = "perp")]
Perpendicular,
Buffalo,
Phoenix,
Lambda,
#[value(alias = "cmulti")]
ComplexMultibrot,
}
#[derive(Copy, Clone, Debug, ValueEnum)]
pub enum RenderingKindArg {
Classic,
Shadow,
#[value(alias = "3d")]
Dimension3,
}
impl From<KindArg> for FractalKind {
fn from(k: KindArg) -> Self {
match k {
KindArg::Mandelbrot => FractalKind::Mandelbrot,
KindArg::BurningShip => FractalKind::BurningShip,
KindArg::Tricorn => FractalKind::Tricorn,
KindArg::Multibrot => FractalKind::Multibrot,
KindArg::Celtic => FractalKind::Celtic,
KindArg::Perpendicular => FractalKind::Perpendicular,
KindArg::Buffalo => FractalKind::Buffalo,
KindArg::Phoenix => FractalKind::Phoenix,
KindArg::Lambda => FractalKind::Lambda,
KindArg::ComplexMultibrot => FractalKind::ComplexMultibrot,
}
}
}
-412
View File
@@ -1,412 +0,0 @@
//! Buddhabrot / Nebulabrot rendering: a Monte-Carlo orbit-density histogram,
//! accumulated progressively across frames by a compute pass and tone-mapped
//! to colour by a fragment pass. See `shaders/buddhabrot.wgsl` for the "why"
//! this is a separate pipeline from the escape-time perturbation renderer.
use std::collections::HashMap;
use eframe::egui_wgpu::{self, wgpu};
/// Random samples dispatched per accumulating frame. Chosen so a frame stays
/// interactive on a modest GPU even when most samples run the full `b_cap`
/// (e.g. the view sits entirely inside the set, so nothing escapes).
const SAMPLES_PER_DISPATCH: u32 = 150_000;
const WORKGROUP_SIZE: u32 = 64;
/// GPU-side parameters for both the accumulate (compute) and tonemap
/// (fragment) passes. Layout must match `Uniforms` in `buddhabrot.wgsl`.
#[repr(C)]
#[derive(Copy, Clone, PartialEq, bytemuck::Pod, bytemuck::Zeroable)]
pub struct BuddhabrotUniforms {
pub center: [f32; 2],
pub half_height: f32,
pub aspect: f32,
pub phoenix_p: [f32; 2],
pub lambda_l: [f32; 2],
pub bailout_sq: f32,
/// Iteration formula (`FractalKind::shader_id`); `KIND_LAMBDA` samples z0
/// instead of c (see the shader's doc comment).
pub kind: u32,
/// Exponent for the Multibrot kind.
pub power: u32,
/// Nested escape-iteration caps (r_cap <= g_cap <= b_cap) that bucket an
/// orbit's points into the R/G/B histogram planes.
pub r_cap: u32,
pub g_cap: u32,
pub b_cap: u32,
/// RNG nonce, bumped every dispatch so each frame samples fresh points.
pub seed: u32,
pub samples_this_dispatch: u32,
/// Tonemap brightness multiplier (user-controlled).
pub exposure: f32,
pub width: u32,
pub height: u32,
/// Running total of samples accumulated into the current histogram
/// (across all dispatches since the last reset); normalizes brightness.
pub total_samples: f32,
/// Tonemap colour style: 0 = classic (R/G/B = raw caps), 1 = nebula
/// (yellow core, blue halo), 2 = grayscale. Display-only, like `exposure`
/// — excluded from `ContentKey` so changing it doesn't reset accumulation.
pub palette: u32,
/// Padding so `complex_power` (a vec2, 8-byte aligned in the shader)
/// starts on an 8-byte boundary.
pub _pad0: u32,
/// Complex exponent for the Complex Multibrot kind; ignored by other kinds.
pub complex_power: [f32; 2],
}
/// The subset of `BuddhabrotUniforms` that determines the *content* of the
/// histogram (as opposed to `exposure`, a display-only rescale). A change in
/// any of these invalidates the accumulated histogram.
#[derive(Copy, Clone, PartialEq)]
struct ContentKey {
center: [f32; 2],
half_height: f32,
aspect: f32,
phoenix_p: [f32; 2],
lambda_l: [f32; 2],
bailout_sq: f32,
kind: u32,
power: u32,
complex_power: [f32; 2],
r_cap: u32,
g_cap: u32,
b_cap: u32,
}
impl From<&BuddhabrotUniforms> for ContentKey {
fn from(u: &BuddhabrotUniforms) -> Self {
Self {
center: u.center,
half_height: u.half_height,
aspect: u.aspect,
phoenix_p: u.phoenix_p,
lambda_l: u.lambda_l,
bailout_sq: u.bailout_sq,
kind: u.kind,
power: u.power,
complex_power: u.complex_power,
r_cap: u.r_cap,
g_cap: u.g_cap,
b_cap: u.b_cap,
}
}
}
/// The histogram buffer and its two bind groups, sized to the widget.
struct Histogram {
buffer: wgpu::Buffer,
compute_bind_group: wgpu::BindGroup,
tonemap_bind_group: wgpu::BindGroup,
width: u32,
height: u32,
}
pub struct BuddhabrotRenderer {
shader: wgpu::ShaderModule,
compute_pipeline_layout: wgpu::PipelineLayout,
/// Accumulation pipelines, specialized per fractal kind (the shader's
/// `override KIND`, so `advance()` has no per-step kind branches) and
/// built lazily on first use.
compute_pipelines: HashMap<u32, wgpu::ComputePipeline>,
compute_bind_group_layout: wgpu::BindGroupLayout,
tonemap_pipeline: wgpu::RenderPipeline,
tonemap_bind_group_layout: wgpu::BindGroupLayout,
uniform_buffer: wgpu::Buffer,
histogram: Option<Histogram>,
/// What the current histogram's content was last accumulated for; a
/// mismatch clears the histogram and restarts accumulation.
last_content: Option<ContentKey>,
/// Running sample count since the last reset (mirrors what was written
/// into `total_samples`, since the callback doesn't own that state).
total_samples: f32,
seed: u32,
}
impl BuddhabrotRenderer {
pub fn new(device: &wgpu::Device, target_format: wgpu::TextureFormat) -> Self {
let shader = device.create_shader_module(wgpu::ShaderModuleDescriptor {
label: Some("buddhabrot"),
source: wgpu::ShaderSource::Wgsl(
concat!(
include_str!("../shaders/common.wgsl"),
include_str!("../shaders/buddhabrot.wgsl"),
)
.into(),
),
});
let uniform_buffer = device.create_buffer(&wgpu::BufferDescriptor {
label: Some("buddhabrot uniforms"),
size: std::mem::size_of::<BuddhabrotUniforms>() as u64,
usage: wgpu::BufferUsages::UNIFORM | wgpu::BufferUsages::COPY_DST,
mapped_at_creation: false,
});
let compute_bind_group_layout =
device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor {
label: Some("buddhabrot compute bind group layout"),
entries: &[
wgpu::BindGroupLayoutEntry {
binding: 0,
visibility: wgpu::ShaderStages::COMPUTE,
ty: wgpu::BindingType::Buffer {
ty: wgpu::BufferBindingType::Uniform,
has_dynamic_offset: false,
min_binding_size: None,
},
count: None,
},
wgpu::BindGroupLayoutEntry {
binding: 1,
visibility: wgpu::ShaderStages::COMPUTE,
ty: wgpu::BindingType::Buffer {
ty: wgpu::BufferBindingType::Storage { read_only: false },
has_dynamic_offset: false,
min_binding_size: None,
},
count: None,
},
],
});
let compute_pipeline_layout =
device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor {
label: Some("buddhabrot compute pipeline layout"),
bind_group_layouts: &[Some(&compute_bind_group_layout)],
immediate_size: 0,
});
let tonemap_bind_group_layout =
device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor {
label: Some("buddhabrot tonemap bind group layout"),
entries: &[
wgpu::BindGroupLayoutEntry {
binding: 0,
visibility: wgpu::ShaderStages::FRAGMENT,
ty: wgpu::BindingType::Buffer {
ty: wgpu::BufferBindingType::Uniform,
has_dynamic_offset: false,
min_binding_size: None,
},
count: None,
},
wgpu::BindGroupLayoutEntry {
binding: 2,
visibility: wgpu::ShaderStages::FRAGMENT,
ty: wgpu::BindingType::Buffer {
ty: wgpu::BufferBindingType::Storage { read_only: true },
has_dynamic_offset: false,
min_binding_size: None,
},
count: None,
},
],
});
let tonemap_pipeline_layout =
device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor {
label: Some("buddhabrot tonemap pipeline layout"),
bind_group_layouts: &[Some(&tonemap_bind_group_layout)],
immediate_size: 0,
});
let tonemap_pipeline = device.create_render_pipeline(&wgpu::RenderPipelineDescriptor {
label: Some("buddhabrot tonemap pipeline"),
layout: Some(&tonemap_pipeline_layout),
vertex: wgpu::VertexState {
module: &shader,
entry_point: Some("vs_main"),
buffers: &[],
compilation_options: Default::default(),
},
fragment: Some(wgpu::FragmentState {
module: &shader,
entry_point: Some("fs_tonemap"),
targets: &[Some(wgpu::ColorTargetState {
format: target_format,
blend: None,
write_mask: wgpu::ColorWrites::ALL,
})],
compilation_options: Default::default(),
}),
primitive: wgpu::PrimitiveState::default(),
depth_stencil: None,
multisample: wgpu::MultisampleState::default(),
multiview_mask: None,
cache: None,
});
Self {
shader,
compute_pipeline_layout,
compute_pipelines: HashMap::new(),
compute_bind_group_layout,
tonemap_pipeline,
tonemap_bind_group_layout,
uniform_buffer,
histogram: None,
last_content: None,
total_samples: 0.0,
seed: 0,
}
}
/// The accumulation pipeline for `kind`, built on first use.
fn compute_pipeline(&mut self, device: &wgpu::Device, kind: u32) -> &wgpu::ComputePipeline {
self.compute_pipelines.entry(kind).or_insert_with(|| {
device.create_compute_pipeline(&wgpu::ComputePipelineDescriptor {
label: Some("buddhabrot compute pipeline"),
layout: Some(&self.compute_pipeline_layout),
module: &self.shader,
entry_point: Some("cs_main"),
compilation_options: wgpu::PipelineCompilationOptions {
constants: &[("KIND", kind as f64)],
..Default::default()
},
cache: None,
})
})
}
/// Ensure the histogram buffer exists at `width`×`height`, recreating (and
/// resetting accumulation) on a size change.
fn ensure_histogram(&mut self, device: &wgpu::Device, width: u32, height: u32) {
if let Some(h) = &self.histogram
&& h.width == width
&& h.height == height
{
return;
}
let plane = (width as u64) * (height as u64);
let buffer = device.create_buffer(&wgpu::BufferDescriptor {
label: Some("buddhabrot histogram"),
size: plane * 3 * std::mem::size_of::<u32>() as u64,
usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_DST,
mapped_at_creation: false,
});
let compute_bind_group = device.create_bind_group(&wgpu::BindGroupDescriptor {
label: Some("buddhabrot compute bind group"),
layout: &self.compute_bind_group_layout,
entries: &[
wgpu::BindGroupEntry {
binding: 0,
resource: self.uniform_buffer.as_entire_binding(),
},
wgpu::BindGroupEntry {
binding: 1,
resource: buffer.as_entire_binding(),
},
],
});
let tonemap_bind_group = device.create_bind_group(&wgpu::BindGroupDescriptor {
label: Some("buddhabrot tonemap bind group"),
layout: &self.tonemap_bind_group_layout,
entries: &[
wgpu::BindGroupEntry {
binding: 0,
resource: self.uniform_buffer.as_entire_binding(),
},
wgpu::BindGroupEntry {
binding: 2,
resource: buffer.as_entire_binding(),
},
],
});
self.histogram = Some(Histogram {
buffer,
compute_bind_group,
tonemap_bind_group,
width,
height,
});
// New (zero-initialized) buffer: accumulation starts fresh.
self.last_content = None;
self.total_samples = 0.0;
}
}
/// Per-frame paint callback. `accumulate` controls whether a new batch of
/// samples is dispatched this frame (a content change always forces one
/// dispatch regardless, so a parameter/view change is never left blank).
pub struct BuddhabrotCallback {
pub uniforms: BuddhabrotUniforms,
pub accumulate: bool,
/// Widget size in physical pixels — the histogram resolution.
pub size_px: [u32; 2],
}
impl egui_wgpu::CallbackTrait for BuddhabrotCallback {
fn prepare(
&self,
device: &wgpu::Device,
queue: &wgpu::Queue,
_screen_descriptor: &egui_wgpu::ScreenDescriptor,
egui_encoder: &mut wgpu::CommandEncoder,
resources: &mut egui_wgpu::CallbackResources,
) -> Vec<wgpu::CommandBuffer> {
let Some(renderer) = resources.get_mut::<BuddhabrotRenderer>() else {
return Vec::new();
};
let width = self.size_px[0].max(1);
let height = self.size_px[1].max(1);
renderer.ensure_histogram(device, width, height);
let pipeline = renderer
.compute_pipeline(device, self.uniforms.kind)
.clone();
let content = ContentKey::from(&self.uniforms);
let content_changed = renderer.last_content != Some(content);
let should_dispatch = content_changed || self.accumulate;
if let Some(histogram) = &renderer.histogram {
if content_changed {
egui_encoder.clear_buffer(&histogram.buffer, 0, None);
renderer.total_samples = 0.0;
renderer.last_content = Some(content);
}
let mut uniforms = self.uniforms;
uniforms.width = width;
uniforms.height = height;
if should_dispatch {
renderer.seed = renderer.seed.wrapping_add(1);
renderer.total_samples += SAMPLES_PER_DISPATCH as f32;
uniforms.seed = renderer.seed;
uniforms.samples_this_dispatch = SAMPLES_PER_DISPATCH;
} else {
uniforms.samples_this_dispatch = 0;
}
uniforms.total_samples = renderer.total_samples;
queue.write_buffer(&renderer.uniform_buffer, 0, bytemuck::bytes_of(&uniforms));
if should_dispatch {
let mut pass = egui_encoder.begin_compute_pass(&wgpu::ComputePassDescriptor {
label: Some("buddhabrot accumulate pass"),
timestamp_writes: None,
});
pass.set_pipeline(&pipeline);
pass.set_bind_group(0, &histogram.compute_bind_group, &[]);
let workgroups = SAMPLES_PER_DISPATCH.div_ceil(WORKGROUP_SIZE);
pass.dispatch_workgroups(workgroups, 1, 1);
}
}
Vec::new()
}
fn paint(
&self,
_info: egui::PaintCallbackInfo,
render_pass: &mut wgpu::RenderPass<'static>,
resources: &egui_wgpu::CallbackResources,
) {
if let Some(renderer) = resources.get::<BuddhabrotRenderer>()
&& let Some(histogram) = &renderer.histogram
{
render_pass.set_pipeline(&renderer.tonemap_pipeline);
render_pass.set_bind_group(0, &histogram.tonemap_bind_group, &[]);
render_pass.draw(0..3, 0..1);
}
}
}
-165
View File
@@ -1,165 +0,0 @@
//! `FractalKind`: the enum selecting which iteration formula is in use, plus
//! everything that only needs to switch on it (UI label/description/formula
//! text, share-link tag, default parameter-plane view). The CPU/GPU orbit
//! math itself lives in `reference.rs` (CPU reference orbit) and
//! `shaders/mandelbrot.wgsl` (GPU perturbation delta) since both must also
//! stay in sync with `common.wgsl`'s `KIND_*` constants.
/// The iteration formula. Must be kept in sync with `advance_delta` and the
/// `KIND_*` constants in the shader.
#[repr(u8)]
#[derive(Clone, Copy, PartialEq, Eq, Debug)]
pub enum FractalKind {
/// `z -> z^2 + c`.
Mandelbrot = 0,
/// `z -> (|Re z| + i|Im z|)^2 + c`.
BurningShip = 1,
/// `z -> conj(z)^2 + c` (the Mandelbar).
Tricorn = 2,
/// `z -> z^power + c` (power >= 2).
Multibrot = 3,
/// `z -> |Re(z^2)| + i·Im(z^2) + c` (abs on the real output of the square).
Celtic = 4,
/// `z -> (x^2 - y^2) - 2·x·|y|·i + c` (abs on the imaginary input).
Perpendicular = 5,
/// `z -> |Re(z^2)| - |Im(z^2)|·i + c` (abs on both outputs).
Buffalo = 6,
/// `z -> z^2 + c + p·z_{n-1}` (two-term recurrence; `p` is `phoenix_p`).
Phoenix = 7,
/// `z -> lambda·z(1 - z)` (logistic map).
Lambda = 8,
/// `z -> z^power + c`, where `power` is a complex constant (the
/// `complex_power` argument), via the principal branch `z^p = exp(p·ln z)`.
ComplexMultibrot = 9,
}
impl FractalKind {
/// Every kind, in declaration/discriminant order. Sized arrays keyed by
/// `kind as usize` (`JULIA_PRESETS`, `SET_PRESETS`) must have one slot per
/// entry here.
pub const ALL: [FractalKind; 10] = [
FractalKind::Mandelbrot,
FractalKind::BurningShip,
FractalKind::Tricorn,
FractalKind::Multibrot,
FractalKind::Celtic,
FractalKind::Perpendicular,
FractalKind::Buffalo,
FractalKind::Phoenix,
FractalKind::Lambda,
FractalKind::ComplexMultibrot,
];
pub fn description(&self) -> &'static str {
match self {
FractalKind::Mandelbrot => {
"The Mandelbrot set is the most famous fractal set, obtained with the simplest escape-time formula. This set represents all Julia fractals: each points of the Mandelbrot set is related to a specific Julia fractal."
}
FractalKind::BurningShip => {
"A variation of the famous Mandelbrot set, using absolute values on the real and imaginary part of each iterations."
}
FractalKind::Tricorn => {
"The Tricorn set is obtained using the same formula as the Mandelbrot set, taking the complex conjugate of the previous iteration."
}
FractalKind::Multibrot => {
"Multibrot use the same formula as the Mandelbrot set, with a bigger exposant."
}
FractalKind::Celtic => "",
FractalKind::Perpendicular => "",
FractalKind::Buffalo => "",
FractalKind::Phoenix => "",
FractalKind::Lambda => "",
FractalKind::ComplexMultibrot => {
"Like Multibrot, but the exponent itself is a complex number instead of a plain integer, via z^p = exp(p·ln z)."
}
}
}
/// UI label for this kind (combo box / info panel heading).
pub fn label(&self) -> &'static str {
match self {
FractalKind::Mandelbrot => "Mandelbrot",
FractalKind::BurningShip => "Burning Ship",
FractalKind::Tricorn => "Tricorn",
FractalKind::Multibrot => "Multibrot",
FractalKind::Celtic => "Celtic",
FractalKind::Perpendicular => "Perpendicular",
FractalKind::Buffalo => "Buffalo",
FractalKind::Phoenix => "Phoenix",
FractalKind::Lambda => "Lambda",
FractalKind::ComplexMultibrot => "Complex Multibrot",
}
}
/// The iteration formula in human-readable notation (mirrors the doc
/// comments on the variants above). `power` is only used by Multibrot;
/// `complex_power` only by Complex Multibrot.
pub fn formula(&self, power: u32, complex_power: (f64, f64)) -> String {
match self {
FractalKind::Mandelbrot => "z = z² + c".to_string(),
FractalKind::BurningShip => "z = (|Re(z)| + i|Im(z)|)² + c".to_string(),
FractalKind::Tricorn => "z = conj(z)² + c".to_string(),
FractalKind::Multibrot => format!("z = z^{power} + c"),
FractalKind::Celtic => "z = |Re(z²)| + i·Im(z²) + c".to_string(),
FractalKind::Perpendicular => "z = (x² − y²) − 2x|y|i + c".to_string(),
FractalKind::Buffalo => "z = |Re(z²)| − i|Im(z²)| + c".to_string(),
FractalKind::Phoenix => "z = z² + c + p·z_prev".to_string(),
FractalKind::Lambda => "z = λ·z(1 − z)".to_string(),
FractalKind::ComplexMultibrot => {
format!("z = z^({:.3}{:+.3}i) + c", complex_power.0, complex_power.1)
}
}
}
/// Short tag used to identify this kind in a share-link fragment.
pub fn share_tag(&self) -> &'static str {
match self {
FractalKind::Mandelbrot => "mandel",
FractalKind::BurningShip => "burning",
FractalKind::Tricorn => "tricorn",
FractalKind::Multibrot => "multi",
FractalKind::Celtic => "celtic",
FractalKind::Perpendicular => "perp",
FractalKind::Buffalo => "buffalo",
FractalKind::Phoenix => "phoenix",
FractalKind::Lambda => "lambda",
FractalKind::ComplexMultibrot => "cmulti",
}
}
/// Inverse of `share_tag`; unknown tags fall back to `None` so the caller
/// can decide the default (matches historical share-link behavior).
pub fn from_share_tag(tag: &str) -> Option<FractalKind> {
Some(match tag {
"mandel" => FractalKind::Mandelbrot,
"burning" => FractalKind::BurningShip,
"tricorn" => FractalKind::Tricorn,
"multi" => FractalKind::Multibrot,
"celtic" => FractalKind::Celtic,
"perp" => FractalKind::Perpendicular,
"buffalo" => FractalKind::Buffalo,
"phoenix" => FractalKind::Phoenix,
"lambda" => FractalKind::Lambda,
"cmulti" => FractalKind::ComplexMultibrot,
_ => return None,
})
}
/// Default parameter-plane (Mandelbrot-mode) view for this kind, as
/// `(center_re, center_im, half_height)`. The Julia (dynamical) plane
/// doesn't vary by kind, so it isn't covered here.
pub fn default_set_view(&self) -> (f64, f64, f64) {
match self {
FractalKind::Mandelbrot => (-0.5, 0.0, 1.25),
FractalKind::BurningShip => (-0.5, -0.5, 1.3),
FractalKind::Tricorn => (-0.25, 0.0, 1.7),
FractalKind::Multibrot => (0.0, 0.0, 1.5),
FractalKind::Celtic => (-0.5, 0.0, 1.6),
FractalKind::Perpendicular => (-0.5, 0.0, 1.5),
FractalKind::Buffalo => (-0.5, 0.5, 1.5),
FractalKind::Phoenix => (-0.5, 0.0, 1.5),
FractalKind::Lambda => (-0.5, 0.0, 2.4),
FractalKind::ComplexMultibrot => (0.0, 0.0, 1.5),
}
}
}
+5 -10
View File
@@ -1,18 +1,13 @@
//! GPU fractal rendering: wgpu pipeline, uniforms, reference orbit, and the
//! egui paint callback.
pub mod buddhabrot;
pub mod kind;
pub mod reference;
pub mod renderer;
pub mod share;
pub use buddhabrot::{BuddhabrotCallback, BuddhabrotRenderer, BuddhabrotUniforms};
pub use kind::FractalKind;
pub use reference::{compute_reference, compute_set_reference};
#[cfg(target_arch = "wasm32")]
pub use renderer::encode_png_with_progress;
#[cfg(not(target_arch = "wasm32"))]
pub use renderer::export_to_png_blocking;
pub use renderer::{ExportRender, FractalCallback, FractalRenderer, MAX_REF_POINTS, Uniforms};
pub use reference::{Bla, FractalKind, build_bla_table, compute_reference, compute_set_reference};
pub use renderer::{
ExportRender, FractalCallback, FractalRenderer, MAX_REF_POINTS, Uniforms,
encode_png_with_progress,
};
pub use share::ShareState;
+295 -501
View File
@@ -10,34 +10,42 @@
//! * Mandelbrot-set: `z0 = 0`, `c = view center` (the c-plane point per pixel).
//! * Julia-set: `z0 = view center`, `c = fractal constant` (fixed per view).
use super::kind::FractalKind;
use crate::view::{Big, big_from_f64};
use crate::view::Big;
/// The iteration formula. Must be kept in sync with `advance_delta` and the
/// `KIND_*` constants in the shader.
#[derive(Clone, Copy, PartialEq, Eq, Debug)]
pub enum FractalKind {
/// `z -> z^2 + c`.
Mandelbrot,
/// `z -> (|Re z| + i|Im z|)^2 + c`.
BurningShip,
/// `z -> conj(z)^2 + c` (the Mandelbar).
Tricorn,
/// `z -> z^power + c` (power >= 2).
Multibrot,
}
impl FractalKind {
/// Integer id matching the shader's `KIND_*` constants.
pub fn shader_id(self) -> u32 {
match self {
FractalKind::Mandelbrot => 0,
FractalKind::BurningShip => 1,
FractalKind::Tricorn => 2,
FractalKind::Multibrot => 3,
}
}
}
/// Reference orbit escapes once |Z|^2 exceeds this. Kept larger than the pixel
/// bailout so pixels escaping alongside the reference can still reach their
/// bailout before the stored orbit runs out.
const REFERENCE_ESCAPE_SQ: f64 = 1.0e10;
/// Up to this working precision (bits) the orbit is iterated in plain `f64`
/// instead of `FBig` — orders of magnitude faster, which matters most on the
/// web (where the reference is computed inline on the UI thread).
///
/// `precision_for` asks for `zoom_bits + 48` guard bits, but the GPU only
/// consumes the orbit as f32 deltas, so two things actually matter:
/// * Each f64 step's rounding (~1e-16 relative) acts like a tiny local error
/// in the pixel orbits too (perturbation reproduces whatever orbit it's
/// given), far below the f32 delta noise — the orbit only has to be a
/// consistent orbit, not the exact one.
/// * The reference center gets rounded to f64 (<= ~2.2e-16 absolute for
/// |c| <= 2), which shifts the image. At 80 bits (zoom_bits <= 32, i.e.
/// half-height >= ~2.3e-10) a pixel is >= ~5e-13 wide, so that shift stays
/// below 0.1% of a pixel.
const F64_MAX_PRECISION: usize = 80;
/// Compute the reference orbit `Z_0..Z_{len-1}` where `Z_0 = z0` and
/// `Z_{n+1} = f(Z_n, c)` for the given `kind` (and `power`, for Multibrot), up
/// to `max_iter` steps at `precision` bits. Each entry is `[re, im]` in f32.
#[allow(clippy::too_many_arguments)]
pub fn compute_reference(
z0_re: &Big,
z0_im: &Big,
@@ -47,143 +55,12 @@ pub fn compute_reference(
precision: usize,
kind: FractalKind,
power: u32,
phoenix_p: (f64, f64),
lambda_l: (f64, f64),
complex_power: (f64, f64),
) -> Vec<[f32; 2]> {
if precision <= F64_MAX_PRECISION {
return compute_reference_f64(
(z0_re.to_f64().value(), z0_im.to_f64().value()),
(c_re.to_f64().value(), c_im.to_f64().value()),
max_iter,
kind,
power,
phoenix_p,
lambda_l,
complex_power,
);
}
compute_reference_big(
z0_re,
z0_im,
c_re,
c_im,
max_iter,
precision,
kind,
power,
phoenix_p,
lambda_l,
complex_power,
)
}
/// [`compute_reference`]'s fast path for shallow views (see
/// [`F64_MAX_PRECISION`]): the same per-kind formulas in plain `f64`.
#[allow(clippy::too_many_arguments)]
fn compute_reference_f64(
z0: (f64, f64),
c: (f64, f64),
max_iter: u32,
kind: FractalKind,
power: u32,
phoenix_p: (f64, f64),
lambda_l: (f64, f64),
complex_power: (f64, f64),
) -> Vec<[f32; 2]> {
let (cr, ci) = c;
let (mut zr, mut zi) = z0;
// Previous iterate, for the Phoenix two-term recurrence (Y_{-1} = 0).
let (mut zr_prev, mut zi_prev) = (0.0f64, 0.0f64);
let (pr, pi) = phoenix_p;
let (lr, li) = lambda_l;
let mut points: Vec<[f32; 2]> = Vec::with_capacity(max_iter as usize + 1);
for _ in 0..=max_iter {
points.push([zr as f32, zi as f32]);
if zr * zr + zi * zi > REFERENCE_ESCAPE_SQ {
break;
}
let (new_zr, new_zi) = match kind {
FractalKind::Mandelbrot => ((zr + zi) * (zr - zi) + cr, 2.0 * zr * zi + ci),
FractalKind::BurningShip => (zr * zr - zi * zi + cr, (2.0 * zr * zi).abs() + ci),
FractalKind::Tricorn => (zr * zr - zi * zi + cr, ci - 2.0 * zr * zi),
FractalKind::Multibrot => {
let (mut rr, mut ri) = (1.0f64, 0.0f64);
for _ in 0..power.max(2) {
(rr, ri) = (rr * zr - ri * zi, rr * zi + ri * zr);
}
(rr + cr, ri + ci)
}
FractalKind::Celtic => ((zr * zr - zi * zi).abs() + cr, 2.0 * zr * zi + ci),
FractalKind::Perpendicular => (zr * zr - zi * zi + cr, ci - 2.0 * zr * zi.abs()),
FractalKind::Buffalo => ((zr * zr - zi * zi).abs() + cr, ci - (2.0 * zr * zi).abs()),
FractalKind::Phoenix => (
zr * zr - zi * zi + cr + (pr * zr_prev - pi * zi_prev),
2.0 * zr * zi + ci + (pr * zi_prev + pi * zr_prev),
),
FractalKind::Lambda => {
// λ·z(1 - z).
let (re2, im2) = (1.0 - zr, -zi);
let (lzr, lzi) = (lr * zr - li * zi, lr * zi + li * zr);
(lzr * re2 - lzi * im2, re2 * lzi + lzr * im2)
}
FractalKind::ComplexMultibrot => {
let (pr, pi) = complex_pow_complex_f64(zr, zi, complex_power.0, complex_power.1);
(pr + cr, pi + ci)
}
};
(zr_prev, zi_prev) = (zr, zi);
(zr, zi) = (new_zr, new_zi);
}
points
}
/// `f64` twin of [`complex_pow_complex`] (principal branch, `0^p = 0`).
fn complex_pow_complex_f64(zr: f64, zi: f64, pr: f64, pi: f64) -> (f64, f64) {
if zr == 0.0 && zi == 0.0 {
return (0.0, 0.0);
}
let ln_r = 0.5 * (zr * zr + zi * zi).ln();
let theta = zi.atan2(zr);
let mag = (pr * ln_r - pi * theta).exp();
let (sin_a, cos_a) = (pr * theta + pi * ln_r).sin_cos();
(mag * cos_a, mag * sin_a)
}
/// [`compute_reference`] at arbitrary precision (`FBig`), for deep views.
#[allow(clippy::too_many_arguments)]
fn compute_reference_big(
z0_re: &Big,
z0_im: &Big,
c_re: &Big,
c_im: &Big,
max_iter: u32,
precision: usize,
kind: FractalKind,
power: u32,
phoenix_p: (f64, f64),
lambda_l: (f64, f64),
complex_power: (f64, f64),
) -> Vec<[f32; 2]> {
let cr = c_re.clone().with_precision(precision).value();
let ci = c_im.clone().with_precision(precision).value();
let mut zr = z0_re.clone().with_precision(precision).value();
let mut zi = z0_im.clone().with_precision(precision).value();
// Previous iterate, for the Phoenix two-term recurrence (Y_{-1} = 0).
let mut zr_prev = big_zero(precision);
let mut zi_prev = big_zero(precision);
// Phoenix distortion constant `p` (a small fixed complex number).
let pr = big_from_f64(phoenix_p.0, precision);
let pi = big_from_f64(phoenix_p.1, precision);
// Lambda distortion constant `l` (a small fixed complex number).
let lr = big_from_f64(lambda_l.0, precision);
let li = big_from_f64(lambda_l.1, precision);
// Complex Multibrot exponent (a fixed complex number).
let cpow_re = big_from_f64(complex_power.0, precision);
let cpow_im = big_from_f64(complex_power.1, precision);
let mut points: Vec<[f32; 2]> = Vec::with_capacity(max_iter as usize + 1);
@@ -199,9 +76,8 @@ fn compute_reference_big(
let (new_zr, new_zi) = match kind {
FractalKind::Mandelbrot => {
// Z^2 = (zr^2 - zi^2) + (2 zr zi) i, with zr^2 - zi^2 as
// (zr + zi)(zr - zi): one multiply instead of two squares.
let re = (&zr + &zi) * (&zr - &zi) + &cr;
// Z^2 = (zr^2 - zi^2) + (2 zr zi) i.
let re = &zr.sqr() - &zi.sqr() + &cr;
let im = ((&zr * &zi) << 1) + &ci; // << 1 is exact ×2 in base 2
(re, im)
}
@@ -221,53 +97,8 @@ fn compute_reference_big(
let (pr, pi) = complex_pow(&zr, &zi, power.max(2), precision);
(pr + &cr, pi + &ci)
}
FractalKind::Celtic => {
// |Re(z^2)| + i·Im(z^2): abs the real output of the square.
let re = big_abs(&zr.sqr() - &zi.sqr()) + &cr;
let im = ((&zr * &zi) << 1) + &ci;
(re, im)
}
FractalKind::Perpendicular => {
// (x^2 - y^2) - 2·x·|y| i: abs the imaginary input.
let re = &zr.sqr() - &zi.sqr() + &cr;
let im = if zi.to_f64().value() < 0.0 {
&ci + ((&zr * &zi) << 1)
} else {
&ci - ((&zr * &zi) << 1)
};
(re, im)
}
FractalKind::Buffalo => {
// |Re(z^2)| - |Im(z^2)| i: abs both outputs.
let re = big_abs(&zr.sqr() - &zi.sqr()) + &cr;
let im = &ci - big_abs((&zr * &zi) << 1);
(re, im)
}
FractalKind::Phoenix => {
// z^2 + c + p·z_{n-1}.
let re2 = &zr.sqr() - &zi.sqr();
let im2 = (&zr * &zi) << 1;
let pzr = &pr * &zr_prev - &pi * &zi_prev;
let pzi = &pr * &zi_prev + &pi * &zr_prev;
(re2 + &cr + pzr, im2 + &ci + pzi)
}
FractalKind::Lambda => {
// λ·z(1 - z): logistic map.
let re2 = 1 - &zr;
let im2 = -&zi;
let lzr = &lr * &zr - &li * &zi;
let lzi = &lr * &zi + &li * &zr;
(&lzr * &re2 - &lzi * &im2, re2 * lzi + lzr * im2)
}
FractalKind::ComplexMultibrot => {
let (pr, pi) = complex_pow_complex(&zr, &zi, &cpow_re, &cpow_im, precision);
(pr + &cr, pi + &ci)
}
};
// Shift the previous iterate (only the Phoenix arm reads it).
zr_prev = zr;
zi_prev = zi;
zr = new_zr.with_precision(precision).value();
zi = new_zi.with_precision(precision).value();
}
@@ -299,35 +130,8 @@ fn complex_pow(zr: &Big, zi: &Big, power: u32, precision: usize) -> (Big, Big) {
(rr, ri)
}
/// `true` if `x` is (numerically) zero. The f64 check is exact for a true
/// zero; only matters here to special-case `ln(0)`.
fn is_big_zero(x: &Big) -> bool {
x.to_f64().value() == 0.0
}
/// `(zr + i zi)^(pr + i pi)` for a complex exponent, via the principal branch
/// `z^p = exp(p·ln z)` where `ln z = ln|z| + i·arg(z)`. Used by
/// `ComplexMultibrot`; must be kept in sync with the shader's `cpow`.
/// `z = 0` is special-cased to `0` (the formula's `ln(0)` would otherwise
/// panic; this is the correct limit for the `Re(p) > 0` region the UI
/// exposes).
fn complex_pow_complex(zr: &Big, zi: &Big, pr: &Big, pi: &Big, precision: usize) -> (Big, Big) {
if is_big_zero(zr) && is_big_zero(zi) {
return (big_zero(precision), big_zero(precision));
}
let r2 = &zr.sqr() + &zi.sqr();
let ln_r = r2.ln() >> 1; // 0.5 * ln(r2) = ln(sqrt(r2)); exact halving.
let theta = zi.atan2(zr);
let exp_re = (pr * &ln_r - pi * &theta).with_precision(precision).value();
let exp_im = (pr * &theta + pi * &ln_r).with_precision(precision).value();
let mag = exp_re.exp();
let (sin_a, cos_a) = exp_im.sin_cos();
(&mag * &cos_a, &mag * &sin_a)
}
/// Convenience: parameter-plane ("Mandelbrot-set") reference (`z0 = 0`,
/// `c = center`) for any `kind`.
#[allow(clippy::too_many_arguments)]
pub fn compute_set_reference(
center_re: &Big,
center_im: &Big,
@@ -335,26 +139,157 @@ pub fn compute_set_reference(
precision: usize,
kind: FractalKind,
power: u32,
phoenix_p: (f64, f64),
lambda_l: (f64, f64),
complex_power: (f64, f64),
) -> Vec<[f32; 2]> {
let zero = big_zero(precision);
compute_reference(
&zero,
&zero,
center_re,
center_im,
max_iter,
precision,
kind,
power,
phoenix_p,
lambda_l,
complex_power,
&zero, &zero, center_re, center_im, max_iter, precision, kind, power,
)
}
// ---------------------------------------------------------------------------
// Bivariate Linear Approximation (BLA)
//
// Deep in a zoom the per-pixel delta stays far smaller than the reference, so
// the nonlinear `e^2` term of the perturbation step is negligible and the step
// is effectively linear: `e -> A e + B dc`. BLA precomputes, for runs of
// iterations, the composed linear coefficients `(A, B)` plus a validity radius
// `r` (the largest `|e|` for which dropping `e^2` stays within tolerance). A
// pixel can then skip a whole run in one multiply whenever `|e| < r`.
//
// Runs are merged pairwise into a binary tree of levels: level `k` holds BLAs of
// length `2^k` starting at multiples of `2^k`. The GPU walks levels high→low to
// take the longest valid skip at the current index. Both sides recompute the
// per-level counts from `ref_len` (count[0] = ref_len-1, count[k] = count[k-1]/2)
// so no offset table needs to travel to the GPU — only this flat array does.
//
// Only the holomorphic square map (Mandelbrot/Julia, `A = 2 Z`, `B = 1`) is
// supported; other kinds fall back to per-iteration stepping on the GPU.
// ---------------------------------------------------------------------------
/// One merged linear step: `e_{n+l} = A e_n + B dc`, valid while `|e_n| < r`.
/// Laid out to match the WGSL `Bla` struct (two `vec2<f32>`, then `f32`, `u32`;
/// 24-byte std430 stride).
#[repr(C)]
#[derive(Clone, Copy, bytemuck::Pod, bytemuck::Zeroable)]
pub struct Bla {
pub a: [f32; 2],
pub b: [f32; 2],
pub r: f32,
pub l: u32,
}
/// Relative budget for the dropped nonlinear term, chosen near the f32 orbit
/// storage noise floor so BLA adds no visible error over plain perturbation.
const BLA_EPS: f64 = 1.0e-6;
/// f64 working form of a BLA (merges accumulate in f64, stored as f32).
#[derive(Clone, Copy)]
struct BlaF {
ar: f64,
ai: f64,
br: f64,
bi: f64,
r: f64,
l: u32,
}
impl BlaF {
fn to_bla(self) -> Bla {
Bla {
a: [self.ar as f32, self.ai as f32],
b: [self.br as f32, self.bi as f32],
r: self.r as f32,
l: self.l,
}
}
}
/// Merge two consecutive BLAs (`x` then `y`) into one covering both runs.
fn merge_bla(x: BlaF, y: BlaF, dc_max: f64) -> BlaF {
// A = Ay Ax ; B = Ay Bx + By (complex).
let ar = y.ar * x.ar - y.ai * x.ai;
let ai = y.ar * x.ai + y.ai * x.ar;
let br = y.ar * x.br - y.ai * x.bi + y.br;
let bi = y.ar * x.bi + y.ai * x.br + y.bi;
// Valid if |e| < rx (so x holds) and |Ax e + Bx dc| < ry (so y holds):
// |e| < (ry - |Bx| dc_max) / |Ax|.
let ax_mag = (x.ar * x.ar + x.ai * x.ai).sqrt();
let bx_mag = (x.br * x.br + x.bi * x.bi).sqrt();
let r_y = if ax_mag > 0.0 {
((y.r - bx_mag * dc_max) / ax_mag).max(0.0)
} else {
x.r
};
BlaF {
ar,
ai,
br,
bi,
r: x.r.min(r_y),
l: x.l + y.l,
}
}
/// Build the BLA table for a square-map reference orbit. `dc_max` is an upper
/// bound on any pixel's `|dc|` in the current view (used to size merge radii).
/// Returns a flat array with levels concatenated (level 0 first). Empty if the
/// orbit is too short to skip.
pub fn build_bla_table(points: &[[f32; 2]], dc_max: f64) -> Vec<Bla> {
let m = points.len();
if m < 2 {
return Vec::new();
}
// Level 0: one single step from each index i (uses Z_i). A = 2 Z, B = 1.
// Dropping e^2 is within tolerance while |e| < BLA_EPS |Z| (since the kept
// linear term is |2 Z e|).
let mut level0: Vec<BlaF> = Vec::with_capacity(m - 1);
for z in &points[..m - 1] {
let (zr, zi) = (z[0] as f64, z[1] as f64);
let zmag = (zr * zr + zi * zi).sqrt();
level0.push(BlaF {
ar: 2.0 * zr,
ai: 2.0 * zi,
br: 1.0,
bi: 0.0,
r: BLA_EPS * zmag,
l: 1,
});
}
let mut levels: Vec<Vec<BlaF>> = vec![level0];
while levels.last().unwrap().len() >= 2 {
let prev = levels.last().unwrap();
let mut next = Vec::with_capacity(prev.len() / 2);
let mut i = 0;
while i + 1 < prev.len() {
next.push(merge_bla(prev[i], prev[i + 1], dc_max));
i += 2;
}
levels.push(next);
}
let mut flat = Vec::with_capacity(levels.iter().map(|l| l.len()).sum());
for level in &levels {
flat.extend(level.iter().map(|b| b.to_bla()));
}
flat
}
/// Start offset of BLA level `lv` within the flat table, computed from the orbit
/// length exactly as the shader does (`count[0] = ref_len-1`, halving each
/// level). Kept here so the traversal test mirrors the GPU indexing.
#[cfg(test)]
fn bla_level_start(ref_len: usize, lv: u32) -> (usize, usize) {
let mut start = 0usize;
let mut count = ref_len - 1;
for _ in 0..lv {
start += count;
count /= 2;
}
(start, count)
}
#[cfg(test)]
mod tests {
use super::*;
@@ -365,17 +300,7 @@ mod tests {
fn reference_matches_naive_f64() {
let cr = Big::try_from(-0.75_f64).unwrap();
let ci = Big::try_from(0.1_f64).unwrap();
let points = compute_set_reference(
&cr,
&ci,
60,
200,
FractalKind::Mandelbrot,
2,
(0.0, 0.0),
(0.0, 0.0),
(0.0, 0.0),
);
let points = compute_set_reference(&cr, &ci, 60, 200, FractalKind::Mandelbrot, 2);
// Independent naive f64 orbit.
let (c_re, c_im) = (-0.75_f64, 0.1_f64);
@@ -385,14 +310,8 @@ mod tests {
// significant figures.
let tol_re = 1e-4 * (1.0 + zr.abs());
let tol_im = 1e-4 * (1.0 + zi.abs());
assert!(
(point[0] as f64 - zr).abs() < tol_re,
"re mismatch: {point:?} vs {zr}"
);
assert!(
(point[1] as f64 - zi).abs() < tol_im,
"im mismatch: {point:?} vs {zi}"
);
assert!((point[0] as f64 - zr).abs() < tol_re, "re mismatch: {point:?} vs {zr}");
assert!((point[1] as f64 - zi).abs() < tol_im, "im mismatch: {point:?} vs {zi}");
let nzr = zr * zr - zi * zi + c_re;
let nzi = 2.0 * zr * zi + c_im;
zr = nzr;
@@ -400,60 +319,12 @@ mod tests {
}
}
/// The f64 fast path (shallow views) must produce the same orbit as the
/// arbitrary-precision path, for every kind, in both planes.
#[test]
fn f64_fast_path_matches_big() {
let bits_fast = F64_MAX_PRECISION;
let bits_big = F64_MAX_PRECISION + 64;
for kind in FractalKind::ALL {
for julia in [false, true] {
let run = |bits: usize| {
let (a, b) = (big_from_f64(-0.3, bits), big_from_f64(0.2, bits));
let (jr, ji) = (big_from_f64(-0.4, bits), big_from_f64(0.55, bits));
let args = (60, bits, kind, 3, (0.1, -0.2), (0.9, 0.3), (2.3, 0.4));
if julia {
compute_reference(
&a, &b, &jr, &ji, args.0, args.1, args.2, args.3, args.4, args.5,
args.6,
)
} else {
compute_set_reference(
&a, &b, args.0, args.1, args.2, args.3, args.4, args.5, args.6,
)
}
};
let (fast, big) = (run(bits_fast), run(bits_big));
assert_eq!(fast.len(), big.len(), "{kind:?} julia={julia}: length");
for (i, (f, b)) in fast.iter().zip(&big).enumerate() {
for k in 0..2 {
let tol = 1e-5 * (1.0 + b[k].abs());
assert!(
(f[k] - b[k]).abs() <= tol,
"{kind:?} julia={julia}: point {i} {f:?} vs {b:?}"
);
}
}
}
}
}
/// A point inside the main cardioid never escapes: full-length orbit.
#[test]
fn interior_orbit_runs_full_length() {
let cr = Big::try_from(-0.2_f64).unwrap();
let ci = Big::try_from(0.0_f64).unwrap();
let points = compute_set_reference(
&cr,
&ci,
500,
120,
FractalKind::Mandelbrot,
2,
(0.0, 0.0),
(0.0, 0.0),
(0.0, 0.0),
);
let points = compute_set_reference(&cr, &ci, 500, 120, FractalKind::Mandelbrot, 2);
assert_eq!(points.len(), 501, "interior orbit should not escape");
}
@@ -462,17 +333,7 @@ mod tests {
fn burning_ship_reference_matches_naive_f64() {
let cr = Big::try_from(-1.75_f64).unwrap();
let ci = Big::try_from(-0.03_f64).unwrap();
let points = compute_set_reference(
&cr,
&ci,
60,
200,
FractalKind::BurningShip,
2,
(0.0, 0.0),
(0.0, 0.0),
(0.0, 0.0),
);
let points = compute_set_reference(&cr, &ci, 60, 200, FractalKind::BurningShip, 2);
let (c_re, c_im) = (-1.75_f64, -0.03_f64);
let (mut zr, mut zi) = (0.0_f64, 0.0_f64);
@@ -492,17 +353,7 @@ mod tests {
fn multibrot3_reference_matches_naive_f64() {
let cr = Big::try_from(0.3_f64).unwrap();
let ci = Big::try_from(0.2_f64).unwrap();
let points = compute_set_reference(
&cr,
&ci,
60,
200,
FractalKind::Multibrot,
3,
(0.0, 0.0),
(0.0, 0.0),
(0.0, 0.0),
);
let points = compute_set_reference(&cr, &ci, 60, 200, FractalKind::Multibrot, 3);
let (c_re, c_im) = (0.3_f64, 0.2_f64);
let (mut zr, mut zi) = (0.0_f64, 0.0_f64);
@@ -535,9 +386,6 @@ mod tests {
200,
FractalKind::Mandelbrot,
2,
(0.0, 0.0),
(0.0, 0.0),
(0.0, 0.0),
);
let (mut zr, mut zi) = (0.15_f64, -0.1_f64);
@@ -553,179 +401,125 @@ mod tests {
}
}
/// Celtic reference matches a naive f64 iteration: real = |x^2 - y^2| + cr.
/// Step-by-step perturbation (drops nothing): `e_{n+1} = 2 Z_n e_n + e_n^2 + dc`,
/// using the stored f32 orbit as `Z_n`. Returns `e` after `target_n` steps.
fn advance_naive(points: &[[f32; 2]], dc: (f64, f64), target_n: usize) -> (f64, f64) {
let (mut er, mut ei) = (0.0f64, 0.0f64);
for z in &points[..target_n] {
let (zr, zi) = (z[0] as f64, z[1] as f64);
let tr = 2.0 * (zr * er - zi * ei);
let ti = 2.0 * (zr * ei + zi * er);
let sr = er * er - ei * ei;
let si = 2.0 * er * ei;
er = tr + sr + dc.0;
ei = ti + si + dc.1;
}
(er, ei)
}
/// Advance `e` to exactly `target_n` steps using the BLA table — the same
/// walk the shader performs (longest valid skip first, else a full step),
/// but never skipping past `target_n`. Returns `(e, took_a_multi_step_skip)`.
fn advance_bla(
points: &[[f32; 2]],
table: &[Bla],
dc: (f64, f64),
target_n: usize,
) -> ((f64, f64), bool) {
let ref_len = points.len();
let (mut er, mut ei) = (0.0f64, 0.0f64);
let mut n = 0usize;
let mut skipped = false;
while n < target_n {
let emag2 = er * er + ei * ei;
let mut applied = false;
// Highest level worth trying is bounded by how far we may advance.
let span = target_n - n;
let max_lv = (usize::BITS - 1 - span.leading_zeros()) as i64; // floor(log2(span))
let mut lv = max_lv;
while lv >= 0 {
let lvu = lv as u32;
let step = 1usize << lvu;
if n % step == 0 && n + step <= target_n {
let (start, count) = bla_level_start(ref_len, lvu);
let idx = n >> lvu;
if idx < count {
let b = table[start + idx];
let r = b.r as f64;
if emag2 < r * r {
let (ar, ai) = (b.a[0] as f64, b.a[1] as f64);
let (br, bi) = (b.b[0] as f64, b.b[1] as f64);
let ner = ar * er - ai * ei + br * dc.0 - bi * dc.1;
let nei = ar * ei + ai * er + br * dc.1 + bi * dc.0;
er = ner;
ei = nei;
n += step;
skipped |= step > 1;
applied = true;
break;
}
}
}
lv -= 1;
}
if !applied {
let (zr, zi) = (points[n][0] as f64, points[n][1] as f64);
let tr = 2.0 * (zr * er - zi * ei);
let ti = 2.0 * (zr * ei + zi * er);
let sr = er * er - ei * ei;
let si = 2.0 * er * ei;
er = tr + sr + dc.0;
ei = ti + si + dc.1;
n += 1;
}
}
((er, ei), skipped)
}
/// The BLA walk must reproduce step-by-step perturbation at deep zoom (where
/// the delta is tiny and skips actually fire).
#[test]
fn celtic_reference_matches_naive_f64() {
let cr = Big::try_from(-0.6_f64).unwrap();
let ci = Big::try_from(0.4_f64).unwrap();
let points = compute_set_reference(
&cr,
&ci,
60,
200,
FractalKind::Celtic,
2,
(0.0, 0.0),
(0.0, 0.0),
(0.0, 0.0),
);
let (c_re, c_im) = (-0.6_f64, 0.4_f64);
let (mut zr, mut zi) = (0.0_f64, 0.0_f64);
for point in &points {
let tol = 1e-4 * (1.0 + zr.abs().max(zi.abs()));
assert!((point[0] as f64 - zr).abs() < tol, "re: {point:?} vs {zr}");
assert!((point[1] as f64 - zi).abs() < tol, "im: {point:?} vs {zi}");
let nzr = (zr * zr - zi * zi).abs() + c_re;
let nzi = 2.0 * zr * zi + c_im;
zr = nzr;
zi = nzi;
}
}
/// Perpendicular reference matches a naive f64 iteration:
/// real = x^2 - y^2 + cr, imag = -2·x·|y| + ci.
#[test]
fn perpendicular_reference_matches_naive_f64() {
let cr = Big::try_from(-0.7_f64).unwrap();
let ci = Big::try_from(-0.2_f64).unwrap();
let points = compute_set_reference(
&cr,
&ci,
60,
200,
FractalKind::Perpendicular,
2,
(0.0, 0.0),
(0.0, 0.0),
(0.0, 0.0),
);
let (c_re, c_im) = (-0.7_f64, -0.2_f64);
let (mut zr, mut zi) = (0.0_f64, 0.0_f64);
for point in &points {
let tol = 1e-4 * (1.0 + zr.abs().max(zi.abs()));
assert!((point[0] as f64 - zr).abs() < tol, "re: {point:?} vs {zr}");
assert!((point[1] as f64 - zi).abs() < tol, "im: {point:?} vs {zi}");
let nzr = zr * zr - zi * zi + c_re;
let nzi = -2.0 * zr * zi.abs() + c_im;
zr = nzr;
zi = nzi;
}
}
/// Buffalo reference matches a naive f64 iteration:
/// real = |x^2 - y^2| + cr, imag = -|2·x·y| + ci.
#[test]
fn buffalo_reference_matches_naive_f64() {
let cr = Big::try_from(-1.2_f64).unwrap();
let ci = Big::try_from(-0.35_f64).unwrap();
let points = compute_set_reference(
&cr,
&ci,
60,
200,
FractalKind::Buffalo,
2,
(0.0, 0.0),
(0.0, 0.0),
(0.0, 0.0),
);
let (c_re, c_im) = (-1.2_f64, -0.35_f64);
let (mut zr, mut zi) = (0.0_f64, 0.0_f64);
for point in &points {
let tol = 1e-4 * (1.0 + zr.abs().max(zi.abs()));
assert!((point[0] as f64 - zr).abs() < tol, "re: {point:?} vs {zr}");
assert!((point[1] as f64 - zi).abs() < tol, "im: {point:?} vs {zi}");
let nzr = (zr * zr - zi * zi).abs() + c_re;
let nzi = -(2.0 * zr * zi).abs() + c_im;
zr = nzr;
zi = nzi;
}
}
/// Phoenix reference matches a naive f64 two-term iteration
/// `z_{n+1} = z_n^2 + c + p·z_{n-1}` (z_0 = 0, z_{-1} = 0).
#[test]
fn phoenix_reference_matches_naive_f64() {
let cr = Big::try_from(0.5667_f64).unwrap();
fn bla_matches_step_by_step() {
// An interior center (never escapes), so `e` stays bounded and we can
// iterate the full orbit; zoomed so |dc| ~ 1e-9 (deep enough for big
// skips). Correctness of the walk is independent of which orbit we pick.
let cr = Big::try_from(-0.5_f64).unwrap();
let ci = Big::try_from(0.0_f64).unwrap();
let p = (-0.5_f64, 0.0_f64);
let points = compute_set_reference(
&cr,
&ci,
60,
200,
FractalKind::Phoenix,
2,
p,
(0.0, 0.0),
let points = compute_set_reference(&cr, &ci, 800, 160, FractalKind::Mandelbrot, 2);
assert!(points.len() > 64, "need a long orbit to exercise BLA levels");
let half_height = 1.0e-9_f64;
let aspect = 1.5_f64;
let dc_max = half_height * (1.0 + aspect * aspect).sqrt();
let table = build_bla_table(&points, dc_max);
assert!(!table.is_empty());
let target_n = points.len() - 1;
// A few pixel offsets across the view (all within dc_max).
let offsets = [
(0.0, 0.0),
(0.6e-9, -0.4e-9),
(-0.9e-9, 0.3e-9),
(0.2e-9, 0.8e-9),
];
let mut any_skip = false;
for dc in offsets {
let naive = advance_naive(&points, dc, target_n);
let (bla, skipped) = advance_bla(&points, &table, dc, target_n);
any_skip |= skipped;
let ref_mag = (naive.0 * naive.0 + naive.1 * naive.1).sqrt();
let err = ((bla.0 - naive.0).powi(2) + (bla.1 - naive.1).powi(2)).sqrt();
// Only the dropped e^2 differs; must stay near the BLA_EPS budget.
let tol = 1e-4 * ref_mag + 1e-15;
assert!(
err <= tol,
"dc={dc:?}: BLA {bla:?} vs naive {naive:?} (err {err:e} > tol {tol:e})"
);
let (c_re, c_im) = (0.5667_f64, 0.0_f64);
let (mut zr, mut zi) = (0.0_f64, 0.0_f64);
let (mut pr, mut pi) = (0.0_f64, 0.0_f64); // previous iterate
for point in &points {
let tol = 1e-4 * (1.0 + zr.abs().max(zi.abs()));
assert!((point[0] as f64 - zr).abs() < tol, "re: {point:?} vs {zr}");
assert!((point[1] as f64 - zi).abs() < tol, "im: {point:?} vs {zi}");
// p·z_{n-1} = (p.0 + i p.1)(pr + i pi).
let pzr = p.0 * pr - p.1 * pi;
let pzi = p.0 * pi + p.1 * pr;
let nzr = zr * zr - zi * zi + c_re + pzr;
let nzi = 2.0 * zr * zi + c_im + pzi;
pr = zr;
pi = zi;
zr = nzr;
zi = nzi;
}
}
/// Complex Multibrot (power 2.5 + 0.3i) reference matches a naive f64
/// iteration of `z^p = exp(p·ln z)`.
#[test]
fn complex_multibrot_reference_matches_naive_f64() {
let cr = Big::try_from(0.1_f64).unwrap();
let ci = Big::try_from(-0.2_f64).unwrap();
let power = (2.5_f64, 0.3_f64);
let points = compute_set_reference(
&cr,
&ci,
60,
200,
FractalKind::ComplexMultibrot,
2,
(0.0, 0.0),
(0.0, 0.0),
power,
);
// Naive f64 complex power via z^p = exp(p * ln z), ln z = ln|z| + i*arg(z).
fn naive_cpow(zr: f64, zi: f64, pr: f64, pi: f64) -> (f64, f64) {
if zr == 0.0 && zi == 0.0 {
return (0.0, 0.0);
}
let ln_r = 0.5 * (zr * zr + zi * zi).ln();
let theta = zi.atan2(zr);
let exp_re = pr * ln_r - pi * theta;
let exp_im = pr * theta + pi * ln_r;
let mag = exp_re.exp();
(mag * exp_im.cos(), mag * exp_im.sin())
}
let (c_re, c_im) = (0.1_f64, -0.2_f64);
let (mut zr, mut zi) = (0.0_f64, 0.0_f64);
for point in &points {
let tol = 1e-4 * (1.0 + zr.abs().max(zi.abs()));
assert!((point[0] as f64 - zr).abs() < tol, "re: {point:?} vs {zr}");
assert!((point[1] as f64 - zi).abs() < tol, "im: {point:?} vs {zi}");
let (pr, pi) = naive_cpow(zr, zi, power.0, power.1);
let nzr = pr + c_re;
let nzi = pi + c_im;
zr = nzr;
zi = nzi;
}
assert!(any_skip, "BLA never took a multi-step skip — test is not exercising it");
}
}
+180 -719
View File
File diff suppressed because it is too large Load Diff
+16 -39
View File
@@ -21,40 +21,31 @@ pub struct ShareState {
pub half_height: f64,
pub iterations: u32,
pub julia_c: (f64, f64),
/// Distortion constant for the Phoenix kind (ignored by others).
pub phoenix_p: (f64, f64),
/// Distortion constant for the Lambda kind (ignored by others).
pub lambda_l: (f64, f64),
/// Complex exponent for the Complex Multibrot kind (ignored by others).
pub complex_power: (f64, f64),
pub color_scale: f32,
pub color_offset: f32,
/// Palette index (`palette_id` in the shader).
pub palette: u32,
/// Shadow palette index (`shadow_palette_id` in the shader).
pub shadow_palette: u32,
}
impl ShareState {
pub fn encode(&self) -> String {
let mut s = String::new();
s.push_str(if self.julia { "m=j" } else { "m=m" });
s.push_str(&format!("&f={}", self.kind.share_tag()));
s.push_str(&format!("&f={}", match self.kind {
FractalKind::Mandelbrot => "mandel",
FractalKind::BurningShip => "burning",
FractalKind::Multibrot => "multi",
FractalKind::Tricorn => "tricorn",
}));
s.push_str(&format!("&pw={}", self.power));
s.push_str(&format!(
"&re={}&im={}&hh={}&it={}",
self.center_re, self.center_im, self.half_height, self.iterations
));
s.push_str(&format!("&jr={}&ji={}", self.julia_c.0, self.julia_c.1));
s.push_str(&format!("&px={}&py={}", self.phoenix_p.0, self.phoenix_p.1));
s.push_str(&format!("&lx={}&ly={}", self.lambda_l.0, self.lambda_l.1));
s.push_str(&format!(
"&cpr={}&cpi={}",
self.complex_power.0, self.complex_power.1
));
s.push_str(&format!(
"&cs={}&co={}&pal={}&spal={}",
self.color_scale, self.color_offset, self.palette, self.shadow_palette
"&cs={}&co={}&pal={}",
self.color_scale, self.color_offset, self.palette
));
s
}
@@ -72,7 +63,13 @@ impl ShareState {
julia: map.get("m").map(|m| *m == "j").unwrap_or(false),
kind: map
.get("f")
.and_then(|f| FractalKind::from_share_tag(f))
.map(|f| match *f {
"mandel" => FractalKind::Mandelbrot,
"multi" => FractalKind::Multibrot,
"burning" => FractalKind::BurningShip,
"tricorn" => FractalKind::Tricorn,
_ => FractalKind::Mandelbrot
})
.unwrap_or(FractalKind::Mandelbrot),
power: map.get("pw").and_then(|s| s.parse().ok()).unwrap_or(2),
center_re: (*map.get("re")?).to_string(),
@@ -83,22 +80,9 @@ impl ShareState {
map.get("jr").and_then(|s| s.parse().ok()).unwrap_or(-0.8),
map.get("ji").and_then(|s| s.parse().ok()).unwrap_or(0.156),
),
phoenix_p: (
map.get("px").and_then(|s| s.parse().ok()).unwrap_or(-0.5),
map.get("py").and_then(|s| s.parse().ok()).unwrap_or(0.0),
),
lambda_l: (
map.get("lx").and_then(|s| s.parse().ok()).unwrap_or(-0.5),
map.get("ly").and_then(|s| s.parse().ok()).unwrap_or(0.0),
),
complex_power: (
map.get("cpr").and_then(|s| s.parse().ok()).unwrap_or(2.0),
map.get("cpi").and_then(|s| s.parse().ok()).unwrap_or(0.0),
),
color_scale: map.get("cs").and_then(|s| s.parse().ok()).unwrap_or(0.02),
color_offset: map.get("co").and_then(|s| s.parse().ok()).unwrap_or(0.0),
palette: map.get("pal").and_then(|s| s.parse().ok()).unwrap_or(0),
shadow_palette: map.get("spal").and_then(|s| s.parse().ok()).unwrap_or(0),
})
}
}
@@ -111,20 +95,16 @@ mod tests {
fn round_trip() {
let s = ShareState {
julia: true,
kind: FractalKind::Phoenix,
kind: FractalKind::Multibrot,
power: 5,
center_re: "-0.743643887037158704752191506114774".into(),
center_im: "0.131825904205311970493132056385139".into(),
half_height: 1.5e-20,
iterations: 4000,
julia_c: (-0.123, 0.745),
phoenix_p: (-0.5, 0.1),
lambda_l: (-0.5, 0.0),
complex_power: (2.5, 0.3),
color_scale: 0.02,
color_offset: 0.25,
palette: 3,
shadow_palette: 1,
};
let d = ShareState::decode(&s.encode()).unwrap();
assert_eq!(d.julia, s.julia);
@@ -135,10 +115,7 @@ mod tests {
assert_eq!(d.half_height, s.half_height);
assert_eq!(d.iterations, s.iterations);
assert_eq!(d.julia_c, s.julia_c);
assert_eq!(d.phoenix_p, s.phoenix_p);
assert_eq!(d.complex_power, s.complex_power);
assert_eq!(d.palette, s.palette);
assert_eq!(d.shadow_palette, s.shadow_palette);
}
#[test]
-245
View File
@@ -1,245 +0,0 @@
// Headless PNG rendering: parse the CLI, build the exact same view/state the
// windowed app would from it, then render straight to a file. No window, no
// event loop, no worker-thread debounce (nothing to debounce for a one-shot
// render); it just creates its own wgpu device, computes the reference orbit
// once, and renders through the same `ExportRender` path the "Export PNG"
// button uses.
use eframe::egui_wgpu::wgpu;
use crate::app::{FractalApp, unix_timestamp};
use crate::cli::Cli;
use crate::fractal::{ExportRender, FractalRenderer, ShareState, export_to_png_blocking};
use crate::view::{
ViewState, big_from_decimal_str, interpolate_f64, interpolate_view, parse_view_spec,
precision_for,
};
/// Cap on the output image dimension (px), to stay within GPU texture limits.
const MAX_DIM: u32 = 8192 * 16;
pub fn run(cli: Cli) -> Result<(), String> {
if cli.buddhabrot {
return Err("headless mode doesn't support --buddhabrot yet".into());
}
let width = cli.width.clamp(16, MAX_DIM);
let height = cli.height.clamp(16, MAX_DIM);
// These drive the animation path below; grab them before `apply_cli`
// consumes `cli` to build the start state.
let to_view = cli.to_view.clone();
let to_share = cli.to_share.clone();
let to_iterations = cli.to_iterations;
let frames_arg = cli.frames;
let fps = cli.fps;
let duration = cli.duration;
let linear = cli.linear;
let export_path = cli.export_path.clone();
let mut app = FractalApp::default_state();
app.apply_cli(cli);
if to_view.is_some() || to_share.is_some() {
return run_animation(
app,
to_view,
to_share,
to_iterations,
frames_arg,
fps,
duration,
linear,
width,
height,
export_path,
);
}
let export_path = export_path.unwrap_or_else(|| format!("fractal-{}.png", unix_timestamp()));
eprintln!("computing reference orbit…");
app.compute_reference_blocking();
let (device, queue) = pollster::block_on(request_device())?;
let format = wgpu::TextureFormat::Bgra8Unorm;
let renderer = FractalRenderer::new(&device, format);
let uniforms = app.make_uniforms(width as f64 / height as f64);
let (pipeline, bind_group_layout, format) = renderer.export_handles(&device, &uniforms);
let er = ExportRender::new(
&device,
&queue,
pipeline,
&bind_group_layout,
format,
width,
height,
uniforms,
app.reference_points(),
app.lights(),
);
eprintln!("rendering {width}×{height}…");
let png = export_to_png_blocking(&device, &queue, &er, |phase, fraction| {
eprint!("\r{phase} {:>3.0}%", fraction * 100.0);
});
eprintln!();
std::fs::write(&export_path, &png).map_err(|e| format!("save failed: {e}"))?;
println!("saved {export_path} ({width}×{height})");
Ok(())
}
/// Render a sequence of frames sweeping the camera from the app's current
/// (start) view to an end view, for feeding into ffmpeg. Everything other
/// than the view (kind, colors, iteration cap policy, ...) stays fixed at
/// whatever `apply_cli` set up for the start; only the camera moves.
#[allow(clippy::too_many_arguments)]
fn run_animation(
mut app: FractalApp,
to_view: Option<String>,
to_share: Option<String>,
mut to_iterations: Option<u32>,
frames_arg: Option<u32>,
fps: f64,
duration: Option<f64>,
linear: bool,
width: u32,
height: u32,
export_path: Option<String>,
) -> Result<(), String> {
let frames = match frames_arg {
Some(n) => n,
None => {
let dur = duration.ok_or("animation needs --frames, or --duration (with --fps)")?;
((fps * dur).round() as u32).max(2)
}
};
if frames < 2 {
return Err("animation needs at least 2 frames".into());
}
let (to, to_iterations_share) =
parse_animation_target(to_view.as_deref(), to_share.as_deref())?;
if to_iterations.is_none()
&& let Some(to_iterations_share) = to_iterations_share
{
to_iterations = Some(to_iterations_share);
}
let from = app.view_state().clone();
let from_iterations = app.max_iterations();
if to_iterations.is_none() {
// Iteration count auto-scales with zoom depth per frame, the same way it
// does while zooming interactively — no need to interpolate it by hand.
app.set_auto_iterations(true);
}
let out_dir = export_path.unwrap_or_else(|| format!("frames-{}", unix_timestamp()));
std::fs::create_dir_all(&out_dir).map_err(|e| format!("failed to create {out_dir}: {e}"))?;
let (device, queue) = pollster::block_on(request_device())?;
let format = wgpu::TextureFormat::Bgra8Unorm;
let renderer = FractalRenderer::new(&device, format);
// Only the camera animates, so the shader specialization (kind, Julia,
// DE) is the same for every frame.
let (pipeline, bind_group_layout, format) =
renderer.export_handles(&device, &app.make_uniforms(width as f64 / height as f64));
for i in 0..frames {
let raw_t = i as f64 / (frames - 1) as f64;
let t = if linear { raw_t } else { smoothstep(raw_t) };
if let Some(to) = to_iterations {
app.set_max_iterations(
interpolate_f64(from_iterations as f64, to as f64, t).round() as u32,
);
}
app.set_view(interpolate_view(&from, &to, t));
eprintln!("[{:>4}/{frames}] computing reference orbit…", i + 1);
app.compute_reference_blocking();
let uniforms = app.make_uniforms(width as f64 / height as f64);
let er = ExportRender::new(
&device,
&queue,
pipeline.clone(),
&bind_group_layout,
format,
width,
height,
uniforms,
app.reference_points(),
app.lights(),
);
let png = export_to_png_blocking(&device, &queue, &er, |phase, fraction| {
eprint!(
"\r[{:>4}/{frames}] {phase} {:>3.0}%",
i + 1,
fraction * 100.0
);
});
eprintln!();
let path = format!("{out_dir}/frame-{:05}.png", i + 1);
std::fs::write(&path, &png).map_err(|e| format!("save failed: {e}"))?;
}
println!("saved {frames} frames to {out_dir}/ ({width}×{height})");
println!(
"tip: ffmpeg -framerate {fps} -i {out_dir}/frame-%05d.png -c:v libx264 -pix_fmt yuv420p out.mp4"
);
Ok(())
}
/// Parse `--to-view`/`--to-share` (exactly one must be set) into the end
/// view of an animation. Only position/zoom/iterations are pulled from a
/// share fragment — the rest of its state (kind, colors, ...) is ignored, so
/// pasting a link from the app doesn't unexpectedly change the fractal kind
/// mid-animation.
fn parse_animation_target(
to_view: Option<&str>,
to_share: Option<&str>,
) -> Result<(ViewState, Option<u32>), String> {
if let Some(spec) = to_view {
return parse_view_spec(spec).ok_or_else(|| format!("invalid --to-view spec: {spec}"));
}
let frag = to_share.expect("run_animation only called with one of to_view/to_share set");
let state =
ShareState::decode(frag).ok_or_else(|| format!("invalid --to-share fragment: {frag}"))?;
let bits = precision_for(state.half_height);
let re =
big_from_decimal_str(&state.center_re, bits).ok_or("invalid --to-share center (re)")?;
let im =
big_from_decimal_str(&state.center_im, bits).ok_or("invalid --to-share center (im)")?;
Ok((
ViewState::with_center(re, im, state.half_height),
Some(state.iterations),
))
}
/// Ease-in/ease-out pacing: slow at both ends, fast through the middle.
fn smoothstep(t: f64) -> f64 {
t * t * (3.0 - 2.0 * t)
}
/// Set up a wgpu device with no surface/window attached, matching the limits
/// `main::wgpu_options` requests for the windowed app (the fractal fragment
/// shader needs storage buffers, which downlevel/WebGL-style limits disallow).
async fn request_device() -> Result<(wgpu::Device, wgpu::Queue), String> {
let instance = wgpu::Instance::default();
let adapter = instance
.request_adapter(&wgpu::RequestAdapterOptions::default())
.await
.map_err(|e| format!("no compatible GPU adapter: {e}"))?;
adapter
.request_device(&wgpu::DeviceDescriptor {
label: Some("headless fractal device"),
required_features: wgpu::Features::empty(),
required_limits: adapter.limits(),
..Default::default()
})
.await
.map_err(|e| format!("failed to create device: {e}"))
}
-94
View File
@@ -1,94 +0,0 @@
use std::f32::consts::PI;
use bytemuck::{Pod, Zeroable};
use egui::{Color32, Ui};
/// Maximum number of simultaneous lights.
pub const MAX_LIGHT_COUNT: usize = 16;
#[derive(Clone, Copy, PartialEq, Zeroable, Pod)]
#[repr(C)]
pub struct Light {
pub azimuth: f32,
pub altitude: f32,
pub color: Color32,
pub _pad: u32,
}
impl Default for Light {
fn default() -> Self {
Self {
azimuth: PI / 4.,
altitude: PI / 4.,
color: Color32::WHITE,
_pad: 0,
}
}
}
impl Light {
pub fn widget(&mut self, ui: &mut Ui) -> bool {
let formater = |v, _| format!("{}°", ((v * 180. / std::f64::consts::PI) as u32));
let parser = |s: &str| {
s.parse::<u32>()
.ok()
.map(|x| x as f64 * std::f64::consts::PI / 180.)
};
ui.horizontal(|ui| {
let del = ui.button("-").clicked();
ui.label("color:");
ui.color_edit_button_srgba(&mut self.color);
ui.label("θ:");
ui.add(
egui::DragValue::new(&mut self.azimuth)
.range(0.0..=PI * 2.)
.custom_formatter(formater)
.custom_parser(parser)
.speed(0.02),
);
ui.label("φ:");
ui.add(
egui::DragValue::new(&mut self.altitude)
.range(0.0..=PI / 2.)
.custom_formatter(formater)
.custom_parser(parser)
.speed(0.02),
);
del
})
.inner
}
}
/// GPU-side light, matching WGSL `Light` in `iterate_uniforms.wgsl`: the unit
/// direction toward the light (precomputed from azimuth/altitude so the
/// shader does no per-pixel trig) plus the packed RGBA colour, whose alpha is
/// the intensity. 16 bytes, so `array<Light, 16>` has a uniform-legal stride.
#[derive(Clone, Copy, PartialEq, Zeroable, Pod, Default)]
#[repr(C)]
pub struct GpuLight {
pub dir: [f32; 3],
pub color: Color32,
}
/// The light buffer's contents: the UI lights with a non-zero colour (the
/// only ones that contribute, and the ones the filmic white point counts),
/// packed to the front, plus how many there are (`Uniforms::light_count`).
pub fn gpu_lights(lights: &[Light]) -> ([GpuLight; MAX_LIGHT_COUNT], u32) {
let mut out = [GpuLight::default(); MAX_LIGHT_COUNT];
let mut n = 0;
for l in lights.iter().filter(|l| l.color != Color32::TRANSPARENT) {
if n == MAX_LIGHT_COUNT {
break;
}
let (sa, ca) = l.altitude.sin_cos();
let (sz, cz) = l.azimuth.sin_cos();
out[n] = GpuLight {
dir: [cz * ca, sz * ca, sa],
color: l.color,
};
n += 1;
}
(out, n as u32)
}
+6 -26
View File
@@ -1,7 +1,3 @@
// Without these, rust fails to infer Send/Sync trait impls
// Probably caused by the new trait solver
#![recursion_limit = "256"]
// Fractal Explorer — Rust + wgpu + egui + WGSL deep-zoom Mandelbrot.
//
// A single binary drives both native and web (WASM/WebGPU) builds; the two
@@ -9,15 +5,9 @@
// and calls the wasm `main`, which boots eframe onto the page's <canvas>.
mod app;
mod camera;
mod fractal;
mod lights;
mod view;
#[cfg(not(target_arch = "wasm32"))]
mod cli;
#[cfg(not(target_arch = "wasm32"))]
mod headless;
#[cfg(not(target_arch = "wasm32"))]
mod worker;
@@ -35,12 +25,13 @@ fn wgpu_options() -> eframe::egui_wgpu::WgpuConfiguration {
let mut options = eframe::egui_wgpu::WgpuConfiguration::default();
if let WgpuSetup::CreateNew(setup) = &mut options.wgpu_setup {
setup.device_descriptor =
std::sync::Arc::new(|adapter: &wgpu::Adapter| wgpu::DeviceDescriptor {
setup.device_descriptor = std::sync::Arc::new(|adapter: &wgpu::Adapter| {
wgpu::DeviceDescriptor {
label: Some("fractal wgpu device"),
required_features: wgpu::Features::empty(),
required_limits: adapter.limits(),
..Default::default()
}
});
#[cfg(target_arch = "wasm32")]
{
@@ -52,24 +43,11 @@ fn wgpu_options() -> eframe::egui_wgpu::WgpuConfiguration {
#[cfg(not(target_arch = "wasm32"))]
fn main() -> eframe::Result {
use clap::Parser as _;
env_logger::builder()
.filter_level(log::LevelFilter::Info)
.parse_default_env()
.init();
let cli = cli::Cli::parse();
if cli.headless {
return match headless::run(cli) {
Ok(()) => Ok(()),
Err(e) => {
eprintln!("error: {e}");
std::process::exit(1);
}
};
}
let native_options = eframe::NativeOptions {
renderer: eframe::Renderer::Wgpu,
wgpu_options: wgpu_options(),
@@ -123,7 +101,9 @@ fn main() {
match result {
Ok(_) => loading.remove(),
Err(e) => {
loading.set_inner_html(&format!("<p>The app has crashed.</br>{e:?}</p>"));
loading.set_inner_html(
"<p>The app has crashed. See the developer console for details.</p>",
);
log::error!("failed to start eframe: {e:?}");
}
}
+6 -1
View File
@@ -13,7 +13,12 @@ struct VsOut {
@vertex
fn vs_main(@builtin(vertex_index) idx: u32) -> VsOut {
let p = fullscreen_triangle_pos(idx);
var verts = array<vec2<f32>, 3>(
vec2<f32>(-1.0, -1.0),
vec2<f32>(3.0, -1.0),
vec2<f32>(-1.0, 3.0),
);
let p = verts[idx];
var out: VsOut;
out.pos = vec4<f32>(p, 0.0, 1.0);
// Map NDC to texture UV. v is flipped so the cache's top row (rendered at
-289
View File
@@ -1,289 +0,0 @@
// Buddhabrot / Nebulabrot rendering: a Monte-Carlo density histogram of
// escaping orbits, accumulated progressively across frames by a compute pass,
// then tone-mapped to colour by a fragment pass every frame.
//
// This does NOT use the deep-zoom perturbation/reference-orbit machinery in
// mandelbrot.wgsl: Buddhabrot's structure is a global Monte-Carlo property of
// the whole basin (a random sample's orbit scatters across the *whole* image,
// not just its own pixel), so the "gather" per-pixel model doesn't apply, and
// deep zoom isn't meaningful for it the way it is for the escape-time set.
// Samples are iterated directly in f32 from the current view's bounds.
//
// Sampling convention: for KIND_LAMBDA the formula z -> l*z*(1-z) has no `c`
// term at all (l is a fixed distortion constant, not a per-sample parameter),
// so the randomly sampled point instead seeds z0 (a "Julia-Buddhabrot" over
// z0 with l fixed). Every other kind samples c with z0 = 0, matching its
// ordinary parameter plane.
//
// A sample's orbit is only plotted if it escapes within b_cap iterations (the
// classic Buddhabrot rule: only escaping orbits are drawn). Its points are
// then splat into up to three histogram channels by cap (r_cap <= g_cap <=
// b_cap): fast-escaping (common) orbits light all three channels (bright),
// slow-escaping (rare) orbits only light the b_cap channel — the classic
// Nebulabrot false-colour split.
//
// Two-pass iteration avoids needing a per-thread orbit buffer sized to
// max_iter: the first pass just finds the escape iteration (if any); the
// second replays the same orbit from scratch, splatting each point.
struct Uniforms {
center: vec2<f32>,
half_height: f32,
aspect: f32,
phoenix_p: vec2<f32>,
lambda_l: vec2<f32>,
bailout_sq: f32,
kind: u32,
power: u32,
r_cap: u32,
g_cap: u32,
b_cap: u32,
seed: u32,
samples_this_dispatch: u32,
exposure: f32,
width: u32,
height: u32,
total_samples: f32,
// Tonemap colour style: 0 = classic (R/G/B = raw caps), 1 = nebula
// (yellow core, blue halo), 2 = grayscale.
palette: u32,
// Padding so `complex_power` (a vec2, 8-byte aligned) starts on an
// 8-byte boundary. NOT vec3<u32> — that type aligns to 16 bytes in WGSL
// (unlike Rust's `[u32; 3]`, which aligns to 4), which silently added 32
// bytes instead of 16 and mismatched the Rust struct's size (a wgpu
// validation error at dispatch time: "size 96 where the shader expects
// 112").
_pad0: u32,
// Complex exponent for the Complex Multibrot kind; unused by other kinds.
complex_power: vec2<f32>,
};
// Fractal kind, as a pipeline-overridable constant (set per compute pipeline
// from `u.kind`, see `BuddhabrotRenderer::compute_pipeline`): every kind
// branch in the iteration loop folds away at pipeline creation. Read this,
// never `u.kind`.
override KIND: u32 = 0u;
const PALETTE_NEBULA: u32 = 0u;
const PALETTE_YELLOW: u32 = 1u;
const PALETTE_GRAYSCALE: u32 = 2u;
@group(0) @binding(0) var<uniform> u: Uniforms;
// Compute pass: read-write atomic histogram (3 planes of width*height, R/G/B).
@group(0) @binding(1) var<storage, read_write> histogram: array<atomic<u32>>;
// Tonemap pass: read-only plain view of the same buffer.
@group(0) @binding(2) var<storage, read> tm_histogram: array<u32>;
// --- RNG: a small, fast integer hash (WGSL has no native RNG). ---
fn hash_u32(x: u32) -> u32 {
var h = x;
h = h ^ (h >> 16u);
h = h * 0x7feb352du;
h = h ^ (h >> 15u);
h = h * 0x846ca68bu;
h = h ^ (h >> 16u);
return h;
}
fn rand01(seed: u32) -> f32 {
return f32(hash_u32(seed)) * (1.0 / 4294967295.0);
}
fn complex_pow(z: vec2<f32>, p: u32) -> vec2<f32> {
var r = vec2<f32>(1.0, 0.0);
for (var i: u32 = 0u; i < p; i = i + 1u) {
r = cmul(r, z);
}
return r;
}
// One iteration step z_n -> z_{n+1} for the current kind. `zp` is the
// previous iterate (z_{n-1}), used only by the Phoenix two-term recurrence.
// Must match `FractalKind` in reference.rs (the direct, non-perturbative form
// of the same formulas).
fn advance(z: vec2<f32>, zp: vec2<f32>, c: vec2<f32>) -> vec2<f32> {
if KIND == KIND_BURNING_SHIP {
return vec2<f32>(z.x * z.x - z.y * z.y, 2.0 * abs(z.x * z.y)) + c;
} else if KIND == KIND_TRICORN {
return vec2<f32>(z.x * z.x - z.y * z.y, -2.0 * z.x * z.y) + c;
} else if KIND == KIND_MULTIBROT {
return complex_pow(z, clamp(u.power, 2u, 8u)) + c;
} else if KIND == KIND_CELTIC {
return vec2<f32>(abs(z.x * z.x - z.y * z.y), 2.0 * z.x * z.y) + c;
} else if KIND == KIND_PERPENDICULAR {
return vec2<f32>(z.x * z.x - z.y * z.y, -2.0 * z.x * abs(z.y)) + c;
} else if KIND == KIND_BUFFALO {
return vec2<f32>(abs(z.x * z.x - z.y * z.y), -abs(2.0 * z.x * z.y)) + c;
} else if KIND == KIND_PHOENIX {
let sq = vec2<f32>(z.x * z.x - z.y * z.y, 2.0 * z.x * z.y);
return sq + c + cmul(u.phoenix_p, zp);
} else if KIND == KIND_LAMBDA {
// l * z * (1 - z); c is unused (see file doc comment above).
return cmul(u.lambda_l, cmul(z, vec2<f32>(1.0 - z.x, -z.y)));
} else if KIND == KIND_COMPLEX_MULTIBROT {
return cpow(z, u.complex_power) + c;
}
return vec2<f32>(z.x * z.x - z.y * z.y, 2.0 * z.x * z.y) + c; // Mandelbrot
}
// Map a complex-plane point to a flat pixel index, or -1 if outside the
// current viewport (the sampling region and the display region are the same).
//
// This must be the exact inverse of how `view.rs::pan_pixels`/`zoom_at_pixel`
// relate screen pixels to world points (those are the confirmed-correct,
// user-tested ground truth — NOT the shader-comment-derived convention tried
// here previously, which was wrong: dragging/zooming treat +y screen exactly
// like +x, no flip, so screen-down means im *increasing*, not decreasing).
fn pixel_index(p: vec2<f32>) -> i32 {
let half_w = u.half_height * u.aspect;
let uu = (p.x - u.center.x) / half_w * 0.5 + 0.5;
let vv = 0.5 + (p.y - u.center.y) / u.half_height * 0.5;
if uu < 0.0 || uu >= 1.0 || vv < 0.0 || vv >= 1.0 {
return -1;
}
let px = i32(uu * f32(u.width));
let py = i32(vv * f32(u.height));
return py * i32(u.width) + px;
}
// Splat one visited orbit point into the R/G/B histogram planes it qualifies
// for by the orbit's total escape iteration `n` (nested caps: a fast escape
// lights all three; only a slow, rare one lights just the blue plane).
fn splat(p: vec2<f32>, n: u32) {
let idx = pixel_index(p);
if idx < 0 {
return;
}
let plane = i32(u.width) * i32(u.height);
if n <= u.b_cap {
atomicAdd(&histogram[idx + 2 * plane], 1u);
}
if n <= u.g_cap {
atomicAdd(&histogram[idx + plane], 1u);
}
if n <= u.r_cap {
atomicAdd(&histogram[idx], 1u);
}
}
@compute @workgroup_size(64)
fn cs_main(@builtin(global_invocation_id) gid: vec3<u32>) {
if gid.x >= u.samples_this_dispatch {
return;
}
let base = hash_u32(gid.x ^ (u.seed * 0x9e3779b9u));
let rx = rand01(base);
let ry = rand01(hash_u32(base ^ 0x68bc21ebu));
let half_w = u.half_height * u.aspect;
let sample = vec2<f32>(
u.center.x + (rx * 2.0 - 1.0) * half_w,
u.center.y + (ry * 2.0 - 1.0) * u.half_height,
);
var c = sample;
var z0 = vec2<f32>(0.0, 0.0);
if KIND == KIND_LAMBDA {
c = vec2<f32>(0.0, 0.0); // unused by the Lambda step
z0 = sample;
}
// First pass: just find the escape iteration (if any).
var zp = vec2<f32>(0.0, 0.0);
var z = z0;
var n: u32 = 0u;
var escaped = false;
loop {
if dot(z, z) > u.bailout_sq {
escaped = true;
break;
}
if n >= u.b_cap {
break;
}
let next = advance(z, zp, c);
zp = z;
z = next;
n = n + 1u;
}
if !escaped || n == 0u {
return;
}
// Second pass: replay the same orbit, splatting each visited point.
// z0 itself is not splat: it's the same fixed point (0,0), or the sample
// itself for Lambda, for every orbit — plotting it would just spike the
// origin instead of showing the orbit's actual shape.
zp = vec2<f32>(0.0, 0.0);
z = z0;
for (var i: u32 = 0u; i < n; i = i + 1u) {
let next = advance(z, zp, c);
zp = z;
z = next;
splat(z, n);
}
}
// --- Tonemap: histogram counts -> colour, drawn as a fullscreen triangle. ---
@vertex
fn vs_main(@builtin(vertex_index) idx: u32) -> @builtin(position) vec4<f32> {
return vec4<f32>(fullscreen_triangle_pos(idx), 0.0, 1.0);
}
@fragment
fn fs_tonemap(@builtin(position) pos: vec4<f32>) -> @location(0) vec4<f32> {
let x = i32(pos.x);
let y = i32(pos.y);
if x < 0 || y < 0 || x >= i32(u.width) || y >= i32(u.height) {
return vec4<f32>(0.0, 0.0, 0.0, 1.0);
}
let idx = y * i32(u.width) + x;
let plane = i32(u.width) * i32(u.height);
let r = f32(tm_histogram[idx]);
let g = f32(tm_histogram[idx + plane]);
let b = f32(tm_histogram[idx + 2 * plane]);
// Normalize by the *average* density (total samples / pixel count) rather
// than total samples alone, so the scale stays sane across widget sizes
// and sample-dispatch rates. Buddhabrot density is extremely peaked (the
// brightest pixels run tens of times the average), so the compressive
// exponential tonemap only needs a small fraction of the average to reach
// full brightness at those peaks; 0.05 is a hand-tuned starting point,
// the exposure slider covers the rest.
let avg_density = max(u.total_samples / f32(u.width * u.height), 1.0e-6);
let scale = u.exposure * 0.05 / avg_density;
// Per-cap brightness, each already compressed to [0,1]. Nested caps mean
// r <= g <= b pointwise (every orbit counted in a smaller cap is also
// counted in every larger one), so fb alone is the full escaping-orbit
// density and fr picks out just the common, fast-escaping ones.
let fr = 1.0 - exp(-r * scale);
let fg = 1.0 - exp(-g * scale);
let fb = 1.0 - exp(-b * scale);
var col: vec3<f32>;
if u.palette == PALETTE_YELLOW {
// fr is *not* a good stand-alone brightness signal: with c sampled
// uniformly over the whole viewport, nearly every sample outside the
// set escapes within a handful of iterations and splats a couple of
// points near itself, so fr is a near-uniform wash across the entire
// image (not concentrated near the boundary the way fb is) — adding
// it directly (tried first, both raw and gamma-lifted) drags that
// wash up to full brightness and floods the background with solid
// colour. Instead use it as a *multiplicative* warm (yellow) tint on
// top of fb's brightness, so it only shows up where fb is already
// bright (i.e. real near-boundary density) and stays near-zero across
// the background (fb ≈ 0 there, so warmth * fb ≈ 0 regardless of fr).
col = vec3<f32>(
fb + fb * fr * 1.3,
fb + fb * fr * 0.6,
fb,
);
} else if u.palette == PALETTE_GRAYSCALE {
// fb is the full escaping-orbit density (the cumulative superset);
// reuse it directly as a single luminance channel.
col = vec3<f32>(fb, fb, fb);
} else {
col = vec3<f32>(fr, fg, fb); // classic: raw per-cap R/G/B
}
return vec4<f32>(clamp(col, vec3<f32>(0.0), vec3<f32>(1.0)), 1.0);
}
-171
View File
@@ -1,171 +0,0 @@
// Colourise pass: map the iteration pass's per-pixel escape data (from
// `mandelbrot.wgsl`'s `fs_data`) through the palette. This is the only
// color-dependent step, so changing the palette / colour scale / offset (e.g.
// colour cycling) re-runs just this cheap pass — the expensive perturbation
// iteration in the data texture is reused untouched.
//
// The data texture holds, per texel: R = ci (palette parameter), G = DE
// darkening factor, B = interior fraction (for boundary anti-aliasing). It is
// the same resolution as this pass's target, so we read it with `textureLoad`
// at the fragment's integer pixel coordinate (nearest — iteration data must not
// be linearly filtered across escape boundaries).
@group(0) @binding(0) var<uniform> u: Uniforms;
@group(0) @binding(1) var data_tex: texture_2d<f32>;
@group(0) @binding(2) var<uniform> lights: array<Light, 16>;
@vertex
fn vs_main(@builtin(vertex_index) idx: u32) -> @builtin(position) vec4<f32> {
return vec4<f32>(fullscreen_triangle_pos(idx), 0.0, 1.0);
}
fn shadow_fragment(pos: vec2<f32>) -> vec4<f32> {
let x = i32(pos.x);
let y = i32(pos.y);
let size = textureDimensions(data_tex);
let here = textureLoad(data_tex, vec2<i32>(x, y), 0);
if here.b != 0. {
return vec4<f32>(shadow_interior_color(), 1.0);
}
// Forward differences, except on the last column/row where x+1 / y+1
// is off the texture: fall back to a backward difference, mirrored
// (h0 + (h0 - h[-1])) so the slope keeps the sign normal_from_heights
// expects — plugging h[-1] in directly would flip the normal there.
let h0 = here.g;
var h1: f32;
if x + 1 < i32(size.x) {
h1 = textureLoad(data_tex, vec2<i32>(x + 1, y), 0).g;
} else {
h1 = 2.0 * h0 - textureLoad(data_tex, vec2<i32>(x - 1, y), 0).g;
}
var h2: f32;
if y + 1 < i32(size.y) {
h2 = textureLoad(data_tex, vec2<i32>(x, y + 1), 0).g;
} else {
h2 = 2.0 * h0 - textureLoad(data_tex, vec2<i32>(x, y - 1), 0).g;
}
let normal = normal_from_heights(h0, h1, h2);
return vec4<f32>(shadow_color(normal, here.r), 1.0);
}
@fragment
fn fs_main(@builtin(position) pos: vec4<f32>) -> @location(0) vec4<f32> {
if u.shadow == 2u {
return ray_marching(pos);
} else if u.shadow == 1u {
return shadow_fragment(pos.xy);
} else {
let d = textureLoad(data_tex, vec2<i32>(i32(pos.x), i32(pos.y)), 0);
let ci = d.r;
let de = d.g;
let interior_frac = d.b;
var col = classic_color(ci, de);
// Anti-alias the set boundary: fade toward black by the fraction of the
// pixel's sub-samples that landed in the interior.
col = col * (1.0 - interior_frac);
return vec4<f32>(col, 1.0);
}
}
// Colour of rays that miss the fractal's footprint.
const RAY_MISS: vec4<f32> = vec4<f32>(1.0, 0.0, 0.0, 1.0);
// Per-frame constants of the raymarch, computed once per pixel in
// `ray_marching` rather than on each of the up-to-100 `sdf` steps.
struct MarchConsts {
size: vec2<f32>,
// (size.x / aspect_ratio, size.y): world xy -> texel scale.
to_texel: vec2<f32>,
size_i: vec2<i32>,
inv_size_y: f32,
};
fn sdf(pos: vec3<f32>, k: MarchConsts) -> f32 {
let texture_pos_f32 = pos.xy * k.to_texel;
let texture_pos = clamp(vec2<i32>(texture_pos_f32), vec2<i32>(0, 0), k.size_i - vec2<i32>(1, 1));
let to_texture = max(-min(texture_pos_f32, vec2(0.)), max(texture_pos_f32 - k.size, vec2(0.)));
let dist_to_texture = length(to_texture) * k.inv_size_y;
let px = textureLoad(data_tex, texture_pos, 0);
let de = (px.g * k.inv_size_y) * 0.5;
// Height is measured toward -z, the side the camera sits on (it looks
// along +z), so the terrain is solid on +z: interior plateau at z = 0,
// exterior sloping away from the camera as `de` grows.
let signed_z = -pos.z;
let z = max(signed_z, 0.);
var d: f32;
if px.b != 0. {
d = z;
} else {
d = min(sqrt(z * z + de * de), signed_z + 1. - exp(-de * 5.));
}
// Outside the texture footprint, `d` is the distance from the clamped
// point q on the footprint's edge. The terrain lies over the (convex)
// footprint, so |p - x|² ≥ |q - x|² + |p - q|² for every terrain point x:
// combine in quadrature (not by adding, which overshoots). p can't be in
// the solid out here, so a negative `d` counts as 0.
if dist_to_texture > 0. {
let d_pos = max(d, 0.);
return sqrt(d_pos * d_pos + dist_to_texture * dist_to_texture);
}
return d;
}
fn ray_marching(pos: vec4<f32>) -> vec4<f32> {
let size_i = vec2<i32>(textureDimensions(data_tex));
let size = vec2<f32>(size_i);
let aspect_ratio = u.screen_dim.x / u.screen_dim.y;
let k = MarchConsts(size, vec2<f32>(size.x / aspect_ratio, size.y), size_i, 1.0 / size.y);
let in_texture = vec2<f32>(
(pos.x / size.x) * 2. - 1.,
(pos.y / size.y) * 2. - 1.,
);
var world_pos = u.camera_inv_proj * vec4<f32>(in_texture, 0., 1.0);
let ray_origin = world_pos.xyz;
let ray_dir = u.camera_direction;
// Start where the ray crosses z = 0, the topmost possible surface (the
// camera pitch is clamped short of ±90°, so ray_dir.z > 0).
let start = ray_origin - ray_dir * (ray_origin.z / ray_dir.z);
// The terrain only exists over the footprint x in [0, aspect],
// y in [0, 1]: clip the ray's xy to it up front, so rays that miss it cost
// nothing and the rest start marching at its edge. A huge finite 1/d on
// an axis the ray doesn't move along (top-down, during the 2D <-> 3D
// transition) keeps the slab maths finite.
let inv = select(1.0 / ray_dir.xy, vec2<f32>(1e30), abs(ray_dir.xy) < vec2<f32>(1e-20));
let ta = -start.xy * inv;
let tb = (vec2<f32>(aspect_ratio, 1.0) - start.xy) * inv;
let t_leave = min(max(ta.x, tb.x), max(ta.y, tb.y));
var t = max(max(min(ta.x, tb.x), min(ta.y, tb.y)), 0.0);
if t >= t_leave {
return RAY_MISS;
}
// About 1/50 of a texel at typical sizes: tighter only adds steps
// without visibly moving the hit.
let dist_threshold = 0.00001;
var hit = false;
for (var i = 0u; i < 100u; i++) {
let dist = sdf(start + t * ray_dir, k);
if dist < dist_threshold {
hit = true;
break;
}
t += dist;
// Past the footprint's far edge: nothing left to hit.
if t >= t_leave {
break;
}
}
// Shade outside the loop, so its registers don't weigh on the march.
if !hit {
return RAY_MISS;
}
return shadow_fragment((start.xy + t * ray_dir.xy) * k.to_texel);
}
-51
View File
@@ -1,51 +0,0 @@
// Shared helpers, concatenated into every shader at build time via
// `concat!`/`include_str!` (see renderer.rs / buddhabrot.rs). Keep this file
// free of anything that differs between pipelines (e.g. a `Uniforms` struct —
// mandelbrot/colorize and buddhabrot each have their own shape) since every
// shader gets the whole thing spliced in.
// Fullscreen triangle vertex position: one triangle that covers the whole
// viewport (cheaper than a quad's two), shared by every full-screen vertex
// shader in this project.
fn fullscreen_triangle_pos(idx: u32) -> vec2<f32> {
var verts = array<vec2<f32>, 3>(
vec2<f32>(-1.0, -1.0),
vec2<f32>(3.0, -1.0),
vec2<f32>(-1.0, 3.0),
);
return verts[idx];
}
// Complex multiply.
fn cmul(a: vec2<f32>, b: vec2<f32>) -> vec2<f32> {
return vec2<f32>(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
}
// z^p for a complex exponent p, via the principal branch z^p = exp(p * ln z),
// ln z = ln|z| + i*arg(z). z = 0 maps to 0 (the correct limit for the
// Re(p) > 0 region the UI exposes; ln(0) would otherwise be -inf).
fn cpow(z: vec2<f32>, p: vec2<f32>) -> vec2<f32> {
let r2 = dot(z, z);
if r2 < 1e-30 {
return vec2<f32>(0.0, 0.0);
}
let ln_r = 0.5 * log(r2);
let theta = atan2(z.y, z.x);
let mag = exp(p.x * ln_r - p.y * theta);
let ang = p.x * theta + p.y * ln_r;
return mag * vec2<f32>(cos(ang), sin(ang));
}
// Iteration formula selector, shared by the perturbation (mandelbrot.wgsl)
// and direct (buddhabrot.wgsl) iteration paths. Must match `FractalKind` in
// reference.rs.
const KIND_MANDELBROT: u32 = 0u;
const KIND_BURNING_SHIP: u32 = 1u;
const KIND_TRICORN: u32 = 2u;
const KIND_MULTIBROT: u32 = 3u;
const KIND_CELTIC: u32 = 4u;
const KIND_PERPENDICULAR: u32 = 5u;
const KIND_BUFFALO: u32 = 6u;
const KIND_PHOENIX: u32 = 7u;
const KIND_LAMBDA: u32 = 8u;
const KIND_COMPLEX_MULTIBROT: u32 = 9u;
-188
View File
@@ -1,188 +0,0 @@
// Shared by mandelbrot.wgsl (writes the per-pixel data texture) and
// colorize.wgsl (reads it): the iteration pass and the colour remap pass
// must agree on both the uniform layout and the palette function.
// Must match the Rust `Uniforms` struct in renderer.rs field-for-field,
// including padding.
struct Uniforms {
span: vec2<f32>,
max_iter: u32,
ref_len: u32,
color_offset: f32,
color_scale: f32,
bailout_sq: f32,
is_julia: u32,
palette_id: u32,
shadow_palette_id: u32,
aa_level: u32,
// Iteration formula (see the KIND_* constants in common.wgsl).
kind: u32,
// Exponent for the Multibrot kind.
power: u32,
dc_offset: vec2<f32>,
// Distortion constant p for the Phoenix map (z^2 + c + p*z_{n-1}); unused
// by other kinds. Placed by dc_offset so both vec2s stay 8-byte aligned.
phoenix_p: vec2<f32>,
// Distortion constant l for the Lambda map (l*z(1 - z_{n-1})); unused
// by other kinds.
lambda_l: vec2<f32>,
// Complex exponent for the Complex Multibrot kind (z^power + c); unused
// by other kinds.
complex_power: vec2<f32>,
// 0 = escape-time coloring, 1 = distance-estimation shading.
de_coloring: u32,
// 0 = classic colors, 1 = shadows, 2 = 3D raymarching rendering
shadow: u32,
// camera direction vector
camera_direction: vec3<f32>,
// Number of live entries at the start of `lights` (fills the vec3's tail
// padding slot).
light_count: u32,
// inverse of the camera's view-projection matrix, for reconstructing a
// world-space ray origin per pixel in the raymarcher
camera_inv_proj: mat4x4<f32>,
// Screen dimensions
screen_dim: vec2<f32>,
// Complex binomial coefficients C(complex_power, k) for k = 1..16, two per
// vec4 (k odd in .xy, k even in .zw), for the Complex Multibrot delta
// series. Precomputed on the CPU since they only depend on the power.
cm_coef: array<vec4<f32>, 8>,
};
// Smooth cyclic palettes (Inigo Quilez cosine palettes), selected by id.
fn palette(id: u32, t: f32) -> vec3<f32> {
if id == 4u {
return vec3<f32>(t, t, t); // grayscale
}
let a = vec3<f32>(0.5, 0.5, 0.5);
let b = vec3<f32>(0.5, 0.5, 0.5);
var c = vec3<f32>(1.0, 1.0, 1.0);
var d = vec3<f32>(0.00, 0.10, 0.20); // 0: amber / blue
if id == 1u {
d = vec3<f32>(0.00, 0.33, 0.67); // rainbow
} else if id == 2u {
d = vec3<f32>(0.30, 0.20, 0.20); // warm ember
} else if id == 3u {
c = vec3<f32>(1.0, 1.0, 0.5);
d = vec3<f32>(0.80, 0.90, 0.30); // lime / magenta
}
return a + b * cos(6.28318530718 * (c * t + d));
}
// Classic (non-shadow) escape colouring: palette lookup at the smoothed
// iteration count `ci`, darkened by the distance-estimate factor `de`
// (sqrt-compressed so the darkening falls off more gently near the
// boundary). Shared by the colourise pass's classic branch (colorize.wgsl,
// applied to an already-averaged data texel) and the PNG-export pass
// (mandelbrot.wgsl's `fs_color`, applied per sub-sample pre-AA) — the two
// places a fully escaped point is turned into a final pixel colour.
fn classic_color(ci: f32, de: f32) -> vec3<f32> {
let t = fract(ci * u.color_scale + u.color_offset);
return palette(u.palette_id, t) * sqrt(de);
}
// A single directional light, built on the CPU from the UI's light list
// (`GpuLight` in lights.rs): `dir` is the unit direction toward the light
// (precomputed from azimuth/altitude so the shader does no trig), `color` a
// packed RGBA8 whose alpha doubles as intensity. Only the first
// `u.light_count` entries are live, all with a non-zero colour. Each shader
// that binds a `lights: array<Light, 16>` uniform (colorize.wgsl,
// mandelbrot.wgsl's export shadow path) uses this same layout.
struct Light {
dir: vec3<f32>,
color: u32,
};
// Lambertian term for a unit `light` direction.
fn compute_light(normal: vec3<f32>, light: vec3<f32>) -> vec3<f32> {
return vec3<f32>(max(0., dot(normal, light)));
}
fn uncharted2tonemap(x: vec3<f32>) -> vec3<f32> {
let A = 0.15; // Shoulder strength
let B = 0.50; // Linear strength
let C = 0.10; // Linear angle
let D = 0.20; // Toe strength
let E = 0.02; // Toe numerator / shoarder angle/etc.
let F = 0.30; // Toe denominator
return ((x * (A * x + C * B) + D * E) / (x * (A * x + B) + D * F)) - E / F;
}
fn filmic(color: vec3<f32>, white_point: f32) -> vec3<f32> {
let exposure_bias = 2.0;
let curr = uncharted2tonemap(color * exposure_bias);
// Valeur blanche maximale de référence
let white_scale = vec3(1.0) / uncharted2tonemap(vec3(white_point));
return curr * white_scale;
}
fn s(color: vec3<f32>, k: f32, c: f32) -> vec3<f32> {
return 1. / (1. + exp(-k * (color - c)));
}
fn contrast(color: vec3<f32>, k: f32, c: f32) -> vec3<f32> {
let color_c = s(color, k, c);
return (color_c - s(vec3<f32>(0), k, c)) / (s(vec3<f32>(1), k, c) - s(vec3<f32>(0), k, c));
}
// Surface normal from three height samples (`h0` at the pixel, `h1` one pixel
// to the right, `h2` one pixel down), treating DE as a height field. Only the
// differences matter, so callers don't need to pass pixel coordinates — a
// texture-backed caller (colorize.wgsl) and a live-sampled caller
// (mandelbrot.wgsl's export shadow path) can share this.
fn normal_from_heights(h0: f32, h1: f32, h2: f32) -> vec3<f32> {
let d0 = vec3<f32>(0.0, 0.0, h0);
let d1 = vec3<f32>(1.0, 0.0, h1);
let d2 = vec3<f32>(0.0, 1.0, h2);
return normalize(cross(d1 - d0, d2 - d0));
}
// Shade a DE-derived surface normal per `u.shadow_palette_id`: 0 = grayscale
// key light, 1 = red/blue two-tone, 2 = the user's custom `lights` list,
// 3 = the classic escape-time palette at `ci` (the smoothed iteration count),
// lit by the grayscale key light.
// Shared by the interactive shadow pass (colorize.wgsl) and the PNG-export
// shadow path (mandelbrot.wgsl's `fs_color`), which must render identically.
fn shadow_color(normal: vec3<f32>, ci: f32) -> vec3<f32> {
var color: vec3<f32>;
if u.shadow_palette_id == 0u {
color = compute_light(normal, vec3<f32>(0.57735027, 0.57735027, 0.57735027)) + vec3<f32>(0.58, 0.85, 1.) * 0.2;
color = filmic(color, 2.5);
color = contrast(color, 4., 0.67);
} else if u.shadow_palette_id == 1u {
color = compute_light(normal, vec3<f32>(0., 0.70710678, 0.70710678)) * vec3<f32>(1., 0.5, 0.5) + compute_light(normal, vec3<f32>(0.70710678, 0., 0.70710678)) * vec3<f32>(0.5, 1., 1.);
color = filmic(color, 4.2);
} else if u.shadow_palette_id == 3u {
// No DE darkening as in `classic_color`: in shadow/3D modes the DE
// is a height (clamped to 1000, not 1), and the lighting already
// shows the relief.
let t = fract(ci * u.color_scale + u.color_offset);
let ambient = 0.25;
let light = compute_light(normal, vec3<f32>(0.57735027, 0.57735027, 0.57735027));
color = palette(u.palette_id, t) * (ambient + (1.0 - ambient) * light);
} else {
color = vec3<f32>(0);
let light_count = min(u.light_count, 16u);
for (var i = 0u; i < light_count; i++) {
let light_color = unpack4x8unorm(lights[i].color);
color += compute_light(normal, lights[i].dir) * light_color.xyz * light_color.a;
}
color = filmic(color, 1. + f32(light_count));
}
return color;
}
// Colour of an interior (non-escaped) pixel in shadow/3D modes: black under
// the classic palette, like classic 2D mode, otherwise a dark gray plateau.
fn shadow_interior_color() -> vec3<f32> {
if u.shadow_palette_id == 3u {
return vec3<f32>(0.0);
}
return vec3<f32>(0.1);
}
+208 -455
View File
@@ -12,25 +12,46 @@
// the reference index to 0 and carry the full value as the new delta (valid
// because X_0 = 0).
struct Uniforms {
span: vec2<f32>,
max_iter: u32,
ref_len: u32,
color_offset: f32,
color_scale: f32,
bailout_sq: f32,
is_julia: u32,
palette_id: u32,
aa_level: u32,
// Iteration formula: 0 Mandelbrot, 1 Burning Ship, 2 Tricorn, 3 Multibrot.
kind: u32,
// Exponent for the Multibrot kind.
power: u32,
dc_offset: vec2<f32>,
// 0 = escape-time coloring, 1 = distance-estimation shading.
de_coloring: u32,
// 0 = per-iteration stepping, 1 = BLA iteration-skipping (square map only).
use_bla: u32,
};
// One merged linear step from the BLA table: e_{n+l} = A e + B dc, valid while
// |e| < r. Matches the Rust `Bla` struct (24-byte std430 stride).
struct Bla {
a: vec2<f32>,
b: vec2<f32>,
r: f32,
l: u32,
};
const KIND_MANDELBROT: u32 = 0u;
const KIND_BURNING_SHIP: u32 = 1u;
const KIND_TRICORN: u32 = 2u;
const KIND_MULTIBROT: u32 = 3u;
@group(0) @binding(0) var<uniform> u: Uniforms;
@group(0) @binding(1) var<storage, read> ref_orbit: array<vec2<f32>>;
// Only read by `fs_color`'s shadow branch (custom-lights palette); the
// iteration pass (`fs_data`) never touches it.
@group(0) @binding(2) var<uniform> lights: array<Light, 16>;
// Only read by the adaptive-AA refine pass (`fs_refine`): the 1-sample-per-
// pixel data texture written by `fs_data`, which decides where to supersample.
@group(1) @binding(0) var coarse_tex: texture_2d<f32>;
// Pipeline-overridable specialization constants, set per pipeline from the
// uniforms' `kind` / `is_julia` / `de_coloring` (see `PipelineKey` in
// renderer.rs). Every per-iteration branch on them folds away at pipeline
// creation, so the hot loop only contains the current kind's math instead of
// testing all of them on every step. The matching uniform fields are still
// uploaded (the layout is shared with colorize.wgsl) but this shader must read
// these constants, never `u.kind` / `u.is_julia` / `u.de_coloring`.
override KIND: u32 = 0u;
override IS_JULIA: bool = false;
override DE: bool = false;
// BLA table, levels concatenated (level 0 first). Per-level counts are recomputed
// from `ref_len` exactly as the CPU builder laid them out.
@group(0) @binding(2) var<storage, read> bla_table: array<Bla>;
struct VsOut {
@builtin(position) pos: vec4<f32>,
@@ -40,7 +61,12 @@ struct VsOut {
@vertex
fn vs_main(@builtin(vertex_index) idx: u32) -> VsOut {
let ndc = fullscreen_triangle_pos(idx);
var verts = array<vec2<f32>, 3>(
vec2<f32>(-1.0, -1.0),
vec2<f32>(3.0, -1.0),
vec2<f32>(-1.0, 3.0),
);
let ndc = verts[idx];
var out: VsOut;
out.pos = vec4<f32>(ndc, 0.0, 1.0);
// Flip y so +imaginary points up the screen.
@@ -48,100 +74,60 @@ fn vs_main(@builtin(vertex_index) idx: u32) -> VsOut {
return out;
}
// Complex multiply.
fn cmul(a: vec2<f32>, b: vec2<f32>) -> vec2<f32> {
return vec2<f32>(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
}
// Complex conjugate.
fn conj(a: vec2<f32>) -> vec2<f32> {
return vec2<f32>(a.x, -a.y);
}
// Complex division a / b.
fn cdiv(a: vec2<f32>, b: vec2<f32>) -> vec2<f32> {
let d = dot(b, b);
return vec2<f32>(a.x * b.x + a.y * b.y, a.y * b.x - a.x * b.y) / d;
}
// |c + d| - |c|, evaluated exactly (no catastrophic cancellation even when the
// sum crosses zero). This is what makes the Burning Ship delta correct through
// the sign flips that happen all along the axes, where the ship's detail lives.
fn diffabs(c: f32, d: f32) -> f32 {
let cd = c + d;
if c >= 0.0 {
if (c >= 0.0) {
return select(-(2.0 * c + d), d, cd >= 0.0);
}
return select(-d, 2.0 * c + d, cd > 0.0);
}
// Perturbation delta for z -> z^p: (Z+e)^p - Z^p = e * sum_{k=0}^{p-1} (Z+e)^k Z^{p-1-k}.
// The large z^p term is never formed (that would cancel catastrophically), and
// the sum is evaluated Horner-style (s <- s*(Z+e) + Z^j) so it needs neither a
// table of powers (a dynamically indexed local array spills to slow memory on
// most GPUs) nor binomial coefficients. Forming Z+e rounds e away when it's
// tiny, but that only perturbs `s` by a relative f32 epsilon, and the result
// is `e * s`, so the delta keeps full relative precision.
// Binomial coefficient C(n, k) as f32 (exact for the small powers we use).
fn binom(n: u32, k: u32) -> f32 {
var num = 1.0;
var den = 1.0;
for (var i: u32 = 0u; i < k; i = i + 1u) {
num = num * f32(n - i);
den = den * f32(i + 1u);
}
return num / den;
}
// Perturbation delta for z -> z^p: sum_{k=1}^{p} C(p,k) Z^{p-k} e^k. Expanded so
// the large z^p term is never formed (that would cancel catastrophically).
fn multibrot_delta(z: vec2<f32>, e: vec2<f32>, p: u32) -> vec2<f32> {
let y = z + e;
var s = vec2<f32>(1.0, 0.0);
var zj = vec2<f32>(1.0, 0.0);
for (var j: u32 = 1u; j < p; j = j + 1u) {
zj = cmul(zj, z); // Z^j
s = cmul(s, y) + zj;
var zp: array<vec2<f32>, 9>; // Z^0 .. Z^8
zp[0] = vec2<f32>(1.0, 0.0);
for (var j: u32 = 1u; j <= p; j = j + 1u) {
zp[j] = cmul(zp[j - 1u], z);
}
return cmul(e, s);
}
// Maximum number of terms in `complex_multibrot_delta`'s series (matches the
// `cm_coef` uniform array: 8 vec4s = 16 complex coefficients). Truncation, not
// exactness: unlike `multibrot_delta` (a finite sum for an integer power), a
// complex power has no finite expansion, so this converges rather than
// terminates. Fine as long as perturbation's usual invariant (|e| << |z|,
// kept true by rebasing) holds, since each extra term is O(w^k) smaller.
const COMPLEX_MULTIBROT_TERMS: u32 = 16u;
// Complex binomial coefficient C(p, k), k in 1..=16, precomputed on the CPU
// (they depend only on p; see `complex_binomials` in app.rs).
fn cm_coef(k: u32) -> vec2<f32> {
let v = u.cm_coef[(k - 1u) / 2u];
return select(v.xy, v.zw, (k & 1u) == 0u);
}
// Perturbation delta for z -> z^p with a complex p: (Z+e)^p - Z^p.
//
// When |e| << |Z| (the common case: it's the whole reason perturbation
// works), forming Z+e directly would round e away in f32, so instead expand
// = Z^p * ((1+w)^p - 1), w = e/Z, as a Taylor series in w: (1+w)^p - 1 =
// sum_{k=1}^N C(p,k) w^k. The series stops as soon as the next w^k is
// negligible against the running sum (below f32 precision) — at deep zoom w
// is tiny, so that's typically after 2-3 terms instead of all 16.
//
// Right after a rebase (or near a reference point close to zero, where w is
// singular), e is *not* small relative to Z — that's normal perturbation
// dynamics, not a deep-zoom edge case — and the series above would diverge.
// But forming Z+e directly is numerically safe exactly there (e isn't many
// orders of magnitude smaller than Z), so fall back to a plain subtraction.
fn complex_multibrot_delta(z: vec2<f32>, e: vec2<f32>, p: vec2<f32>) -> vec2<f32> {
// |w|^2 = |e|^2 / |Z|^2; inf or nan (Z ~ 0, or both ~ 0) correctly fails
// the `< 0.25` test below and falls through to the direct branch.
let w2 = dot(e, e) / dot(z, z);
if w2 < 0.25 {
let w = cdiv(e, z);
var wk = w; // w^1
var acc = vec2<f32>(0.0, 0.0);
for (var k: u32 = 1u; k <= COMPLEX_MULTIBROT_TERMS; k = k + 1u) {
acc = acc + cmul(cm_coef(k), wk);
wk = cmul(wk, w);
if dot(wk, wk) < 1e-18 * dot(acc, acc) {
break;
var ek = vec2<f32>(1.0, 0.0); // e^0
for (var k: u32 = 1u; k <= p; k = k + 1u) {
ek = cmul(ek, e); // e^k
acc = acc + binom(p, k) * cmul(zp[p - k], ek);
}
}
return cmul(cpow(z, p), acc);
}
return cpow(z + e, p) - cpow(z, p);
return acc;
}
// One perturbation step of the current fractal's delta: e -> f(Z+e) - f(Z),
// where `z` is the reference orbit value X_m. `step_add` (dc) is added by the
// caller. Must match `FractalKind` on the CPU side.
fn advance_delta(z: vec2<f32>, e: vec2<f32>) -> vec2<f32> {
if KIND == KIND_BURNING_SHIP {
if (u.kind == KIND_BURNING_SHIP) {
// (|x| + i|y|)^2 has real part x^2 - y^2 (an ordinary square delta) and
// imaginary part 2|x y|. The imaginary delta is 2(|x y| - |X Y|); diffabs
// computes it exactly, even where the product x y changes sign — which the
@@ -150,37 +136,14 @@ fn advance_delta(z: vec2<f32>, e: vec2<f32>) -> vec2<f32> {
let base = 2.0 * cmul(z, e) + cmul(e, e);
let dp = z.x * e.y + z.y * e.x + e.x * e.y;
return vec2<f32>(base.x, 2.0 * diffabs(z.x * z.y, dp));
} else if KIND == KIND_TRICORN {
} else if (u.kind == KIND_TRICORN) {
let cz = conj(z);
let ce = conj(e);
return 2.0 * cmul(cz, ce) + cmul(ce, ce);
} else if KIND == KIND_MULTIBROT {
} else if (u.kind == KIND_MULTIBROT) {
return multibrot_delta(z, e, clamp(u.power, 2u, 8u));
} else if KIND == KIND_CELTIC {
// z^2 delta split: sq.x = delta of Re(z^2), sq.y = delta of Im(z^2).
// Celtic abs the real output, so |Re(z^2)| delta = diffabs(Re(Z^2), sq.x).
let sq = 2.0 * cmul(z, e) + cmul(e, e);
return vec2<f32>(diffabs(z.x * z.x - z.y * z.y, sq.x), sq.y);
} else if KIND == KIND_BUFFALO {
// Abs both outputs: real |Re(z^2)|, imag -|Im(z^2)| (Im(Z^2) = 2 X Y).
let sq = 2.0 * cmul(z, e) + cmul(e, e);
return vec2<f32>(diffabs(z.x * z.x - z.y * z.y, sq.x),
-diffabs(2.0 * z.x * z.y, sq.y));
} else if KIND == KIND_PERPENDICULAR {
// real x^2 - y^2 (ordinary square delta), imag -2 x |y|.
// d(-2 x |y|) = -2[ X·(|Y+ey|-|Y|) + ex·|Y+ey| ]; diffabs gives |Y+ey|-|Y|.
let sq = 2.0 * cmul(z, e) + cmul(e, e);
let da = diffabs(z.y, e.y); // |Y + ey| - |Y|
let abs_yf = abs(z.y) + da; // |Y + ey|
return vec2<f32>(sq.x, -2.0 * (z.x * da + e.x * abs_yf));
} else if KIND == KIND_LAMBDA {
// Lambda map: z^{n+1} = λ·z·(1-z). Delta: e = λ·e·(1-2z-e).
let one_minus_2z_minus_e = vec2<f32>(1.0 - 2.0 * z.x - e.x, -2.0 * z.y - e.y);
return cmul(u.lambda_l, cmul(e, one_minus_2z_minus_e));
} else if KIND == KIND_COMPLEX_MULTIBROT {
return complex_multibrot_delta(z, e, u.complex_power);
}
return 2.0 * cmul(z, e) + cmul(e, e); // Mandelbrot (and Phoenix square part)
return 2.0 * cmul(z, e) + cmul(e, e); // Mandelbrot
}
// Derivative f'(Z) of the iteration map at the full value Z, used to propagate
@@ -189,427 +152,217 @@ fn advance_delta(z: vec2<f32>, e: vec2<f32>) -> vec2<f32> {
// Burning Ship / Tricorn we use |f'| ~ |2Z|, which keeps the DE magnitude close
// enough to de-speckle filaments.
fn fprime(z: vec2<f32>) -> vec2<f32> {
if KIND == KIND_MULTIBROT {
if (u.kind == KIND_MULTIBROT) {
let p = clamp(u.power, 2u, 8u);
var zk = z; // Z^1
for (var k: u32 = 2u; k < p; k = k + 1u) {
var zk = vec2<f32>(1.0, 0.0); // Z^0
for (var k: u32 = 1u; k < p; k = k + 1u) {
zk = cmul(zk, z); // -> Z^{p-1}
}
return f32(p) * zk;
} else if KIND == KIND_LAMBDA {
// Lambda: f'(z) = λ·(1-2z).
return cmul(u.lambda_l, vec2<f32>(1.0 - 2.0 * z.x, -2.0 * z.y));
} else if KIND == KIND_COMPLEX_MULTIBROT {
// f'(z) = p * z^(p-1).
return cmul(u.complex_power, cpow(z, u.complex_power - vec2<f32>(1.0, 0.0)));
}
return 2.0 * z;
}
// Periodicity (interior) detection, Brent-style: the full orbit value is
// saved at iterations PERIOD_FIRST_CHECK, 2x that, 4x ..., and every later
// iterate is compared against the last saved one. Returning within
// PERIOD_EPS2 (relative, squared) means the orbit has closed a cycle.
//
// A close return alone isn't trusted. A pixel just *outside* the set (at a
// minibrot's edge, or a cusp) can shadow a cycle for thousands of iterations
// before escaping. So three safeguards apply, tuned against an f64 simulation
// of this exact algorithm and f64 ground truth on cusp, bulb-contact,
// minibrot-edge and deep seahorse views:
// * Multiplier: |product of f'(z)|^2 over the steps since the save must be
// < PERIOD_MAX_MULT2, so the cycle it closed is clearly attracting. Plain
// "< 1" let near-parabolic exterior points (|multiplier| ~ 1) through at
// cusps; the margin fixes that.
// * Confirmation: the contracting return must happen in
// PERIOD_CONFIRMATIONS consecutive windows, each twice as long as the
// last. Exterior orbits passing near the critical point can look strongly
// contracting for one window (seen: flagged at iteration 245, escaped at
// 2275). A second, longer window rules that out.
// * Tolerance: PERIOD_EPS2 is relative and near f32 precision.
// Every kind here except Phoenix is
// (piecewise) conformal, so |f'| from `fprime` is the exact local scale
// factor, including the abs-folding kinds, whose folds are isometries.
// Phoenix's two-term map would need a 2x2 Jacobian, so it's excluded. So is
// Complex Multibrot without DE, where `fprime` would add a second `cpow`
// (log/atan2/exp) per step for a check that rarely fires on its views.
const PERIOD_FIRST_CHECK: u32 = 16u;
const PERIOD_EPS2: f32 = 1e-12;
const PERIOD_MAX_MULT2: f32 = 0.25;
const PERIOD_CONFIRMATIONS: u32 = 2u;
// Whether `iterate_sample` runs periodicity detection for this kind (folds to
// a constant per pipeline).
fn periodic_enabled() -> bool {
if KIND == KIND_PHOENIX {
return false;
// Smooth cyclic palettes (Inigo Quilez cosine palettes), selected by id.
fn palette(id: u32, t: f32) -> vec3<f32> {
if (id == 4u) {
return vec3<f32>(t, t, t); // grayscale
}
if KIND == KIND_COMPLEX_MULTIBROT && !DE {
return false;
let a = vec3<f32>(0.5, 0.5, 0.5);
let b = vec3<f32>(0.5, 0.5, 0.5);
var c = vec3<f32>(1.0, 1.0, 1.0);
var d = vec3<f32>(0.00, 0.10, 0.20); // 0: amber / blue
if (id == 1u) {
d = vec3<f32>(0.00, 0.33, 0.67); // rainbow
} else if (id == 2u) {
d = vec3<f32>(0.30, 0.20, 0.20); // warm ember
} else if (id == 3u) {
c = vec3<f32>(1.0, 1.0, 0.5);
d = vec3<f32>(0.80, 0.90, 0.30); // lime / magenta
}
return true;
return a + b * cos(6.28318530718 * (c * t + d));
}
// Escape data for one sample: `ci` is the (color-independent) palette parameter,
// `de` the distance-estimate darkening factor in [0,1], `escaped` false for the
// interior of the set. Splitting iteration from coloring lets a colour change be
// remapped cheaply (see the colourise pass) without re-iterating.
struct Sample {
ci: f32,
de: f32,
escaped: bool,
// Result of a BLA lookup at an orbit index.
struct Hop {
found: bool,
a: vec2<f32>,
b: vec2<f32>,
l: u32,
};
// Perturbation iterate a single sample. `offset` is the per-pixel offset in
// complex units. For Mandelbrot it is the c-plane offset added every step (delta
// starts at 0); for Julia it is the z-plane offset that seeds the initial delta
// (c is fixed, so nothing is added per step).
fn iterate_sample(offset: vec2<f32>, px: f32) -> Sample {
// Loop invariants, read once instead of on every iteration.
let max_iter = u.max_iter;
let bailout_sq = u.bailout_sq;
let ref_len = u.ref_len;
// Longest valid BLA skip starting at orbit index `n` for a delta of squared
// magnitude `emag2`. Walks levels low->high, recomputing each level's flat-array
// start and count from `ref_len` (count[0] = ref_len-1, halving each level).
// Radii shrink with level and alignment is monotonic, so the valid levels form a
// prefix: we keep the last valid one and stop at the first that fails.
fn bla_find(n: u32, emag2: f32) -> Hop {
var res: Hop;
res.found = false;
var start: u32 = 0u;
var count: u32 = u.ref_len - 1u;
var lv: u32 = 0u;
loop {
if (count == 0u) { break; }
let step = 1u << lv;
if ((n & (step - 1u)) != 0u) { break; } // n not aligned to this level
let idx = n >> lv;
if (idx >= count) { break; } // run would exceed the orbit
let b = bla_table[start + idx];
if (emag2 >= b.r * b.r) { break; } // delta too large: not valid
res.found = true;
res.a = b.a;
res.b = b.b;
res.l = b.l;
start = start + count;
count = count / 2u;
lv = lv + 1u;
}
return res;
}
// Perturbation iterate + color a single sample. `offset` is the per-pixel
// offset in complex units. For Mandelbrot it is the c-plane offset added every
// step (delta starts at 0); for Julia it is the z-plane offset that seeds the
// initial delta (c is fixed, so nothing is added per step). Interior pixels
// return black.
fn shade(offset: vec2<f32>, px: f32) -> vec3<f32> {
let z0 = ref_orbit[0]; // reference start (0 for Mandelbrot, center for Julia)
// Main cardioid / period-2 bulb bypass: those points never escape, so skip
// iterating them (they'd otherwise all burn the full max_iter). `offset` is
// relative to the reference center; the absolute c is recovered from the
// orbit itself, since X_1 = X_0^2 + C_ref = C_ref. That's only f32-accurate,
// so skip the test once a pixel is smaller than that error (deep zoom),
// where it could misclassify pixels right at the boundary.
if KIND == KIND_MANDELBROT && !IS_JULIA && ref_len > 1u && px > 1e-6 {
let c = ref_orbit[1] + offset;
let xq = c.x - 0.25;
let q = xq * xq + c.y * c.y;
let in_cardioid = q * (q + xq) <= 0.25 * c.y * c.y;
let xb = c.x + 1.0;
let in_bulb = xb * xb + c.y * c.y <= 0.0625;
if in_cardioid || in_bulb {
return Sample(0.0, 1.0, false); // interior of the set
}
}
// Set plane: delta starts at 0 and gains dc every step. Julia: the offset
// seeds the delta and nothing is added per step.
var step_add = offset;
var e = vec2<f32>(0.0, 0.0);
// Orbit derivative for distance estimation, pre-multiplied by the pixel
// size `px`. For the set plane it is px·d/dc (starts at 0, gains +px each
// step); for Julia it is px·d/dz0 (starts at px). The raw derivative grows
// like 1/px, so unscaled its square overflows f32 at deep zoom (~1e-12),
// which zeroed DE along iteration bands; scaled, it stays ~|z|ln|z| / DE
// in pixels at any depth.
var dzs = vec2<f32>(0.0, 0.0);
if IS_JULIA {
// Orbit derivative for distance estimation. For the set plane it is d/dc
// (starts at 0, gains +1 each step); for Julia it is d/dz0 (starts at 1).
var dz = vec2<f32>(0.0, 0.0);
var dz_seed = vec2<f32>(1.0, 0.0);
if (u.is_julia != 0u) {
step_add = vec2<f32>(0.0, 0.0);
e = offset;
dzs = vec2<f32>(px, 0.0);
dz = vec2<f32>(1.0, 0.0);
dz_seed = vec2<f32>(0.0, 0.0);
}
// Previous-iterate state for the Phoenix two-term recurrence (delta of
// y_{n-1}, and its scaled derivative for DE). Both start at 0 (y_{-1} = 0).
var e_prev = vec2<f32>(0.0, 0.0);
var dzs_prev = vec2<f32>(0.0, 0.0);
var m: u32 = 0u; // reference index; invariant: y_n = xm + e, xm = X[m]
var m: u32 = 0u; // reference index; invariant: y_n = X[m] + e
var n: u32 = 0u; // total iteration count
var xm = z0; // X[m], carried so each step loads the orbit once
var z = xm + e; // full value y_n, kept for coloring
var z2 = dot(z, z);
var z = vec2<f32>(0.0, 0.0); // full value y_n, kept for coloring
var escaped = false;
// Periodicity detection (see PERIOD_FIRST_CHECK): last saved orbit value,
// |f'|^2 product of the steps since it was saved, next save iteration.
let periodic = periodic_enabled();
// Plus whether this window already had a contracting return, and how many
// consecutive windows have.
var z_saved = z;
var mult2 = 1.0;
var check_at = PERIOD_FIRST_CHECK;
var period_hit = false;
var period_streak = 0u;
loop {
if z2 > bailout_sq {
let xm = ref_orbit[m];
z = xm + e;
let z2 = dot(z, z);
if (z2 > u.bailout_sq) {
escaped = true;
break;
}
if n >= max_iter {
if (n >= u.max_iter) {
break; // interior
}
// Advance one step, or skip a whole run via BLA when the delta is small
// enough (square map only; other kinds keep `use_bla == 0`). The orbit
// derivative for DE follows the same linear map: over a run it advances
// by the run's own (A, B) coefficients, matching the per-step recurrence.
var did_skip = false;
if (u.use_bla != 0u) {
let hop = bla_find(m, dot(e, e));
if (hop.found) {
if (u.de_coloring != 0u) {
dz = cmul(hop.a, dz) + cmul(hop.b, dz_seed);
}
e = cmul(hop.a, e) + cmul(hop.b, step_add);
m = m + hop.l;
n = n + hop.l;
did_skip = true;
}
}
if (!did_skip) {
// Propagate the derivative of the full orbit (unaffected by rebasing,
// which only re-expresses the same value). Only when DE is enabled.
// Phoenix's two-term map adds p·dz_{n-1} and carries the previous dz.
// f'(z) of this step, shared by DE and the periodicity multiplier.
var fp = vec2<f32>(0.0, 0.0);
if DE || periodic {
fp = fprime(z);
if (u.de_coloring != 0u) {
dz = cmul(fprime(z), dz) + dz_seed;
}
if periodic {
mult2 = mult2 * dot(fp, fp);
}
if DE {
var dzs_new = cmul(fp, dzs);
if !IS_JULIA {
dzs_new.x = dzs_new.x + px;
}
if KIND == KIND_PHOENIX {
dzs_new = dzs_new + cmul(u.phoenix_p, dzs_prev);
dzs_prev = dzs;
}
dzs = dzs_new;
}
// Advance the delta by this fractal's formula (+ dc for the set plane).
// Phoenix additionally adds p·e_{n-1} and carries the previous delta.
let e_old = e;
let z_old = z;
e = advance_delta(xm, e) + step_add;
if KIND == KIND_PHOENIX {
e = e + cmul(u.phoenix_p, e_prev);
e_prev = e_old;
}
m = m + 1u;
n = n + 1u;
}
// Keep the reference index valid and the delta small.
if m >= ref_len {
if (m >= u.ref_len) {
// Reference exhausted: any pixel that followed it this far has
// effectively escaped (interior pixels rebase before reaching here).
z = xm + e;
z = ref_orbit[u.ref_len - 1u] + e;
escaped = true;
break;
}
xm = ref_orbit[m];
z = xm + e;
z2 = dot(z, z);
if z2 < dot(e, e) {
let y = ref_orbit[m] + e;
if (dot(y, y) < dot(e, e)) {
// Rebase to index 0: carry the full value as the new delta. Valid
// because y_n = X[0] + (y_n - X[0]); for Mandelbrot X[0]=0. The
// full value `z` (and `z2`) is unchanged by the re-expression.
// Phoenix: after rebasing the implied previous reference is Y[-1]=0,
// so the previous delta becomes the full previous value y_{n-1}.
if KIND == KIND_PHOENIX {
e_prev = z_old;
}
e = z - z0;
xm = z0;
// because y_n = X[0] + (y_n - X[0]); for Mandelbrot X[0]=0.
e = y - z0;
m = 0u;
}
if periodic {
// Closed an attracting cycle in enough consecutive windows:
// interior (see PERIOD_FIRST_CHECK).
let d = z - z_saved;
if !period_hit && mult2 < PERIOD_MAX_MULT2 && dot(d, d) <= PERIOD_EPS2 * z2 {
period_hit = true;
period_streak = period_streak + 1u;
if period_streak >= PERIOD_CONFIRMATIONS {
break;
}
}
if n == check_at {
if !period_hit {
period_streak = 0u;
}
period_hit = false;
z_saved = z;
mult2 = 1.0;
check_at = check_at * 2u;
}
}
}
if !escaped {
return Sample(0.0, 1.0, false); // interior of the set
if (!escaped) {
return vec3<f32>(0.0, 0.0, 0.0); // interior of the set
}
// Both escape formulas below (smooth count, DE) assume |f(z)| ~ |z|^2 near
// escape. Lambda's λz(1-z) + c ~ -λz^2 adds a factor |λ| per step, which
// made both jump at every band boundary (contour lines in shadow/3D).
// w = -λz conjugates it to an exact w^2 + C, so measure |w| and |dw|.
var dzs_esc = dzs;
if KIND == KIND_LAMBDA {
z = cmul(u.lambda_l, z);
dzs_esc = cmul(u.lambda_l, dzs);
}
z2 = dot(z, z);
let z2 = dot(z, z);
// Continuous (smooth) iteration count.
let log_zn = 0.5 * log(max(z2, 1.0));
let nu = log2(log_zn * INV_LN2);
let nu = log2(log_zn / log(2.0));
let smooth_i = f32(n) + 1.0 - nu;
// sqrt compresses the huge iteration counts of deep zooms so the palette
// varies smoothly instead of aliasing into speckle.
let ci = sqrt(max(smooth_i, 0.0));
let t = fract(ci * u.color_scale + u.color_offset);
var col = palette(u.palette_id, t);
var de = 1.0;
if DE {
// Exterior distance estimate |z|·ln|z| / |dz|, already in pixels since
// `dzs` = px·dz. We darken toward the boundary (< ~1 px away) so
// filaments stay crisp instead of aliasing into speckle. If |dzs|
// overflowed (far sub-pixel from the set), de -> 0 and the boundary
// simply reads as dark, which is the correct limit.
if (u.de_coloring != 0u) {
// Exterior distance estimate (complex-plane units): |z|·ln|z| / |dz|.
// Divided by the pixel footprint it becomes a distance in pixels; we
// darken toward the boundary (< ~1 px away) so filaments stay crisp
// instead of aliasing into speckle. If |dz| overflowed, de -> 0 and the
// boundary simply reads as dark, which is the correct limit.
let zmag = sqrt(max(z2, 1.0));
let dzmag = sqrt(max(dot(dzs_esc, dzs_esc), 1e-30));
let max_de = select(1.0, 1000.0, u.shadow != 0u);
de = clamp(zmag * log(zmag) / dzmag, 0.0, max_de);
let dzmag = sqrt(max(dot(dz, dz), 1e-20));
let de = zmag * log(zmag) / dzmag;
let de_px = de / max(px, 1e-30);
col = col * clamp(de_px, 0.0, 1.0);
}
return Sample(ci, de, true);
return col;
}
// 1 / ln(2), for the smooth iteration count's log2(ln|z| / ln 2).
const INV_LN2: f32 = 1.4426950408889634;
// Map a sample's escape data through the palette (+ DE darkening). This is the
// only color-dependent step, so it can be redone without re-iterating. Interior
// samples are black.
fn color_sample(s: Sample) -> vec3<f32> {
if !s.escaped {
return vec3<f32>(0.0, 0.0, 0.0);
}
return classic_color(s.ci, s.de);
}
// Supersampled escape data at one point: average (ci, DE factor) over an
// `aa`×`aa` grid's escaped sub-samples, plus the fraction that landed in the
// interior. Shared by `fs_data` (1 sample), `fs_refine` (the AA grid, only on
// pixels that need it) and `fs_color`'s shadow branch (used both at the pixel
// and at its two neighbours, to build a DE height field without a texture
// round-trip).
fn aggregate_sample(base: vec2<f32>, dx: vec2<f32>, dy: vec2<f32>, px: f32, aa: u32) -> vec3<f32> {
let inv = 1.0 / f32(aa);
var ci_sum = 0.0;
var de_sum = 0.0;
var escaped_n = 0u;
for (var sy: u32 = 0u; sy < aa; sy = sy + 1u) {
for (var sx: u32 = 0u; sx < aa; sx = sx + 1u) {
let jx = (f32(sx) + 0.5) * inv - 0.5;
let jy = (f32(sy) + 0.5) * inv - 0.5;
let s = iterate_sample(base + jx * dx + jy * dy, px);
if s.escaped {
ci_sum = ci_sum + s.ci;
de_sum = de_sum + s.de;
escaped_n = escaped_n + 1u;
}
}
}
let total = f32(aa * aa);
let ci_avg = select(0.0, ci_sum / f32(escaped_n), escaped_n > 0u);
let de_avg = select(1.0, de_sum / f32(escaped_n), escaped_n > 0u);
let interior_frac = 1.0 - f32(escaped_n) / total;
return vec3<f32>(ci_avg, de_avg, interior_frac);
}
// Pixel footprint in complex units, |(|dx| + |dy|)|. Not `length()` directly:
// that squares its argument, and below ~1e-19 per pixel (half-height ~1e-16,
// sooner for the 2x-resolution 3D texture) the square drops under f32's
// smallest normal and flushes to 0, making px = 0 and DE meaningless.
// Normalizing by the largest component first keeps the square near 1.
fn pixel_size(dx: vec2<f32>, dy: vec2<f32>) -> f32 {
let a = abs(dx) + abs(dy);
let m = max(a.x, a.y);
if m == 0.0 {
return 0.0;
}
return m * length(a / m);
}
// Iteration pass: write per-pixel escape data (color-independent) so a colour
// change is remapped by the cheap colourise pass without re-iterating.
// R = ci (palette parameter), G = DE factor, B = interior fraction (for AA).
// Always one sample per pixel: anti-aliasing is added afterwards, only where
// it matters, by `fs_refine`.
@fragment
fn fs_data(in: VsOut) -> @location(0) vec4<f32> {
fn fs_main(in: VsOut) -> @location(0) vec4<f32> {
let base = in.centered * u.span + u.dc_offset;
// Screen-space complex-units-per-pixel. Derivatives must be evaluated in
// uniform control flow, so take them here; used to place sub-pixel AA
// samples and to convert the distance estimate into pixels.
let dx = dpdx(base);
let dy = dpdy(base);
let px = pixel_size(dx, dy);
let px = length(abs(dx) + abs(dy)); // ~ complex units per pixel (footprint)
return vec4<f32>(aggregate_sample(base, dx, dy, px, 1u), 1.0);
}
// Adaptive-AA thresholds for `fs_refine`: a pixel is supersampled only if a
// 4-neighbour's 1-spp sample differs from its own by more than this. `ci`
// steps are palette-phase steps of `ci * color_scale` (color_scale <= 1 in the
// UI), so 0.02 keeps anything visibly banded; DE is compared relative to its
// own magnitude (it's in pixels, up to 1000 for shadow/3D height fields).
const AA_CI_EPS: f32 = 0.02;
const AA_DE_EPS: f32 = 0.1;
fn aa_differs(c: vec4<f32>, n: vec4<f32>) -> bool {
if c.b != n.b {
return true; // interior / exterior boundary
}
if c.b != 0.0 {
return false; // both interior: uniformly black
}
return abs(n.r - c.r) > AA_CI_EPS || abs(n.g - c.g) > AA_DE_EPS * max(c.g, 0.1);
}
// Adaptive anti-aliasing pass (only run when AA is on): reads `fs_data`'s
// 1-spp texture and re-iterates the full AA grid only for pixels whose
// neighbourhood isn't smooth (set boundary, filaments, palette discontinuities).
// Everywhere else the centre sample already equals the grid average to within
// the thresholds above, so it's copied — which skips the AA cost entirely for
// the interior (the most expensive pixels, each burning max_iter) and for the
// smooth exterior.
@fragment
fn fs_refine(in: VsOut) -> @location(0) vec4<f32> {
// Derivatives first, while control flow is still uniform.
let base = in.centered * u.span + u.dc_offset;
let dx = dpdx(base);
let dy = dpdy(base);
let px = pixel_size(dx, dy);
let p = vec2<i32>(in.pos.xy);
let hi = vec2<i32>(textureDimensions(coarse_tex)) - vec2<i32>(1, 1);
let c = textureLoad(coarse_tex, p, 0);
let l = textureLoad(coarse_tex, max(p - vec2<i32>(1, 0), vec2<i32>(0, 0)), 0);
let r = textureLoad(coarse_tex, min(p + vec2<i32>(1, 0), hi), 0);
let t = textureLoad(coarse_tex, max(p - vec2<i32>(0, 1), vec2<i32>(0, 0)), 0);
let b = textureLoad(coarse_tex, min(p + vec2<i32>(0, 1), hi), 0);
if aa_differs(c, l) || aa_differs(c, r) || aa_differs(c, t) || aa_differs(c, b) {
return vec4<f32>(aggregate_sample(base, dx, dy, px, max(u.aa_level, 1u)), 1.0);
}
return c;
}
// Combined iterate + colour in a single pass, for PNG export (which never needs
// incremental recolouring). The interactive path uses fs_data (+ fs_refine) +
// the colourise pass so colour changes skip iteration. Export always runs the
// full AA grid on every pixel, for maximum quality.
@fragment
fn fs_color(in: VsOut) -> @location(0) vec4<f32> {
let base = in.centered * u.span + u.dc_offset;
let dx = dpdx(base);
let dy = dpdy(base);
let px = pixel_size(dx, dy);
let aa = max(u.aa_level, 1u);
if u.shadow != 0u {
// No data texture to sample neighbours from (this pass never runs
// one), so build the same DE height field colorize.wgsl reads from
// the texture by aggregating live, at the pixel and its two
// neighbours a `dx`/`dy` step away.
let here = aggregate_sample(base, dx, dy, px, aa);
if here.z != 0.0 {
return vec4<f32>(shadow_interior_color(), 1.0);
}
let right = aggregate_sample(base + dx, dx, dy, px, aa);
let down = aggregate_sample(base + dy, dx, dy, px, aa);
let normal = normal_from_heights(here.y, right.y, down.y);
return vec4<f32>(shadow_color(normal, here.x), 1.0);
if (aa <= 1u) {
return vec4<f32>(shade(base, px), 1.0);
}
let inv = 1.0 / f32(aa);
var acc = vec3<f32>(0.0, 0.0, 0.0);
let inv = 1.0 / f32(aa);
for (var sy: u32 = 0u; sy < aa; sy = sy + 1u) {
for (var sx: u32 = 0u; sx < aa; sx = sx + 1u) {
// Sample centers evenly spread across the pixel, jitter in (-0.5, 0.5).
let jx = (f32(sx) + 0.5) * inv - 0.5;
let jy = (f32(sy) + 0.5) * inv - 0.5;
acc = acc + color_sample(iterate_sample(base + jx * dx + jy * dy, px));
acc = acc + shade(base + jx * dx + jy * dy, px);
}
}
return vec4<f32>(acc / f32(aa * aa), 1.0);
+2 -148
View File
@@ -58,11 +58,6 @@ impl ViewState {
DEFAULT_HALF_HEIGHT / self.half_height
}
/// Current zoom level.
pub fn zoom(&self) -> f64 {
self.half_height
}
/// Bits of precision the center currently needs for this zoom level.
pub fn precision_bits(&self) -> usize {
precision_for(self.half_height)
@@ -87,7 +82,7 @@ impl ViewState {
let bits = self.precision_bits();
// Grab-and-drag: moving the mouse right shows content to the left.
self.center_re = &self.center_re - &big_from_f64(dx * cpp, bits);
self.center_im = &self.center_im - &big_from_f64(dy * cpp, bits);
self.center_im = &self.center_im - &big_from_f64(dy * cpp, bits); // y-down -> imag-up
}
/// Zoom by `factor` (<1 zooms in) keeping the complex point currently under
@@ -102,7 +97,7 @@ impl ViewState {
// off * cpp * (1 - factor). (Derivation: new_c = fixed + (c-fixed)*f.)
let k = cpp * (1.0 - factor);
self.center_re = &self.center_re + &big_from_f64(off_x * k, bits);
self.center_im = &self.center_im + &big_from_f64(off_y * k, bits);
self.center_im = &self.center_im + &big_from_f64(off_y * k, bits); // y flip
self.half_height *= factor;
}
@@ -125,86 +120,6 @@ pub fn big_from_decimal_str(s: &str, bits: usize) -> Option<Big> {
Some(dec.with_base_and_precision::<2>(bits.max(53)).value())
}
/// Parse a "re,im,half_height[,iterations]" spec (re/im decimal, parsed at
/// full precision) into a view and an optional iteration count. Shared by
/// `FractalApp::apply_view_spec` (the `--view` CLI flag) and headless
/// animation's `--to-view`.
pub fn parse_view_spec(spec: &str) -> Option<(ViewState, Option<u32>)> {
let parts: Vec<&str> = spec.split(',').collect();
if parts.len() < 3 {
return None;
}
let half_height = parts[2].trim().parse::<f64>().ok()?;
if !(half_height > 0.0 && half_height.is_finite()) {
return None;
}
let bits = precision_for(half_height);
let re = big_from_decimal_str(parts[0], bits)?;
let im = big_from_decimal_str(parts[1], bits)?;
let iterations = parts.get(3).and_then(|s| s.trim().parse::<u32>().ok());
Some((ViewState::with_center(re, im, half_height), iterations))
}
/// Parse a half_height spec. Shared by
/// `FractalApp::apply_half_height_spec` (the `--zoom` CLI flag) and headless
/// animation's `--to-zoom`.
pub fn parse_half_height_spec(spec: &str) -> Option<f64> {
let half_height = spec.trim().parse::<f64>().ok()?;
if !(half_height > 0.0 && half_height.is_finite()) {
return None;
}
Some(half_height)
}
/// Parse a "re,im" spec (re/im decimal, parsed at
/// full precision) into a view. Shared by
/// `FractalApp::apply_re_im_spec` (the `--position` CLI flag) and headless
/// animation's `--to-position`.
pub fn parse_re_im_spec(spec: &str, bits: usize) -> Option<(Big, Big)> {
let parts: Vec<&str> = spec.split(',').collect();
if parts.len() != 2 {
return None;
}
let re = big_from_decimal_str(parts[0], bits)?;
let im = big_from_decimal_str(parts[1], bits)?;
Some((re, im))
}
/// Interpolate between two views for an animation frame, `t` in `[0, 1]`.
/// The half-height interpolates geometrically (log-linear), since zoom depth
/// spans many decades and a linear sweep would crawl at the start and blow
/// past the target at the end. The center has to shrink its offset from the
/// target at that *same* geometric rate: blending it linearly in `t` instead
/// barely moves it while the view is still huge (early frames), so the
/// target stays effectively off-screen — offset/half_height ratio blows up —
/// for nearly the whole animation, and only lands on `to`'s center in the
/// literal last frame where `t == 1` forces an exact match. `g(t)` below
/// tracks the same `q^t` decay used for `half_height` (keeping the
/// offset/half_height ratio roughly constant, i.e. the target's on-screen
/// position steady) but is shifted so it lands on exactly 1 at `t = 0` and
/// exactly 0 at `t = 1`.
pub fn interpolate_view(from: &ViewState, to: &ViewState, t: f64) -> ViewState {
let q = to.half_height / from.half_height;
let half_height = from.half_height * q.powf(t);
let bits = precision_for(half_height);
let g = if (q - 1.0).abs() < 1e-12 {
1.0 - t
} else {
(q.powf(t) - q) / (1.0 - q)
};
let g_big = big_from_f64(g, bits);
let re0 = from.center_re.clone().with_precision(bits).value();
let im0 = from.center_im.clone().with_precision(bits).value();
let re1 = to.center_re.clone().with_precision(bits).value();
let im1 = to.center_im.clone().with_precision(bits).value();
let center_re = &re1 + &(&(&re0 - &re1) * &g_big);
let center_im = &im1 + &(&(&im0 - &im1) * &g_big);
ViewState::with_center(center_re, center_im, half_height)
}
pub fn interpolate_f64(from: f64, to: f64, t: f64) -> f64 {
from + (to - from) * t
}
/// Render a `Big` as a decimal string with `sig_digits` significant digits.
pub fn big_to_decimal_str(x: &Big, sig_digits: usize) -> String {
let dec = x
@@ -234,64 +149,3 @@ pub fn big_from_f64(x: f64, bits: usize) -> Big {
.with_precision(bits)
.value()
}
#[cfg(test)]
mod tests {
use super::*;
fn re_im_f64(v: &ViewState) -> (f64, f64) {
let re: f64 = v.center_re.to_decimal().value().to_f64().value();
let im: f64 = v.center_im.to_decimal().value().to_f64().value();
(re, im)
}
#[test]
fn interpolate_view_hits_exact_endpoints() {
let bits = precision_for(1.0);
let from = ViewState::with_center(big_from_f64(-0.5, bits), big_from_f64(0.0, bits), 1.5);
let to = ViewState::with_center(
big_from_f64(-0.7515, precision_for(1e-20)),
big_from_f64(0.1013, precision_for(1e-20)),
1e-20,
);
let start = interpolate_view(&from, &to, 0.0);
assert_eq!(re_im_f64(&start), re_im_f64(&from));
assert_eq!(start.half_height, from.half_height);
let end = interpolate_view(&from, &to, 1.0);
assert_eq!(re_im_f64(&end), re_im_f64(&to));
assert_eq!(end.half_height, to.half_height);
}
/// Regression test: a deep zoom's center used to be blended linearly in
/// `t` while `half_height` shrank geometrically, so partway through the
/// animation the offset from the target would already be far larger than
/// the (tiny, geometrically-shrunk) view — the target only snapped into
/// frame on the very last frame. The offset/half_height ratio should
/// instead stay roughly bounded throughout.
#[test]
fn interpolate_view_keeps_target_offset_bounded() {
let bits = precision_for(1.0);
let from = ViewState::with_center(big_from_f64(-0.5, bits), big_from_f64(0.0, bits), 1.5);
let to = ViewState::with_center(
big_from_f64(-0.7515, precision_for(1e-20)),
big_from_f64(0.1013, precision_for(1e-20)),
1e-20,
);
let (to_re, to_im) = re_im_f64(&to);
for i in 1..10 {
let t = i as f64 / 10.0;
let mid = interpolate_view(&from, &to, t);
let (re, im) = re_im_f64(&mid);
let offset = ((re - to_re).powi(2) + (im - to_im).powi(2)).sqrt();
let ratio = offset / mid.half_height;
assert!(
ratio < 10.0,
"t={t}: offset/half_height ratio {ratio} blew up (offset={offset}, half_height={})",
mid.half_height
);
}
}
}
+12 -13
View File
@@ -9,7 +9,7 @@
use std::sync::mpsc::{Receiver, Sender, TryRecvError, channel};
use std::thread;
use crate::fractal::{FractalKind, compute_reference, compute_set_reference};
use crate::fractal::{Bla, FractalKind, build_bla_table, compute_reference, compute_set_reference};
use crate::view::{Big, big_from_f64};
pub struct RefRequest {
@@ -22,12 +22,10 @@ pub struct RefRequest {
pub precision: usize,
pub kind: FractalKind,
pub power: u32,
/// Distortion constant for the Phoenix map (ignored by other kinds).
pub phoenix_p: (f64, f64),
/// Distortion constant for the Lambda map (ignored by other kinds).
pub lambda_l: (f64, f64),
/// Complex exponent for the Complex Multibrot kind (ignored by other kinds).
pub complex_power: (f64, f64),
/// Build the BLA iteration-skip table for this orbit (square map only).
pub build_bla: bool,
/// Upper bound on any pixel's `|dc|`, sizing the BLA merge radii.
pub dc_max: f64,
}
pub struct RefResult {
@@ -35,6 +33,7 @@ pub struct RefResult {
pub center_im: Big,
pub half_height: f64,
pub points: Vec<[f32; 2]>,
pub bla: Vec<Bla>,
}
pub struct RefWorker {
@@ -62,12 +61,18 @@ impl RefWorker {
}
let points = compute(&req);
let bla = if req.build_bla {
build_bla_table(&points, req.dc_max)
} else {
Vec::new()
};
if res_tx
.send(RefResult {
center_re: req.center_re,
center_im: req.center_im,
half_height: req.half_height,
points,
bla,
})
.is_err()
{
@@ -107,9 +112,6 @@ fn compute(req: &RefRequest) -> Vec<[f32; 2]> {
req.precision,
req.kind,
req.power,
req.phoenix_p,
req.lambda_l,
req.complex_power,
)
} else {
compute_set_reference(
@@ -119,9 +121,6 @@ fn compute(req: &RefRequest) -> Vec<[f32; 2]> {
req.precision,
req.kind,
req.power,
req.phoenix_p,
req.lambda_l,
req.complex_power,
)
}
}
+6 -125
View File
@@ -3,7 +3,7 @@
//! shader with the same `naga` version wgpu uses — catching shader errors
//! without needing a GPU or a display.
fn validate(name: &str, src: &str) -> (naga::Module, naga::valid::ModuleInfo) {
fn validate(name: &str, src: &str) {
let module = match naga::front::wgsl::parse_str(src) {
Ok(m) => m,
Err(e) => panic!("{name}: WGSL parse error:\n{}", e.emit_to_string(src)),
@@ -12,139 +12,20 @@ fn validate(name: &str, src: &str) -> (naga::Module, naga::valid::ModuleInfo) {
naga::valid::ValidationFlags::all(),
naga::valid::Capabilities::all(),
);
match validator.validate(&module) {
Ok(info) => (module, info),
Err(e) => panic!("{name}: WGSL validation error:\n{}", e.emit_to_string(src)),
if let Err(e) = validator.validate(&module) {
panic!("{name}: WGSL validation error:\n{}", e.emit_to_string(src));
}
}
/// Number of fractal kinds, i.e. the `const KIND_*` declarations in
/// common.wgsl (one per `FractalKind` variant, values 0..N).
fn kind_count() -> u32 {
let n = include_str!("../src/shaders/common.wgsl")
.lines()
.filter(|l| l.starts_with("const KIND_"))
.count() as u32;
assert!(n >= 10, "found only {n} KIND_* constants in common.wgsl");
n
}
/// Specialize `module`'s `override`s with `constants` for `entry_point` (as
/// wgpu does at pipeline creation) and compile the result to SPIR-V, so a
/// shader that only breaks once a particular override value folds a branch
/// in or out is still caught.
fn specialize(
name: &str,
module: &naga::Module,
info: &naga::valid::ModuleInfo,
stage: naga::ShaderStage,
entry_point: &str,
constants: &[(&str, f64)],
) {
let mut pc = naga::back::PipelineConstants::default();
for (k, v) in constants {
pc.insert((*k).to_string(), *v);
}
let (module, info) = naga::back::pipeline_constants::process_overrides(
module,
info,
Some((stage, entry_point)),
&pc,
)
.unwrap_or_else(|e| panic!("{name} {entry_point} {constants:?}: override error: {e:?}"));
let pipeline = naga::back::spv::PipelineOptions {
shader_stage: stage,
entry_point: entry_point.to_string(),
};
naga::back::spv::write_vec(
&module,
&info,
&naga::back::spv::Options::default(),
Some(&pipeline),
)
.unwrap_or_else(|e| panic!("{name} {entry_point} {constants:?}: SPIR-V error: {e:?}"));
}
const MANDELBROT_SRC: &str = concat!(
include_str!("../src/shaders/common.wgsl"),
include_str!("../src/shaders/iterate_uniforms.wgsl"),
include_str!("../src/shaders/mandelbrot.wgsl"),
);
#[test]
fn mandelbrot_shader_is_valid() {
validate("mandelbrot.wgsl", MANDELBROT_SRC);
}
/// Every specialization renderer.rs can build (`PipelineKey`: kind × Julia ×
/// DE), for every fragment entry point.
#[test]
fn mandelbrot_shader_specializations_compile() {
let (module, info) = validate("mandelbrot.wgsl", MANDELBROT_SRC);
for kind in 0..kind_count() {
for julia in [0.0, 1.0] {
for de in [0.0, 1.0] {
let constants = [("KIND", kind as f64), ("IS_JULIA", julia), ("DE", de)];
for entry in ["fs_data", "fs_refine", "fs_color"] {
specialize(
"mandelbrot.wgsl",
&module,
&info,
naga::ShaderStage::Fragment,
entry,
&constants,
);
}
}
}
}
}
#[test]
fn colorize_shader_is_valid() {
validate(
"colorize.wgsl",
concat!(
include_str!("../src/shaders/common.wgsl"),
include_str!("../src/shaders/iterate_uniforms.wgsl"),
include_str!("../src/shaders/colorize.wgsl"),
),
"mandelbrot.wgsl",
include_str!("../src/shaders/mandelbrot.wgsl"),
);
}
#[test]
fn blit_shader_is_valid() {
validate(
"blit.wgsl",
concat!(
include_str!("../src/shaders/common.wgsl"),
include_str!("../src/shaders/blit.wgsl"),
),
);
}
const BUDDHABROT_SRC: &str = concat!(
include_str!("../src/shaders/common.wgsl"),
include_str!("../src/shaders/buddhabrot.wgsl"),
);
#[test]
fn buddhabrot_shader_is_valid() {
validate("buddhabrot.wgsl", BUDDHABROT_SRC);
}
/// Every per-kind accumulation pipeline buddhabrot.rs can build.
#[test]
fn buddhabrot_shader_specializations_compile() {
let (module, info) = validate("buddhabrot.wgsl", BUDDHABROT_SRC);
for kind in 0..kind_count() {
specialize(
"buddhabrot.wgsl",
&module,
&info,
naga::ShaderStage::Compute,
"cs_main",
&[("KIND", kind as f64)],
);
}
validate("blit.wgsl", include_str!("../src/shaders/blit.wgsl"));
}