From fb5edb43cf013c45b312ac53251b2ee990528896 Mon Sep 17 00:00:00 2001 From: Donovan Hutchence Date: Sat, 1 Aug 2026 00:14:43 +0100 Subject: [PATCH 1/5] latest --- src/lib/readers/read-lod.ts | 6 + src/lib/writers/write-lod.ts | 323 ++++++++++++++++++++++++++++++++++- test/write-lod.test.mjs | 99 +++++++++++ 3 files changed, 426 insertions(+), 2 deletions(-) diff --git a/src/lib/readers/read-lod.ts b/src/lib/readers/read-lod.ts index 2f5ce630..384f2309 100644 --- a/src/lib/readers/read-lod.ts +++ b/src/lib/readers/read-lod.ts @@ -13,6 +13,7 @@ type LodReference = { type LodNode = { children?: LodNode[]; lods?: Record; + errors?: number[]; }; type LodMeta = { @@ -65,6 +66,11 @@ const collectFilesByLod = (meta: LodMeta, filename: string): Map typeof error !== 'number' || !Number.isFinite(error) || error < 0))) { + throw new Error(`Invalid lod-meta.json errors: ${filename}`); + } for (const [key, ref] of Object.entries(node.lods ?? {})) { const lod = Number(key); if (!Number.isInteger(lod) || lod < 0 || lod >= meta.lodLevels || diff --git a/src/lib/writers/write-lod.ts b/src/lib/writers/write-lod.ts index 31cffcae..83ca1426 100644 --- a/src/lib/writers/write-lod.ts +++ b/src/lib/writers/write-lod.ts @@ -5,6 +5,8 @@ import { logWrittenFile } from './utils'; import { writeSogSource } from './write-sog.js'; import { type ChunkDataPool, type ChunkSource, type ReadRequest, createChunkDataPool } from '../chunk'; import { Column, DataTable } from '../data-table'; +import { det3, sigmoid, type SplatView } from '../decimate/moment-match'; +import { buildCostCache, type CostCache } from '../decimate-uniform/edge-cost-cpu'; import { type FileSystem } from '../io/write'; import { bakeTransform, permuteSource, sortMortonColumns } from '../ops'; import { BTreeNode, BTree } from '../spatial'; @@ -27,6 +29,7 @@ type MetaNode = { bound: Aabb; children?: MetaNode[]; lods?: { [key: number]: MetaLod }; + errors?: number[]; }; type LodMeta = { @@ -241,6 +244,321 @@ const binIndices = (parent: BTreeNode, lodOf: (index: number) => number): Map { + const { min, max } = node.aabb; + const dx = x < min[0] ? min[0] - x : x > max[0] ? x - max[0] : 0; + const dy = y < min[1] ? min[1] - y : y > max[1] ? y - max[1] : 0; + const dz = z < min[2] ? min[2] - z : z > max[2] ? z - max[2] : 0; + return dx * dx + dy * dy + dz * dz; +}; + +const findNearest = ( + root: BTreeNode, slim: SlimColumns, lodOf: (index: number) => number, + queries: Uint32Array, targetLod: number, k: number +): Int32Array => { + const result = new Int32Array(queries.length * k); + result.fill(-1); + const distances = new Float64Array(k); + + for (let q = 0; q < queries.length; q++) { + distances.fill(Infinity); + const g = queries[q]; + const x = slim.x[g], y = slim.y[g], z = slim.z[g]; + const output = q * k; + + const visit = (node: BTreeNode): void => { + if (distanceToAabbSq(node, x, y, z) > distances[k - 1]) return; + + if (node.indices) { + for (let i = 0; i < node.indices.length; i++) { + const candidate = node.indices[i]; + if (lodOf(candidate) !== targetLod) continue; + const dx = slim.x[candidate] - x; + const dy = slim.y[candidate] - y; + const dz = slim.z[candidate] - z; + const distance = dx * dx + dy * dy + dz * dz; + if (distance >= distances[k - 1]) continue; + + let insert = k - 1; + while (insert > 0 && distance < distances[insert - 1]) { + distances[insert] = distances[insert - 1]; + result[output + insert] = result[output + insert - 1]; + insert--; + } + distances[insert] = distance; + result[output + insert] = candidate; + } + return; + } + + const left = node.left!; + const right = node.right!; + if (distanceToAabbSq(left, x, y, z) <= distanceToAabbSq(right, x, y, z)) { + visit(left); + visit(right); + } else { + visit(right); + visit(left); + } + }; + visit(root); + } + + return result; +}; + +const uniqueGlobals = (queries: Uint32Array, neighbours: Int32Array): Uint32Array => { + const combined = new Uint32Array(queries.length + neighbours.length); + combined.set(queries); + let count = queries.length; + for (let i = 0; i < neighbours.length; i++) { + if (neighbours[i] >= 0) combined[count++] = neighbours[i]; + } + const sorted = combined.subarray(0, count); + sorted.sort(); + let unique = 0; + for (let i = 0; i < sorted.length; i++) { + if (i === 0 || sorted[i] !== sorted[i - 1]) sorted[unique++] = sorted[i]; + } + return sorted.slice(0, unique); +}; + +const indexOfSorted = (values: Uint32Array, value: number): number => { + let lo = 0; + let hi = values.length - 1; + while (lo <= hi) { + const mid = (lo + hi) >>> 1; + const candidate = values[mid]; + if (candidate < value) lo = mid + 1; + else if (candidate > value) hi = mid - 1; + else return mid; + } + return -1; +}; + +const gatherView = async ( + source: ChunkSource, pool: ChunkDataPool, slim: SlimColumns, + globals: Uint32Array, lod: number, base: number +): Promise => { + const { layouts } = source.meta; + const colorDim = layouts.color!.stride >> 2; + const view: SplatView = { + pos: new Float32Array(globals.length * 3), + geo: new Float32Array(globals.length * 8), + color: new Float32Array(globals.length * colorDim), + colorDim + }; + const local = new Uint32Array(globals.length); + for (let i = 0; i < globals.length; i++) { + const g = globals[i]; + local[i] = g - base; + view.pos[i * 3] = slim.x[g]; + view.pos[i * 3 + 1] = slim.y[g]; + view.pos[i * 3 + 2] = slim.z[g]; + } + + for (let off = 0; off < local.length; off += pool.chunkSize) { + const count = Math.min(pool.chunkSize, local.length - off); + const geo = pool.acquire('geometric', layouts.geometric!, count); + const color = pool.acquire('color', layouts.color!, count); + await source.read({ indices: local, indexOffset: off, count, lod, geometric: geo, color }); + view.geo.set(new Float32Array(geo.data, 0, count * 8), off * 8); + view.color.set(new Float32Array(color.data, 0, count * colorDim), off * colorDim); + geo.release(); + color.release(); + } + + return view; +}; + +const combineViews = (a: SplatView, b: SplatView): SplatView => { + const countA = a.pos.length / 3; + const countB = b.pos.length / 3; + const result: SplatView = { + pos: new Float32Array((countA + countB) * 3), + geo: new Float32Array((countA + countB) * 8), + color: new Float32Array((countA + countB) * a.colorDim), + colorDim: a.colorDim + }; + result.pos.set(a.pos); result.pos.set(b.pos, a.pos.length); + result.geo.set(a.geo); result.geo.set(b.geo, a.geo.length); + result.color.set(a.color); result.color.set(b.color, a.color.length); + return result; +}; + +const PI_1_5 = Math.PI ** 1.5; +const TWO_PI_1_5 = (2 * Math.PI) ** 1.5; + +// Per-splat terms of the field-L2: the gaussian's opacity, sqrt(det sigma), and +// its self inner product = alpha^2 * pi^1.5 * sqrt(det sigma). +type ErrorCache = { + alpha: Float64Array; + sqrtDet: Float64Array; + self: Float64Array; +}; + +const buildErrorCache = (view: SplatView, cache: CostCache): ErrorCache => { + const n = view.geo.length / 8; + const alpha = new Float64Array(n); + const sqrtDet = new Float64Array(n); + const self = new Float64Array(n); + + for (let i = 0; i < n; i++) { + const a = sigmoid(view.geo[i * 8 + 7]); + const root = Math.sqrt(Math.max(det3(cache.sigma, i * 9), 1e-300)); + alpha[i] = a; + sqrtDet[i] = root; + self[i] = a * a * PI_1_5 * root; + } + + return { alpha, sqrtDet, self }; +}; + +const sumScratch = new Float64Array(9); + +/** + * Approximation error between two splats: the relative field-L2 of the density + * each one paints plus the stored-SH L2. + * + * Writing f(x) = alpha * exp(-0.5 (x-mu)^T sigma^-1 (x-mu)), every inner product + * is closed form, so + * + * ||f_i - f_j||^2 / (||f_i||^2 + ||f_j||^2) + * + * is exact: zero only for identical splats, monotone in their separation, and + * bounded by 1 once they no longer overlap. The decimator's edge cost is not + * usable here — as a one-sample MC estimate of a KL it carries ~1.2 nats of + * noise, far more than the error a well-decimated leaf actually has, and it is + * blind to opacity (alpha only sets the merge weights there). + * + * @param view - Splat columns. + * @param cache - Per-splat cache from `buildCostCache` (supplies sigma). + * @param errorCache - Per-splat cache from {@link buildErrorCache}. + * @param i - First splat (view row). + * @param j - Second splat (view row). + * @returns The error, zero when the two splats are identical. + */ +const splatError = ( + view: SplatView, cache: CostCache, errorCache: ErrorCache, i: number, j: number +): number => { + const s = sumScratch; + const i9 = 9 * i, j9 = 9 * j; + for (let a = 0; a < 9; a++) s[a] = cache.sigma[i9 + a] + cache.sigma[j9 + a]; + const detS = Math.max(det3(s, 0), 1e-300); + + // d^T (sigma_i + sigma_j)^-1 d, via the adjugate (symmetric, so no transpose) + const { pos } = view; + const dx = pos[3 * i] - pos[3 * j]; + const dy = pos[3 * i + 1] - pos[3 * j + 1]; + const dz = pos[3 * i + 2] - pos[3 * j + 2]; + const ax = (s[4] * s[8] - s[5] * s[7]) * dx + (s[2] * s[7] - s[1] * s[8]) * dy + (s[1] * s[5] - s[2] * s[4]) * dz; + const ay = (s[5] * s[6] - s[3] * s[8]) * dx + (s[0] * s[8] - s[2] * s[6]) * dy + (s[2] * s[3] - s[0] * s[5]) * dz; + const az = (s[3] * s[7] - s[4] * s[6]) * dx + (s[1] * s[6] - s[0] * s[7]) * dy + (s[0] * s[4] - s[1] * s[3]) * dz; + const quad = (dx * ax + dy * ay + dz * az) / detS; + + const cross = errorCache.alpha[i] * errorCache.alpha[j] * TWO_PI_1_5 * + errorCache.sqrtDet[i] * errorCache.sqrtDet[j] / Math.sqrt(detS) * Math.exp(-0.5 * quad); + const total = errorCache.self[i] + errorCache.self[j]; + const geoError = Math.max(0, total - 2 * cross) / Math.max(total, 1e-300); + + const { color, colorDim } = view; + let colorError = 0; + for (let c = 0; c < colorDim; c++) { + const d = color[i * colorDim + c] - color[j * colorDim + c]; + colorError += d * d; + } + + return geoError + colorError; +}; + +const calcErrors = async ( + source: ChunkSource, pool: ChunkDataPool, slim: SlimColumns, root: BTreeNode, + bins: Map, lodOf: (index: number) => number, cum: number[], numLods: number +): Promise => { + // Symmetric Chamfer-style error against the finest representation present in + // this leaf, plus the alpha-mass a level fails to carry. Center-space KNN + // supplies a small candidate set (including candidates across leaf + // boundaries); the final match uses the closed-form {@link splatError}. This + // compares arbitrary input LODs, so it does not depend on how those LODs + // were produced. + const errors = new Array(numLods).fill(0); + const referenceLod = Math.min(...bins.keys()); + const referenceQueries = bins.get(referenceLod)!; + const k = 4; + + for (const [lod, targetQueries] of bins) { + if (lod === referenceLod) continue; + + const forward = findNearest(root, slim, lodOf, referenceQueries, lod, k); + const reverse = findNearest(root, slim, lodOf, targetQueries, referenceLod, k); + const referenceGlobals = uniqueGlobals(referenceQueries, reverse); + const targetGlobals = uniqueGlobals(targetQueries, forward); + const referenceView = await gatherView(source, pool, slim, referenceGlobals, referenceLod, cum[referenceLod]); + const targetView = await gatherView(source, pool, slim, targetGlobals, lod, cum[lod]); + const view = combineViews(referenceView, targetView); + const cache = buildCostCache(view); + const errorCache = buildErrorCache(view, cache); + const targetOffset = referenceGlobals.length; + + // `mass` (alpha * ellipsoid area) is a splat's share of what the leaf + // paints, so it weights both the per-splat mean and the coverage term — + // a hair-thin splat and one spanning the whole leaf are not equals. + const directional = ( + queries: Uint32Array, neighbours: Int32Array, + queryGlobals: Uint32Array, neighbourGlobals: Uint32Array, + queryOffset: number, neighbourOffset: number + ): number => { + let total = 0; + let weight = 0; + for (let i = 0; i < queries.length; i++) { + const queryRow = queryOffset + indexOfSorted(queryGlobals, queries[i]); + let best = Infinity; + for (let j = 0; j < k; j++) { + const neighbour = neighbours[i * k + j]; + if (neighbour < 0) continue; + const neighbourRow = neighbourOffset + indexOfSorted(neighbourGlobals, neighbour); + best = Math.min(best, splatError(view, cache, errorCache, queryRow, neighbourRow)); + } + if (best === Infinity) continue; + total += cache.mass[queryRow] * best; + weight += cache.mass[queryRow]; + } + return weight > 0 ? total / weight : 0; + }; + + const massOf = (queries: Uint32Array, globals: Uint32Array, offset: number): number => { + let mass = 0; + for (let i = 0; i < queries.length; i++) { + mass += cache.mass[offset + indexOfSorted(globals, queries[i])]; + } + return mass; + }; + + // Nearest-neighbour matching cannot see thinning: drop every second + // splat of an overlapping group and each survivor still has a + // near-identical partner, while the alpha the group accumulates halves. + // Sky gaps are exactly that, so compare the mass the levels carry. + const referenceMass = massOf(referenceQueries, referenceGlobals, 0); + const targetMass = massOf(targetQueries, targetGlobals, targetOffset); + const coverageError = Math.max(0, 1 - targetMass / Math.max(referenceMass, 1e-300)); + + const forwardError = directional(referenceQueries, forward, referenceGlobals, targetGlobals, 0, targetOffset); + const reverseError = directional(targetQueries, reverse, targetGlobals, referenceGlobals, targetOffset, 0); + errors[lod] = 0.5 * (forwardError + reverseError) + coverageError; + } + + // The engine's budget balancer drops a level from its frontier whenever a + // coarser one claims less error, which pins the node to that coarser level + // at any budget. Keep the table monotone so a level can only ever be + // dropped for being genuinely no better. + let previous = 0; + for (const lod of [...bins.keys()].sort((a, b) => a - b)) { + previous = errors[lod] = Math.max(errors[lod], previous); + } + + return errors; +}; + /** * Read positions out of a multi-LOD source into flat per-gaussian arrays — one * sequential pass across every structural LOD (LOD 0 first, then 1, …, laid out @@ -424,10 +742,11 @@ const writeLodSource = async (options: WriteLodSourceOptions, fs: FileSystem) => lodLevels = Math.max(lodLevels, lodValue + 1); } - // bound over the leaf's gaussians, gathered per structural LOD. + // Bound and approximation errors over the leaf's full structural LOD data. const bound = await calcBound(mainSource, pool, bins, cum, n => chunkingBar.tick(n)); + const errors = await calcErrors(mainSource, pool, slim, bTree!.root, bins, lodOf, cum, numLods); - return { bound, lods }; + return { bound, lods, errors }; }; let tree: MetaNode; diff --git a/test/write-lod.test.mjs b/test/write-lod.test.mjs index d4b750e3..3d565818 100644 --- a/test/write-lod.test.mjs +++ b/test/write-lod.test.mjs @@ -79,6 +79,53 @@ const makeTable = (n) => { ], Transform.PLY); }; +// A table from explicit per-splat values, for exercising the LOD error metric. +// Defaults put every splat at the origin with log-scale -3 (sigma ~0.0498) and +// opacity logit 0 (alpha 0.5). +const makeSplatTable = (splats) => { + const col = (key, fallback) => new Float32Array(splats.map(s => s[key] ?? fallback)); + return new DataTable([ + new Column('x', col('x', 0)), + new Column('y', col('y', 0)), + new Column('z', col('z', 0)), + new Column('rot_0', col('rot_0', 1)), + new Column('rot_1', col('rot_1', 0)), + new Column('rot_2', col('rot_2', 0)), + new Column('rot_3', col('rot_3', 0)), + new Column('scale_0', col('scale', -3)), + new Column('scale_1', col('scale', -3)), + new Column('scale_2', col('scale', -3)), + new Column('f_dc_0', col('f_dc', 0)), + new Column('f_dc_1', col('f_dc', 0)), + new Column('f_dc_2', col('f_dc', 0)), + new Column('opacity', col('opacity', 0)) + ], Transform.PLY); +}; + +// The error table of a scene whose levels are given splat-by-splat. +const writeErrors = async (levels) => { + const fs = new MemoryFileSystem(); + const sources = levels.map(splats => dataTableToChunkSource(makeSplatTable(splats), 1 << 20)); + await writeLodSource({ + filename: '/scene/lod-meta.json', + mainSource: sources.length === 1 ? sources[0] : stackLods(sources), + envSource: null, + iterations: 1, + chunkCount: 1, + chunkExtent: 16 + }, fs); + const meta = JSON.parse(new TextDecoder().decode(fs.results.get('/scene/lod-meta.json'))); + return meta.tree.errors; +}; + +const makeShTable = (restValue) => { + const table = makeTable(1); + for (let i = 0; i < 9; i++) { + table.addColumn(new Column(`f_rest_${i}`, new Float32Array([restValue]))); + } + return table; +}; + // Build a structural multi-LOD source from per-level row counts (each level a // single-LOD resident source; stacked when there is more than one level). const makeSource = (levelCounts) => { @@ -139,11 +186,63 @@ describe('writeLodSource: lod-meta.json contract', function () { { offset: meta.tree.lods['1'].offset, count: meta.tree.lods['1'].count }, { offset: 0, count: 2 } ); + assert.strictEqual(meta.tree.errors.length, 2); + assert.strictEqual(meta.tree.errors[0], 0); + assert.ok(Number.isFinite(meta.tree.errors[1]) && meta.tree.errors[1] >= 0); assert.ok(fs.results.has('/scene/0_0/meta.json')); assert.ok(fs.results.has('/scene/1_0/meta.json')); }); + it('includes stored spherical harmonics in the LOD error', async function () { + const fs = new MemoryFileSystem(); + await writeLodSource({ + filename: '/scene/lod-meta.json', + mainSource: stackLods([ + dataTableToChunkSource(makeShTable(0), 1 << 20), + dataTableToChunkSource(makeShTable(2), 1 << 20) + ]), + envSource: null, + iterations: 1, + chunkCount: 1, + chunkExtent: 16 + }, fs); + + const meta = JSON.parse(new TextDecoder().decode(fs.results.get('/scene/lod-meta.json'))); + assert.ok(meta.tree.errors[1] > 0); + }); + + it('reports no error for a level identical to the finest', async function () { + assert.deepStrictEqual(await writeErrors([[{}], [{}]]), [0, 0]); + }); + + it('resolves a sub-sigma displacement', async function () { + // half a sigma apart: 1 - exp(-d^2/4) = 0.0606 for the relative field-L2 + const sigma = Math.exp(-3); + const errors = await writeErrors([[{}], [{ x: 0.5 * sigma }]]); + assert.ok(errors[1] > 0.05 && errors[1] < 0.07, `expected ~0.0606, got ${errors[1]}`); + }); + + it('penalises an opacity drop at identical geometry', async function () { + const errors = await writeErrors([[{ opacity: 2 }], [{ opacity: -2 }]]); + assert.ok(errors[1] > 0.5, `expected a large error, got ${errors[1]}`); + }); + + it('penalises thinning even when the survivors are identical', async function () { + // two coincident splats decimated to one: every nearest-neighbour match + // is exact, but the level carries half the alpha mass + const errors = await writeErrors([[{}, {}], [{}]]); + assert.ok(Math.abs(errors[1] - 0.5) < 1e-6, `expected 0.5, got ${errors[1]}`); + }); + + it('keeps the error table monotone across levels', async function () { + // level 2 matches a level-0 splat exactly while level 1 sits between + // both, so the raw errors would rank the coarser level as the better one + const errors = await writeErrors([[{ x: 0 }, { x: 1 }], [{ x: 0.5 }], [{ x: 0 }]]); + assert.ok(errors[1] > 0, `expected a non-zero error, got ${errors[1]}`); + assert.ok(errors[2] >= errors[1], `expected monotone errors, got ${errors}`); + }); + it('references the environment SOG when environment splats are present', async function () { const { fs, meta } = await writeScene([3], 2); assert.strictEqual(meta.environment, 'env/meta.json'); From a3009d0090c2eb29a073675408aa7d04883292c4 Mon Sep 17 00:00:00 2001 From: Donovan Hutchence Date: Sat, 1 Aug 2026 21:01:39 +0100 Subject: [PATCH 2/5] latest --- src/lib/writers/write-lod.ts | 368 +++++++++++++++++++++++------------ 1 file changed, 247 insertions(+), 121 deletions(-) diff --git a/src/lib/writers/write-lod.ts b/src/lib/writers/write-lod.ts index 83ca1426..75114ce2 100644 --- a/src/lib/writers/write-lod.ts +++ b/src/lib/writers/write-lod.ts @@ -5,8 +5,9 @@ import { logWrittenFile } from './utils'; import { writeSogSource } from './write-sog.js'; import { type ChunkDataPool, type ChunkSource, type ReadRequest, createChunkDataPool } from '../chunk'; import { Column, DataTable } from '../data-table'; -import { det3, sigmoid, type SplatView } from '../decimate/moment-match'; -import { buildCostCache, type CostCache } from '../decimate-uniform/edge-cost-cpu'; +import { + EPS_COV, det3, ellipsoidArea, quatToRotmat, sigmaFromRotVar, sigmoid, type SplatView +} from '../decimate/moment-match'; import { type FileSystem } from '../io/write'; import { bakeTransform, permuteSource, sortMortonColumns } from '../ops'; import { BTreeNode, BTree } from '../spatial'; @@ -252,13 +253,40 @@ const distanceToAabbSq = (node: BTreeNode, x: number, y: number, z: number): num return dx * dx + dy * dy + dz * dz; }; +/** + * The k nearest candidates of each of several structural LODs to every query, in + * center space — one traversal per query, serving every LOD. + * + * `ranges` holds `[lo, hi)` per LOD in flat analysis index space: LOD is + * structural, so a level's gaussians are exactly one contiguous index range and + * the level test is two comparisons rather than a lookup. Sharing the traversal + * is what makes the error pass affordable — every splat of the reference level is + * matched against every coarser level, so searching them separately would repeat + * the same descent once per level. + * + * A node is descended while any level can still improve on it, and a leaf's + * candidates are examined for a level only when that level's own k-th distance + * says they could, so each level's result is identical to a search of its own. + * Each node's distance is computed once, by its parent, and passed in. + * + * @param root - Root of the tree to search. + * @param slim - Resident position columns. + * @param queries - Flat analysis indices to search from. + * @param ranges - Flat `[lo, hi)` pairs, one per LOD to match against. + * @param k - Neighbours per query per LOD. + * @returns Per LOD, `k` flat indices per query, nearest first, `-1` for none. + */ const findNearest = ( - root: BTreeNode, slim: SlimColumns, lodOf: (index: number) => number, - queries: Uint32Array, targetLod: number, k: number -): Int32Array => { - const result = new Int32Array(queries.length * k); - result.fill(-1); - const distances = new Float64Array(k); + root: BTreeNode, slim: SlimColumns, queries: Uint32Array, ranges: Int32Array, k: number +): Int32Array[] => { + const targets = ranges.length >> 1; + const results: Int32Array[] = []; + for (let t = 0; t < targets; t++) { + const result = new Int32Array(queries.length * k); + result.fill(-1); + results.push(result); + } + const distances = new Float64Array(targets * k); for (let q = 0; q < queries.length; q++) { distances.fill(Infinity); @@ -266,53 +294,78 @@ const findNearest = ( const x = slim.x[g], y = slim.y[g], z = slim.z[g]; const output = q * k; - const visit = (node: BTreeNode): void => { - if (distanceToAabbSq(node, x, y, z) > distances[k - 1]) return; - - if (node.indices) { - for (let i = 0; i < node.indices.length; i++) { - const candidate = node.indices[i]; - if (lodOf(candidate) !== targetLod) continue; - const dx = slim.x[candidate] - x; - const dy = slim.y[candidate] - y; - const dz = slim.z[candidate] - z; - const distance = dx * dx + dy * dy + dz * dz; - if (distance >= distances[k - 1]) continue; - - let insert = k - 1; - while (insert > 0 && distance < distances[insert - 1]) { - distances[insert] = distances[insert - 1]; - result[output + insert] = result[output + insert - 1]; - insert--; + // the walk can stop where no level can improve: the largest k-th distance + let radius = Infinity; + + const visit = (node: BTreeNode, distance: number): void => { + if (distance > radius) return; + + const { indices } = node; + if (indices) { + for (let t = 0; t < targets; t++) { + const base = t * k; + const worst = base + k - 1; + if (distance > distances[worst]) continue; + + const lo = ranges[2 * t]; + const hi = ranges[2 * t + 1]; + const result = results[t]; + for (let i = 0; i < indices.length; i++) { + const candidate = indices[i]; + if (candidate < lo || candidate >= hi) continue; + const dx = slim.x[candidate] - x; + const dy = slim.y[candidate] - y; + const dz = slim.z[candidate] - z; + const d = dx * dx + dy * dy + dz * dz; + if (d >= distances[worst]) continue; + + let insert = k - 1; + while (insert > 0 && d < distances[base + insert - 1]) { + distances[base + insert] = distances[base + insert - 1]; + result[output + insert] = result[output + insert - 1]; + insert--; + } + distances[base + insert] = d; + result[output + insert] = candidate; } - distances[insert] = distance; - result[output + insert] = candidate; + } + + radius = 0; + for (let t = 0; t < targets; t++) { + const d = distances[t * k + k - 1]; + if (d > radius) radius = d; } return; } const left = node.left!; const right = node.right!; - if (distanceToAabbSq(left, x, y, z) <= distanceToAabbSq(right, x, y, z)) { - visit(left); - visit(right); + const dl = distanceToAabbSq(left, x, y, z); + const dr = distanceToAabbSq(right, x, y, z); + if (dl <= dr) { + visit(left, dl); + visit(right, dr); } else { - visit(right); - visit(left); + visit(right, dr); + visit(left, dl); } }; - visit(root); + visit(root, distanceToAabbSq(root, x, y, z)); } - return result; + return results; }; -const uniqueGlobals = (queries: Uint32Array, neighbours: Int32Array): Uint32Array => { - const combined = new Uint32Array(queries.length + neighbours.length); +const uniqueGlobals = (queries: Uint32Array, neighbourSets: Int32Array[]): Uint32Array => { + let capacity = queries.length; + for (const neighbours of neighbourSets) capacity += neighbours.length; + const combined = new Uint32Array(capacity); combined.set(queries); let count = queries.length; - for (let i = 0; i < neighbours.length; i++) { - if (neighbours[i] >= 0) combined[count++] = neighbours[i]; + for (const neighbours of neighbourSets) { + for (let i = 0; i < neighbours.length; i++) { + if (neighbours[i] >= 0) combined[count++] = neighbours[i]; + } } const sorted = combined.subarray(0, count); sorted.sort(); @@ -371,47 +424,80 @@ const gatherView = async ( return view; }; -const combineViews = (a: SplatView, b: SplatView): SplatView => { - const countA = a.pos.length / 3; - const countB = b.pos.length / 3; +const concatViews = (views: SplatView[]): SplatView => { + const { colorDim } = views[0]; + let count = 0; + for (const v of views) count += v.pos.length / 3; const result: SplatView = { - pos: new Float32Array((countA + countB) * 3), - geo: new Float32Array((countA + countB) * 8), - color: new Float32Array((countA + countB) * a.colorDim), - colorDim: a.colorDim + pos: new Float32Array(count * 3), + geo: new Float32Array(count * 8), + color: new Float32Array(count * colorDim), + colorDim }; - result.pos.set(a.pos); result.pos.set(b.pos, a.pos.length); - result.geo.set(a.geo); result.geo.set(b.geo, a.geo.length); - result.color.set(a.color); result.color.set(b.color, a.color.length); + let row = 0; + for (const v of views) { + result.pos.set(v.pos, row * 3); + result.geo.set(v.geo, row * 8); + result.color.set(v.color, row * colorDim); + row += v.pos.length / 3; + } return result; }; const PI_1_5 = Math.PI ** 1.5; const TWO_PI_1_5 = (2 * Math.PI) ** 1.5; -// Per-splat terms of the field-L2: the gaussian's opacity, sqrt(det sigma), and -// its self inner product = alpha^2 * pi^1.5 * sqrt(det sigma). +/** + * Per-splat terms of the field-L2: the covariance, the footprint mass used to + * weight it, the gaussian's opacity, sqrt(det sigma), and its self inner product + * = alpha^2 * pi^1.5 * sqrt(det sigma). + * + * `sigma` and `mass` are what the decimator's `buildCostCache` computes (same + * math, same f32 storage, so the same values); the rest of that cache serves the + * MC-KL edge cost, which this path does not evaluate, so it is not built here. + */ type ErrorCache = { + sigma: Float32Array; + mass: Float32Array; alpha: Float64Array; sqrtDet: Float64Array; self: Float64Array; }; -const buildErrorCache = (view: SplatView, cache: CostCache): ErrorCache => { - const n = view.geo.length / 8; +const buildErrorCache = (view: SplatView): ErrorCache => { + const { geo } = view; + const n = geo.length / 8; + const sigma = new Float32Array(n * 9); + const mass = new Float32Array(n); const alpha = new Float64Array(n); const sqrtDet = new Float64Array(n); const self = new Float64Array(n); + const rot = new Float32Array(9); for (let i = 0; i < n; i++) { - const a = sigmoid(view.geo[i * 8 + 7]); - const root = Math.sqrt(Math.max(det3(cache.sigma, i * 9), 1e-300)); - alpha[i] = a; + const i8 = 8 * i; + const i9 = 9 * i; + + const linAlpha = sigmoid(geo[i8 + 7]); + const sx = Math.max(Math.exp(geo[i8 + 4]), 1e-12); + const sy = Math.max(Math.exp(geo[i8 + 5]), 1e-12); + const sz = Math.max(Math.exp(geo[i8 + 6]), 1e-12); + + let qw = geo[i8], qx = geo[i8 + 1], qy = geo[i8 + 2], qz = geo[i8 + 3]; + const invq = 1 / Math.max(Math.hypot(qw, qx, qy, qz), 1e-12); + qw *= invq; qx *= invq; qy *= invq; qz *= invq; + + quatToRotmat(qw, qx, qy, qz, rot, 0); + sigmaFromRotVar(rot, 0, sx * sx + EPS_COV, sy * sy + EPS_COV, sz * sz + EPS_COV, sigma, i9); + + const root = Math.sqrt(Math.max(det3(sigma, i9), 1e-300)); + alpha[i] = linAlpha; sqrtDet[i] = root; - self[i] = a * a * PI_1_5 * root; + self[i] = linAlpha * linAlpha * PI_1_5 * root; + mass[i] = linAlpha * ellipsoidArea(sx, sy, sz) + 1e-12; } - return { alpha, sqrtDet, self }; + return { sigma, mass, alpha, sqrtDet, self }; }; const sumScratch = new Float64Array(9); @@ -431,15 +517,19 @@ const sumScratch = new Float64Array(9); * noise, far more than the error a well-decimated leaf actually has, and it is * blind to opacity (alpha only sets the merge weights there). * + * Both terms are non-negative, so a pair whose geometric term alone already + * reaches `cutoff` cannot beat it: `Infinity` is returned without summing the SH + * coefficients, which is the bulk of the arithmetic at 3 SH bands. + * * @param view - Splat columns. - * @param cache - Per-splat cache from `buildCostCache` (supplies sigma). - * @param errorCache - Per-splat cache from {@link buildErrorCache}. + * @param cache - Per-splat cache from {@link buildErrorCache}. * @param i - First splat (view row). * @param j - Second splat (view row). + * @param cutoff - Error to beat; pass `Infinity` for the unconditional error. * @returns The error, zero when the two splats are identical. */ const splatError = ( - view: SplatView, cache: CostCache, errorCache: ErrorCache, i: number, j: number + view: SplatView, cache: ErrorCache, i: number, j: number, cutoff: number ): number => { const s = sumScratch; const i9 = 9 * i, j9 = 9 * j; @@ -456,10 +546,11 @@ const splatError = ( const az = (s[3] * s[7] - s[4] * s[6]) * dx + (s[1] * s[6] - s[0] * s[7]) * dy + (s[0] * s[4] - s[1] * s[3]) * dz; const quad = (dx * ax + dy * ay + dz * az) / detS; - const cross = errorCache.alpha[i] * errorCache.alpha[j] * TWO_PI_1_5 * - errorCache.sqrtDet[i] * errorCache.sqrtDet[j] / Math.sqrt(detS) * Math.exp(-0.5 * quad); - const total = errorCache.self[i] + errorCache.self[j]; + const cross = cache.alpha[i] * cache.alpha[j] * TWO_PI_1_5 * + cache.sqrtDet[i] * cache.sqrtDet[j] / Math.sqrt(detS) * Math.exp(-0.5 * quad); + const total = cache.self[i] + cache.self[j]; const geoError = Math.max(0, total - 2 * cross) / Math.max(total, 1e-300); + if (geoError >= cutoff) return Infinity; const { color, colorDim } = view; let colorError = 0; @@ -473,7 +564,7 @@ const splatError = ( const calcErrors = async ( source: ChunkSource, pool: ChunkDataPool, slim: SlimColumns, root: BTreeNode, - bins: Map, lodOf: (index: number) => number, cum: number[], numLods: number + bins: Map, cum: number[], numLods: number ): Promise => { // Symmetric Chamfer-style error against the finest representation present in // this leaf, plus the alpha-mass a level fails to carry. Center-space KNN @@ -486,66 +577,101 @@ const calcErrors = async ( const referenceQueries = bins.get(referenceLod)!; const k = 4; - for (const [lod, targetQueries] of bins) { - if (lod === referenceLod) continue; - - const forward = findNearest(root, slim, lodOf, referenceQueries, lod, k); - const reverse = findNearest(root, slim, lodOf, targetQueries, referenceLod, k); - const referenceGlobals = uniqueGlobals(referenceQueries, reverse); - const targetGlobals = uniqueGlobals(targetQueries, forward); - const referenceView = await gatherView(source, pool, slim, referenceGlobals, referenceLod, cum[referenceLod]); - const targetView = await gatherView(source, pool, slim, targetGlobals, lod, cum[lod]); - const view = combineViews(referenceView, targetView); - const cache = buildCostCache(view); - const errorCache = buildErrorCache(view, cache); - const targetOffset = referenceGlobals.length; - - // `mass` (alpha * ellipsoid area) is a splat's share of what the leaf - // paints, so it weights both the per-splat mean and the coverage term — - // a hair-thin splat and one spanning the whole leaf are not equals. - const directional = ( - queries: Uint32Array, neighbours: Int32Array, - queryGlobals: Uint32Array, neighbourGlobals: Uint32Array, - queryOffset: number, neighbourOffset: number - ): number => { - let total = 0; - let weight = 0; - for (let i = 0; i < queries.length; i++) { - const queryRow = queryOffset + indexOfSorted(queryGlobals, queries[i]); - let best = Infinity; - for (let j = 0; j < k; j++) { - const neighbour = neighbours[i * k + j]; - if (neighbour < 0) continue; - const neighbourRow = neighbourOffset + indexOfSorted(neighbourGlobals, neighbour); - best = Math.min(best, splatError(view, cache, errorCache, queryRow, neighbourRow)); - } - if (best === Infinity) continue; - total += cache.mass[queryRow] * best; - weight += cache.mass[queryRow]; - } - return weight > 0 ? total / weight : 0; - }; + const targetLods = [...bins.keys()].filter(lod => lod !== referenceLod).sort((a, b) => a - b); + if (targetLods.length === 0) return errors; + + // Forward (reference -> each level) shares one traversal per reference splat; + // reverse (each level -> reference) is a separate query set per level. + const targetRanges = new Int32Array(targetLods.length * 2); + targetLods.forEach((lod, i) => { + targetRanges[2 * i] = cum[lod]; + targetRanges[2 * i + 1] = cum[lod + 1]; + }); + const referenceRange = Int32Array.of(cum[referenceLod], cum[referenceLod + 1]); + + const forwards = findNearest(root, slim, referenceQueries, targetRanges, k); + const reverses = targetLods.map(lod => findNearest(root, slim, bins.get(lod)!, referenceRange, k)[0]); + + // One gathered section per level, holding that level's own splats plus every + // splat another level matched into it, concatenated into a single view. All + // levels are handled in one pass so the reference level — the largest, and a + // participant in every pair — is read and cached once instead of once per + // level pair. + const order = [referenceLod, ...targetLods]; + const globals: Uint32Array[] = new Array(numLods); + globals[referenceLod] = uniqueGlobals(referenceQueries, reverses); + targetLods.forEach((lod, i) => { + globals[lod] = uniqueGlobals(bins.get(lod)!, [forwards[i]]); + }); + + const offsets = new Int32Array(numLods); + const sections: SplatView[] = []; + let rows = 0; + for (const lod of order) { + offsets[lod] = rows; + sections.push(await gatherView(source, pool, slim, globals[lod], lod, cum[lod])); + rows += globals[lod].length; + } - const massOf = (queries: Uint32Array, globals: Uint32Array, offset: number): number => { - let mass = 0; - for (let i = 0; i < queries.length; i++) { - mass += cache.mass[offset + indexOfSorted(globals, queries[i])]; + const view = concatViews(sections); + const cache = buildErrorCache(view); + + // Rows of a level's own splats, in query order. The reference level's are + // reused by every pair. + const rowsOf = (queries: Uint32Array, lod: number): Int32Array => { + const result = new Int32Array(queries.length); + const offset = offsets[lod]; + const g = globals[lod]; + for (let i = 0; i < queries.length; i++) result[i] = offset + indexOfSorted(g, queries[i]); + return result; + }; + const referenceRows = rowsOf(referenceQueries, referenceLod); + + // `mass` (alpha * ellipsoid area) is a splat's share of what the leaf + // paints, so it weights both the per-splat mean and the coverage term — + // a hair-thin splat and one spanning the whole leaf are not equals. + const directional = ( + queryRows: Int32Array, neighbours: Int32Array, neighbourLod: number + ): number => { + let total = 0; + let weight = 0; + const neighbourOffset = offsets[neighbourLod]; + const neighbourGlobals = globals[neighbourLod]; + for (let i = 0; i < queryRows.length; i++) { + const queryRow = queryRows[i]; + let best = Infinity; + for (let j = 0; j < k; j++) { + const neighbour = neighbours[i * k + j]; + if (neighbour < 0) continue; + const neighbourRow = neighbourOffset + indexOfSorted(neighbourGlobals, neighbour); + const error = splatError(view, cache, queryRow, neighbourRow, best); + if (error < best) best = error; } - return mass; - }; + if (best === Infinity) continue; + total += cache.mass[queryRow] * best; + weight += cache.mass[queryRow]; + } + return weight > 0 ? total / weight : 0; + }; - // Nearest-neighbour matching cannot see thinning: drop every second - // splat of an overlapping group and each survivor still has a - // near-identical partner, while the alpha the group accumulates halves. - // Sky gaps are exactly that, so compare the mass the levels carry. - const referenceMass = massOf(referenceQueries, referenceGlobals, 0); - const targetMass = massOf(targetQueries, targetGlobals, targetOffset); - const coverageError = Math.max(0, 1 - targetMass / Math.max(referenceMass, 1e-300)); + const massOf = (queryRows: Int32Array): number => { + let mass = 0; + for (let i = 0; i < queryRows.length; i++) mass += cache.mass[queryRows[i]]; + return mass; + }; - const forwardError = directional(referenceQueries, forward, referenceGlobals, targetGlobals, 0, targetOffset); - const reverseError = directional(targetQueries, reverse, targetGlobals, referenceGlobals, targetOffset, 0); + // Nearest-neighbour matching cannot see thinning: drop every second + // splat of an overlapping group and each survivor still has a + // near-identical partner, while the alpha the group accumulates halves. + // Sky gaps are exactly that, so compare the mass the levels carry. + const referenceMass = massOf(referenceRows); + targetLods.forEach((lod, i) => { + const targetRows = rowsOf(bins.get(lod)!, lod); + const coverageError = Math.max(0, 1 - massOf(targetRows) / Math.max(referenceMass, 1e-300)); + const forwardError = directional(referenceRows, forwards[i], lod); + const reverseError = directional(targetRows, reverses[i], referenceLod); errors[lod] = 0.5 * (forwardError + reverseError) + coverageError; - } + }); // The engine's budget balancer drops a level from its frontier whenever a // coarser one claims less error, which pins the node to that coarser level @@ -744,7 +870,7 @@ const writeLodSource = async (options: WriteLodSourceOptions, fs: FileSystem) => // Bound and approximation errors over the leaf's full structural LOD data. const bound = await calcBound(mainSource, pool, bins, cum, n => chunkingBar.tick(n)); - const errors = await calcErrors(mainSource, pool, slim, bTree!.root, bins, lodOf, cum, numLods); + const errors = await calcErrors(mainSource, pool, slim, bTree!.root, bins, cum, numLods); return { bound, lods, errors }; }; From 14cdf1be58440814968a43fd82ada97f4d65c4be Mon Sep 17 00:00:00 2001 From: Donovan Hutchence Date: Mon, 10 Aug 2026 11:22:34 +0200 Subject: [PATCH 3/5] latest --- src/lib/writers/write-lod.ts | 169 ++++++++++++++++++++++------------- test/write-lod.test.mjs | 72 ++++++++++++++- 2 files changed, 176 insertions(+), 65 deletions(-) diff --git a/src/lib/writers/write-lod.ts b/src/lib/writers/write-lod.ts index 75114ce2..e148a8a5 100644 --- a/src/lib/writers/write-lod.ts +++ b/src/lib/writers/write-lod.ts @@ -37,10 +37,20 @@ type LodMeta = { version: number; asset: { generator: string; + /** Gaussians per file unit the partition aimed for (`--lod-chunk-count` × 1024). */ + chunkGaussians: number; + /** Largest leaf extent the partition allowed, in world units (`--lod-chunk-extent`). */ + chunkExtent: number; }; count: number; counts: number[]; lodLevels: number; + /** + * Whether every leaf carries an `errors` table. Declared here so a consumer + * can pick its LOD allocation strategy up front instead of searching the tree + * for the field (the engine's budget balancer needs exactly this answer). + */ + lodErrors: boolean; environment?: string; filenames: string[]; tree: MetaNode; @@ -287,69 +297,81 @@ const findNearest = ( results.push(result); } const distances = new Float64Array(targets * k); + const targetK = new Int32Array(targets); + for (let t = 0; t < targets; t++) { + targetK[t] = Math.min(k, ranges[2 * t + 1] - ranges[2 * t]); + } - for (let q = 0; q < queries.length; q++) { - distances.fill(Infinity); - const g = queries[q]; - const x = slim.x[g], y = slim.y[g], z = slim.z[g]; - const output = q * k; - - // the walk can stop where no level can improve: the largest k-th distance - let radius = Infinity; - - const visit = (node: BTreeNode, distance: number): void => { - if (distance > radius) return; - - const { indices } = node; - if (indices) { - for (let t = 0; t < targets; t++) { - const base = t * k; - const worst = base + k - 1; - if (distance > distances[worst]) continue; - - const lo = ranges[2 * t]; - const hi = ranges[2 * t + 1]; - const result = results[t]; - for (let i = 0; i < indices.length; i++) { - const candidate = indices[i]; - if (candidate < lo || candidate >= hi) continue; - const dx = slim.x[candidate] - x; - const dy = slim.y[candidate] - y; - const dz = slim.z[candidate] - z; - const d = dx * dx + dy * dy + dz * dz; - if (d >= distances[worst]) continue; - - let insert = k - 1; - while (insert > 0 && d < distances[base + insert - 1]) { - distances[base + insert] = distances[base + insert - 1]; - result[output + insert] = result[output + insert - 1]; - insert--; - } - distances[base + insert] = d; - result[output + insert] = candidate; + const px = slim.x, py = slim.y, pz = slim.z; + + // Per-query state lives out here so the traversal is one function compiled + // once, not a fresh closure allocated per query. + let x = 0, y = 0, z = 0; + let output = 0; + + // the walk can stop where no level can improve: the largest k-th distance + let radius = Infinity; + + const visit = (node: BTreeNode, distance: number): void => { + if (distance > radius) return; + + const { indices } = node; + if (indices) { + for (let t = 0; t < targets; t++) { + const base = t * k; + const worst = base + targetK[t] - 1; + if (distance > distances[worst]) continue; + + const lo = ranges[2 * t]; + const hi = ranges[2 * t + 1]; + const result = results[t]; + for (let i = 0; i < indices.length; i++) { + const candidate = indices[i]; + if (candidate < lo || candidate >= hi) continue; + const dx = px[candidate] - x; + const dy = py[candidate] - y; + const dz = pz[candidate] - z; + const d = dx * dx + dy * dy + dz * dz; + if (d >= distances[worst]) continue; + + let insert = targetK[t] - 1; + while (insert > 0 && d < distances[base + insert - 1]) { + distances[base + insert] = distances[base + insert - 1]; + result[output + insert] = result[output + insert - 1]; + insert--; } + distances[base + insert] = d; + result[output + insert] = candidate; } - - radius = 0; - for (let t = 0; t < targets; t++) { - const d = distances[t * k + k - 1]; - if (d > radius) radius = d; - } - return; } - const left = node.left!; - const right = node.right!; - const dl = distanceToAabbSq(left, x, y, z); - const dr = distanceToAabbSq(right, x, y, z); - if (dl <= dr) { - visit(left, dl); - visit(right, dr); - } else { - visit(right, dr); - visit(left, dl); + radius = 0; + for (let t = 0; t < targets; t++) { + const d = distances[t * k + targetK[t] - 1]; + if (d > radius) radius = d; } - }; + return; + } + + const left = node.left!; + const right = node.right!; + const dl = distanceToAabbSq(left, x, y, z); + const dr = distanceToAabbSq(right, x, y, z); + if (dl <= dr) { + visit(left, dl); + visit(right, dr); + } else { + visit(right, dr); + visit(left, dl); + } + }; + + for (let q = 0; q < queries.length; q++) { + distances.fill(Infinity); + const g = queries[q]; + x = px[g]; y = py[g]; z = pz[g]; + output = q * k; + radius = Infinity; visit(root, distanceToAabbSq(root, x, y, z)); } @@ -673,13 +695,23 @@ const calcErrors = async ( errors[lod] = 0.5 * (forwardError + reverseError) + coverageError; }); - // The engine's budget balancer drops a level from its frontier whenever a - // coarser one claims less error, which pins the node to that coarser level - // at any budget. Keep the table monotone so a level can only ever be - // dropped for being genuinely no better. + // Keep the table monotone across levels. This does not keep a finer level on + // the engine's Pareto frontier — that domination test accepts an equal error, + // so a cheaper coarser level displaces the finer one either way — but the + // budget balancer ranks transitions by error reduction per splat *across* + // nodes, so a coarser level must never advertise less error than the finer + // one it stands in for. + // + // A NaN scale or opacity in the input poisons that splat's mass, and a NaN + // mass reaches the coverage term (unlike a NaN error, which the nearest + // match discards). Left alone it serializes as `null`, and a single + // non-finite error drops the engine back to distance-based allocation for the + // whole scene, so treat a level whose error is not a number as no better than + // the one before it. let previous = 0; for (const lod of [...bins.keys()].sort((a, b) => a - b)) { - previous = errors[lod] = Math.max(errors[lod], previous); + const error = errors[lod]; + previous = errors[lod] = Number.isFinite(error) ? Math.max(error, previous) : previous; } return errors; @@ -882,6 +914,12 @@ const writeLodSource = async (options: WriteLodSourceOptions, fs: FileSystem) => chunkingBar.end(); } + const trimErrors = (node: MetaNode): void => { + if (node.errors) node.errors.length = lodLevels; + for (const child of node.children ?? []) trimErrors(child); + }; + trimErrors(tree); + // The kd-tree is dead once the partition is built (lodFiles holds its own // index copies): release its N×4B index buffer and node AABBs before the // unit writes, where peak memory lives. @@ -898,11 +936,14 @@ const writeLodSource = async (options: WriteLodSourceOptions, fs: FileSystem) => const meta: LodMeta = { version: 1, asset: { - generator: `splat-transform v${version}` + generator: `splat-transform v${version}`, + chunkGaussians: binSize, + chunkExtent: binDim }, count: counts.reduce((acc, curr) => acc + curr, 0), counts, lodLevels, + lodErrors: true, ...(hasEnv ? { environment: 'env/meta.json' } : {}), filenames, tree @@ -1021,4 +1062,4 @@ const writeLodSource = async (options: WriteLodSourceOptions, fs: FileSystem) => writingGroup.end(); }; -export { positionsFromSlim, writeLodSource, type WriteLodSourceOptions }; +export { findNearest, positionsFromSlim, writeLodSource, type WriteLodSourceOptions }; diff --git a/test/write-lod.test.mjs b/test/write-lod.test.mjs index 3d565818..bbe3b035 100644 --- a/test/write-lod.test.mjs +++ b/test/write-lod.test.mjs @@ -21,7 +21,7 @@ import { bakeTransform, mapSource, stackLods } from '../src/lib/ops/index.js'; import { readPly } from '../src/lib/readers/read-ply.js'; import { collectFilesByLod, readLodEnvironmentSource } from '../src/lib/readers/read-lod.js'; import { createChunkDataPool } from '../src/lib/chunk/index.js'; -import { positionsFromSlim, writeLodSource } from '../src/lib/writers/write-lod.js'; +import { findNearest, positionsFromSlim, writeLodSource } from '../src/lib/writers/write-lod.js'; import { version } from '../src/lib/version.js'; import { encodePlyBinary } from './helpers/test-utils.mjs'; @@ -169,9 +169,14 @@ describe('writeLodSource: lod-meta.json contract', function () { assert.strictEqual(meta.version, 1); assert.strictEqual(meta.asset.generator, `splat-transform v${version}`); + // the partition parameters the caller asked for: chunkCount is in units of + // 1024 gaussians, chunkExtent in world units + assert.strictEqual(meta.asset.chunkGaussians, 1024); + assert.strictEqual(meta.asset.chunkExtent, 16); assert.strictEqual(meta.count, 5); assert.deepStrictEqual(meta.counts, [3, 2]); assert.strictEqual(meta.lodLevels, 2); + assert.strictEqual(meta.lodErrors, true, 'error tables are declared in the header'); assert.ok(!('environment' in meta), 'environment omitted when there are no environment splats'); assert.deepStrictEqual([...meta.filenames].sort(), ['0_0/meta.json', '1_0/meta.json']); @@ -194,6 +199,45 @@ describe('writeLodSource: lod-meta.json contract', function () { assert.ok(fs.results.has('/scene/1_0/meta.json')); }); + it('matches errors to lodLevels when trailing structural LODs are empty', async function () { + const { meta } = await writeScene([1, 0], 0); + assert.strictEqual(meta.lodLevels, 1); + assert.deepStrictEqual(meta.tree.errors, [0]); + assert.doesNotThrow(() => collectFilesByLod(meta, '/scene/lod-meta.json')); + }); + + it('prunes KNN traversal when a target LOD contains fewer than k splats', function () { + let farLeafVisits = 0; + const nearLeaf = { + count: 2, + aabb: { min: [0, 0, 0], max: [0, 0, 0] }, + indices: Uint32Array.of(0, 2) + }; + const farLeaf = { + count: 1, + aabb: { min: [1000, 0, 0], max: [1000, 0, 0] }, + get indices() { + farLeafVisits++; + return Uint32Array.of(1); + } + }; + const root = { + count: 3, + aabb: { min: [0, 0, 0], max: [1000, 0, 0] }, + left: nearLeaf, + right: farLeaf + }; + const slim = { + x: Float32Array.of(0, 1000, 0), + y: new Float32Array(3), + z: new Float32Array(3) + }; + + const [result] = findNearest(root, slim, Uint32Array.of(0), Int32Array.of(2, 3), 4); + assert.deepStrictEqual(result, Int32Array.of(2, -1, -1, -1)); + assert.strictEqual(farLeafVisits, 0); + }); + it('includes stored spherical harmonics in the LOD error', async function () { const fs = new MemoryFileSystem(); await writeLodSource({ @@ -235,6 +279,32 @@ describe('writeLodSource: lod-meta.json contract', function () { assert.ok(Math.abs(errors[1] - 0.5) < 1e-6, `expected 0.5, got ${errors[1]}`); }); + it('keeps the error table readable when a NaN input poisons a level', async function () { + // A NaN scale makes that splat's mass NaN, which reaches the coverage + // term; left alone it serializes as `null` and the reader rejects it. + const fs = new MemoryFileSystem(); + await writeLodSource({ + filename: '/scene/lod-meta.json', + mainSource: stackLods([ + dataTableToChunkSource(makeSplatTable([{}, { scale: NaN }]), 1 << 20), + dataTableToChunkSource(makeSplatTable([{}]), 1 << 20) + ]), + envSource: null, + iterations: 1, + chunkCount: 1, + chunkExtent: 16 + }, fs); + + const text = new TextDecoder().decode(fs.results.get('/scene/lod-meta.json')); + assert.ok(!text.includes('null'), 'no null slipped into the meta'); + const meta = JSON.parse(text); + assert.ok( + meta.tree.errors.every(error => Number.isFinite(error) && error >= 0), + `expected finite non-negative errors, got ${meta.tree.errors}` + ); + assert.doesNotThrow(() => collectFilesByLod(meta, '/scene/lod-meta.json')); + }); + it('keeps the error table monotone across levels', async function () { // level 2 matches a level-0 splat exactly while level 1 sits between // both, so the raw errors would rank the coarser level as the better one From 8cea42cc63a6732f77562e85359ea4f38242c705 Mon Sep 17 00:00:00 2001 From: Donovan Hutchence Date: Mon, 10 Aug 2026 12:00:13 +0200 Subject: [PATCH 4/5] latest --- src/lib/writers/write-lod.ts | 71 ++++++++++++++++++++++++++++++------ test/write-lod.test.mjs | 47 +++++++++++++----------- 2 files changed, 85 insertions(+), 33 deletions(-) diff --git a/src/lib/writers/write-lod.ts b/src/lib/writers/write-lod.ts index e148a8a5..a360bb1b 100644 --- a/src/lib/writers/write-lod.ts +++ b/src/lib/writers/write-lod.ts @@ -165,6 +165,57 @@ const accumulateBound = ( } }; +/** + * Reject a batch of gaussians whose geometry is not finite. The bounds pass reads + * every gaussian's geometric record once, so this is the one place the LOD writer + * sees the whole scene: everything downstream — the error metric above all, whose + * footprint mass runs through `sigmoid`/`exp` — may then assume finite input + * rather than each stage carrying its own opinion about invalid data. + * + * The rules mirror `filterNaNRows`, including its two deliberate exceptions + * (`scale_*` may be `-Infinity`, `opacity` may be `+Infinity`, both harmless + * here), so anything `--filter-nan` keeps is accepted. + * + * @param pos - Packed xyz for the batch. + * @param geo - Packed 8-float geometric records for the batch. + * @param count - Gaussians in the batch. + * @param rows - Rows local to `lod`, indexed from `offset`, for error messages. + * @param offset - Index of the batch's first row within `rows`. + * @param lod - Structural LOD the batch was read from. + */ +const assertFiniteGeometry = ( + pos: Float32Array, geo: Float32Array, count: number, + rows: Uint32Array, offset: number, lod: number +): void => { + const reject = (i: number, what: string) => { + throw new Error( + `LOD ${lod} gaussian ${rows[offset + i]} has ${what}; ` + + 'run --filter-nan to drop invalid gaussians before writing LODs' + ); + }; + + for (let i = 0; i < count; i++) { + const p = i * 3; + if (!isFinite(pos[p]) || !isFinite(pos[p + 1]) || !isFinite(pos[p + 2])) { + reject(i, 'a non-finite position'); + } + + const o = i * 8; + if (!isFinite(geo[o]) || !isFinite(geo[o + 1]) || !isFinite(geo[o + 2]) || !isFinite(geo[o + 3])) { + reject(i, 'a non-finite rotation'); + } + if (geo[o] === 0 && geo[o + 1] === 0 && geo[o + 2] === 0 && geo[o + 3] === 0) { + reject(i, 'a zero-norm rotation'); + } + for (let e = 4; e <= 6; e++) { + const v = geo[o + e]; + if (!isFinite(v) && v !== -Infinity) reject(i, 'a non-finite scale'); + } + const opacity = geo[o + 7]; + if (!isFinite(opacity) && opacity !== Infinity) reject(i, 'a non-finite opacity'); + } +}; + // Per-leaf ellipsoid AABB, computed per structural LOD. Positions are resident, // but rotation/scale are gathered from the source by index so the geometric layer // is never wholly resident — the bounds-pass analog of the per-unit heavy gather. @@ -190,11 +241,15 @@ const calcBound = async ( const pos = pool.acquire('position', layouts.position!, count); const geo = pool.acquire('geometric', layouts.geometric!, count); await source.read({ indices: local, indexOffset: off, count, lod: lodValue, position: pos, geometric: geo }); + // position is full-stride packed xyz — read the pool buffer in place + // rather than copying it out per batch + const posBatch = new Float32Array(pos.data, 0, count * 3); + assertFiniteGeometry( + posBatch, new Float32Array(geo.data, 0, count * 8), count, local, off, lodValue + ); accumulateBound( min, max, - // position is full-stride packed xyz — read the pool buffer - // in place rather than copying it out per batch - new Float32Array(pos.data, 0, count * 3), + posBatch, geo.field('rotation') as Float32Array, geo.field('scale') as Float32Array, count @@ -701,17 +756,9 @@ const calcErrors = async ( // budget balancer ranks transitions by error reduction per splat *across* // nodes, so a coarser level must never advertise less error than the finer // one it stands in for. - // - // A NaN scale or opacity in the input poisons that splat's mass, and a NaN - // mass reaches the coverage term (unlike a NaN error, which the nearest - // match discards). Left alone it serializes as `null`, and a single - // non-finite error drops the engine back to distance-based allocation for the - // whole scene, so treat a level whose error is not a number as no better than - // the one before it. let previous = 0; for (const lod of [...bins.keys()].sort((a, b) => a - b)) { - const error = errors[lod]; - previous = errors[lod] = Number.isFinite(error) ? Math.max(error, previous) : previous; + previous = errors[lod] = Math.max(errors[lod], previous); } return errors; diff --git a/test/write-lod.test.mjs b/test/write-lod.test.mjs index bbe3b035..eaf1ede3 100644 --- a/test/write-lod.test.mjs +++ b/test/write-lod.test.mjs @@ -279,30 +279,35 @@ describe('writeLodSource: lod-meta.json contract', function () { assert.ok(Math.abs(errors[1] - 0.5) < 1e-6, `expected 0.5, got ${errors[1]}`); }); - it('keeps the error table readable when a NaN input poisons a level', async function () { - // A NaN scale makes that splat's mass NaN, which reaches the coverage - // term; left alone it serializes as `null` and the reader rejects it. - const fs = new MemoryFileSystem(); - await writeLodSource({ - filename: '/scene/lod-meta.json', - mainSource: stackLods([ - dataTableToChunkSource(makeSplatTable([{}, { scale: NaN }]), 1 << 20), - dataTableToChunkSource(makeSplatTable([{}]), 1 << 20) - ]), - envSource: null, - iterations: 1, - chunkCount: 1, - chunkExtent: 16 - }, fs); + // Non-finite geometry is rejected up front rather than tolerated: a NaN scale + // or opacity would otherwise poison that splat's footprint mass, and the error + // table would quietly claim a coarse level costs nothing. + const rejects = [ + ['a NaN scale', { scale: NaN }, /non-finite scale/], + ['a NaN opacity', { opacity: NaN }, /non-finite opacity/], + ['a NaN position', { x: NaN }, /non-finite position/], + ['a NaN rotation', { rot_0: NaN }, /non-finite rotation/], + ['a zero-norm rotation', { rot_0: 0 }, /zero-norm rotation/] + ]; + + for (const [label, splat, expected] of rejects) { + it(`refuses to write LODs for input with ${label}`, async function () { + await assert.rejects(() => writeErrors([[{}, splat], [{}]]), (err) => { + assert.match(err.message, expected); + assert.match(err.message, /--filter-nan/); + return true; + }); + }); + } - const text = new TextDecoder().decode(fs.results.get('/scene/lod-meta.json')); - assert.ok(!text.includes('null'), 'no null slipped into the meta'); - const meta = JSON.parse(text); + it('accepts the non-finite values --filter-nan deliberately keeps', async function () { + // a flat splat (scale -Inf) and a fully opaque one (opacity +Inf) survive + // filterNaN, so the writer must not reject them + const errors = await writeErrors([[{}, { scale: -Infinity }, { opacity: Infinity }], [{}]]); assert.ok( - meta.tree.errors.every(error => Number.isFinite(error) && error >= 0), - `expected finite non-negative errors, got ${meta.tree.errors}` + errors.every(error => Number.isFinite(error) && error >= 0), + `expected finite non-negative errors, got ${errors}` ); - assert.doesNotThrow(() => collectFilesByLod(meta, '/scene/lod-meta.json')); }); it('keeps the error table monotone across levels', async function () { From 9206d0d9082d8265ff09d68543e297c1f5f76285 Mon Sep 17 00:00:00 2001 From: Donovan Hutchence Date: Mon, 10 Aug 2026 15:25:32 +0200 Subject: [PATCH 5/5] latest --- src/lib/writers/write-lod.ts | 52 +++++++++++++++++++++++++++++------- test/write-lod.test.mjs | 7 ++++- 2 files changed, 49 insertions(+), 10 deletions(-) diff --git a/src/lib/writers/write-lod.ts b/src/lib/writers/write-lod.ts index a360bb1b..8f2b6ebe 100644 --- a/src/lib/writers/write-lod.ts +++ b/src/lib/writers/write-lod.ts @@ -165,12 +165,18 @@ const accumulateBound = ( } }; +const invalidGaussian = (lod: number, row: number, what: string): Error => new Error( + `LOD ${lod} gaussian ${row} has ${what}; ` + + 'run --filter-nan to drop invalid gaussians before writing LODs' +); + /** * Reject a batch of gaussians whose geometry is not finite. The bounds pass reads - * every gaussian's geometric record once, so this is the one place the LOD writer - * sees the whole scene: everything downstream — the error metric above all, whose - * footprint mass runs through `sigmoid`/`exp` — may then assume finite input - * rather than each stage carrying its own opinion about invalid data. + * every gaussian's geometric record once, so this covers the whole scene: + * everything downstream — the error metric above all, whose footprint mass runs + * through `sigmoid`/`exp` — may then assume finite input rather than each stage + * carrying its own opinion about invalid data. {@link assertFiniteColor} does the + * same for the layer this pass does not read. * * The rules mirror `filterNaNRows`, including its two deliberate exceptions * (`scale_*` may be `-Infinity`, `opacity` may be `+Infinity`, both harmless @@ -188,10 +194,7 @@ const assertFiniteGeometry = ( rows: Uint32Array, offset: number, lod: number ): void => { const reject = (i: number, what: string) => { - throw new Error( - `LOD ${lod} gaussian ${rows[offset + i]} has ${what}; ` + - 'run --filter-nan to drop invalid gaussians before writing LODs' - ); + throw invalidGaussian(lod, rows[offset + i], what); }; for (let i = 0; i < count; i++) { @@ -466,6 +469,35 @@ const indexOfSorted = (values: Uint32Array, value: number): number => { return -1; }; +/** + * Reject a batch of gaussians whose color or stored SH is not finite — the + * companion to {@link assertFiniteGeometry} for the layer the bounds pass does not + * read. A single non-finite coefficient makes {@link splatError} return `NaN` for + * every pair the splat takes part in, and `directional` discards `NaN` matches, so + * a level would report *less* error the more of it is broken. Coverage is again + * the whole scene: every gaussian of a leaf is gathered as its own level's query. + * + * @param color - Packed color/SH records for the batch. + * @param count - Gaussians in the batch. + * @param colorDim - Floats per gaussian (3 + stored SH). + * @param rows - Rows local to `lod`, indexed from `offset`, for error messages. + * @param offset - Index of the batch's first row within `rows`. + * @param lod - Structural LOD the batch was read from. + */ +const assertFiniteColor = ( + color: Float32Array, count: number, colorDim: number, + rows: Uint32Array, offset: number, lod: number +): void => { + for (let i = 0; i < count; i++) { + const base = i * colorDim; + for (let c = 0; c < colorDim; c++) { + if (!isFinite(color[base + c])) { + throw invalidGaussian(lod, rows[offset + i], 'a non-finite color or SH coefficient'); + } + } + } +}; + const gatherView = async ( source: ChunkSource, pool: ChunkDataPool, slim: SlimColumns, globals: Uint32Array, lod: number, base: number @@ -492,8 +524,10 @@ const gatherView = async ( const geo = pool.acquire('geometric', layouts.geometric!, count); const color = pool.acquire('color', layouts.color!, count); await source.read({ indices: local, indexOffset: off, count, lod, geometric: geo, color }); + const colorBatch = new Float32Array(color.data, 0, count * colorDim); + assertFiniteColor(colorBatch, count, colorDim, local, off, lod); view.geo.set(new Float32Array(geo.data, 0, count * 8), off * 8); - view.color.set(new Float32Array(color.data, 0, count * colorDim), off * colorDim); + view.color.set(colorBatch, off * colorDim); geo.release(); color.release(); } diff --git a/test/write-lod.test.mjs b/test/write-lod.test.mjs index eaf1ede3..0e4b8f26 100644 --- a/test/write-lod.test.mjs +++ b/test/write-lod.test.mjs @@ -287,7 +287,12 @@ describe('writeLodSource: lod-meta.json contract', function () { ['a NaN opacity', { opacity: NaN }, /non-finite opacity/], ['a NaN position', { x: NaN }, /non-finite position/], ['a NaN rotation', { rot_0: NaN }, /non-finite rotation/], - ['a zero-norm rotation', { rot_0: 0 }, /zero-norm rotation/] + ['a zero-norm rotation', { rot_0: 0 }, /zero-norm rotation/], + // a NaN colour makes splatError NaN for every pair the splat is in, and + // directional discards NaN matches — so the more of a level is broken, the + // less error it reports. Left unchecked, two displaced levels whose colours + // are NaN report error 0 and the engine drops the finer one. + ['a NaN color', { f_dc: NaN }, /non-finite color or SH/] ]; for (const [label, splat, expected] of rejects) {