Compare commits
1
Commits
9cc5a80650
..
BLA
| Author | SHA1 | Date | |
|---|---|---|---|
|
|
ea2a23999d |
@@ -1,6 +1,3 @@
|
|||||||
/target
|
/target
|
||||||
Cargo.lock
|
Cargo.lock
|
||||||
dist
|
dist
|
||||||
frames*
|
|
||||||
out.mp4
|
|
||||||
|
|
||||||
|
|||||||
@@ -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
@@ -8,17 +8,14 @@ bytemuck = { version = "1.25.2", features = ["derive"] }
|
|||||||
dashu-float = "0.6.0"
|
dashu-float = "0.6.0"
|
||||||
eframe = { version = "0.36.2", default-features = false, features = ["wgpu", "default_fonts", "x11", "wayland", "accesskit"] }
|
eframe = { version = "0.36.2", default-features = false, features = ["wgpu", "default_fonts", "x11", "wayland", "accesskit"] }
|
||||||
egui = "0.36.2"
|
egui = "0.36.2"
|
||||||
glam = "0.33.8"
|
futures-channel = { version = "0.3.34", default-features = false, features = ["alloc", "std"] }
|
||||||
log = "0.4.34"
|
log = "0.4.34"
|
||||||
png = "0.18.1"
|
png = "0.18.1"
|
||||||
|
|
||||||
[target.'cfg(not(target_arch = "wasm32"))'.dependencies]
|
[target.'cfg(not(target_arch = "wasm32"))'.dependencies]
|
||||||
env_logger = "0.11.11"
|
env_logger = "0.11.11"
|
||||||
clap = { version = "4.5.51", features = ["derive"] }
|
|
||||||
pollster = "1.0.1"
|
|
||||||
|
|
||||||
[target.'cfg(target_arch = "wasm32")'.dependencies]
|
[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_error_panic_hook = "0.1.7"
|
||||||
console_log = "1.1.0"
|
console_log = "1.1.0"
|
||||||
js-sys = "0.3.105"
|
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).
|
# Release: optimize hard (fractal math is hot).
|
||||||
[profile.release]
|
[profile.release]
|
||||||
opt-level = 3
|
opt-level = 3
|
||||||
# codegen-units = 1
|
|
||||||
debug = true
|
|
||||||
|
|
||||||
# Dev: keep our own crate debuggable, but optimize dependencies (dashu, wgpu,
|
# Dev: keep our own crate debuggable, but optimize dependencies (dashu, wgpu,
|
||||||
# egui) so the explorer is actually interactive during development.
|
# egui) so the explorer is actually interactive during development.
|
||||||
@@ -41,4 +36,4 @@ opt-level = 1
|
|||||||
opt-level = 3
|
opt-level = 3
|
||||||
|
|
||||||
[dev-dependencies]
|
[dev-dependencies]
|
||||||
naga = { version = "30", features = ["wgsl-in", "spv-out"] }
|
naga = { version = "30", features = ["wgsl-in"] }
|
||||||
|
|||||||
@@ -20,7 +20,6 @@ wasm-bindgen \
|
|||||||
target/wasm32-unknown-unknown/release/mandelbrot.wasm
|
target/wasm32-unknown-unknown/release/mandelbrot.wasm
|
||||||
|
|
||||||
cp index.html "$OUT/index.html"
|
cp index.html "$OUT/index.html"
|
||||||
cp favicon.ico "$OUT/favicon.ico"
|
|
||||||
|
|
||||||
echo "==> done: $OUT/ (index.html, mandelbrot.js, mandelbrot_bg.wasm)"
|
echo "==> done: $OUT/ (index.html, mandelbrot.js, mandelbrot_bg.wasm)"
|
||||||
echo " serve: python3 -m http.server -d $OUT 8080"
|
echo " serve: python3 -m http.server -d $OUT 8080"
|
||||||
|
|||||||
BIN
Binary file not shown.
|
Before Width: | Height: | Size: 422 KiB |
+256
-1555
File diff suppressed because it is too large
Load Diff
@@ -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
@@ -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,
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
@@ -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);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
@@ -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
@@ -1,18 +1,13 @@
|
|||||||
//! GPU fractal rendering: wgpu pipeline, uniforms, reference orbit, and the
|
//! GPU fractal rendering: wgpu pipeline, uniforms, reference orbit, and the
|
||||||
//! egui paint callback.
|
//! egui paint callback.
|
||||||
|
|
||||||
pub mod buddhabrot;
|
|
||||||
pub mod kind;
|
|
||||||
pub mod reference;
|
pub mod reference;
|
||||||
pub mod renderer;
|
pub mod renderer;
|
||||||
pub mod share;
|
pub mod share;
|
||||||
|
|
||||||
pub use buddhabrot::{BuddhabrotCallback, BuddhabrotRenderer, BuddhabrotUniforms};
|
pub use reference::{Bla, FractalKind, build_bla_table, compute_reference, compute_set_reference};
|
||||||
pub use kind::FractalKind;
|
pub use renderer::{
|
||||||
pub use reference::{compute_reference, compute_set_reference};
|
ExportRender, FractalCallback, FractalRenderer, MAX_REF_POINTS, Uniforms,
|
||||||
#[cfg(target_arch = "wasm32")]
|
encode_png_with_progress,
|
||||||
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 share::ShareState;
|
pub use share::ShareState;
|
||||||
|
|||||||
+295
-501
@@ -10,34 +10,42 @@
|
|||||||
//! * Mandelbrot-set: `z0 = 0`, `c = view center` (the c-plane point per pixel).
|
//! * Mandelbrot-set: `z0 = 0`, `c = view center` (the c-plane point per pixel).
|
||||||
//! * Julia-set: `z0 = view center`, `c = fractal constant` (fixed per view).
|
//! * Julia-set: `z0 = view center`, `c = fractal constant` (fixed per view).
|
||||||
|
|
||||||
use super::kind::FractalKind;
|
use crate::view::Big;
|
||||||
use crate::view::{Big, big_from_f64};
|
|
||||||
|
/// 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
|
/// 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 so pixels escaping alongside the reference can still reach their
|
||||||
/// bailout before the stored orbit runs out.
|
/// bailout before the stored orbit runs out.
|
||||||
const REFERENCE_ESCAPE_SQ: f64 = 1.0e10;
|
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
|
/// 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
|
/// `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.
|
/// to `max_iter` steps at `precision` bits. Each entry is `[re, im]` in f32.
|
||||||
#[allow(clippy::too_many_arguments)]
|
|
||||||
pub fn compute_reference(
|
pub fn compute_reference(
|
||||||
z0_re: &Big,
|
z0_re: &Big,
|
||||||
z0_im: &Big,
|
z0_im: &Big,
|
||||||
@@ -47,143 +55,12 @@ pub fn compute_reference(
|
|||||||
precision: usize,
|
precision: usize,
|
||||||
kind: FractalKind,
|
kind: FractalKind,
|
||||||
power: u32,
|
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]> {
|
) -> Vec<[f32; 2]> {
|
||||||
let cr = c_re.clone().with_precision(precision).value();
|
let cr = c_re.clone().with_precision(precision).value();
|
||||||
let ci = c_im.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 zr = z0_re.clone().with_precision(precision).value();
|
||||||
let mut zi = z0_im.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);
|
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 {
|
let (new_zr, new_zi) = match kind {
|
||||||
FractalKind::Mandelbrot => {
|
FractalKind::Mandelbrot => {
|
||||||
// Z^2 = (zr^2 - zi^2) + (2 zr zi) i, with zr^2 - zi^2 as
|
// Z^2 = (zr^2 - zi^2) + (2 zr zi) i.
|
||||||
// (zr + zi)(zr - zi): one multiply instead of two squares.
|
let re = &zr.sqr() - &zi.sqr() + &cr;
|
||||||
let re = (&zr + &zi) * (&zr - &zi) + &cr;
|
|
||||||
let im = ((&zr * &zi) << 1) + &ci; // << 1 is exact ×2 in base 2
|
let im = ((&zr * &zi) << 1) + &ci; // << 1 is exact ×2 in base 2
|
||||||
(re, im)
|
(re, im)
|
||||||
}
|
}
|
||||||
@@ -221,53 +97,8 @@ fn compute_reference_big(
|
|||||||
let (pr, pi) = complex_pow(&zr, &zi, power.max(2), precision);
|
let (pr, pi) = complex_pow(&zr, &zi, power.max(2), precision);
|
||||||
(pr + &cr, pi + &ci)
|
(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();
|
zr = new_zr.with_precision(precision).value();
|
||||||
zi = new_zi.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)
|
(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`,
|
/// Convenience: parameter-plane ("Mandelbrot-set") reference (`z0 = 0`,
|
||||||
/// `c = center`) for any `kind`.
|
/// `c = center`) for any `kind`.
|
||||||
#[allow(clippy::too_many_arguments)]
|
|
||||||
pub fn compute_set_reference(
|
pub fn compute_set_reference(
|
||||||
center_re: &Big,
|
center_re: &Big,
|
||||||
center_im: &Big,
|
center_im: &Big,
|
||||||
@@ -335,26 +139,157 @@ pub fn compute_set_reference(
|
|||||||
precision: usize,
|
precision: usize,
|
||||||
kind: FractalKind,
|
kind: FractalKind,
|
||||||
power: u32,
|
power: u32,
|
||||||
phoenix_p: (f64, f64),
|
|
||||||
lambda_l: (f64, f64),
|
|
||||||
complex_power: (f64, f64),
|
|
||||||
) -> Vec<[f32; 2]> {
|
) -> Vec<[f32; 2]> {
|
||||||
let zero = big_zero(precision);
|
let zero = big_zero(precision);
|
||||||
compute_reference(
|
compute_reference(
|
||||||
&zero,
|
&zero, &zero, center_re, center_im, max_iter, precision, kind, power,
|
||||||
&zero,
|
|
||||||
center_re,
|
|
||||||
center_im,
|
|
||||||
max_iter,
|
|
||||||
precision,
|
|
||||||
kind,
|
|
||||||
power,
|
|
||||||
phoenix_p,
|
|
||||||
lambda_l,
|
|
||||||
complex_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)]
|
#[cfg(test)]
|
||||||
mod tests {
|
mod tests {
|
||||||
use super::*;
|
use super::*;
|
||||||
@@ -365,17 +300,7 @@ mod tests {
|
|||||||
fn reference_matches_naive_f64() {
|
fn reference_matches_naive_f64() {
|
||||||
let cr = Big::try_from(-0.75_f64).unwrap();
|
let cr = Big::try_from(-0.75_f64).unwrap();
|
||||||
let ci = Big::try_from(0.1_f64).unwrap();
|
let ci = Big::try_from(0.1_f64).unwrap();
|
||||||
let points = compute_set_reference(
|
let points = compute_set_reference(&cr, &ci, 60, 200, FractalKind::Mandelbrot, 2);
|
||||||
&cr,
|
|
||||||
&ci,
|
|
||||||
60,
|
|
||||||
200,
|
|
||||||
FractalKind::Mandelbrot,
|
|
||||||
2,
|
|
||||||
(0.0, 0.0),
|
|
||||||
(0.0, 0.0),
|
|
||||||
(0.0, 0.0),
|
|
||||||
);
|
|
||||||
|
|
||||||
// Independent naive f64 orbit.
|
// Independent naive f64 orbit.
|
||||||
let (c_re, c_im) = (-0.75_f64, 0.1_f64);
|
let (c_re, c_im) = (-0.75_f64, 0.1_f64);
|
||||||
@@ -385,14 +310,8 @@ mod tests {
|
|||||||
// significant figures.
|
// significant figures.
|
||||||
let tol_re = 1e-4 * (1.0 + zr.abs());
|
let tol_re = 1e-4 * (1.0 + zr.abs());
|
||||||
let tol_im = 1e-4 * (1.0 + zi.abs());
|
let tol_im = 1e-4 * (1.0 + zi.abs());
|
||||||
assert!(
|
assert!((point[0] as f64 - zr).abs() < tol_re, "re mismatch: {point:?} vs {zr}");
|
||||||
(point[0] as f64 - zr).abs() < tol_re,
|
assert!((point[1] as f64 - zi).abs() < tol_im, "im mismatch: {point:?} vs {zi}");
|
||||||
"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 nzr = zr * zr - zi * zi + c_re;
|
||||||
let nzi = 2.0 * zr * zi + c_im;
|
let nzi = 2.0 * zr * zi + c_im;
|
||||||
zr = nzr;
|
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.
|
/// A point inside the main cardioid never escapes: full-length orbit.
|
||||||
#[test]
|
#[test]
|
||||||
fn interior_orbit_runs_full_length() {
|
fn interior_orbit_runs_full_length() {
|
||||||
let cr = Big::try_from(-0.2_f64).unwrap();
|
let cr = Big::try_from(-0.2_f64).unwrap();
|
||||||
let ci = Big::try_from(0.0_f64).unwrap();
|
let ci = Big::try_from(0.0_f64).unwrap();
|
||||||
let points = compute_set_reference(
|
let points = compute_set_reference(&cr, &ci, 500, 120, FractalKind::Mandelbrot, 2);
|
||||||
&cr,
|
|
||||||
&ci,
|
|
||||||
500,
|
|
||||||
120,
|
|
||||||
FractalKind::Mandelbrot,
|
|
||||||
2,
|
|
||||||
(0.0, 0.0),
|
|
||||||
(0.0, 0.0),
|
|
||||||
(0.0, 0.0),
|
|
||||||
);
|
|
||||||
assert_eq!(points.len(), 501, "interior orbit should not escape");
|
assert_eq!(points.len(), 501, "interior orbit should not escape");
|
||||||
}
|
}
|
||||||
|
|
||||||
@@ -462,17 +333,7 @@ mod tests {
|
|||||||
fn burning_ship_reference_matches_naive_f64() {
|
fn burning_ship_reference_matches_naive_f64() {
|
||||||
let cr = Big::try_from(-1.75_f64).unwrap();
|
let cr = Big::try_from(-1.75_f64).unwrap();
|
||||||
let ci = Big::try_from(-0.03_f64).unwrap();
|
let ci = Big::try_from(-0.03_f64).unwrap();
|
||||||
let points = compute_set_reference(
|
let points = compute_set_reference(&cr, &ci, 60, 200, FractalKind::BurningShip, 2);
|
||||||
&cr,
|
|
||||||
&ci,
|
|
||||||
60,
|
|
||||||
200,
|
|
||||||
FractalKind::BurningShip,
|
|
||||||
2,
|
|
||||||
(0.0, 0.0),
|
|
||||||
(0.0, 0.0),
|
|
||||||
(0.0, 0.0),
|
|
||||||
);
|
|
||||||
|
|
||||||
let (c_re, c_im) = (-1.75_f64, -0.03_f64);
|
let (c_re, c_im) = (-1.75_f64, -0.03_f64);
|
||||||
let (mut zr, mut zi) = (0.0_f64, 0.0_f64);
|
let (mut zr, mut zi) = (0.0_f64, 0.0_f64);
|
||||||
@@ -492,17 +353,7 @@ mod tests {
|
|||||||
fn multibrot3_reference_matches_naive_f64() {
|
fn multibrot3_reference_matches_naive_f64() {
|
||||||
let cr = Big::try_from(0.3_f64).unwrap();
|
let cr = Big::try_from(0.3_f64).unwrap();
|
||||||
let ci = Big::try_from(0.2_f64).unwrap();
|
let ci = Big::try_from(0.2_f64).unwrap();
|
||||||
let points = compute_set_reference(
|
let points = compute_set_reference(&cr, &ci, 60, 200, FractalKind::Multibrot, 3);
|
||||||
&cr,
|
|
||||||
&ci,
|
|
||||||
60,
|
|
||||||
200,
|
|
||||||
FractalKind::Multibrot,
|
|
||||||
3,
|
|
||||||
(0.0, 0.0),
|
|
||||||
(0.0, 0.0),
|
|
||||||
(0.0, 0.0),
|
|
||||||
);
|
|
||||||
|
|
||||||
let (c_re, c_im) = (0.3_f64, 0.2_f64);
|
let (c_re, c_im) = (0.3_f64, 0.2_f64);
|
||||||
let (mut zr, mut zi) = (0.0_f64, 0.0_f64);
|
let (mut zr, mut zi) = (0.0_f64, 0.0_f64);
|
||||||
@@ -535,9 +386,6 @@ mod tests {
|
|||||||
200,
|
200,
|
||||||
FractalKind::Mandelbrot,
|
FractalKind::Mandelbrot,
|
||||||
2,
|
2,
|
||||||
(0.0, 0.0),
|
|
||||||
(0.0, 0.0),
|
|
||||||
(0.0, 0.0),
|
|
||||||
);
|
);
|
||||||
|
|
||||||
let (mut zr, mut zi) = (0.15_f64, -0.1_f64);
|
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]
|
#[test]
|
||||||
fn celtic_reference_matches_naive_f64() {
|
fn bla_matches_step_by_step() {
|
||||||
let cr = Big::try_from(-0.6_f64).unwrap();
|
// An interior center (never escapes), so `e` stays bounded and we can
|
||||||
let ci = Big::try_from(0.4_f64).unwrap();
|
// iterate the full orbit; zoomed so |dc| ~ 1e-9 (deep enough for big
|
||||||
let points = compute_set_reference(
|
// skips). Correctness of the walk is independent of which orbit we pick.
|
||||||
&cr,
|
let cr = Big::try_from(-0.5_f64).unwrap();
|
||||||
&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();
|
|
||||||
let ci = Big::try_from(0.0_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, 800, 160, FractalKind::Mandelbrot, 2);
|
||||||
let points = compute_set_reference(
|
assert!(points.len() > 64, "need a long orbit to exercise BLA levels");
|
||||||
&cr,
|
|
||||||
&ci,
|
let half_height = 1.0e-9_f64;
|
||||||
60,
|
let aspect = 1.5_f64;
|
||||||
200,
|
let dc_max = half_height * (1.0 + aspect * aspect).sqrt();
|
||||||
FractalKind::Phoenix,
|
let table = build_bla_table(&points, dc_max);
|
||||||
2,
|
assert!(!table.is_empty());
|
||||||
p,
|
|
||||||
(0.0, 0.0),
|
let target_n = points.len() - 1;
|
||||||
|
// A few pixel offsets across the view (all within dc_max).
|
||||||
|
let offsets = [
|
||||||
(0.0, 0.0),
|
(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;
|
|
||||||
}
|
}
|
||||||
}
|
assert!(any_skip, "BLA never took a multi-step skip — test is not exercising it");
|
||||||
|
|
||||||
/// 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;
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|||||||
+180
-719
File diff suppressed because it is too large
Load Diff
+16
-39
@@ -21,40 +21,31 @@ pub struct ShareState {
|
|||||||
pub half_height: f64,
|
pub half_height: f64,
|
||||||
pub iterations: u32,
|
pub iterations: u32,
|
||||||
pub julia_c: (f64, f64),
|
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_scale: f32,
|
||||||
pub color_offset: f32,
|
pub color_offset: f32,
|
||||||
/// Palette index (`palette_id` in the shader).
|
/// Palette index (`palette_id` in the shader).
|
||||||
pub palette: u32,
|
pub palette: u32,
|
||||||
/// Shadow palette index (`shadow_palette_id` in the shader).
|
|
||||||
pub shadow_palette: u32,
|
|
||||||
}
|
}
|
||||||
|
|
||||||
impl ShareState {
|
impl ShareState {
|
||||||
pub fn encode(&self) -> String {
|
pub fn encode(&self) -> String {
|
||||||
let mut s = String::new();
|
let mut s = String::new();
|
||||||
s.push_str(if self.julia { "m=j" } else { "m=m" });
|
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!("&pw={}", self.power));
|
||||||
s.push_str(&format!(
|
s.push_str(&format!(
|
||||||
"&re={}&im={}&hh={}&it={}",
|
"&re={}&im={}&hh={}&it={}",
|
||||||
self.center_re, self.center_im, self.half_height, self.iterations
|
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!("&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!(
|
s.push_str(&format!(
|
||||||
"&cpr={}&cpi={}",
|
"&cs={}&co={}&pal={}",
|
||||||
self.complex_power.0, self.complex_power.1
|
self.color_scale, self.color_offset, self.palette
|
||||||
));
|
|
||||||
s.push_str(&format!(
|
|
||||||
"&cs={}&co={}&pal={}&spal={}",
|
|
||||||
self.color_scale, self.color_offset, self.palette, self.shadow_palette
|
|
||||||
));
|
));
|
||||||
s
|
s
|
||||||
}
|
}
|
||||||
@@ -72,7 +63,13 @@ impl ShareState {
|
|||||||
julia: map.get("m").map(|m| *m == "j").unwrap_or(false),
|
julia: map.get("m").map(|m| *m == "j").unwrap_or(false),
|
||||||
kind: map
|
kind: map
|
||||||
.get("f")
|
.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),
|
.unwrap_or(FractalKind::Mandelbrot),
|
||||||
power: map.get("pw").and_then(|s| s.parse().ok()).unwrap_or(2),
|
power: map.get("pw").and_then(|s| s.parse().ok()).unwrap_or(2),
|
||||||
center_re: (*map.get("re")?).to_string(),
|
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("jr").and_then(|s| s.parse().ok()).unwrap_or(-0.8),
|
||||||
map.get("ji").and_then(|s| s.parse().ok()).unwrap_or(0.156),
|
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_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),
|
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),
|
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() {
|
fn round_trip() {
|
||||||
let s = ShareState {
|
let s = ShareState {
|
||||||
julia: true,
|
julia: true,
|
||||||
kind: FractalKind::Phoenix,
|
kind: FractalKind::Multibrot,
|
||||||
power: 5,
|
power: 5,
|
||||||
center_re: "-0.743643887037158704752191506114774".into(),
|
center_re: "-0.743643887037158704752191506114774".into(),
|
||||||
center_im: "0.131825904205311970493132056385139".into(),
|
center_im: "0.131825904205311970493132056385139".into(),
|
||||||
half_height: 1.5e-20,
|
half_height: 1.5e-20,
|
||||||
iterations: 4000,
|
iterations: 4000,
|
||||||
julia_c: (-0.123, 0.745),
|
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_scale: 0.02,
|
||||||
color_offset: 0.25,
|
color_offset: 0.25,
|
||||||
palette: 3,
|
palette: 3,
|
||||||
shadow_palette: 1,
|
|
||||||
};
|
};
|
||||||
let d = ShareState::decode(&s.encode()).unwrap();
|
let d = ShareState::decode(&s.encode()).unwrap();
|
||||||
assert_eq!(d.julia, s.julia);
|
assert_eq!(d.julia, s.julia);
|
||||||
@@ -135,10 +115,7 @@ mod tests {
|
|||||||
assert_eq!(d.half_height, s.half_height);
|
assert_eq!(d.half_height, s.half_height);
|
||||||
assert_eq!(d.iterations, s.iterations);
|
assert_eq!(d.iterations, s.iterations);
|
||||||
assert_eq!(d.julia_c, s.julia_c);
|
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.palette, s.palette);
|
||||||
assert_eq!(d.shadow_palette, s.shadow_palette);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
#[test]
|
#[test]
|
||||||
|
|||||||
-245
@@ -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}"))
|
|
||||||
}
|
|
||||||
@@ -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
@@ -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.
|
// Fractal Explorer — Rust + wgpu + egui + WGSL deep-zoom Mandelbrot.
|
||||||
//
|
//
|
||||||
// A single binary drives both native and web (WASM/WebGPU) builds; the two
|
// 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>.
|
// and calls the wasm `main`, which boots eframe onto the page's <canvas>.
|
||||||
|
|
||||||
mod app;
|
mod app;
|
||||||
mod camera;
|
|
||||||
mod fractal;
|
mod fractal;
|
||||||
mod lights;
|
|
||||||
mod view;
|
mod view;
|
||||||
|
|
||||||
#[cfg(not(target_arch = "wasm32"))]
|
|
||||||
mod cli;
|
|
||||||
#[cfg(not(target_arch = "wasm32"))]
|
|
||||||
mod headless;
|
|
||||||
#[cfg(not(target_arch = "wasm32"))]
|
#[cfg(not(target_arch = "wasm32"))]
|
||||||
mod worker;
|
mod worker;
|
||||||
|
|
||||||
@@ -35,12 +25,13 @@ fn wgpu_options() -> eframe::egui_wgpu::WgpuConfiguration {
|
|||||||
|
|
||||||
let mut options = eframe::egui_wgpu::WgpuConfiguration::default();
|
let mut options = eframe::egui_wgpu::WgpuConfiguration::default();
|
||||||
if let WgpuSetup::CreateNew(setup) = &mut options.wgpu_setup {
|
if let WgpuSetup::CreateNew(setup) = &mut options.wgpu_setup {
|
||||||
setup.device_descriptor =
|
setup.device_descriptor = std::sync::Arc::new(|adapter: &wgpu::Adapter| {
|
||||||
std::sync::Arc::new(|adapter: &wgpu::Adapter| wgpu::DeviceDescriptor {
|
wgpu::DeviceDescriptor {
|
||||||
label: Some("fractal wgpu device"),
|
label: Some("fractal wgpu device"),
|
||||||
required_features: wgpu::Features::empty(),
|
required_features: wgpu::Features::empty(),
|
||||||
required_limits: adapter.limits(),
|
required_limits: adapter.limits(),
|
||||||
..Default::default()
|
..Default::default()
|
||||||
|
}
|
||||||
});
|
});
|
||||||
#[cfg(target_arch = "wasm32")]
|
#[cfg(target_arch = "wasm32")]
|
||||||
{
|
{
|
||||||
@@ -52,24 +43,11 @@ fn wgpu_options() -> eframe::egui_wgpu::WgpuConfiguration {
|
|||||||
|
|
||||||
#[cfg(not(target_arch = "wasm32"))]
|
#[cfg(not(target_arch = "wasm32"))]
|
||||||
fn main() -> eframe::Result {
|
fn main() -> eframe::Result {
|
||||||
use clap::Parser as _;
|
|
||||||
|
|
||||||
env_logger::builder()
|
env_logger::builder()
|
||||||
.filter_level(log::LevelFilter::Info)
|
.filter_level(log::LevelFilter::Info)
|
||||||
.parse_default_env()
|
.parse_default_env()
|
||||||
.init();
|
.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 {
|
let native_options = eframe::NativeOptions {
|
||||||
renderer: eframe::Renderer::Wgpu,
|
renderer: eframe::Renderer::Wgpu,
|
||||||
wgpu_options: wgpu_options(),
|
wgpu_options: wgpu_options(),
|
||||||
@@ -123,7 +101,9 @@ fn main() {
|
|||||||
match result {
|
match result {
|
||||||
Ok(_) => loading.remove(),
|
Ok(_) => loading.remove(),
|
||||||
Err(e) => {
|
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:?}");
|
log::error!("failed to start eframe: {e:?}");
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|||||||
@@ -13,7 +13,12 @@ struct VsOut {
|
|||||||
|
|
||||||
@vertex
|
@vertex
|
||||||
fn vs_main(@builtin(vertex_index) idx: u32) -> VsOut {
|
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;
|
var out: VsOut;
|
||||||
out.pos = vec4<f32>(p, 0.0, 1.0);
|
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
|
// Map NDC to texture UV. v is flipped so the cache's top row (rendered at
|
||||||
|
|||||||
@@ -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);
|
|
||||||
}
|
|
||||||
@@ -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);
|
|
||||||
}
|
|
||||||
@@ -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;
|
|
||||||
@@ -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
@@ -12,25 +12,46 @@
|
|||||||
// the reference index to 0 and carry the full value as the new delta (valid
|
// the reference index to 0 and carry the full value as the new delta (valid
|
||||||
// because X_0 = 0).
|
// 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(0) var<uniform> u: Uniforms;
|
||||||
@group(0) @binding(1) var<storage, read> ref_orbit: array<vec2<f32>>;
|
@group(0) @binding(1) var<storage, read> ref_orbit: array<vec2<f32>>;
|
||||||
// Only read by `fs_color`'s shadow branch (custom-lights palette); the
|
// BLA table, levels concatenated (level 0 first). Per-level counts are recomputed
|
||||||
// iteration pass (`fs_data`) never touches it.
|
// from `ref_len` exactly as the CPU builder laid them out.
|
||||||
@group(0) @binding(2) var<uniform> lights: array<Light, 16>;
|
@group(0) @binding(2) var<storage, read> bla_table: array<Bla>;
|
||||||
// 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;
|
|
||||||
|
|
||||||
struct VsOut {
|
struct VsOut {
|
||||||
@builtin(position) pos: vec4<f32>,
|
@builtin(position) pos: vec4<f32>,
|
||||||
@@ -40,7 +61,12 @@ struct VsOut {
|
|||||||
|
|
||||||
@vertex
|
@vertex
|
||||||
fn vs_main(@builtin(vertex_index) idx: u32) -> VsOut {
|
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;
|
var out: VsOut;
|
||||||
out.pos = vec4<f32>(ndc, 0.0, 1.0);
|
out.pos = vec4<f32>(ndc, 0.0, 1.0);
|
||||||
// Flip y so +imaginary points up the screen.
|
// Flip y so +imaginary points up the screen.
|
||||||
@@ -48,100 +74,60 @@ fn vs_main(@builtin(vertex_index) idx: u32) -> VsOut {
|
|||||||
return out;
|
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.
|
// Complex conjugate.
|
||||||
fn conj(a: vec2<f32>) -> vec2<f32> {
|
fn conj(a: vec2<f32>) -> vec2<f32> {
|
||||||
return vec2<f32>(a.x, -a.y);
|
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
|
// |c + d| - |c|, evaluated exactly (no catastrophic cancellation even when the
|
||||||
// sum crosses zero). This is what makes the Burning Ship delta correct through
|
// 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.
|
// the sign flips that happen all along the axes, where the ship's detail lives.
|
||||||
fn diffabs(c: f32, d: f32) -> f32 {
|
fn diffabs(c: f32, d: f32) -> f32 {
|
||||||
let cd = c + d;
|
let cd = c + d;
|
||||||
if c >= 0.0 {
|
if (c >= 0.0) {
|
||||||
return select(-(2.0 * c + d), d, cd >= 0.0);
|
return select(-(2.0 * c + d), d, cd >= 0.0);
|
||||||
}
|
}
|
||||||
return select(-d, 2.0 * c + 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}.
|
// Binomial coefficient C(n, k) as f32 (exact for the small powers we use).
|
||||||
// The large z^p term is never formed (that would cancel catastrophically), and
|
fn binom(n: u32, k: u32) -> f32 {
|
||||||
// the sum is evaluated Horner-style (s <- s*(Z+e) + Z^j) so it needs neither a
|
var num = 1.0;
|
||||||
// table of powers (a dynamically indexed local array spills to slow memory on
|
var den = 1.0;
|
||||||
// most GPUs) nor binomial coefficients. Forming Z+e rounds e away when it's
|
for (var i: u32 = 0u; i < k; i = i + 1u) {
|
||||||
// tiny, but that only perturbs `s` by a relative f32 epsilon, and the result
|
num = num * f32(n - i);
|
||||||
// is `e * s`, so the delta keeps full relative precision.
|
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> {
|
fn multibrot_delta(z: vec2<f32>, e: vec2<f32>, p: u32) -> vec2<f32> {
|
||||||
let y = z + e;
|
var zp: array<vec2<f32>, 9>; // Z^0 .. Z^8
|
||||||
var s = vec2<f32>(1.0, 0.0);
|
zp[0] = vec2<f32>(1.0, 0.0);
|
||||||
var zj = vec2<f32>(1.0, 0.0);
|
for (var j: u32 = 1u; j <= p; j = j + 1u) {
|
||||||
for (var j: u32 = 1u; j < p; j = j + 1u) {
|
zp[j] = cmul(zp[j - 1u], z);
|
||||||
zj = cmul(zj, z); // Z^j
|
|
||||||
s = cmul(s, y) + zj;
|
|
||||||
}
|
}
|
||||||
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);
|
var acc = vec2<f32>(0.0, 0.0);
|
||||||
for (var k: u32 = 1u; k <= COMPLEX_MULTIBROT_TERMS; k = k + 1u) {
|
var ek = vec2<f32>(1.0, 0.0); // e^0
|
||||||
acc = acc + cmul(cm_coef(k), wk);
|
for (var k: u32 = 1u; k <= p; k = k + 1u) {
|
||||||
wk = cmul(wk, w);
|
ek = cmul(ek, e); // e^k
|
||||||
if dot(wk, wk) < 1e-18 * dot(acc, acc) {
|
acc = acc + binom(p, k) * cmul(zp[p - k], ek);
|
||||||
break;
|
|
||||||
}
|
}
|
||||||
}
|
return acc;
|
||||||
return cmul(cpow(z, p), acc);
|
|
||||||
}
|
|
||||||
return cpow(z + e, p) - cpow(z, p);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
// One perturbation step of the current fractal's delta: e -> f(Z+e) - f(Z),
|
// 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
|
// where `z` is the reference orbit value X_m. `step_add` (dc) is added by the
|
||||||
// caller. Must match `FractalKind` on the CPU side.
|
// caller. Must match `FractalKind` on the CPU side.
|
||||||
fn advance_delta(z: vec2<f32>, e: vec2<f32>) -> vec2<f32> {
|
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
|
// (|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
|
// 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
|
// 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 base = 2.0 * cmul(z, e) + cmul(e, e);
|
||||||
let dp = z.x * e.y + z.y * e.x + e.x * e.y;
|
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));
|
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 cz = conj(z);
|
||||||
let ce = conj(e);
|
let ce = conj(e);
|
||||||
return 2.0 * cmul(cz, ce) + cmul(ce, ce);
|
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));
|
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
|
// 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
|
// Burning Ship / Tricorn we use |f'| ~ |2Z|, which keeps the DE magnitude close
|
||||||
// enough to de-speckle filaments.
|
// enough to de-speckle filaments.
|
||||||
fn fprime(z: vec2<f32>) -> vec2<f32> {
|
fn fprime(z: vec2<f32>) -> vec2<f32> {
|
||||||
if KIND == KIND_MULTIBROT {
|
if (u.kind == KIND_MULTIBROT) {
|
||||||
let p = clamp(u.power, 2u, 8u);
|
let p = clamp(u.power, 2u, 8u);
|
||||||
var zk = z; // Z^1
|
var zk = vec2<f32>(1.0, 0.0); // Z^0
|
||||||
for (var k: u32 = 2u; k < p; k = k + 1u) {
|
for (var k: u32 = 1u; k < p; k = k + 1u) {
|
||||||
zk = cmul(zk, z); // -> Z^{p-1}
|
zk = cmul(zk, z); // -> Z^{p-1}
|
||||||
}
|
}
|
||||||
return f32(p) * zk;
|
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;
|
return 2.0 * z;
|
||||||
}
|
}
|
||||||
|
|
||||||
// Periodicity (interior) detection, Brent-style: the full orbit value is
|
// Smooth cyclic palettes (Inigo Quilez cosine palettes), selected by id.
|
||||||
// saved at iterations PERIOD_FIRST_CHECK, 2x that, 4x ..., and every later
|
fn palette(id: u32, t: f32) -> vec3<f32> {
|
||||||
// iterate is compared against the last saved one. Returning within
|
if (id == 4u) {
|
||||||
// PERIOD_EPS2 (relative, squared) means the orbit has closed a cycle.
|
return vec3<f32>(t, t, t); // grayscale
|
||||||
//
|
|
||||||
// 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;
|
|
||||||
}
|
}
|
||||||
if KIND == KIND_COMPLEX_MULTIBROT && !DE {
|
let a = vec3<f32>(0.5, 0.5, 0.5);
|
||||||
return false;
|
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,
|
// Result of a BLA lookup at an orbit index.
|
||||||
// `de` the distance-estimate darkening factor in [0,1], `escaped` false for the
|
struct Hop {
|
||||||
// interior of the set. Splitting iteration from coloring lets a colour change be
|
found: bool,
|
||||||
// remapped cheaply (see the colourise pass) without re-iterating.
|
a: vec2<f32>,
|
||||||
struct Sample {
|
b: vec2<f32>,
|
||||||
ci: f32,
|
l: u32,
|
||||||
de: f32,
|
|
||||||
escaped: bool,
|
|
||||||
};
|
};
|
||||||
|
|
||||||
// Perturbation iterate a single sample. `offset` is the per-pixel offset in
|
// Longest valid BLA skip starting at orbit index `n` for a delta of squared
|
||||||
// complex units. For Mandelbrot it is the c-plane offset added every step (delta
|
// magnitude `emag2`. Walks levels low->high, recomputing each level's flat-array
|
||||||
// starts at 0); for Julia it is the z-plane offset that seeds the initial delta
|
// start and count from `ref_len` (count[0] = ref_len-1, halving each level).
|
||||||
// (c is fixed, so nothing is added per step).
|
// Radii shrink with level and alignment is monotonic, so the valid levels form a
|
||||||
fn iterate_sample(offset: vec2<f32>, px: f32) -> Sample {
|
// prefix: we keep the last valid one and stop at the first that fails.
|
||||||
// Loop invariants, read once instead of on every iteration.
|
fn bla_find(n: u32, emag2: f32) -> Hop {
|
||||||
let max_iter = u.max_iter;
|
var res: Hop;
|
||||||
let bailout_sq = u.bailout_sq;
|
res.found = false;
|
||||||
let ref_len = u.ref_len;
|
|
||||||
|
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)
|
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 step_add = offset;
|
||||||
var e = vec2<f32>(0.0, 0.0);
|
var e = vec2<f32>(0.0, 0.0);
|
||||||
// Orbit derivative for distance estimation, pre-multiplied by the pixel
|
// Orbit derivative for distance estimation. For the set plane it is d/dc
|
||||||
// size `px`. For the set plane it is px·d/dc (starts at 0, gains +px each
|
// (starts at 0, gains +1 each step); for Julia it is d/dz0 (starts at 1).
|
||||||
// step); for Julia it is px·d/dz0 (starts at px). The raw derivative grows
|
var dz = vec2<f32>(0.0, 0.0);
|
||||||
// like 1/px, so unscaled its square overflows f32 at deep zoom (~1e-12),
|
var dz_seed = vec2<f32>(1.0, 0.0);
|
||||||
// which zeroed DE along iteration bands; scaled, it stays ~|z|ln|z| / DE
|
if (u.is_julia != 0u) {
|
||||||
// in pixels at any depth.
|
|
||||||
var dzs = vec2<f32>(0.0, 0.0);
|
|
||||||
if IS_JULIA {
|
|
||||||
step_add = vec2<f32>(0.0, 0.0);
|
step_add = vec2<f32>(0.0, 0.0);
|
||||||
e = offset;
|
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 n: u32 = 0u; // total iteration count
|
||||||
var xm = z0; // X[m], carried so each step loads the orbit once
|
var z = vec2<f32>(0.0, 0.0); // full value y_n, kept for coloring
|
||||||
var z = xm + e; // full value y_n, kept for coloring
|
|
||||||
var z2 = dot(z, z);
|
|
||||||
var escaped = false;
|
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 {
|
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;
|
escaped = true;
|
||||||
break;
|
break;
|
||||||
}
|
}
|
||||||
if n >= max_iter {
|
if (n >= u.max_iter) {
|
||||||
break; // interior
|
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,
|
// Propagate the derivative of the full orbit (unaffected by rebasing,
|
||||||
// which only re-expresses the same value). Only when DE is enabled.
|
// 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.
|
if (u.de_coloring != 0u) {
|
||||||
// f'(z) of this step, shared by DE and the periodicity multiplier.
|
dz = cmul(fprime(z), dz) + dz_seed;
|
||||||
var fp = vec2<f32>(0.0, 0.0);
|
|
||||||
if DE || periodic {
|
|
||||||
fp = fprime(z);
|
|
||||||
}
|
}
|
||||||
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).
|
// 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;
|
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;
|
m = m + 1u;
|
||||||
n = n + 1u;
|
n = n + 1u;
|
||||||
|
}
|
||||||
|
|
||||||
// Keep the reference index valid and the delta small.
|
// 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
|
// Reference exhausted: any pixel that followed it this far has
|
||||||
// effectively escaped (interior pixels rebase before reaching here).
|
// effectively escaped (interior pixels rebase before reaching here).
|
||||||
z = xm + e;
|
z = ref_orbit[u.ref_len - 1u] + e;
|
||||||
escaped = true;
|
escaped = true;
|
||||||
break;
|
break;
|
||||||
}
|
}
|
||||||
xm = ref_orbit[m];
|
let y = ref_orbit[m] + e;
|
||||||
z = xm + e;
|
if (dot(y, y) < dot(e, e)) {
|
||||||
z2 = dot(z, z);
|
|
||||||
if z2 < dot(e, e) {
|
|
||||||
// Rebase to index 0: carry the full value as the new delta. Valid
|
// 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
|
// because y_n = X[0] + (y_n - X[0]); for Mandelbrot X[0]=0.
|
||||||
// full value `z` (and `z2`) is unchanged by the re-expression.
|
e = y - z0;
|
||||||
// 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;
|
|
||||||
m = 0u;
|
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 {
|
if (!escaped) {
|
||||||
return Sample(0.0, 1.0, false); // interior of the set
|
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
|
let z2 = dot(z, z);
|
||||||
// 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);
|
|
||||||
|
|
||||||
// Continuous (smooth) iteration count.
|
// Continuous (smooth) iteration count.
|
||||||
let log_zn = 0.5 * log(max(z2, 1.0));
|
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;
|
let smooth_i = f32(n) + 1.0 - nu;
|
||||||
|
|
||||||
// sqrt compresses the huge iteration counts of deep zooms so the palette
|
// sqrt compresses the huge iteration counts of deep zooms so the palette
|
||||||
// varies smoothly instead of aliasing into speckle.
|
// varies smoothly instead of aliasing into speckle.
|
||||||
let ci = sqrt(max(smooth_i, 0.0));
|
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 (u.de_coloring != 0u) {
|
||||||
if DE {
|
// Exterior distance estimate (complex-plane units): |z|·ln|z| / |dz|.
|
||||||
// Exterior distance estimate |z|·ln|z| / |dz|, already in pixels since
|
// Divided by the pixel footprint it becomes a distance in pixels; we
|
||||||
// `dzs` = px·dz. We darken toward the boundary (< ~1 px away) so
|
// darken toward the boundary (< ~1 px away) so filaments stay crisp
|
||||||
// filaments stay crisp instead of aliasing into speckle. If |dzs|
|
// instead of aliasing into speckle. If |dz| overflowed, de -> 0 and the
|
||||||
// overflowed (far sub-pixel from the set), de -> 0 and the boundary
|
// boundary simply reads as dark, which is the correct limit.
|
||||||
// simply reads as dark, which is the correct limit.
|
|
||||||
let zmag = sqrt(max(z2, 1.0));
|
let zmag = sqrt(max(z2, 1.0));
|
||||||
let dzmag = sqrt(max(dot(dzs_esc, dzs_esc), 1e-30));
|
let dzmag = sqrt(max(dot(dz, dz), 1e-20));
|
||||||
let max_de = select(1.0, 1000.0, u.shadow != 0u);
|
let de = zmag * log(zmag) / dzmag;
|
||||||
de = clamp(zmag * log(zmag) / dzmag, 0.0, max_de);
|
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
|
@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;
|
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 dx = dpdx(base);
|
||||||
let dy = dpdy(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);
|
let aa = max(u.aa_level, 1u);
|
||||||
|
if (aa <= 1u) {
|
||||||
if u.shadow != 0u {
|
return vec4<f32>(shade(base, px), 1.0);
|
||||||
// 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);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
let inv = 1.0 / f32(aa);
|
|
||||||
var acc = vec3<f32>(0.0, 0.0, 0.0);
|
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 sy: u32 = 0u; sy < aa; sy = sy + 1u) {
|
||||||
for (var sx: u32 = 0u; sx < aa; sx = sx + 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 jx = (f32(sx) + 0.5) * inv - 0.5;
|
||||||
let jy = (f32(sy) + 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);
|
return vec4<f32>(acc / f32(aa * aa), 1.0);
|
||||||
|
|||||||
+2
-148
@@ -58,11 +58,6 @@ impl ViewState {
|
|||||||
DEFAULT_HALF_HEIGHT / self.half_height
|
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.
|
/// Bits of precision the center currently needs for this zoom level.
|
||||||
pub fn precision_bits(&self) -> usize {
|
pub fn precision_bits(&self) -> usize {
|
||||||
precision_for(self.half_height)
|
precision_for(self.half_height)
|
||||||
@@ -87,7 +82,7 @@ impl ViewState {
|
|||||||
let bits = self.precision_bits();
|
let bits = self.precision_bits();
|
||||||
// Grab-and-drag: moving the mouse right shows content to the left.
|
// 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_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
|
/// 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.)
|
// off * cpp * (1 - factor). (Derivation: new_c = fixed + (c-fixed)*f.)
|
||||||
let k = cpp * (1.0 - factor);
|
let k = cpp * (1.0 - factor);
|
||||||
self.center_re = &self.center_re + &big_from_f64(off_x * k, bits);
|
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;
|
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())
|
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.
|
/// Render a `Big` as a decimal string with `sig_digits` significant digits.
|
||||||
pub fn big_to_decimal_str(x: &Big, sig_digits: usize) -> String {
|
pub fn big_to_decimal_str(x: &Big, sig_digits: usize) -> String {
|
||||||
let dec = x
|
let dec = x
|
||||||
@@ -234,64 +149,3 @@ pub fn big_from_f64(x: f64, bits: usize) -> Big {
|
|||||||
.with_precision(bits)
|
.with_precision(bits)
|
||||||
.value()
|
.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
@@ -9,7 +9,7 @@
|
|||||||
use std::sync::mpsc::{Receiver, Sender, TryRecvError, channel};
|
use std::sync::mpsc::{Receiver, Sender, TryRecvError, channel};
|
||||||
use std::thread;
|
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};
|
use crate::view::{Big, big_from_f64};
|
||||||
|
|
||||||
pub struct RefRequest {
|
pub struct RefRequest {
|
||||||
@@ -22,12 +22,10 @@ pub struct RefRequest {
|
|||||||
pub precision: usize,
|
pub precision: usize,
|
||||||
pub kind: FractalKind,
|
pub kind: FractalKind,
|
||||||
pub power: u32,
|
pub power: u32,
|
||||||
/// Distortion constant for the Phoenix map (ignored by other kinds).
|
/// Build the BLA iteration-skip table for this orbit (square map only).
|
||||||
pub phoenix_p: (f64, f64),
|
pub build_bla: bool,
|
||||||
/// Distortion constant for the Lambda map (ignored by other kinds).
|
/// Upper bound on any pixel's `|dc|`, sizing the BLA merge radii.
|
||||||
pub lambda_l: (f64, f64),
|
pub dc_max: f64,
|
||||||
/// Complex exponent for the Complex Multibrot kind (ignored by other kinds).
|
|
||||||
pub complex_power: (f64, f64),
|
|
||||||
}
|
}
|
||||||
|
|
||||||
pub struct RefResult {
|
pub struct RefResult {
|
||||||
@@ -35,6 +33,7 @@ pub struct RefResult {
|
|||||||
pub center_im: Big,
|
pub center_im: Big,
|
||||||
pub half_height: f64,
|
pub half_height: f64,
|
||||||
pub points: Vec<[f32; 2]>,
|
pub points: Vec<[f32; 2]>,
|
||||||
|
pub bla: Vec<Bla>,
|
||||||
}
|
}
|
||||||
|
|
||||||
pub struct RefWorker {
|
pub struct RefWorker {
|
||||||
@@ -62,12 +61,18 @@ impl RefWorker {
|
|||||||
}
|
}
|
||||||
|
|
||||||
let points = compute(&req);
|
let points = compute(&req);
|
||||||
|
let bla = if req.build_bla {
|
||||||
|
build_bla_table(&points, req.dc_max)
|
||||||
|
} else {
|
||||||
|
Vec::new()
|
||||||
|
};
|
||||||
if res_tx
|
if res_tx
|
||||||
.send(RefResult {
|
.send(RefResult {
|
||||||
center_re: req.center_re,
|
center_re: req.center_re,
|
||||||
center_im: req.center_im,
|
center_im: req.center_im,
|
||||||
half_height: req.half_height,
|
half_height: req.half_height,
|
||||||
points,
|
points,
|
||||||
|
bla,
|
||||||
})
|
})
|
||||||
.is_err()
|
.is_err()
|
||||||
{
|
{
|
||||||
@@ -107,9 +112,6 @@ fn compute(req: &RefRequest) -> Vec<[f32; 2]> {
|
|||||||
req.precision,
|
req.precision,
|
||||||
req.kind,
|
req.kind,
|
||||||
req.power,
|
req.power,
|
||||||
req.phoenix_p,
|
|
||||||
req.lambda_l,
|
|
||||||
req.complex_power,
|
|
||||||
)
|
)
|
||||||
} else {
|
} else {
|
||||||
compute_set_reference(
|
compute_set_reference(
|
||||||
@@ -119,9 +121,6 @@ fn compute(req: &RefRequest) -> Vec<[f32; 2]> {
|
|||||||
req.precision,
|
req.precision,
|
||||||
req.kind,
|
req.kind,
|
||||||
req.power,
|
req.power,
|
||||||
req.phoenix_p,
|
|
||||||
req.lambda_l,
|
|
||||||
req.complex_power,
|
|
||||||
)
|
)
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|||||||
+7
-126
@@ -3,7 +3,7 @@
|
|||||||
//! shader with the same `naga` version wgpu uses — catching shader errors
|
//! shader with the same `naga` version wgpu uses — catching shader errors
|
||||||
//! without needing a GPU or a display.
|
//! 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) {
|
let module = match naga::front::wgsl::parse_str(src) {
|
||||||
Ok(m) => m,
|
Ok(m) => m,
|
||||||
Err(e) => panic!("{name}: WGSL parse error:\n{}", e.emit_to_string(src)),
|
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::ValidationFlags::all(),
|
||||||
naga::valid::Capabilities::all(),
|
naga::valid::Capabilities::all(),
|
||||||
);
|
);
|
||||||
match validator.validate(&module) {
|
if let Err(e) = validator.validate(&module) {
|
||||||
Ok(info) => (module, info),
|
panic!("{name}: WGSL validation error:\n{}", e.emit_to_string(src));
|
||||||
Err(e) => 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]
|
#[test]
|
||||||
fn colorize_shader_is_valid() {
|
fn mandelbrot_shader_is_valid() {
|
||||||
validate(
|
validate(
|
||||||
"colorize.wgsl",
|
"mandelbrot.wgsl",
|
||||||
concat!(
|
include_str!("../src/shaders/mandelbrot.wgsl"),
|
||||||
include_str!("../src/shaders/common.wgsl"),
|
|
||||||
include_str!("../src/shaders/iterate_uniforms.wgsl"),
|
|
||||||
include_str!("../src/shaders/colorize.wgsl"),
|
|
||||||
),
|
|
||||||
);
|
);
|
||||||
}
|
}
|
||||||
|
|
||||||
#[test]
|
#[test]
|
||||||
fn blit_shader_is_valid() {
|
fn blit_shader_is_valid() {
|
||||||
validate(
|
validate("blit.wgsl", include_str!("../src/shaders/blit.wgsl"));
|
||||||
"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)],
|
|
||||||
);
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
|||||||
Reference in New Issue
Block a user