ponder.ooo
馃悕馃悕馃悕
1
2struct Uniforms /* buffer 0 0 */ {
3 extent: vec2u,
4 center: vec2f,
5 zoom: f32,
6 c: vec2f, // 0 to 1 = 0.85,0
7 iterations: u32, // hard 0 to 500 = 100
8 twist: f32, // -pi to pi = 0
9 phi: f32, // 0 to tau = 0
10 psi: f32, // 0 to 12 = 1
11 squoosh: vec2f, // 0 to 10 = 1,1
12 escape_distance: f32 // 0 to 10 = 2
13}
14
15@group(0) @binding(0) var<uniform> uniforms: Uniforms;
16@group(0) @binding(1) var output_texture: texture_storage_2d<rgba8unorm, write>;
17
18fn complex_mag_sq(z: vec2<f32>) -> f32 {
19 return z.x * z.x + z.y * z.y;
20}
21
22fn complex_mag(z: vec2<f32>) -> f32 {
23 return sqrt(complex_mag_sq(z));
24}
25
26fn complex_angle(z: vec2<f32>) -> f32 {
27 return atan2(z.y, z.x);
28}
29
30fn projective_shift(x: vec2<f32>, phi: f32, psi: f32) -> vec2<f32> {
31 let x_mag = complex_mag(x);
32 let x_angle = complex_angle(x);
33 let angle_diff = x_angle - phi;
34 let new_mag = x_mag + cos(angle_diff);
35 return vec2<f32>(new_mag * cos(x_angle * psi), new_mag * sin(x_angle * psi));
36}
37
38fn iterate_polar(
39 x: vec2<f32>,
40 phi: f32, psi: f32,
41 c: vec2f,
42 cosTwist: f32, sinTwist: f32,
43 squoosh: vec2f
44 ) -> vec2<f32> {
45 let shifted = projective_shift(x, phi, psi);
46 let halfEbb = -c * 0.5;
47 let mid = shifted + halfEbb;
48 let mids = vec2<f32>(mid.x / squoosh.x, mid.y / squoosh.y);
49 let rotated = vec2<f32>(
50 mids.x * cosTwist - mids.y * sinTwist,
51 mids.x * sinTwist + mids.y * cosTwist
52 );
53 return rotated + halfEbb;
54}
55
56fn iterate_cartesian(z: vec2<f32>, phi: f32, c: vec2f) -> vec2<f32> {
57 let phi_hat = vec2<f32>(cos(phi), sin(phi));
58 let z_mag_sq = complex_mag_sq(z);
59
60 // Avoid division by zero
61 if (z_mag_sq < 1e-10) {
62 return -c;
63 }
64
65 let z_dot_phi = z.x * phi_hat.x + z.y * phi_hat.y;
66 let projection_scale = z_dot_phi / z_mag_sq;
67
68 // z' = z + (z 路 蠁虃) / |z|虏 * z - c
69 let result = z + projection_scale * z;
70 return result - c;
71}
72
73fn pixel_to_complex(px: u32, py: u32) -> vec2<f32> {
74 let aspect = f32(uniforms.extent.x) / f32(uniforms.extent.y);
75 let scale = 4.0 / uniforms.zoom;
76 let half_width = f32(uniforms.extent.x) * 0.5;
77 let half_height = f32(uniforms.extent.y) * 0.5;
78 let x = (f32(px) - half_width) * scale * aspect / f32(uniforms.extent.x) + uniforms.center.x;
79 let y = (f32(py) - half_height) * scale / f32(uniforms.extent.y) + uniforms.center.y;
80 return vec2<f32>(x, y);
81}
82
83fn pixel_to_complex_alt(px: u32, py: u32) -> vec2<f32> {
84 let aspect = f32(uniforms.extent.x) / f32(uniforms.extent.y);
85 let scale = 4.0 / uniforms.zoom;
86 let half_width = f32(uniforms.extent.x) * 0.5;
87 let half_height = f32(uniforms.extent.y) * 0.5;
88
89 // map y-axis to magnitude, x-axis to angle
90 let angle = (f32(px) - half_width) * scale * aspect / f32(uniforms.extent.x) + uniforms.center.x;
91 let mag = (f32(py) - half_height) * scale / (2.0 * f32(uniforms.extent.y)) + uniforms.center.y;
92
93 return vec2<f32>(mag * cos(angle), mag * sin(angle));
94}
95
96fn pixel_to_complex_inverted(px: u32, py: u32) -> vec2<f32> {
97 let aspect = f32(uniforms.extent.x) / f32(uniforms.extent.y);
98 let scale = 4.0 / uniforms.zoom;
99 let half_width = f32(uniforms.extent.x) * 0.5;
100 let half_height = f32(uniforms.extent.y) * 0.5;
101
102 let x = (f32(px) - half_width) * scale * aspect / f32(uniforms.extent.x) + uniforms.center.x;
103 let y = (f32(py) - half_height) * scale / f32(uniforms.extent.y) + uniforms.center.y;
104
105 let z = vec2<f32>(x, y);
106 let mag_sq = x * x + y * y;
107
108 // handle singularity at origin
109 if (mag_sq < 1e-10) {
110 return vec2<f32>(1e10, 0.0); // or whatever you want infinity to map to
111 }
112
113 return vec2<f32>(x / mag_sq, -y / mag_sq);
114}
115
116fn pixel_to_complex_inverted_d(px: u32, py: u32) -> vec2<f32> {
117 let aspect = f32(uniforms.extent.x) / f32(uniforms.extent.y);
118 let scale = 4.0 / uniforms.zoom;
119 let half_width = f32(uniforms.extent.x) * 0.5;
120 let half_height = f32(uniforms.extent.y) * 0.5;
121
122 // get position relative to center
123 let x = (f32(px) - half_width) * scale * aspect / f32(uniforms.extent.x);
124 let y = (f32(py) - half_height) * scale / f32(uniforms.extent.y);
125
126 let mag_sq = x * x + y * y;
127 if (mag_sq < 1e-10) {
128 return vec2<f32>(uniforms.center.x + 1e10, uniforms.center.y);
129 }
130
131 // invert then add center back
132 return vec2<f32>(
133 uniforms.center.x + x / mag_sq,
134 uniforms.center.y - y / mag_sq
135 );
136}
137
138fn pixel_to_complex_inverted_b(px: u32, py: u32) -> vec2<f32> {
139 let aspect = f32(uniforms.extent.x) / f32(uniforms.extent.y);
140 let scale = 4.0 / uniforms.zoom;
141 let half_width = f32(uniforms.extent.x) * 0.5;
142 let half_height = f32(uniforms.extent.y) * 0.5;
143
144 let x = (f32(px) - half_width) * scale * aspect / f32(uniforms.extent.x) + uniforms.center.x;
145 let y = (f32(py) - half_height) * scale / f32(uniforms.extent.y) + uniforms.center.y;
146
147 // translate by (c/2, d/2)
148 let tx = x + uniforms.c.x * 0.5;
149 let ty = y + uniforms.c.y * 0.5;
150
151 let mag_sq = tx * tx + ty * ty;
152
153 if (mag_sq < 1e-10) {
154 return vec2<f32>(1e10, 0.0);
155 }
156
157 // invert, then translate back
158 return vec2<f32>(tx / mag_sq - uniforms.c.x * 0.5, -ty / mag_sq - uniforms.c.y * 0.5);
159}
160
161fn pixel_to_complex_inverted_c(px: u32, py: u32) -> vec2<f32> {
162 let aspect = f32(uniforms.extent.x) / f32(uniforms.extent.y);
163 let scale = 4.0 / uniforms.zoom;
164 let half_width = f32(uniforms.extent.x) * 0.5;
165 let half_height = f32(uniforms.extent.y) * 0.5;
166
167 let x = (f32(px) - half_width) * scale * aspect / f32(uniforms.extent.x) + uniforms.center.x;
168 let y = (f32(py) - half_height) * scale / f32(uniforms.extent.y) + uniforms.center.y;
169
170 // translate by (c/2, d/2), then scale down by 0.5
171 let tx = (x + uniforms.c.x * 0.5) * 2.0;
172 let ty = (y + uniforms.c.y * 0.5) * 2.0;
173
174 let mag_sq = tx * tx + ty * ty;
175
176 if (mag_sq < 1e-10) {
177 return vec2<f32>(1e10, 0.0);
178 }
179
180 // invert, scale up by 2x, then translate back
181 return vec2<f32>((tx / mag_sq) * 0.5 - uniforms.c.x * 0.5, (-ty / mag_sq) * 0.5 - uniforms.c.y * 0.5);
182}
183
184fn c_avg(prev: f32, val: f32, n: f32) -> f32 {
185 return prev + (val - prev) / n;
186}
187
188fn hash(seed: u32) -> u32 {
189 var x = seed;
190 x = ((x >> 16u) ^ x) * 0x45d9f3bu;
191 x = ((x >> 16u) ^ x) * 0x45d9f3bu;
192 x = (x >> 16u) ^ x;
193 return x;
194}
195
196fn pcg_hash(x: u32) -> u32 {
197 let state = x * 747796405u + 2891336453u;
198 let word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u;
199 return (word >> 22u) ^ word;
200}
201
202fn random_float(seed: u32) -> f32 {
203 return f32(hash(seed)) / 4294967296.0; // 2^32
204}
205
206@compute @workgroup_size(16, 16)
207fn main(@builtin(global_invocation_id) id: vec3<u32>) {
208 let px = id.x;
209 let py = id.y;
210
211 if (px >= uniforms.extent.x || py >= uniforms.extent.y) {
212 return;
213 }
214
215 var z = ${pixel_mapping}(px, py);
216 let orig_z = z;
217 let n_perturbations = 1u;
218 let escape_threshold = uniforms.escape_distance * uniforms.escape_distance;
219
220 var escaped = false;
221 var iter = 0u;
222 var cavg = f32(0.0);
223 let epsilon = 1.19e-6;
224
225 let transient_skip = u32(f32(uniforms.iterations) * f32(0.95));
226
227 let cosTwist = cos(uniforms.twist);
228 let sinTwist = sin(uniforms.twist);
229
230 for (iter = 0u; iter < uniforms.iterations; iter = iter + 1u) {
231 let mag_sq = complex_mag_sq(z);
232 if (mag_sq > escape_threshold) {
233 escaped = true;
234 //break;
235 }
236 var z_p: array<vec2<f32>, 10>;
237 var df_z_p: array<vec2<f32>, 10>;
238 var abs_df_z_p: array<f32, 10>;
239 var avg_abs_df_z_p = f32(0.0);
240
241 var f_z = iterate_polar(z, uniforms.phi, uniforms.psi, uniforms.c, cosTwist, sinTwist, uniforms.squoosh);
242
243 if (iter >= transient_skip) {
244 if (cavg < -10.0) { continue; }
245 if (cavg > 10.0) { continue; }
246 for (var p = 0u; p < n_perturbations; p = p + 1u) {
247 let rng_seed = pcg_hash(px + pcg_hash(py << 1)) ^ pcg_hash(p);
248 let random_angle = random_float(rng_seed) * 6.283185307;
249 let perturb_direction = vec2<f32>(cos(random_angle), sin(random_angle));
250 z_p[p] = z + perturb_direction * epsilon;
251 df_z_p[p] = iterate_polar(z_p[p], uniforms.phi, uniforms.psi, uniforms.c, cosTwist, sinTwist, uniforms.squoosh) - orig_z;
252 abs_df_z_p[p] = abs(complex_mag(df_z_p[p]));
253 avg_abs_df_z_p += abs_df_z_p[p];
254 }
255
256 avg_abs_df_z_p /= f32(n_perturbations);
257
258 cavg = c_avg(cavg, log(avg_abs_df_z_p), f32(iter + 1 - transient_skip));
259 }
260
261 z = f_z;
262 //z = avg_abs_df_z_p;
263 }
264
265 var color: vec4<f32>;
266 if (escaped) {
267 let escape_speed = 1.0 - pow(f32(iter) / f32(uniforms.iterations), 2.0);
268 color = vec4<f32>(escape_speed, escape_speed, escape_speed, 1.0);
269 } else {
270 let final_mag = complex_mag(z);
271 let final_angle = complex_angle(z);
272 //let scaled_mag = log(final_mag);
273 //let final_x = scaled_mag * (z.x / final_mag);
274 //let final_y = scaled_mag * (z.y / final_mag);
275 //color = vec4<f32>(0.5 + 0.5 * final_x, 0.0, 0.5 + 0.5 * final_y, 1.0);
276 let scaled_mag = pow(final_mag / escape_threshold, 0.5);
277 let scaled_angle = (final_angle + 3.1415926535) / 6.283185307179586;
278 let r = max(0, 2.0 * scaled_angle - 1.0);
279 let b = max(0, 2.0 * (1.0 - scaled_angle) - 1.0);
280 color = vec4<f32>(r, scaled_mag, b, 1.0);
281 }
282
283 //color = vec4<f32>(1.0 - cavg / 10.0, 1.0 - (-cavg / 5.0), 1.0 - (-cavg / 5.0), 1.0);
284
285 textureStore(output_texture, vec2<i32>(i32(px), i32(py)), color);
286}
287