Gaussian splat spherical harmonics support

CHANGES.md

@@ -23,6 +23,7 @@
 - Expand the CustomShader Sample to support real-time modification of CustomShader. [#12702](https://github.com/CesiumGS/cesium/pull/12702)
 - Add wrapR property to Sampler and Texture3D, to support the newly added third dimension wrap.[#12701](https://github.com/CesiumGS/cesium/pull/12701)
 - Added the ability to load a specific changeset for iTwin Mesh Exports using `ITwinData.createTilesetFromIModelId` [#12778](https://github.com/CesiumGS/cesium/issues/12778)
+- Added spherical harmonics support for Gaussian Splats. Supports degrees 1, 2, and 3 in the SPZ format.
 
 #### Deprecated :hourglass_flowing_sand:
 

packages/engine/Source/Scene/GaussianSplat3DTileContent.js

@@ -78,6 +78,20 @@ function GaussianSplat3DTileContent(loader, tileset, tile, resource) {
    * @private
    */
   this._transformed = false;
+
+  /**
+   * The degree of the spherical harmonics used for the Gaussian splats.
+   * @type {number}
+   * @private
+   */
+  this.shDegree = 0;
+
+  /**
+   * The number of spherical harmonic coefficients used for the Gaussian splats.
+   * @type {number}
+   * @private
+   */
+  this.shCoefficientCount = 0;
 }
 
 Object.defineProperties(GaussianSplat3DTileContent.prototype, {
@@ -311,6 +325,144 @@ Object.defineProperties(GaussianSplat3DTileContent.prototype, {
   },
 });
 
+/**
+ * Determine Spherical Harmonics degree and coefficient count from attributes
+ * @param {Array} attributes - The list of attributes.
+ * @returns {Object} An object containing the degree (l) and coefficient (n).
+ */
+function degreeAndCoefFromAttributes(attributes) {
+  const prefix = "_SH_DEGREE_";
+  const shAttributes = attributes.filter((attr) =>
+    attr.name.startsWith(prefix),
+  );
+
+  switch (shAttributes.length) {
+    default:
+    case 0:
+      return { l: 0, n: 0 };
+    case 3:
+      return { l: 1, n: 9 };
+    case 8:
+      return { l: 2, n: 24 };
+    case 15:
+      return { l: 3, n: 45 };
+  }
+}
+
+/**
+ * Converts a 32-bit floating point number to a 16-bit floating point number.
+ * @param {Float} float32 input
+ * @returns {number} Half precision float
+ * @private
+ */
+const buffer = new ArrayBuffer(4);
+const floatView = new Float32Array(buffer);
+const intView = new Uint32Array(buffer);
+
+function float32ToFloat16(float32) {
+  floatView[0] = float32;
+  const bits = intView[0];
+
+  const sign = (bits >> 31) & 0x1;
+  const exponent = (bits >> 23) & 0xff;
+  const mantissa = bits & 0x7fffff;
+
+  let half;
+
+  if (exponent === 0xff) {
+    half = (sign << 15) | (0x1f << 10) | (mantissa ? 0x200 : 0);
+  } else if (exponent === 0) {
+    half = sign << 15;
+  } else {
+    const newExponent = exponent - 127 + 15;
+    if (newExponent >= 31) {
+      half = (sign << 15) | (0x1f << 10);
+    } else if (newExponent <= 0) {
+      half = sign << 15;
+    } else {
+      half = (sign << 15) | (newExponent << 10) | (mantissa >>> 13);
+    }
+  }
+
+  return half;
+}
+
+/**
+ * Extracts the spherical harmonic degree and coefficient from the attribute name.
+ * @param {String} attribute - The attribute name.
+ * @returns {Object} An object containing the degree (l) and coefficient (n).
+ * @private
+ */
+function extractSHDegreeAndCoef(attribute) {
+  const prefix = "_SH_DEGREE_";
+  const separator = "_COEF_";
+
+  const lStart = prefix.length;
+  const coefIndex = attribute.indexOf(separator, lStart);
+
+  const l = parseInt(attribute.slice(lStart, coefIndex), 10);
+  const n = parseInt(attribute.slice(coefIndex + separator.length), 10);
+
+  return { l, n };
+}
+
+/**
+ * Packs spherical harmonic data into half-precision floats.
+ * @param {*} data - The input data to pack.
+ * @param {*} shDegree - The spherical harmonic degree.
+ * @returns {Uint32Array} - The packed data.
+ */
+function packSphericalHarmonicData(tileContent) {
+  const degree = tileContent.shDegree;
+  const coefs = tileContent.shCoefficientCount;
+  const totalLength = tileContent.pointsLength * (coefs * (2 / 3)); //3 packs into 2
+  const packedData = new Uint32Array(totalLength);
+
+  const shAttributes = tileContent.splatPrimitive.attributes.filter((attr) =>
+    attr.name.startsWith("_SH_DEGREE_"),
+  );
+  let stride = 0;
+  const base = [0, 9, 24];
+  switch (degree) {
+    case 1:
+      stride = 9;
+      break;
+    case 2:
+      stride = 24;
+      break;
+    case 3:
+      stride = 45;
+      break;
+  }
+  shAttributes.sort((a, b) => {
+    if (a.name < b.name) {
+      return -1;
+    }
+    if (a.name > b.name) {
+      return 1;
+    }
+
+    return 0;
+  });
+  const packedStride = stride * (2 / 3);
+  for (let i = 0; i < shAttributes.length; i++) {
+    const { l, n } = extractSHDegreeAndCoef(shAttributes[i].name);
+    for (let j = 0; j < tileContent.pointsLength; j++) {
+      //interleave the data
+      const packedBase = (base[l - 1] * 2) / 3;
+      const idx = j * packedStride + packedBase + n * 2;
+      const src = j * 3;
+      packedData[idx] =
+        float32ToFloat16(shAttributes[i].typedArray[src]) |
+        (float32ToFloat16(shAttributes[i].typedArray[src + 1]) << 16);
+      packedData[idx + 1] = float32ToFloat16(
+        shAttributes[i].typedArray[src + 2],
+      );
+    }
+  }
+  return packedData;
+}
+
 /**
  * Creates a new instance of {@link GaussianSplat3DTileContent} from a glTF or glb resource.
  *
@@ -412,6 +564,14 @@ GaussianSplat3DTileContent.prototype.update = function (primitive, frameState) {
       ).typedArray,
     );
 
+    const { l, n } = degreeAndCoefFromAttributes(
+      this.splatPrimitive.attributes,
+    );
+    this.shDegree = l;
+    this.shCoefficientCount = n;
+
+    this._packedShData = packSphericalHarmonicData(this);
+
     return;
   }
 

packages/engine/Source/Scene/GaussianSplatPrimitive.js

@@ -31,10 +31,12 @@ import Cartesian3 from "../Core/Cartesian3.js";
 import Quaternion from "../Core/Quaternion.js";
 import SplitDirection from "./SplitDirection.js";
 import destroyObject from "../Core/destroyObject.js";
+import ContextLimits from "../Renderer/ContextLimits.js";
 
 const scratchMatrix4A = new Matrix4();
 const scratchMatrix4B = new Matrix4();
 const scratchMatrix4C = new Matrix4();
+const scratchMatrix4D = new Matrix4();
 
 const GaussianSplatSortingState = {
   IDLE: 0,
@@ -44,6 +46,25 @@ const GaussianSplatSortingState = {
   ERROR: 4,
 };
 
+function createGaussianSplatSHTexture(context, shData) {
+  const texture = new Texture({
+    context: context,
+    source: {
+      width: shData.width,
+      height: shData.height,
+      arrayBufferView: shData.data,
+    },
+    preMultiplyAlpha: false,
+    skipColorSpaceConversion: true,
+    pixelFormat: PixelFormat.RG_INTEGER,
+    pixelDatatype: PixelDatatype.UNSIGNED_INT,
+    flipY: false,
+    sampler: Sampler.NEAREST,
+  });
+
+  return texture;
+}
+
 function createGaussianSplatTexture(context, splatTextureData) {
   return new Texture({
     context: context,
@@ -148,6 +169,14 @@ function GaussianSplatPrimitive(options) {
    * @see {@link GaussianSplatTextureGenerator}
    */
   this.gaussianSplatTexture = undefined;
+
+  /**
+   * The texture used to store the spherical harmonics coefficients for the Gaussian splats.
+   * @type {undefined|Texture}
+   * @private
+   */
+  this.gaussianSplatSHTexture = undefined;
+
   /**
    * The last width of the Gaussian splat texture.
    * This is used to track changes in the texture size and update the primitive accordingly.
@@ -601,12 +630,65 @@ GaussianSplatPrimitive.buildGSplatDrawCommand = function (
     ShaderDestination.VERTEX,
   );
 
+  shaderBuilder.addUniform("float", "u_shDegree", ShaderDestination.VERTEX);
+
+  shaderBuilder.addUniform("float", "u_splatScale", ShaderDestination.VERTEX);
+
+  shaderBuilder.addUniform(
+    "vec3",
+    "u_cameraPositionWC",
+    ShaderDestination.VERTEX,
+  );
+
+  shaderBuilder.addUniform(
+    "mat3",
+    "u_inverseModelRotation",
+    ShaderDestination.VERTEX,
+  );
+
   const uniformMap = renderResources.uniformMap;
 
   uniformMap.u_splatAttributeTexture = function () {
     return primitive.gaussianSplatTexture;
   };
 
+  if (primitive._shDegree > 0) {
+    shaderBuilder.addDefine(
+      "HAS_SPHERICAL_HARMONICS",
+      "1",
+      ShaderDestination.VERTEX,
+    );
+    shaderBuilder.addUniform(
+      "highp usampler2D",
+      "u_gaussianSplatSHTexture",
+      ShaderDestination.VERTEX,
+    );
+    uniformMap.u_gaussianSplatSHTexture = function () {
+      return primitive.gaussianSplatSHTexture;
+    };
+  }
+  uniformMap.u_shDegree = function () {
+    return primitive._shDegree;
+  };
+
+  uniformMap.u_cameraPositionWC = function () {
+    return Cartesian3.clone(frameState.camera.positionWC);
+  };
+
+  uniformMap.u_inverseModelRotation = function () {
+    const tileset = primitive._tileset;
+    const modelMatrix = Matrix4.multiply(
+      tileset.modelMatrix,
+      Matrix4.fromArray(tileset.root.transform),
+      scratchMatrix4A,
+    );
+    const inverseModelRotation = Matrix4.getRotation(
+      Matrix4.inverse(modelMatrix, scratchMatrix4C),
+      scratchMatrix4D,
+    );
+    return inverseModelRotation;
+  };
+
   uniformMap.u_splitDirection = function () {
     return primitive.splitDirection;
   };
@@ -756,6 +838,7 @@ GaussianSplatPrimitive.prototype.update = function (frameState) {
       this._scales = undefined;
       this._colors = undefined;
       this._indexes = undefined;
+      this._shData = undefined;
       this._needsGaussianSplatTexture = true;
       this._gaussianSplatTexturePending = false;
 
@@ -786,6 +869,31 @@ GaussianSplatPrimitive.prototype.update = function (frameState) {
         return aggregate;
       };
 
+      const aggregateShData = () => {
+        let offset = 0;
+        for (const tile of tiles) {
+          const shData = tile.content._packedShData;
+          if (tile.content.shDegree > 0) {
+            if (!defined(this._shData)) {
+              let coefs;
+              switch (tile.content.shDegree) {
+                case 1:
+                  coefs = 9;
+                  break;
+                case 2:
+                  coefs = 24;
+                  break;
+                case 3:
+                  coefs = 45;
+              }
+              this._shData = new Uint32Array(totalElements * (coefs * (2 / 3)));
+            }
+            this._shData.set(shData, offset);
+            offset += shData.length;
+          }
+        }
+      };
+
       this._positions = aggregateAttributeValues(
         ComponentDatatype.FLOAT,
         (splatPrimitive) =>
@@ -822,6 +930,9 @@ GaussianSplatPrimitive.prototype.update = function (frameState) {
           ),
       );
 
+      aggregateShData();
+      this._shDegree = tiles[0].content.shDegree;
+
       this._numSplats = totalElements;
       this.selectedTileLength = tileset._selectedTiles.length;
     }
@@ -833,6 +944,36 @@ GaussianSplatPrimitive.prototype.update = function (frameState) {
     if (this._needsGaussianSplatTexture) {
       if (!this._gaussianSplatTexturePending) {
         GaussianSplatPrimitive.generateSplatTexture(this, frameState);
+        if (defined(this._shData)) {
+          const oldTex = this.gaussianSplatSHTexture;
+          const width = ContextLimits.maximumTextureSize;
+          const dims = tileset._selectedTiles[0].content.shCoefficientCount / 3;
+          const splatsPerRow = Math.floor(width / dims);
+          const floatsPerRow = splatsPerRow * (dims * 2);
+          const texBuf = new Uint32Array(
+            width * Math.ceil(this._numSplats / splatsPerRow) * 2,
+          );
+
+          let dataIndex = 0;
+          for (let i = 0; dataIndex < this._shData.length; i += width * 2) {
+            texBuf.set(
+              this._shData.subarray(dataIndex, dataIndex + floatsPerRow),
+              i,
+            );
+            dataIndex += floatsPerRow;
+          }
+          this.gaussianSplatSHTexture = createGaussianSplatSHTexture(
+            frameState.context,
+            {
+              data: texBuf,
+              width: width,
+              height: Math.ceil(this._numSplats / splatsPerRow),
+            },
+          );
+          if (defined(oldTex)) {
+            oldTex.destroy();
+          }
+        }
       }
       return;
     }