2 * The scene viewer: the instrument drawn from its own fields.
3 *
4 * A `visual` section names the passes that compute the medium (normally the
5 * same medium shader the precompute runs) and which buffers to draw. The
6 * viewer runs those passes at the current parameter values and raymarches
7 * the results on the GPU: a 'solid' volume (the body mask) is shaded as an
8 * opaque surface where it crosses one half, a 'cloud' volume (a coupling
9 * profile) accumulates as a translucent density. Because the fields are the
10 * medium shader's own output on the actual air grid, the picture is what
11 * the solver sees, grid effects included: coarsen nx and the sound hole
12 * staircases, thin a wall under two cells and it visibly stops sealing.
13 *
14 * Markers (an object's outline, the microphone) are the one declarative
15 * part: crisp annotations whose coordinates are expressions over the scope,
16 * drawn on a 2D overlay with the same camera.
17 */
18import { evalExpr } from '../expr.ts';
19import type { InstrumentFile, VisualSpec } from '../format.ts';
21export interface Camera {
22 /** Orbit angles in radians and distance in metres, about the origin. */
23 yaw: number;
24 pitch: number;
25 dist: number;
26}
28const WG = 64;
30interface Vec3 {
31 x: number;
32 y: number;
33 z: number;
34}
36const basis = (cam: Camera): { eye: Vec3; right: Vec3; up: Vec3; fwd: Vec3 } => {
37 const cp = Math.cos(cam.pitch);
38 const eye = {
39 x: cam.dist * cp * Math.cos(cam.yaw),
40 y: cam.dist * cp * Math.sin(cam.yaw),
41 z: cam.dist * Math.sin(cam.pitch),
42 };
43 const norm = (v: Vec3): Vec3 => {
44 const l = Math.hypot(v.x, v.y, v.z) || 1;
45 return { x: v.x / l, y: v.y / l, z: v.z / l };
46 };
47 const fwd = norm({ x: -eye.x, y: -eye.y, z: -eye.z });
48 const right = norm({ x: -fwd.y, y: fwd.x, z: 0 }); // cross(fwd, z-up)
49 const up = {
50 x: right.y * fwd.z - right.z * fwd.y,
51 y: right.z * fwd.x - right.x * fwd.z,
52 z: right.x * fwd.y - right.y * fwd.x,
53 };
54 return { eye, right, up: { x: -up.x, y: -up.y, z: -up.z }, fwd };
55};
57const TAN_HALF = Math.tan((35 * Math.PI) / 180 / 2);
59/** The raymarch fragment shader, generated for this instrument's volumes.
60 * Exported so the test suite can compile-check it headless. */
61export function raymarchWgsl(visual: VisualSpec): string {
62 const n = visual.volumes.length;
63 const binds = visual.volumes
64 .map((_, i) => `@group(0) @binding(${i + 1}) var<storage, read> f${i}: array<f32>;`)
65 .join('\n');
66 const samplers = visual.volumes
67 .map(
68 (v, i) => `
69fn tri${i}(g: vec3f) -> f32 {
70 let gc = clamp(g, vec3f(0.0), vec3f(f32(NX) - 1.001, f32(NY) - 1.001, f32(NZ) - 1.001));
71 let i0 = vec3u(floor(gc));
72 let fr = gc - floor(gc);
73 let s = i0.x + NX * (i0.y + NY * i0.z);
74 let v000 = f${i}[s];
75 let v100 = f${i}[s + 1u];
76 let v010 = f${i}[s + NX];
77 let v110 = f${i}[s + NX + 1u];
78 let v001 = f${i}[s + NX * NY];
79 let v101 = f${i}[s + NX * NY + 1u];
80 let v011 = f${i}[s + NX * NY + NX];
81 let v111 = f${i}[s + NX * NY + NX + 1u];
82 let a = mix(mix(v000, v100, fr.x), mix(v010, v110, fr.x), fr.y);
83 let b = mix(mix(v001, v101, fr.x), mix(v011, v111, fr.x), fr.y);
84 let v = mix(a, b, fr.z);
85 return ${v.invert ? '1.0 - v' : 'v'};
86}`,
87 )
88 .join('\n');
89 const marchBody = visual.volumes
90 .map((v, i) => {
91 const col = `vec3f(${v.color.join(', ')})`;
92 if (v.style === 'solid') {
93 return `
94 {
95 let v = tri${i}(g);
96 // A mask that never saturates (a wall thinner than the grid can
97 // carry) draws as a haze rather than vanishing: seeing the leak is
98 // the point of drawing the real fields.
99 if (v <= 0.5) {
100 let a = clamp(v, 0.0, 0.5) * 0.25;
101 accum += trans * a * ${col};
102 trans *= 1.0 - a;
103 } else {
104 let e = 1.0;
105 let grad = vec3f(
106 tri${i}(g + vec3f(e, 0.0, 0.0)) - tri${i}(g - vec3f(e, 0.0, 0.0)),
107 tri${i}(g + vec3f(0.0, e, 0.0)) - tri${i}(g - vec3f(0.0, e, 0.0)),
108 tri${i}(g + vec3f(0.0, 0.0, e)) - tri${i}(g - vec3f(0.0, 0.0, e)));
109 var nrm = -grad;
110 if (dot(nrm, nrm) < 1e-12) { nrm = vec3f(0.0, 0.0, 1.0); }
111 nrm = normalize(nrm);
112 let light = normalize(vec3f(0.45, -0.6, 0.75));
113 let shade = 0.35 + 0.65 * max(dot(nrm, light), 0.0);
114 let depth = 1.0 - 0.3 * clamp((t - t0) / (t1 - t0 + 1e-6), 0.0, 1.0);
115 accum += trans * ${col} * shade * depth;
116 trans = 0.0;
117 break;
118 }
119 }`;
120 }
121 return `
122 {
123 let vn = tri${i}(g) * u[6][${i}];
124 let a = 1.0 - exp(-vn * 30.0 / u[4].w);
125 accum += trans * a * ${col};
126 trans *= 1.0 - a;
127 }`;
128 })
129 .join('\n');
131 return `// Raymarch the instrument's medium fields. Uniforms u:
132// 0: eye.xyz, tanHalf 1: right.xyz, aspect 2: up.xyz 3: fwd.xyz
133// 4: halfExtents.xyz, steps 5: (unused) 6: per-cloud value scales
134override NX: u32;
135override NY: u32;
136override NZ: u32;
137@group(0) @binding(0) var<uniform> u: array<vec4f, 7>;
138${binds}
139${samplers}
141struct VOut {
142 @builtin(position) pos: vec4f,
143 @location(0) ndc: vec2f,
144}
146@vertex
147fn vmain(@builtin(vertex_index) vi: u32) -> VOut {
148 let xy = vec2f(f32((vi << 1u) & 2u), f32(vi & 2u)) * 2.0 - 1.0;
149 return VOut(vec4f(xy, 0.0, 1.0), xy);
150}
152@fragment
153fn fmain(v: VOut) -> @location(0) vec4f {
154 let eye = u[0].xyz;
155 let dir = normalize(u[3].xyz + v.ndc.x * u[1].w * u[0].w * u[1].xyz + v.ndc.y * u[0].w * u[2].xyz);
156 let half = u[4].xyz;
157 // Slab intersection with the domain box.
158 let inv = 1.0 / dir;
159 let tA = (-half - eye) * inv;
160 let tB = (half - eye) * inv;
161 let tmin = min(tA, tB);
162 let tmax = max(tA, tB);
163 let t0 = max(max(tmin.x, tmin.y), max(tmin.z, 0.0));
164 let t1 = min(min(tmax.x, tmax.y), tmax.z);
165 let bg = mix(vec3f(0.10, 0.11, 0.13), vec3f(0.16, 0.17, 0.20), 0.5 * (v.ndc.y + 1.0));
166 if (t1 <= t0) { return vec4f(bg, 1.0); }
167 var accum = vec3f(0.0);
168 var trans = 1.0;
169 let steps = u[4].w;
170 let dt = (t1 - t0) / steps;
171 for (var k = 0.0; k < steps; k += 1.0) {
172 let t = t0 + (k + 0.5) * dt;
173 let p = eye + t * dir;
174 let g = (p + half) / (2.0 * half) * vec3f(f32(NX), f32(NY), f32(NZ)) - 0.5;
175${marchBody}
176 if (trans < 0.02) { break; }
177 }
178 return vec4f(accum + trans * bg, 1.0);
179}
180`;
181}
183export class Viewer {
184 #device: GPUDevice;
185 #inst: InstrumentFile;
186 #visual: VisualSpec;
187 #canvas: HTMLCanvasElement;
188 #overlay: HTMLCanvasElement;
189 #ctx: GPUCanvasContext;
190 #format: GPUTextureFormat;
191 #buffers = new Map<string, GPUBuffer>();
192 #uniform: GPUBuffer;
193 #pipeline: GPURenderPipeline | null = null;
194 #bind: GPUBindGroup | null = null;
195 #scope: Record<string, number> = {};
196 #scales: number[] = [1, 1, 1, 1];
197 #dims: [number, number, number] = [1, 1, 1];
198 #half: [number, number, number] = [0.5, 0.25, 0.25];
200 constructor(device: GPUDevice, inst: InstrumentFile, canvas: HTMLCanvasElement, overlay: HTMLCanvasElement) {
201 if (!inst.visual) throw new Error(`instrument '${inst.name}' has no visual section`);
202 this.#device = device;
203 this.#inst = inst;
204 this.#visual = inst.visual;
205 this.#canvas = canvas;
206 this.#overlay = overlay;
207 this.#format = navigator.gpu.getPreferredCanvasFormat();
208 const ctx = canvas.getContext('webgpu');
209 if (!ctx) throw new Error('no webgpu canvas context');
210 this.#ctx = ctx;
211 ctx.configure({ device, format: this.#format, alphaMode: 'opaque' });
212 this.#uniform = device.createBuffer({
213 label: 'view-uniforms',
214 size: 7 * 16,
215 usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST,
216 });
217 }
219 /** Run the visual passes at these values and refresh normalization. */
220 async update(scope: Record<string, number>): Promise<void> {
221 const device = this.#device;
222 const v = this.#visual;
223 this.#scope = scope;
224 this.#dims = v.dims.map((d) => evalExpr(d, scope)) as [number, number, number];
225 this.#half = v.extents.map((e) => evalExpr(e, scope) / 2) as [number, number, number];
227 for (const b of v.buffers) {
228 const bytes = 4 * Math.max(1, evalExpr(b.count, scope));
229 const have = this.#buffers.get(b.name);
230 if (!have || have.size !== bytes) {
231 have?.destroy();
232 this.#buffers.set(
233 b.name,
234 device.createBuffer({
235 label: `view-${b.name}`,
236 size: bytes,
237 usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC,
238 }),
239 );
240 }
241 }
243 device.pushErrorScope('validation');
244 const enc = device.createCommandEncoder({ label: 'view-update' });
245 const cp = enc.beginComputePass();
246 for (const pass of v.setup) {
247 const s = this.#inst.shaders[pass.shader];
248 const constants: Record<string, number> = {};
249 for (const c of s.constants) {
250 if (!(c in scope)) throw new Error(`visual pass '${pass.shader}': constant '${c}' not in scope`);
251 constants[c] = scope[c];
252 }
253 const pipe = await device.createComputePipelineAsync({
254 label: `view-${pass.shader}`,
255 layout: 'auto',
256 compute: {
257 module: device.createShaderModule({ code: s.code, label: `view-${pass.shader}` }),
258 entryPoint: 'main',
259 constants,
260 },
261 });
262 cp.setPipeline(pipe);
263 cp.setBindGroup(
264 0,
265 device.createBindGroup({
266 layout: pipe.getBindGroupLayout(0),
267 entries: pass.buffers.map((name, i) => {
268 const buf = this.#buffers.get(name);
269 if (!buf) throw new Error(`visual pass binds undeclared buffer '${name}'`);
270 return { binding: i, resource: { buffer: buf } };
271 }),
272 }),
273 );
274 const wg = Math.ceil(evalExpr(pass.dispatch, scope) / WG);
275 cp.dispatchWorkgroups(Math.min(wg, 65535), Math.ceil(wg / 65535));
276 }
277 cp.end();
278 device.queue.submit([enc.finish()]);
280 // Normalize each cloud volume to its own peak.
281 for (let i = 0; i < v.volumes.length; i++) {
282 const vol = v.volumes[i];
283 if (vol.style !== 'cloud') {
284 this.#scales[i] = 1;
285 continue;
286 }
287 const src = this.#buffers.get(vol.buffer)!;
288 const rb = device.createBuffer({ size: src.size, usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST });
289 const e2 = device.createCommandEncoder();
290 e2.copyBufferToBuffer(src, 0, rb, 0, src.size);
291 device.queue.submit([e2.finish()]);
292 await rb.mapAsync(GPUMapMode.READ);
293 const data = new Float32Array(rb.getMappedRange());
294 let peak = 0;
295 for (const x of data) peak = Math.max(peak, Math.abs(x));
296 rb.unmap();
297 rb.destroy();
298 this.#scales[i] = peak > 0 ? 1 / peak : 1;
299 }
301 // The raymarch pipeline, rebuilt because the dims are baked in.
302 const module = device.createShaderModule({ code: raymarchWgsl(v), label: 'view-raymarch' });
303 this.#pipeline = await device.createRenderPipelineAsync({
304 label: 'view-raymarch',
305 layout: 'auto',
306 vertex: { module, entryPoint: 'vmain' },
307 fragment: {
308 module,
309 entryPoint: 'fmain',
310 targets: [{ format: this.#format }],
311 constants: { NX: this.#dims[0], NY: this.#dims[1], NZ: this.#dims[2] },
312 },
313 primitive: { topology: 'triangle-list' },
314 });
315 this.#bind = device.createBindGroup({
316 layout: this.#pipeline.getBindGroupLayout(0),
317 entries: [
318 { binding: 0, resource: { buffer: this.#uniform } },
319 ...v.volumes.map((vol, i) => ({ binding: i + 1, resource: { buffer: this.#buffers.get(vol.buffer)! } })),
320 ],
321 });
322 const err = await device.popErrorScope();
323 if (err) throw new Error(`viewer: ${err.message}`);
324 }
326 draw(cam: Camera): void {
327 if (!this.#pipeline || !this.#bind) return;
328 const device = this.#device;
329 const { eye, right, up, fwd } = basis(cam);
330 const w = this.#canvas.width;
331 const h = this.#canvas.height;
332 const steps = 2 * Math.max(...this.#dims);
333 const u = new Float32Array(7 * 4);
334 u.set([eye.x, eye.y, eye.z, TAN_HALF], 0);
335 u.set([right.x, right.y, right.z, w / h], 4);
336 u.set([up.x, up.y, up.z, 0], 8);
337 u.set([fwd.x, fwd.y, fwd.z, 0], 12);
338 u.set([...this.#half, steps], 16);
339 u.set(this.#scales, 24);
340 device.queue.writeBuffer(this.#uniform, 0, u);
341 const enc = device.createCommandEncoder({ label: 'view-draw' });
342 const rp = enc.beginRenderPass({
343 colorAttachments: [
344 {
345 view: this.#ctx.getCurrentTexture().createView(),
346 loadOp: 'clear',
347 storeOp: 'store',
348 clearValue: { r: 0, g: 0, b: 0, a: 1 },
349 },
350 ],
351 });
352 rp.setPipeline(this.#pipeline);
353 rp.setBindGroup(0, this.#bind);
354 rp.draw(3);
355 rp.end();
356 device.queue.submit([enc.finish()]);
357 this.#drawOverlay(cam);
358 }
360 /** Project a world point to overlay pixels; null when behind the eye. */
361 #project(cam: ReturnType<typeof basis>, p: [number, number, number]): [number, number] | null {
362 const d = { x: p[0] - cam.eye.x, y: p[1] - cam.eye.y, z: p[2] - cam.eye.z };
363 const dot = (a: Vec3): number => a.x * d.x + a.y * d.y + a.z * d.z;
364 const z = dot(cam.fwd);
365 if (z <= 1e-4) return null;
366 const w = this.#overlay.width;
367 const h = this.#overlay.height;
368 const x = dot(cam.right) / (z * TAN_HALF * (w / h));
369 const y = dot(cam.up) / (z * TAN_HALF);
370 return [(x * 0.5 + 0.5) * w, (1 - (y * 0.5 + 0.5)) * h];
371 }
373 #drawOverlay(cam: Camera): void {
374 const b = basis(cam);
375 const g = this.#overlay.getContext('2d');
376 if (!g) return;
377 const w = this.#overlay.width;
378 const h = this.#overlay.height;
379 g.clearRect(0, 0, w, h);
380 const seg = (p: [number, number, number], q: [number, number, number]): void => {
381 const a = this.#project(b, p);
382 const c = this.#project(b, q);
383 if (!a || !c) return;
384 g.beginPath();
385 g.moveTo(a[0], a[1]);
386 g.lineTo(c[0], c[1]);
387 g.stroke();
388 };
389 // The domain box, faint.
390 const [hx, hy, hz] = this.#half;
391 g.strokeStyle = 'rgba(255,255,255,0.18)';
392 g.lineWidth = 1;
393 for (const [i, j] of [[0, 1], [0, 2], [1, 3], [2, 3]] as const) {
394 const c = (k: number): [number, number] => [k & 1 ? hx : -hx, k & 2 ? hy : -hy];
395 const [xi, yi] = c(i);
396 const [xj, yj] = c(j);
397 seg([xi, yi, -hz], [xj, yj, -hz]);
398 seg([xi, yi, hz], [xj, yj, hz]);
399 seg([xi, yi, -hz], [xi, yi, hz]);
400 }
401 const ev = (e: [string, string, string]): [number, number, number] =>
402 e.map((x) => evalExpr(x, this.#scope)) as [number, number, number];
403 for (const m of this.#visual.markers) {
404 const col = `rgb(${m.color.map((c) => Math.round(255 * c)).join(',')})`;
405 g.strokeStyle = col;
406 g.fillStyle = col;
407 g.lineWidth = 1.5;
408 let anchor: [number, number] | null = null;
409 if (m.kind === 'point') {
410 const p = this.#project(b, ev(m.at));
411 if (p) {
412 g.beginPath();
413 g.arc(p[0], p[1], 4, 0, 2 * Math.PI);
414 g.fill();
415 anchor = p;
416 }
417 } else if (m.kind === 'line') {
418 const from = ev(m.from);
419 const to = ev(m.to);
420 seg(from, to);
421 anchor = this.#project(b, to);
422 } else {
423 const c0 = ev(m.corner);
424 const du = ev(m.du);
425 const dv = ev(m.dv);
426 const add = (a: [number, number, number], ...rest: [number, number, number][]): [number, number, number] =>
427 rest.reduce((s, r) => [s[0] + r[0], s[1] + r[1], s[2] + r[2]] as [number, number, number], a);
428 const c1 = add(c0, du);
429 const c2 = add(c0, du, dv);
430 const c3 = add(c0, dv);
431 seg(c0, c1);
432 seg(c1, c2);
433 seg(c2, c3);
434 seg(c3, c0);
435 anchor = this.#project(b, c2);
436 }
437 if (m.label && anchor) {
438 g.font = '11px system-ui, sans-serif';
439 g.fillText(m.label, anchor[0] + 6, anchor[1] - 5);
440 }
441 }
442 }
444 destroy(): void {
445 for (const b of this.#buffers.values()) b.destroy();
446 this.#buffers.clear();
447 this.#uniform.destroy();
448 }
449}
451/** Predicted partial frequencies read off a probed slots buffer: one entry
452 * per underdamped root, in Hz with a relative weight, deduplicated. */
453export function slotPartials(
454 slots: Float32Array,
455 layout: { stride: number; zRe: number; zIm: number },
456 dt: number,
457): { freq: number; weight: number }[] {
458 const byFreq = new Map<number, number>();
459 for (let s = 0; s + layout.stride <= slots.length; s += layout.stride) {
460 const zRe = slots[s + layout.zRe];
461 const zIm = slots[s + layout.zIm];
462 const cAbs = Math.hypot(slots[s], slots[s + 1]);
463 if (zIm <= 0 || cAbs === 0) continue;
464 const w = Math.abs(Math.atan2(zIm, zRe)) / dt;
465 const freq = w / (2 * Math.PI);
466 // Weight by rough displacement (the coefficients are acceleration
467 // series, which flatter high modes); this only sets tick opacity.
468 const key = Math.round(freq * 2) / 2;
469 byFreq.set(key, (byFreq.get(key) ?? 0) + cAbs / Math.max(w * w, 1));
470 }
471 const out = [...byFreq.entries()].map(([freq, weight]) => ({ freq, weight }));
472 const max = out.reduce((a, p) => Math.max(a, p.weight), 0);
473 return out
474 .filter((p) => p.weight > 1e-4 * max)
475 .map((p) => ({ freq: p.freq, weight: p.weight / max }))
476 .sort((a, b) => a.freq - b.freq);
477}