From aac44ceb634dd2fd421631189291bc770cef8c53 Mon Sep 17 00:00:00 2001 From: Jozkah <38704522+Jozkah@users.noreply.github.com> Date: Sun, 1 Nov 2026 17:42:26 +0000 Subject: [PATCH] Deviation: nearest-face search on the GPU where WebGPU is available Scans of 100k+ vertices look up each vertex's nearest reference face in a WebGPU compute shader walking the reference's own BVH node buffer. The reading itself (distance, pseudonormal sign, facing limit, direction) is still resolved on the CPU in doubles on that face, so the map matches the CPU search. No WebGPU, an oversized reference or any GPU error falls back to the CPU BVH search unchanged. --- src/core/deviation/deviation.ts | 18 ++- src/core/deviation/gpuNearest.ts | 247 +++++++++++++++++++++++++++++++ src/core/deviation/surface.ts | 18 +++ src/core/meshWorker.ts | 77 ++++++---- tests/gpuDeviation.test.ts | 140 ++++++++++++++++++ 5 files changed, 467 insertions(+), 33 deletions(-) create mode 100644 src/core/deviation/gpuNearest.ts create mode 100644 tests/gpuDeviation.test.ts diff --git a/src/core/deviation/deviation.ts b/src/core/deviation/deviation.ts index 01b5d863..587d57f1 100644 --- a/src/core/deviation/deviation.ts +++ b/src/core/deviation/deviation.ts @@ -75,6 +75,11 @@ export interface DeviationOptions { * may be from facing the way the scan faces there. Null, or no `normals`, * takes the nearest surface whatever it faces. */ maxNormalDeviation?: number | null + /** The nearest reference face for each vertex, found ahead of time — by + * the GPU search in gpuNearest.ts — or -1 where it is to be searched for + * here. Only the face is taken from it: the reading itself is resolved on + * that face in doubles, so it is the one the tree search gives. */ + nearestFaces?: Int32Array | null } /** @@ -103,7 +108,7 @@ export function computeDeviation( surface: NominalSurface, scanPositions: Float32Array, transform: Rigid, - { onProgress, directions, normals, maxNormalDeviation = null }: DeviationOptions = {}, + { onProgress, directions, normals, maxNormalDeviation = null, nearestFaces = null }: DeviationOptions = {}, ): Float32Array { const n = scanPositions.length / 3 const values = new Float32Array(n) @@ -119,7 +124,9 @@ export function computeDeviation( for (let v = 0; v < n; v++) { rigidApply(transform, scanPositions[v * 3], scanPositions[v * 3 + 1], scanPositions[v * 3 + 2], p) let found: boolean - if (minFacing === null) found = surface.closest(p[0], p[1], p[2], hit) + const face = nearestFaces ? nearestFaces[v] : -1 + if (face >= 0) surface.closestOnFace(p[0], p[1], p[2], face, hit) + if (minFacing === null) found = face >= 0 || surface.closest(p[0], p[1], p[2], hit) else { // The scan's own normal, carried into the reference's frame. rigidRotate( @@ -129,8 +136,11 @@ export function computeDeviation( outward * normals![v * 3 + 2], sn, ) - const read = closestFacingWay(surface, p[0], p[1], p[2], sn[0], sn[1], sn[2], minFacing, hit, Infinity, reach) - found = read === 'nearest' || read === 'refaced' + if (face >= 0 && sn[0] * hit.nx + sn[1] * hit.ny + sn[2] * hit.nz >= minFacing) found = true + else { + const read = closestFacingWay(surface, p[0], p[1], p[2], sn[0], sn[1], sn[2], minFacing, hit, Infinity, reach) + found = read === 'nearest' || read === 'refaced' + } } if (!found) { values[v] = NaN diff --git a/src/core/deviation/gpuNearest.ts b/src/core/deviation/gpuNearest.ts new file mode 100644 index 00000000..edd652f3 --- /dev/null +++ b/src/core/deviation/gpuNearest.ts @@ -0,0 +1,247 @@ +// SPDX-License-Identifier: AGPL-3.0-only +// The nearest-face search of the deviation map on the GPU, where WebGPU is +// there to run it. Every scan vertex walks the reference's own tree — the +// node buffer three-mesh-bvh built, uploaded as it stands — in a compute +// shader, and comes back with the face it lies nearest. Only the face: the +// reading itself, distance and sign, is still resolved on the CPU in doubles +// (NominalSurface.closestOnFace), so the map is the one the CPU search gives, +// short of a near tie between two faces the shader's single precision splits +// the other way — and there the two readings agree to that precision anyway. +// Anything missing or failing answers null, and the CPU search runs instead. +// +// Timed in Chrome on a desktop GPU against an icosphere reference, the scan a +// noisy sphere lying on it: 655k vertices on 328k faces, 990 ms on the CPU +// tree and 71 ms here, search and readings together; 2.6M on 1.3M, 5.2 s and +// 0.28 s. A scan standing well off the reference, where the CPU walk opens +// far more of the tree, went from 24 s to 0.11 s. The readings agreed to +// 3e-6 mm throughout. + +import { rigidApply, type Rigid } from './rigid' +import type { NominalSurface } from './surface' + +/** Below this many scan vertices the CPU tree finishes before a GPU is even + * set up, and the map is not worth sending out for. */ +export const GPU_MIN_VERTICES = 100_000 + +/** Vertices per dispatch: each is awaited before the next goes out, so no + * one submission holds the GPU long enough for the driver to reset it, and + * progress has somewhere to come from. */ +const CHUNK = 1 << 18 + +/** A face index no search found: the stack ran out, or nothing was there. */ +const NO_FACE = 0xffffffff + +const SHADER = /* wgsl */ ` +struct Params { count: u32, pad0: u32, pad1: u32, pad2: u32 } +@group(0) @binding(0) var nodes: array; +@group(0) @binding(1) var positions: array; +@group(0) @binding(2) var tris: array; +@group(0) @binding(3) var points: array; +@group(0) @binding(4) var faces: array; +@group(0) @binding(5) var params: Params; + +fn boxDist2(n: u32, p: vec3f) -> f32 { + let lo = vec3f(bitcast(nodes[n]), bitcast(nodes[n + 1u]), bitcast(nodes[n + 2u])); + let hi = vec3f(bitcast(nodes[n + 3u]), bitcast(nodes[n + 4u]), bitcast(nodes[n + 5u])); + let d = max(max(lo - p, p - hi), vec3f(0.0)); + return dot(d, d); +} + +fn vertex(i: u32) -> vec3f { + return vec3f(positions[i * 3u], positions[i * 3u + 1u], positions[i * 3u + 2u]); +} + +// Ericson's closest point on a triangle, as squared distance. +fn triDist2(p: vec3f, a: vec3f, b: vec3f, c: vec3f) -> f32 { + let ab = b - a; let ac = c - a; let ap = p - a; + let d1 = dot(ab, ap); let d2 = dot(ac, ap); + if (d1 <= 0.0 && d2 <= 0.0) { return dot(ap, ap); } + let bp = p - b; + let d3 = dot(ab, bp); let d4 = dot(ac, bp); + if (d3 >= 0.0 && d4 <= d3) { return dot(bp, bp); } + let vc = d1 * d4 - d3 * d2; + if (vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0) { + let q = a + ab * (d1 / (d1 - d3)); let e = p - q; return dot(e, e); + } + let cp = p - c; + let d5 = dot(ab, cp); let d6 = dot(ac, cp); + if (d6 >= 0.0 && d5 <= d6) { return dot(cp, cp); } + let vb = d5 * d2 - d1 * d6; + if (vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0) { + let q = a + ac * (d2 / (d2 - d6)); let e = p - q; return dot(e, e); + } + let va = d3 * d6 - d5 * d4; + if (va <= 0.0 && d4 - d3 >= 0.0 && d5 - d6 >= 0.0) { + let q = b + (c - b) * ((d4 - d3) / ((d4 - d3) + (d5 - d6))); let e = p - q; return dot(e, e); + } + let denom = 1.0 / (va + vb + vc); + let q = a + ab * (vb * denom) + ac * (vc * denom); + let e = p - q; + return dot(e, e); +} + +@compute @workgroup_size(64) +fn main(@builtin(global_invocation_id) gid: vec3u) { + let i = gid.x; + if (i >= params.count) { return; } + let p = vec3f(points[i * 3u], points[i * 3u + 1u], points[i * 3u + 2u]); + var best = 3.4e38; + var bestFace = 0xffffffffu; + var stack: array; + stack[0] = 0u; + var sp = 1u; + loop { + if (sp == 0u) { break; } + sp -= 1u; + let n = stack[sp]; + if (boxDist2(n, p) >= best) { continue; } + let info = nodes[n + 7u]; + if ((info >> 16u) == 0xffffu) { + let first = nodes[n + 6u]; + let end = first + (info & 0xffffu); + for (var t = first; t < end; t += 1u) { + let d = triDist2(p, vertex(tris[t * 3u]), vertex(tris[t * 3u + 1u]), vertex(tris[t * 3u + 2u])); + if (d < best) { best = d; bestFace = t; } + } + } else { + if (sp + 2u > 64u) { bestFace = 0xffffffffu; break; } + let l = n + 8u; + let r = nodes[n + 6u]; + let dl = boxDist2(l, p); + let dr = boxDist2(r, p); + // The nearer child goes on top, to be opened first. + if (dl < dr) { + if (dr < best) { stack[sp] = r; sp += 1u; } + if (dl < best) { stack[sp] = l; sp += 1u; } + } else { + if (dl < best) { stack[sp] = l; sp += 1u; } + if (dr < best) { stack[sp] = r; sp += 1u; } + } + } + } + faces[i] = bestFace; +} +` + +// The few WebGPU names used here, loosely: the DOM library this project +// builds against does not declare them. +type GPUDeviceLike = any + +const STORAGE = 0x80, COPY_DST = 0x08, COPY_SRC = 0x04, UNIFORM = 0x40, MAP_READ = 0x01 + +let device: Promise | null = null +let pipeline: any = null + +function deviceOf(): Promise { + if (device) return device + device = (async () => { + const gpu = (globalThis.navigator as any)?.gpu + if (!gpu) return null + const adapter = await gpu.requestAdapter({ powerPreference: 'high-performance' }) + if (!adapter) return null + const d = await adapter.requestDevice({ + requiredLimits: { + maxStorageBufferBindingSize: adapter.limits.maxStorageBufferBindingSize, + maxBufferSize: adapter.limits.maxBufferSize, + }, + }) + d.lost.then(() => { + device = null + pipeline = null + held = null + }) + pipeline = await d.createComputePipelineAsync({ + layout: 'auto', + compute: { module: d.createShaderModule({ code: SHADER }), entryPoint: 'main' }, + }) + return d + })().catch(() => null) + return device +} + +/** The reference uploaded, kept for the next map against the same one. */ +let held: { surface: NominalSurface; buffers: any[] } | null = null + +function upload(d: GPUDeviceLike, data: ArrayBufferView): any { + const size = Math.max(16, Math.ceil(data.byteLength / 4) * 4) + const buffer = d.createBuffer({ size, usage: STORAGE | COPY_DST }) + d.queue.writeBuffer(buffer, 0, data.buffer, data.byteOffset, data.byteLength) + return buffer +} + +function referenceBuffers(d: GPUDeviceLike, surface: NominalSurface): any[] | null { + if (held?.surface === surface) return held.buffers + const tree = surface.flatTree() + if (!tree) return null + const max = Math.min(d.limits.maxStorageBufferBindingSize, d.limits.maxBufferSize) + if (Math.max(tree.byteLength, surface.positions.byteLength, surface.index.byteLength) > max) return null + if (held) for (const b of held.buffers) b.destroy() + held = { surface, buffers: [upload(d, tree), upload(d, surface.positions), upload(d, surface.index)] } + return held.buffers +} + +/** + * The nearest reference face for every scan vertex, -1 where the shader could + * not say, found on the GPU; null when there is no GPU to ask or it failed, + * for the caller to search on the CPU instead. + */ +export async function gpuNearestFaces( + surface: NominalSurface, + scanPositions: Float32Array, + transform: Rigid, + onProgress?: (fraction: number) => void, +): Promise { + const d = await deviceOf() + if (!d) return null + try { + const reference = referenceBuffers(d, surface) + if (!reference) return null + const n = scanPositions.length / 3 + const faces = new Int32Array(n) + const span = Math.min(CHUNK, n) + const points = d.createBuffer({ size: span * 12, usage: STORAGE | COPY_DST }) + const out = d.createBuffer({ size: span * 4, usage: STORAGE | COPY_SRC }) + const readback = d.createBuffer({ size: span * 4, usage: MAP_READ | COPY_DST }) + const params = d.createBuffer({ size: 16, usage: UNIFORM | COPY_DST }) + const bind = d.createBindGroup({ + layout: pipeline.getBindGroupLayout(0), + entries: [...reference, points, out, params].map((buffer, binding) => ({ binding, resource: { buffer } })), + }) + const local = new Float32Array(span * 3) + const p = new Float64Array(3) + try { + for (let start = 0; start < n; start += span) { + const count = Math.min(span, n - start) + // Into the reference's frame here, where the transform is in doubles. + for (let i = 0; i < count; i++) { + const v = (start + i) * 3 + rigidApply(transform, scanPositions[v], scanPositions[v + 1], scanPositions[v + 2], p) + local[i * 3] = p[0]; local[i * 3 + 1] = p[1]; local[i * 3 + 2] = p[2] + } + d.queue.writeBuffer(points, 0, local.buffer, 0, count * 12) + d.queue.writeBuffer(params, 0, new Uint32Array([count, 0, 0, 0])) + d.pushErrorScope('validation') + const encoder = d.createCommandEncoder() + const pass = encoder.beginComputePass() + pass.setPipeline(pipeline) + pass.setBindGroup(0, bind) + pass.dispatchWorkgroups(Math.ceil(count / 64)) + pass.end() + encoder.copyBufferToBuffer(out, 0, readback, 0, count * 4) + d.queue.submit([encoder.finish()]) + const error = await d.popErrorScope() + if (error) return null + await readback.mapAsync(MAP_READ, 0, count * 4) + const got = new Uint32Array(readback.getMappedRange(0, count * 4)) + for (let i = 0; i < count; i++) faces[start + i] = got[i] === NO_FACE ? -1 : got[i] + readback.unmap() + onProgress?.((start + count) / n) + } + } finally { + for (const b of [points, out, readback, params]) b.destroy() + } + return faces + } catch { + return null + } +} diff --git a/src/core/deviation/surface.ts b/src/core/deviation/surface.ts index 1e536446..256ac272 100644 --- a/src/core/deviation/surface.ts +++ b/src/core/deviation/surface.ts @@ -321,6 +321,24 @@ export class NominalSurface { return true } + /** The signed reading off one triangle already known to be the nearest — + * the face a GPU search found (see gpuNearest.ts) — resolved here in + * doubles, with the same pseudonormal sign `closest` gives. */ + closestOnFace(x: number, y: number, z: number, face: number, out: ClosestHit): void { + this.fillHit(x, y, z, face, out) + } + + /** The tree flattened for a GPU search: three-mesh-bvh's own node buffer, + * 32 bytes a node — six float bounds, then for a leaf the first triangle + * and a 16-bit count flagged 0xFFFF, and for an inner node the word offset + * of its right child, the left following directly. Its triangle order is + * the reordered `index`, so a leaf's triangles are the face ids `closest` + * reports. Null for a tree built in more than one root. */ + flatTree(): Uint32Array | null { + const roots = (this.bvh as unknown as { _roots: ArrayBuffer[] })._roots + return roots.length === 1 ? new Uint32Array(roots[0]) : null + } + /** * The same query, restricted to surface that faces roughly the way the query * point does: the nearest triangle whose outward normal satisfies diff --git a/src/core/meshWorker.ts b/src/core/meshWorker.ts index 320732ec..49e6a8ee 100644 --- a/src/core/meshWorker.ts +++ b/src/core/meshWorker.ts @@ -30,6 +30,7 @@ import { } from './deviation/align' import type { Steps } from './deviation/steps' import { computeDeviation, defaultMaxDistance, suggestRange } from './deviation/deviation' +import { GPU_MIN_VERTICES, gpuNearestFaces } from './deviation/gpuNearest' import { rigidApplyToPoints, rigidInvert, rigidRotateVectors, type Rigid } from './deviation/rigid' import { meshArea, meshLongEdges, readReference, sampleSpacing, subdivideForMap } from './deviation/onReference' import { buildSolidIndex, computeThickness, suggestThicknessScale } from './thickness/thickness' @@ -887,37 +888,55 @@ function handle(msg: Exclude): void { }) return } - try { - let lastPercent = -1 - const directions = new Int8Array(graph.positions.length) - const values = computeDeviation(nominal, graph.positions, msg.transform, { - onProgress: (f) => { - const percent = Math.round(f * 100) - if (percent === lastPercent) return - lastPercent = percent - progress(`Measuring the scan against the reference — ${percent}%…`) - }, - directions, - // Rotated along with the positions whenever an alignment is baked in, - // so they stand in the same frame the transform starts from. - normals: graph.normals, - maxNormalDeviation: msg.facingDeg === null ? null : (msg.facingDeg * Math.PI) / 180, - }) - const suggestedMaxDistance = defaultMaxDistance(graph.bboxDiag) - post( - { - type: 'deviation-ok', - requestId: msg.requestId, - values, + const scan = graph, reference = nominal, transform = msg.transform, requestId = msg.requestId + let lastPercent = -1 + const report = (f: number) => { + const percent = Math.round(f * 100) + if (percent === lastPercent) return + lastPercent = percent + progress(`Measuring the scan against the reference — ${percent}%…`) + } + const deviate = (nearestFaces: Int32Array | null) => { + try { + const directions = new Int8Array(scan.positions.length) + const values = computeDeviation(reference, scan.positions, transform, { + // With the faces in hand the GPU had the first half of the bar. + onProgress: (f) => report(nearestFaces ? 0.5 + f * 0.5 : f), + nearestFaces, directions, - suggestedRange: suggestRange(values, suggestedMaxDistance), - suggestedMaxDistance, - }, - [values.buffer, directions.buffer], - ) - } catch (e) { - post({ type: 'error', requestId: msg.requestId, message: errorText(e) }) + // Rotated along with the positions whenever an alignment is baked in, + // so they stand in the same frame the transform starts from. + normals: scan.normals, + maxNormalDeviation: msg.facingDeg === null ? null : (msg.facingDeg * Math.PI) / 180, + }) + const suggestedMaxDistance = defaultMaxDistance(scan.bboxDiag) + post( + { + type: 'deviation-ok', + requestId, + values, + directions, + suggestedRange: suggestRange(values, suggestedMaxDistance), + suggestedMaxDistance, + }, + [values.buffer, directions.buffer], + ) + } catch (e) { + post({ type: 'error', requestId, message: errorText(e) }) + } } + // A big scan's nearest faces are looked up on the GPU first, where there + // is one, and the CPU tree searched as ever where there is not. The wait + // holds the queue, as an awaiting plugin's does. + if (scan.positions.length / 3 < GPU_MIN_VERTICES) return deviate(null) + holding = true + gpuNearestFaces(reference, scan.positions, transform, (f) => report(f * 0.5)) + .then(deviate) + .finally(() => { + holding = false + drainQueue() + }) + return } } diff --git a/tests/gpuDeviation.test.ts b/tests/gpuDeviation.test.ts new file mode 100644 index 00000000..123c653f --- /dev/null +++ b/tests/gpuDeviation.test.ts @@ -0,0 +1,140 @@ +// SPDX-License-Identifier: AGPL-3.0-only +import { describe, expect, it } from 'vitest' +import { buildMeshGraph } from '../src/core/geometry/buildGraph' +import { computeDeviation } from '../src/core/deviation/deviation' +import { gpuNearestFaces } from '../src/core/deviation/gpuNearest' +import { rigidApply, rigidFromAxisAngle } from '../src/core/deviation/rigid' +import { NominalSurface } from '../src/core/deviation/surface' +import { boxMesh, icosphere } from './helpers' + +/** + * The compute shader's search, line for line in JS and single precision: the + * same walk over the same flattened node buffer, so the layout it reads and + * the faces it settles on are checked where no GPU is to hand. + */ +function shaderNearestFaces(surface: NominalSurface, points: Float32Array): Int32Array { + const nodes = surface.flatTree()! + const f = new Float32Array(nodes.buffer, nodes.byteOffset, nodes.length) + const pos = surface.positions + const tris = surface.index + const n = points.length / 3 + const faces = new Int32Array(n) + const box = (node: number, x: number, y: number, z: number) => { + const dx = Math.max(f[node] - x, x - f[node + 3], 0) + const dy = Math.max(f[node + 1] - y, y - f[node + 4], 0) + const dz = Math.max(f[node + 2] - z, z - f[node + 5], 0) + return Math.fround(dx * dx + dy * dy + dz * dz) + } + const hit = new Float64Array(3) + for (let i = 0; i < n; i++) { + const x = points[i * 3], y = points[i * 3 + 1], z = points[i * 3 + 2] + let best = Infinity + let bestFace = -1 + const stack = [0] + while (stack.length) { + const node = stack.pop()! + if (box(node, x, y, z) >= best) continue + const meta = nodes[node + 7] + if (meta >>> 16 === 0xffff) { + for (let t = nodes[node + 6], end = t + (meta & 0xffff); t < end; t++) { + const a = tris[t * 3], b = tris[t * 3 + 1], c = tris[t * 3 + 2] + triClosest(x, y, z, pos, a, b, c, hit) + const d = Math.fround((hit[0] - x) ** 2 + (hit[1] - y) ** 2 + (hit[2] - z) ** 2) + if (d < best) { best = d; bestFace = t } + } + } else { + const l = node + 8, r = nodes[node + 6] + const dl = box(l, x, y, z), dr = box(r, x, y, z) + const [far, dFar, near, dNear] = dl < dr ? [r, dr, l, dl] : [l, dl, r, dr] + if (dFar < best) stack.push(far) + if (dNear < best) stack.push(near) + } + } + faces[i] = bestFace + } + return faces +} + +/** Closest point on a triangle, by brute force over a fine barycentric + * projection — plain enough to trust as the yardstick here. */ +function triClosest( + x: number, y: number, z: number, + p: Float32Array, a: number, b: number, c: number, + out: Float64Array, +): void { + const ax = p[a * 3], ay = p[a * 3 + 1], az = p[a * 3 + 2] + const abx = p[b * 3] - ax, aby = p[b * 3 + 1] - ay, abz = p[b * 3 + 2] - az + const acx = p[c * 3] - ax, acy = p[c * 3 + 1] - ay, acz = p[c * 3 + 2] - az + const apx = x - ax, apy = y - ay, apz = z - az + // Solve for the plane's barycentrics, then clamp onto the triangle by + // projecting onto whichever edge is nearest when outside. + const d00 = abx * abx + aby * aby + abz * abz + const d01 = abx * acx + aby * acy + abz * acz + const d11 = acx * acx + acy * acy + acz * acz + const d20 = apx * abx + apy * aby + apz * abz + const d21 = apx * acx + apy * acy + apz * acz + const den = d00 * d11 - d01 * d01 + const v = (d11 * d20 - d01 * d21) / den + const w = (d00 * d21 - d01 * d20) / den + if (v >= 0 && w >= 0 && v + w <= 1) { + out[0] = ax + abx * v + acx * w; out[1] = ay + aby * v + acy * w; out[2] = az + abz * v + acz * w + return + } + let best = Infinity + const seg = (sx: number, sy: number, sz: number, ex: number, ey: number, ez: number) => { + const dx = ex - sx, dy = ey - sy, dz = ez - sz + const t = Math.min(1, Math.max(0, ((x - sx) * dx + (y - sy) * dy + (z - sz) * dz) / (dx * dx + dy * dy + dz * dz))) + const qx = sx + dx * t, qy = sy + dy * t, qz = sz + dz * t + const d = (qx - x) ** 2 + (qy - y) ** 2 + (qz - z) ** 2 + if (d < best) { best = d; out[0] = qx; out[1] = qy; out[2] = qz } + } + seg(ax, ay, az, ax + abx, ay + aby, az + abz) + seg(ax, ay, az, ax + acx, ay + acy, az + acz) + seg(ax + abx, ay + aby, az + abz, ax + acx, ay + acy, az + acz) +} + +describe('the GPU deviation search', () => { + const sphere = icosphere(4, 25) + const reference = new NominalSurface(sphere.positions, sphere.indices) + const scan = buildMeshGraph({ kind: 'soup', positions: boxMesh(30, 24, 0.4) }) + const transform = rigidFromAxisAngle([1, 2, 3], 0.4) + const local = new Float32Array(scan.positions.length) + const p = new Float64Array(3) + for (let v = 0; v < local.length; v += 3) { + rigidApply(transform, scan.positions[v], scan.positions[v + 1], scan.positions[v + 2], p) + local.set(p, v) + } + const faces = shaderNearestFaces(reference, local) + + it('finds a face for every point', () => { + expect(reference.flatTree()).not.toBeNull() + expect(faces.every((face) => face >= 0)).toBe(true) + }) + + it('gives the CPU map with the faces it found', () => { + const cpu = computeDeviation(reference, scan.positions, transform) + const gpu = computeDeviation(reference, scan.positions, transform, { nearestFaces: faces }) + for (let i = 0; i < cpu.length; i++) expect(Math.abs(gpu[i] - cpu[i])).toBeLessThan(1e-4) + }) + + it('gives the CPU map under a facing limit too', () => { + const options = { normals: scan.normals, maxNormalDeviation: Math.PI / 2 } + const cpu = computeDeviation(reference, scan.positions, transform, options) + const gpu = computeDeviation(reference, scan.positions, transform, { ...options, nearestFaces: faces }) + for (let i = 0; i < cpu.length; i++) { + if (Number.isNaN(cpu[i])) expect(gpu[i]).toBeNaN() + else expect(Math.abs(gpu[i] - cpu[i])).toBeLessThan(1e-4) + } + }) + + it('searches the rest where a face is missing', () => { + const cpu = computeDeviation(reference, scan.positions, transform) + const some = faces.map((face, i) => (i % 3 === 0 ? -1 : face)) + const mixed = computeDeviation(reference, scan.positions, transform, { nearestFaces: some }) + for (let i = 0; i < cpu.length; i++) expect(Math.abs(mixed[i] - cpu[i])).toBeLessThan(1e-4) + }) + + it('answers null with no WebGPU, for the CPU search to run', async () => { + expect(await gpuNearestFaces(reference, scan.positions, transform)).toBeNull() + }) +})