Implement float32 sorting option (#129)

* Add commented benchmarkSort to compare JS vs Wasm sorting.

* Add 32-bit float sort via SparkViewpooint.sort32, using 2-pass radix-65536 sort. Turn on sort32 by default for examples/editor. Implemented Rust sort and JS sort, updated benchmarking code to include float16 and float32 sort.

* Add sort32 to docs/spark-viewpoint.md.

src/SparkViewpoint.ts

@@ -18,6 +18,7 @@ import {
   dyno,
   dynoBlock,
   dynoConst,
+  floatBitsToUint,
   mul,
   packHalf2x16,
   readPackedSplat,
@@ -117,6 +118,11 @@ export type SparkViewpointOptions = {
    * @default false
    */
   sort360?: boolean;
+  /*
+   * Set this to true to sort with float32 precision with two-pass sort.
+   * @default true
+   */
+  sort32?: boolean;
 };
 
 // A SparkViewpoint is created from and tied to a SparkRenderer, and represents
@@ -149,6 +155,7 @@ export class SparkViewpoint {
   sortCoorient?: boolean;
   depthBias?: number;
   sort360?: boolean;
+  sort32?: boolean;
 
   display: {
     accumulator: SplatAccumulator;
@@ -164,7 +171,8 @@ export class SparkViewpoint {
   } | null = null;
   private sortingCheck = false;
 
-  private readback: Uint16Array = new Uint16Array(0);
+  private readback16: Uint16Array = new Uint16Array(0);
+  private readback32: Uint32Array = new Uint32Array(0);
   private orderingFreelist: FreeList<Uint32Array, number>;
 
   constructor(options: SparkViewpointOptions & { spark: SparkRenderer }) {
@@ -209,6 +217,7 @@ export class SparkViewpoint {
     this.sortCoorient = options.sortCoorient;
     this.depthBias = options.depthBias;
     this.sort360 = options.sort360;
+    this.sort32 = options.sort32;
 
     this.orderingFreelist = new FreeList({
       allocate: (maxSplats) => new Uint32Array(maxSplats),
@@ -557,6 +566,7 @@ export class SparkViewpoint {
       const {
         reader,
         doubleSortReader,
+        sort32Reader,
         dynoSortRadial,
         dynoOrigin,
         dynoDirection,
@@ -564,8 +574,16 @@ export class SparkViewpoint {
         dynoSort360,
         dynoSplats,
       } = SparkViewpoint.makeSorter();
-      const halfMaxSplats = Math.ceil(maxSplats / 2);
-      this.readback = reader.ensureBuffer(halfMaxSplats, this.readback);
+      const sort32 = this.sort32 ?? false;
+      let readback: Uint16Array | Uint32Array;
+      if (sort32) {
+        this.readback32 = reader.ensureBuffer(maxSplats, this.readback32);
+        readback = this.readback32;
+      } else {
+        const halfMaxSplats = Math.ceil(maxSplats / 2);
+        this.readback16 = reader.ensureBuffer(halfMaxSplats, this.readback16);
+        readback = this.readback16;
+      }
 
       const worldToOrigin = accumulator.toWorld.clone().invert();
       const viewToOrigin = viewToWorld.clone().premultiply(worldToOrigin);
@@ -581,25 +599,33 @@ export class SparkViewpoint {
       dynoSort360.value = this.sort360 ?? false;
       dynoSplats.packedSplats = accumulator.splats;
 
+      const sortReader = sort32 ? sort32Reader : doubleSortReader;
+      const count = sort32 ? numSplats : Math.ceil(numSplats / 2);
       await reader.renderReadback({
         renderer: this.spark.renderer,
-        reader: doubleSortReader,
-        count: Math.ceil(numSplats / 2),
-        readback: this.readback,
+        reader: sortReader,
+        count,
+        readback,
       });
 
       const result = (await withWorker(async (worker) => {
-        return worker.call("sortDoubleSplats", {
+        const rpcName = sort32 ? "sort32Splats" : "sortDoubleSplats";
+        return worker.call(rpcName, {
+          maxSplats,
           numSplats,
-          readback: this.readback,
+          readback,
           ordering,
         });
       })) as {
-        readback: Uint16Array;
+        readback: Uint16Array | Uint32Array;
         ordering: Uint32Array;
         activeSplats: number;
       };
-      this.readback = result.readback;
+      if (sort32) {
+        this.readback32 = result.readback as Uint32Array;
+      } else {
+        this.readback16 = result.readback as Uint16Array;
+      }
       ordering = result.ordering;
       activeSplats = result.activeSplats;
     }
@@ -669,6 +695,7 @@ export class SparkViewpoint {
     dynoSplats: DynoPackedSplats;
     reader: Readback;
     doubleSortReader: DynoBlock<{ index: "int" }, { rgba8: "vec4" }>;
+    sort32Reader: DynoBlock<{ index: "int" }, { rgba8: "vec4" }>;
   } | null = null;
 
   private static makeSorter() {
@@ -716,6 +743,28 @@ export class SparkViewpoint {
         },
       );
 
+      const sort32Reader = dynoBlock(
+        { index: "int" },
+        { rgba8: "vec4" },
+        ({ index }) => {
+          if (!index) {
+            throw new Error("No index");
+          }
+          const sortParams = {
+            sortRadial: dynoSortRadial,
+            sortOrigin: dynoOrigin,
+            sortDirection: dynoDirection,
+            sortDepthBias: dynoDepthBias,
+            sort360: dynoSort360,
+          };
+
+          const gsplat = readPackedSplat(dynoSplats, index);
+          const metric = computeSortMetric({ gsplat, ...sortParams });
+          const rgba8 = uintToRgba8(floatBitsToUint(metric));
+          return { rgba8 };
+        },
+      );
+
       SparkViewpoint.dynos = {
         dynoSortRadial,
         dynoOrigin,
@@ -725,6 +774,7 @@ export class SparkViewpoint {
         dynoSplats,
         reader,
         doubleSortReader,
+        sort32Reader,
       };
     }
     return SparkViewpoint.dynos;

src/worker.ts

@@ -1,4 +1,4 @@
-import init_wasm, { sort_splats } from "spark-internal-rs";
+import init_wasm, { sort_splats, sort32_splats } from "spark-internal-rs";
 import type { PcSogsJson, TranscodeSpzInput } from "./SplatLoader";
 import { unpackAntiSplat } from "./antisplat";
 import { WASM_SPLAT_SORT } from "./defines";
@@ -18,6 +18,7 @@ import {
   setPackedSplatQuat,
   setPackedSplatRgb,
   setPackedSplatScales,
+  toHalf,
 } from "./utils";
 
 // WebWorker for Spark's background CPU tasks, such as Gsplat file decoding
@@ -134,11 +135,6 @@ async function onMessage(event: MessageEvent) {
           readback: Uint16Array;
           ordering: Uint32Array;
         };
-        result = {
-          id,
-          readback,
-          ordering,
-        };
         if (WASM_SPLAT_SORT) {
           result = {
             id,
@@ -155,6 +151,31 @@ async function onMessage(event: MessageEvent) {
         }
         break;
       }
+      case "sort32Splats": {
+        const { maxSplats, numSplats, readback, ordering } = args as {
+          maxSplats: number;
+          numSplats: number;
+          readback: Uint32Array;
+          ordering: Uint32Array;
+        };
+        // Benchmark sort
+        // benchmarkSort(numSplats, readback, ordering);
+        if (WASM_SPLAT_SORT) {
+          result = {
+            id,
+            readback,
+            ordering,
+            activeSplats: sort32_splats(numSplats, readback, ordering),
+          };
+        } else {
+          result = {
+            id,
+            readback,
+            ...sort32Splats({ maxSplats, numSplats, readback, ordering }),
+          };
+        }
+        break;
+      }
       case "transcodeSpz": {
         const input = args as TranscodeSpzInput;
         const spzBytes = await transcodeSpz(input);
@@ -171,6 +192,7 @@ async function onMessage(event: MessageEvent) {
     }
   } catch (e) {
     error = e;
+    console.error(error);
   }
 
   // Send the result or error back to the main thread, making sure to transfer any ArrayBuffers
@@ -180,6 +202,82 @@ async function onMessage(event: MessageEvent) {
   );
 }
 
+function benchmarkSort(
+  numSplats: number,
+  readback32: Uint32Array,
+  ordering: Uint32Array,
+) {
+  if (numSplats > 0) {
+    console.log("Running sort benchmark");
+    const readbackF32 = new Float32Array(readback32.buffer);
+    const readback16 = new Uint16Array(readback32.length);
+    for (let i = 0; i < numSplats; ++i) {
+      readback16[i] = toHalf(readbackF32[i]);
+    }
+
+    const WARMUP = 10;
+    for (let i = 0; i < WARMUP; ++i) {
+      const activeSplats = sort_splats(numSplats, readback16, ordering);
+      const activeSplats32 = sort32_splats(numSplats, readback32, ordering);
+      const results = sortDoubleSplats({
+        numSplats,
+        readback: readback16,
+        ordering,
+      });
+      const results32 = sort32Splats({
+        maxSplats: numSplats,
+        numSplats,
+        readback: readback32,
+        ordering,
+      });
+    }
+
+    const TIMING_SAMPLES = 1000;
+    let start: number;
+
+    start = performance.now();
+    for (let i = 0; i < TIMING_SAMPLES; ++i) {
+      const activeSplats = sort_splats(numSplats, readback16, ordering);
+    }
+    const wasmTime = (performance.now() - start) / TIMING_SAMPLES;
+
+    start = performance.now();
+    for (let i = 0; i < TIMING_SAMPLES; ++i) {
+      const results = sortDoubleSplats({
+        numSplats,
+        readback: readback16,
+        ordering,
+      });
+    }
+    const jsTime = (performance.now() - start) / TIMING_SAMPLES;
+
+    console.log(
+      `JS: ${jsTime} ms, WASM: ${wasmTime} ms, numSplats: ${numSplats}`,
+    );
+
+    start = performance.now();
+    for (let i = 0; i < TIMING_SAMPLES; ++i) {
+      const activeSplats32 = sort32_splats(numSplats, readback32, ordering);
+    }
+    const wasm32Time = (performance.now() - start) / TIMING_SAMPLES;
+
+    start = performance.now();
+    for (let i = 0; i < TIMING_SAMPLES; ++i) {
+      const results = sort32Splats({
+        maxSplats: numSplats,
+        numSplats,
+        readback: readback32,
+        ordering,
+      });
+    }
+    const js32Time = (performance.now() - start) / TIMING_SAMPLES;
+
+    console.log(
+      `JS32: ${js32Time} ms, WASM32: ${wasm32Time} ms, numSplats: ${numSplats}`,
+    );
+  }
+}
+
 async function unpackPly({
   packedArray,
   fileBytes,
@@ -308,9 +406,9 @@ function unpackSpz(fileBytes: Uint8Array): {
 }
 
 // Array of buckets for sorting float16 distances with range [0, DEPTH_INFINITY].
-const DEPTH_INFINITY = 0x7c00;
-const DEPTH_SIZE = DEPTH_INFINITY + 1;
-let depthArray: Uint32Array | null = null;
+const DEPTH_INFINITY_F16 = 0x7c00;
+const DEPTH_SIZE_16 = DEPTH_INFINITY_F16 + 1;
+let depthArray16: Uint32Array | null = null;
 
 function sortSplats({
   totalSplats,
@@ -323,10 +421,10 @@ function sortSplats({
   // Sort totalSplats Gsplats, each with 4 bytes of readback, and outputs Uint32Array
   // of indices from most distant to nearest. Each 4 bytes encode a float16 distance
   // and unused high bytes.
-  if (!depthArray) {
-    depthArray = new Uint32Array(DEPTH_SIZE);
+  if (!depthArray16) {
+    depthArray16 = new Uint32Array(DEPTH_SIZE_16);
   }
-  depthArray.fill(0);
+  depthArray16.fill(0);
 
   const readbackUint32 = readback.map((layer) => new Uint32Array(layer.buffer));
   const layerSize = readbackUint32[0].length;
@@ -338,17 +436,17 @@ function sortSplats({
     const layerSplats = Math.min(readbackLayer.length, totalSplats - layerBase);
     for (let i = 0; i < layerSplats; ++i) {
       const pri = readbackLayer[i] & 0x7fff;
-      if (pri < DEPTH_INFINITY) {
-        depthArray[pri] += 1;
+      if (pri < DEPTH_INFINITY_F16) {
+        depthArray16[pri] += 1;
       }
     }
     layerBase += layerSplats;
   }
 
   let activeSplats = 0;
-  for (let j = 0; j < DEPTH_SIZE; ++j) {
-    const nextIndex = activeSplats + depthArray[j];
-    depthArray[j] = activeSplats;
+  for (let j = 0; j < DEPTH_SIZE_16; ++j) {
+    const nextIndex = activeSplats + depthArray16[j];
+    depthArray16[j] = activeSplats;
     activeSplats = nextIndex;
   }
 
@@ -358,16 +456,16 @@ function sortSplats({
     const layerSplats = Math.min(readbackLayer.length, totalSplats - layerBase);
     for (let i = 0; i < layerSplats; ++i) {
       const pri = readbackLayer[i] & 0x7fff;
-      if (pri < DEPTH_INFINITY) {
-        ordering[depthArray[pri]] = layerBase + i;
-        depthArray[pri] += 1;
+      if (pri < DEPTH_INFINITY_F16) {
+        ordering[depthArray16[pri]] = layerBase + i;
+        depthArray16[pri] += 1;
       }
     }
     layerBase += layerSplats;
   }
-  if (depthArray[DEPTH_SIZE - 1] !== activeSplats) {
+  if (depthArray16[DEPTH_SIZE_16 - 1] !== activeSplats) {
     throw new Error(
-      `Expected ${activeSplats} active splats but got ${depthArray[DEPTH_SIZE - 1]}`,
+      `Expected ${activeSplats} active splats but got ${depthArray16[DEPTH_SIZE_16 - 1]}`,
     );
   }
 
@@ -385,16 +483,16 @@ function sortDoubleSplats({
   ordering: Uint32Array;
 } {
   // Ensure depthArray is allocated and zeroed out for our buckets.
-  if (!depthArray) {
-    depthArray = new Uint32Array(DEPTH_SIZE);
+  if (!depthArray16) {
+    depthArray16 = new Uint32Array(DEPTH_SIZE_16);
   }
-  depthArray.fill(0);
+  depthArray16.fill(0);
 
   // Count the number of splats in each bucket (cull Gsplats at infinity).
   for (let i = 0; i < numSplats; ++i) {
     const pri = readback[i];
-    if (pri < DEPTH_INFINITY) {
-      depthArray[pri] += 1;
+    if (pri < DEPTH_INFINITY_F16) {
+      depthArray16[pri] += 1;
     }
   }
 
@@ -402,9 +500,9 @@ function sortDoubleSplats({
   // total number of active (non-infinity) splats, going in reverse order
   // because we want most distant Gsplats to be first in the output array.
   let activeSplats = 0;
-  for (let j = DEPTH_INFINITY - 1; j >= 0; --j) {
-    const nextIndex = activeSplats + depthArray[j];
-    depthArray[j] = activeSplats;
+  for (let j = DEPTH_INFINITY_F16 - 1; j >= 0; --j) {
+    const nextIndex = activeSplats + depthArray16[j];
+    depthArray16[j] = activeSplats;
     activeSplats = nextIndex;
   }
 
@@ -412,16 +510,106 @@ function sortDoubleSplats({
   // bucket order.
   for (let i = 0; i < numSplats; ++i) {
     const pri = readback[i];
-    if (pri < DEPTH_INFINITY) {
-      ordering[depthArray[pri]] = i;
-      depthArray[pri] += 1;
+    if (pri < DEPTH_INFINITY_F16) {
+      ordering[depthArray16[pri]] = i;
+      depthArray16[pri] += 1;
     }
   }
   // Sanity check that the end of the closest bucket is the same as
   // our total count of active splats (not at infinity).
-  if (depthArray[0] !== activeSplats) {
+  if (depthArray16[0] !== activeSplats) {
     throw new Error(
-      `Expected ${activeSplats} active splats but got ${depthArray[0]}`,
+      `Expected ${activeSplats} active splats but got ${depthArray16[0]}`,
+    );
+  }
+
+  return { activeSplats, ordering };
+}
+
+const DEPTH_INFINITY_F32 = 0x7f800000;
+let bucket16: Uint32Array | null = null;
+let scratchSplats: Uint32Array | null = null;
+
+// two-pass radix sort (base 65536) of 32-bit keys in readback,
+// but placing largest values first.
+function sort32Splats({
+  maxSplats,
+  numSplats,
+  readback, // Uint32Array of bit‑patterns
+  ordering, // Uint32Array to fill with sorted indices
+}: {
+  maxSplats: number;
+  numSplats: number;
+  readback: Uint32Array;
+  ordering: Uint32Array;
+}): { activeSplats: number; ordering: Uint32Array } {
+  const BASE = 1 << 16; // 65536
+
+  // allocate once
+  if (!bucket16) {
+    bucket16 = new Uint32Array(BASE);
+  }
+  if (!scratchSplats || scratchSplats.length < maxSplats) {
+    scratchSplats = new Uint32Array(maxSplats);
+  }
+
+  //
+  // ——— Pass #1: bucket by inv(lo 16 bits) ———
+  //
+  bucket16.fill(0);
+  for (let i = 0; i < numSplats; ++i) {
+    const key = readback[i];
+    if (key < DEPTH_INFINITY_F32) {
+      const inv = ~key >>> 0;
+      bucket16[inv & 0xffff] += 1;
+    }
+  }
+  // exclusive prefix‑sum → starting offsets
+  let total = 0;
+  for (let b = 0; b < BASE; ++b) {
+    const c = bucket16[b];
+    bucket16[b] = total;
+    total += c;
+  }
+  const activeSplats = total;
+
+  // scatter into scratch by low bits of inv
+  for (let i = 0; i < numSplats; ++i) {
+    const key = readback[i];
+    if (key < DEPTH_INFINITY_F32) {
+      const inv = ~key >>> 0;
+      scratchSplats[bucket16[inv & 0xffff]++] = i;
+    }
+  }
+
+  //
+  // ——— Pass #2: bucket by inv(hi 16 bits) ———
+  //
+  bucket16.fill(0);
+  for (let k = 0; k < activeSplats; ++k) {
+    const idx = scratchSplats[k];
+    const inv = ~readback[idx] >>> 0;
+    bucket16[inv >>> 16] += 1;
+  }
+  // exclusive prefix‑sum again
+  let sum = 0;
+  for (let b = 0; b < BASE; ++b) {
+    const c = bucket16[b];
+    bucket16[b] = sum;
+    sum += c;
+  }
+
+  // scatter into final ordering by high bits of inv
+  for (let k = 0; k < activeSplats; ++k) {
+    const idx = scratchSplats[k];
+    const inv = ~readback[idx] >>> 0;
+    ordering[bucket16[inv >>> 16]++] = idx;
+  }
+
+  // sanity‑check: the last bucket should have eaten all entries
+  if (bucket16[BASE - 1] !== activeSplats) {
+    throw new Error(
+      `Expected ${activeSplats} active splats but got ${bucket16[BASE - 1]}`,
     );
   }