diff --git a/Sources/CShaderTypes/ShaderTypes.h b/Sources/CShaderTypes/ShaderTypes.h
index 50a547c0..0d7b53eb 100644
--- a/Sources/CShaderTypes/ShaderTypes.h
+++ b/Sources/CShaderTypes/ShaderTypes.h
@@ -445,10 +445,12 @@ typedef enum {
 } ARTextureIndices;
 
 typedef enum{
-    gaussianSplatIndex,
+    gaussianEncodedSplatIndex,
     gaussianUniformIndex,
     gaussianNumberOfSplatsIndex,
     gaussianIndicesIndex,
+    gaussianVisibleIndicesIndex,
+    gaussianVisibleCountIndex,
 }GaussianDepthBufferIndices;
 
 typedef struct{
@@ -470,12 +472,25 @@ typedef struct{
     float _pad2;
 }EncodedGaussianSplat;
 
+// Higher-order SH coefficients are stored in a separate packed half buffer.
+// Layout per splat: R[1...n], G[1...n], B[1...n]. The DC coefficient remains
+// represented by EncodedGaussianSplat.color.
+typedef struct{
+    uint degree;
+    uint coefficientsPerChannel;
+    uint higherOrderCoefficientsPerSplat;
+    uint _pad0;
+}GaussianSHMetadata;
+
 typedef enum{
       gaussianTBDRRenderIndicesIndex = 0,
       gaussianTBDRRenderSplatIndex,
       gaussianTBDRRenderUniformIndex,
       gaussianTBDRRenderViewPortIndex,
       gaussianTBDRRenderReverseZIndex,
+      gaussianTBDRRenderSHIndex,
+      gaussianTBDRRenderSHMetadataIndex,
+      gaussianTBDRRenderLocalCameraIndex,
   }GaussianTBDRRenderBufferIndices;
 
 typedef enum{
diff --git a/Sources/UntoldEngine/ECS/Components.swift b/Sources/UntoldEngine/ECS/Components.swift
index f3949f7c..994418db 100644
--- a/Sources/UntoldEngine/ECS/Components.swift
+++ b/Sources/UntoldEngine/ECS/Components.swift
@@ -8,6 +8,7 @@
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
 
+import CShaderTypes
 import Foundation
 import Metal
 import MetalKit
@@ -89,9 +90,13 @@ public class EntitySceneChannelsComponent: Component {
 }
 
 public class GaussianComponent: Component {
-    var splatData: MTLBuffer?
     var encodedSplatData: MTLBuffer?
+    var sphericalHarmonicsData: MTLBuffer?
+    var sphericalHarmonicsMetadata: GaussianSHMetadata?
     var gaussianSortedIndices: MTLBuffer?
+    var gaussianVisibleIndices: MTLBuffer?
+    var gaussianVisibleCount: MTLBuffer?
+    var visibleSplatCountForRendering: UInt = 0
     public var spaceUniform: [MTLBuffer?] = Array(repeating: nil, count: totalPerMeshUniformBuffers())
     var splatCount: UInt = 0
 
diff --git a/Sources/UntoldEngine/Profiling/GaussianProfiling.swift b/Sources/UntoldEngine/Profiling/GaussianProfiling.swift
new file mode 100644
index 00000000..fead55ed
--- /dev/null
+++ b/Sources/UntoldEngine/Profiling/GaussianProfiling.swift
@@ -0,0 +1,112 @@
+//
+//  GaussianProfiling.swift
+//  UntoldEngine
+//
+// Copyright (C) Untold Engine Studios
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at https://mozilla.org/MPL/2.0/.
+
+import Foundation
+import QuartzCore
+
+struct GaussianProfileTotals {
+    var entityCount: Int = 0
+    var splatCount: Int = 0
+    var drawCallCount: Int = 0
+    var dispatchCount: Int = 0
+    var radixPassCount: Int = 0
+    var encodedBytes: Int = 0
+    var sortedIndexBytes: Int = 0
+    var visibleIndexBytes: Int = 0
+    var visibleCountBytes: Int = 0
+    var sphericalHarmonicsBytes: Int = 0
+    var uniformBytes: Int = 0
+    var scratchBytes: Int = 0
+    var maxSphericalHarmonicsDegree: UInt32 = 0
+    var higherOrderCoefficientsPerSplat: UInt32 = 0
+
+    var totalResidentBytes: Int {
+        encodedBytes + sortedIndexBytes + visibleIndexBytes + visibleCountBytes + sphericalHarmonicsBytes + uniformBytes + scratchBytes
+    }
+
+    mutating func include(component: GaussianComponent) {
+        entityCount += 1
+        splatCount += Int(component.splatCount)
+        encodedBytes += component.encodedSplatData?.length ?? 0
+        sortedIndexBytes += component.gaussianSortedIndices?.length ?? 0
+        visibleIndexBytes += component.gaussianVisibleIndices?.length ?? 0
+        visibleCountBytes += component.gaussianVisibleCount?.length ?? 0
+        if let shBuffer = component.sphericalHarmonicsData {
+            sphericalHarmonicsBytes += shBuffer.length
+        }
+        uniformBytes += component.spaceUniform.reduce(0) { total, buffer in
+            total + (buffer?.length ?? 0)
+        }
+        if let metadata = component.sphericalHarmonicsMetadata {
+            maxSphericalHarmonicsDegree = max(maxSphericalHarmonicsDegree, metadata.degree)
+            higherOrderCoefficientsPerSplat = max(
+                higherOrderCoefficientsPerSplat,
+                metadata.higherOrderCoefficientsPerSplat
+            )
+        }
+    }
+}
+
+@inline(__always)
+func gaussianProfilingStartTime() -> CFTimeInterval? {
+    guard Logger.isEnabled(category: .gaussian) else { return nil }
+    return CACurrentMediaTime()
+}
+
+func logGaussianProfile(
+    stage: String,
+    startTime: CFTimeInterval?,
+    totals: GaussianProfileTotals,
+    extra: String = ""
+) {
+    guard let startTime else { return }
+
+    let elapsedMs = (CACurrentMediaTime() - startTime) * 1000.0
+    let bytesPerSplat = totals.splatCount > 0
+        ? Double(totals.totalResidentBytes) / Double(totals.splatCount)
+        : 0.0
+    let suffix = extra.isEmpty ? "" : " \(extra)"
+
+    Logger.log(
+        message: String(
+            format: "[Gaussian][%@] cpuEncodeMs=%.3f entities=%d splats=%d draws=%d dispatches=%d radixPasses=%d shDegree=%u shRestCoeffsPerSplat=%u memory=%@ bytesPerSplat=%.1f encoded=%@ sorted=%@ visible=%@ sh=%@ uniforms=%@ scratch=%@%@",
+            stage,
+            elapsedMs,
+            totals.entityCount,
+            totals.splatCount,
+            totals.drawCallCount,
+            totals.dispatchCount,
+            totals.radixPassCount,
+            totals.maxSphericalHarmonicsDegree,
+            totals.higherOrderCoefficientsPerSplat,
+            gaussianFormatBytes(totals.totalResidentBytes),
+            bytesPerSplat,
+            gaussianFormatBytes(totals.encodedBytes),
+            gaussianFormatBytes(totals.sortedIndexBytes),
+            gaussianFormatBytes(totals.visibleIndexBytes + totals.visibleCountBytes),
+            gaussianFormatBytes(totals.sphericalHarmonicsBytes),
+            gaussianFormatBytes(totals.uniformBytes),
+            gaussianFormatBytes(totals.scratchBytes),
+            suffix
+        ),
+        category: LogCategory.gaussian.rawValue
+    )
+}
+
+func gaussianFormatBytes(_ bytes: Int) -> String {
+    let value = Double(bytes)
+    if bytes >= 1024 * 1024 {
+        return String(format: "%.2fMiB", value / 1_048_576.0)
+    }
+    if bytes >= 1024 {
+        return String(format: "%.2fKiB", value / 1024.0)
+    }
+    return "\(bytes)B"
+}
diff --git a/Sources/UntoldEngine/Renderer/RenderPasses.swift b/Sources/UntoldEngine/Renderer/RenderPasses.swift
index c6bd8d2a..2d3de2cd 100644
--- a/Sources/UntoldEngine/Renderer/RenderPasses.swift
+++ b/Sources/UntoldEngine/Renderer/RenderPasses.swift
@@ -3993,6 +3993,10 @@ public enum RenderPasses {
     }
 
     public static let gaussianExecution: RenderPassExecution = { commandBuffer in
+        let profileStart = gaussianProfilingStartTime()
+        var profileTotals = GaussianProfileTotals()
+        var activeSplatTotal = 0
+
         let pipelines = PipelineManager.shared.renderPipelinesByType
         guard let initializePipeline = pipelines[.gaussianTBDRInitialize],
               let drawPipeline = pipelines[.gaussianTBDRDraw],
@@ -4046,6 +4050,7 @@ public enum RenderPasses {
         renderEncoder.pushDebugGroup("Initialize")
         renderEncoder.setRenderPipelineState(initializePipelineState)
         renderEncoder.dispatchThreadsPerTile(MTLSize(width: 32, height: 32, depth: 1))
+        profileTotals.dispatchCount += 1
         renderEncoder.popDebugGroup()
 
         renderEncoder.pushDebugGroup("Draw Splats")
@@ -4063,6 +4068,11 @@ public enum RenderPasses {
                 handleError(.noGaussianComponent, entityId)
                 continue
             }
+            profileTotals.include(component: gaussianComponent)
+
+            let activeSplatCount = min(Int(gaussianComponent.visibleSplatCountForRendering), Int(gaussianComponent.splatCount))
+            guard activeSplatCount > 0 else { continue }
+            activeSplatTotal += activeSplatCount
 
             guard gaussianComponent.encodedSplatData != nil else {
                 handleError(.bufferAllocationFailed, "Encoded Gaussian splat buffer")
@@ -4135,10 +4145,37 @@ public enum RenderPasses {
             )
             renderEncoder.setVertexBytes(&renderInfo.viewPort, length: MemoryLayout<simd_float2>.stride, index: Int(gaussianTBDRRenderViewPortIndex.rawValue))
 
+            var shMetadata = gaussianComponent.sphericalHarmonicsMetadata ?? GaussianSHMetadata(
+                degree: 0,
+                coefficientsPerChannel: 0,
+                higherOrderCoefficientsPerSplat: 0,
+                _pad0: 0
+            )
+            renderEncoder.setVertexBuffer(
+                gaussianComponent.sphericalHarmonicsData ?? gaussianComponent.encodedSplatData,
+                offset: 0,
+                index: Int(gaussianTBDRRenderSHIndex.rawValue)
+            )
+            renderEncoder.setVertexBytes(
+                &shMetadata,
+                length: MemoryLayout<GaussianSHMetadata>.stride,
+                index: Int(gaussianTBDRRenderSHMetadataIndex.rawValue)
+            )
+            var localCameraPosition = gaussianLocalCameraPosition(
+                cameraWorldPosition: effectiveCameraPosition,
+                modelMatrix: modelMatrix
+            )
+            renderEncoder.setVertexBytes(
+                &localCameraPosition,
+                length: MemoryLayout<simd_float3>.stride,
+                index: Int(gaussianTBDRRenderLocalCameraIndex.rawValue)
+            )
+
             renderEncoder.drawPrimitivesTracked(type: .triangleStrip,
                                                 vertexStart: 0,
                                                 vertexCount: 4,
-                                                instanceCount: Int(gaussianComponent.splatCount))
+                                                instanceCount: activeSplatCount)
+            profileTotals.drawCallCount += 1
         }
 
         renderEncoder.popDebugGroup()
@@ -4149,11 +4186,18 @@ public enum RenderPasses {
         var reverseZ = renderInfo.reverseZEnabled
         renderEncoder.setFragmentBytes(&reverseZ, length: MemoryLayout<Bool>.stride, index: Int(gaussianTBDRRenderReverseZIndex.rawValue))
         renderEncoder.drawPrimitivesTracked(type: .triangle, vertexStart: 0, vertexCount: 3)
+        profileTotals.drawCallCount += 1
         renderEncoder.popDebugGroup()
 
         renderEncoder.updateFence(renderInfo.fence, after: .fragment)
         renderEncoder.popDebugGroup()
         renderEncoder.endEncoding()
+        logGaussianProfile(
+            stage: "Render",
+            startTime: profileStart,
+            totals: profileTotals,
+            extra: String(format: "activeSplats=%d viewport=%.0fx%.0f tile=32x32", activeSplatTotal, renderInfo.viewPort.x, renderInfo.viewPort.y)
+        )
     }
 
     static func executePostProcess(
diff --git a/Sources/UntoldEngine/Shaders/BitonicSort.metal b/Sources/UntoldEngine/Shaders/BitonicSort.metal
index 03b29cdb..f92b5e0b 100644
--- a/Sources/UntoldEngine/Shaders/BitonicSort.metal
+++ b/Sources/UntoldEngine/Shaders/BitonicSort.metal
@@ -38,17 +38,62 @@ inline float eye_space_depth(const float4x4 modelView, float3 worldPos) {
     return max(-v.z, 0.0f);
 }
 
+kernel void gaussianResetVisibleCount(
+    device atomic_uint *visibleCount [[buffer(gaussianVisibleCountIndex)]],
+    uint index [[thread_position_in_grid]])
+{
+    if (index == 0u) {
+        atomic_store_explicit(visibleCount, 0u, memory_order_relaxed);
+    }
+}
+
+kernel void gaussianFrustumCull(
+    const device EncodedGaussianSplat *splats [[buffer(gaussianEncodedSplatIndex)]],
+    constant Uniforms &uniforms [[buffer(gaussianUniformIndex)]],
+    constant uint &numOfSplats [[buffer(gaussianNumberOfSplatsIndex)]],
+    constant float &clipGuardBand [[buffer(gaussianIndicesIndex)]],
+    device uint *visibleIndices [[buffer(gaussianVisibleIndicesIndex)]],
+    device atomic_uint *visibleCount [[buffer(gaussianVisibleCountIndex)]],
+    uint index [[thread_position_in_grid]])
+{
+    if (index >= numOfSplats) return;
+
+    float4 centerClip = uniforms.projectionMatrix *
+                        uniforms.modelViewMatrix *
+                        float4(splats[index].position, 1.0f);
+
+    if (centerClip.w <= 0.0f) return;
+
+    float limit = max(0.0f, 1.0f + clipGuardBand);
+    float2 ndc = centerClip.xy / centerClip.w;
+    if (abs(ndc.x) > limit || abs(ndc.y) > limit) return;
+    if (centerClip.z < -centerClip.w * clipGuardBand || centerClip.z > centerClip.w * limit) return;
+
+    uint writeIndex = atomic_fetch_add_explicit(visibleCount, 1u, memory_order_relaxed);
+    visibleIndices[writeIndex] = index;
+}
+
 kernel void gaussianDepthKeys(device uint64_t *outKeys              [[buffer(gaussianIndicesIndex)]],
-                              const device GaussianSplat *splats     [[buffer(gaussianSplatIndex)]],
+                              const device EncodedGaussianSplat *splats [[buffer(gaussianEncodedSplatIndex)]],
                               constant uint &numOfSplats [[buffer(gaussianNumberOfSplatsIndex)]],
                               constant Uniforms &uniforms     [[buffer(gaussianUniformIndex)]],
+                              const device uint *visibleIndices [[buffer(gaussianVisibleIndicesIndex)]],
+                              const device uint *visibleCount [[buffer(gaussianVisibleCountIndex)]],
                               uint index                              [[thread_position_in_grid]],
                               uint gridSize                           [[threads_per_grid]])
 {
     if (index >= numOfSplats) return;
 
+    if (index >= visibleCount[0]) {
+        uint64_t invalidPacked = ((uint64_t)0xffffffffu << 32) | (uint64_t)0xffffffffu;
+        outKeys[index] = invalidPacked;
+        return;
+    }
+
+    uint splatIndex = visibleIndices[index];
+
     // 1) Eye-space depth
-    float z = eye_space_depth(uniforms.modelViewMatrix, splats[index].center.xyz);
+    float z = eye_space_depth(uniforms.modelViewMatrix, splats[splatIndex].position);
 
     // 2) Convert to lexicographically sortable unsigned int
     uint keyFrontToBack = float_to_sortable_u32(z);
@@ -57,7 +102,6 @@ kernel void gaussianDepthKeys(device uint64_t *outKeys              [[buffer(gau
     //uint keyBackToFront = 0xffffffffu - keyFrontToBack;
 
     // 3) Pack [key | index] for 64-bit radix sort (key in high bits)
-    uint64_t packed = (uint64_t)keyFrontToBack << 32 | (uint64_t)index;
+    uint64_t packed = (uint64_t)keyFrontToBack << 32 | (uint64_t)splatIndex;
     outKeys[index] = packed;
 }
-
diff --git a/Sources/UntoldEngine/Shaders/Gaussians.metal b/Sources/UntoldEngine/Shaders/Gaussians.metal
index eda94e5c..60f8c7ce 100644
--- a/Sources/UntoldEngine/Shaders/Gaussians.metal
+++ b/Sources/UntoldEngine/Shaders/Gaussians.metal
@@ -15,9 +15,124 @@
 
 using namespace metal;
 
+constant float GAUSSIAN_SH_C1 = 0.4886025119029199f;
+constant float GAUSSIAN_SH_C2[5] = {
+    1.0925484305920792f,
+   -1.0925484305920792f,
+    0.31539156525252005f,
+   -1.0925484305920792f,
+    0.5462742152960396f
+};
+constant float GAUSSIAN_SH_C3[7] = {
+   -0.5900435899266435f,
+    2.890611442640554f,
+   -0.4570457994644658f,
+    0.3731763325901154f,
+   -0.4570457994644658f,
+    1.445305721320277f,
+   -0.5900435899266435f
+};
+
 inline uint unpackIndex(uint64_t packed)  { return (uint)(packed & 0xffffffffu); }
 inline uint unpackDepthKey(uint64_t packed) { return (uint)(packed >> 32); }
 
+float3 loadGaussianSHCoefficient(
+    const device half *coefficients,
+    constant GaussianSHMetadata &metadata,
+    uint splatIndex,
+    uint coefficientIndex)
+{
+    uint perChannel = metadata.coefficientsPerChannel - 1;
+    uint splatBase = splatIndex * metadata.higherOrderCoefficientsPerSplat;
+    uint offset = coefficientIndex - 1;
+    return float3(
+        coefficients[splatBase + offset],
+        coefficients[splatBase + perChannel + offset],
+        coefficients[splatBase + 2 * perChannel + offset]
+    );
+}
+
+float3 evaluateGaussianSphericalHarmonics(
+    float3 baseColor,
+    const device half *coefficients,
+    constant GaussianSHMetadata &metadata,
+    uint splatIndex,
+    float3 direction)
+{
+    if (metadata.degree == 0 || metadata.higherOrderCoefficientsPerSplat == 0) {
+        return baseColor;
+    }
+
+    float lengthSquared = dot(direction, direction);
+    if (lengthSquared <= 1.0e-8f) {
+        return baseColor;
+    }
+    float3 d = direction * rsqrt(lengthSquared);
+    float x = d.x;
+    float y = d.y;
+    float z = d.z;
+    float3 result = baseColor;
+    result += GAUSSIAN_SH_C1 * (
+        -y * loadGaussianSHCoefficient(coefficients, metadata, splatIndex, 1) +
+         z * loadGaussianSHCoefficient(coefficients, metadata, splatIndex, 2) -
+         x * loadGaussianSHCoefficient(coefficients, metadata, splatIndex, 3)
+    );
+
+    if (metadata.degree > 1) {
+        float xx = x * x;
+        float yy = y * y;
+        float zz = z * z;
+        result += GAUSSIAN_SH_C2[0] * x * y * loadGaussianSHCoefficient(coefficients, metadata, splatIndex, 4);
+        result += GAUSSIAN_SH_C2[1] * y * z * loadGaussianSHCoefficient(coefficients, metadata, splatIndex, 5);
+        result += GAUSSIAN_SH_C2[2] * (2.0f * zz - xx - yy) * loadGaussianSHCoefficient(coefficients, metadata, splatIndex, 6);
+        result += GAUSSIAN_SH_C2[3] * x * z * loadGaussianSHCoefficient(coefficients, metadata, splatIndex, 7);
+        result += GAUSSIAN_SH_C2[4] * (xx - yy) * loadGaussianSHCoefficient(coefficients, metadata, splatIndex, 8);
+
+        if (metadata.degree > 2) {
+            result += GAUSSIAN_SH_C3[0] * y * (3.0f * xx - yy) * loadGaussianSHCoefficient(coefficients, metadata, splatIndex, 9);
+            result += GAUSSIAN_SH_C3[1] * x * y * z * loadGaussianSHCoefficient(coefficients, metadata, splatIndex, 10);
+            result += GAUSSIAN_SH_C3[2] * y * (4.0f * zz - xx - yy) * loadGaussianSHCoefficient(coefficients, metadata, splatIndex, 11);
+            result += GAUSSIAN_SH_C3[3] * z * (2.0f * zz - 3.0f * xx - 3.0f * yy) * loadGaussianSHCoefficient(coefficients, metadata, splatIndex, 12);
+            result += GAUSSIAN_SH_C3[4] * x * (4.0f * zz - xx - yy) * loadGaussianSHCoefficient(coefficients, metadata, splatIndex, 13);
+            result += GAUSSIAN_SH_C3[5] * z * (xx - yy) * loadGaussianSHCoefficient(coefficients, metadata, splatIndex, 14);
+            result += GAUSSIAN_SH_C3[6] * x * (xx - 3.0f * yy) * loadGaussianSHCoefficient(coefficients, metadata, splatIndex, 15);
+        }
+    }
+
+    return max(result, float3(0.0f));
+}
+
+// Standard 3DGS coefficients are fitted to normalized image-code values.
+// Decode those display-referred values before writing to Untold's linear HDR
+// Gaussian target. The final output transform will encode them exactly once.
+float3 gaussianSRGBToLinear(float3 color)
+{
+    color = max(color, float3(0.0f));
+    float3 low = color / 12.92f;
+    float3 high = pow((color + 0.055f) / 1.055f, float3(2.4f));
+    return select(high, low, color <= 0.04045f);
+}
+
+// Diagnostic entry point for validating the packed GPU SH contract against
+// the exact evaluator used by the Gaussian vertex shader.
+kernel void gaussianSphericalHarmonicsDiagnostic(
+    const device half *coefficients       [[buffer(0)]],
+    constant GaussianSHMetadata &metadata [[buffer(1)]],
+    constant float4 &baseColor            [[buffer(2)]],
+    constant float4 &direction            [[buffer(3)]],
+    device float4 *result                 [[buffer(4)]])
+{
+    float3 evaluated = evaluateGaussianSphericalHarmonics(
+        baseColor.xyz,
+        coefficients,
+        metadata,
+        0,
+        direction.xyz
+    );
+    result[0] = float4(evaluated, 1.0f);
+    result[1] = float4(gaussianSRGBToLinear(evaluated), 1.0f);
+}
+
 // Project 3D covariance into 2D screen-pixel space
 float3 computeCov2D(float4      splatCenter,
                     float3x3    cov3D,
@@ -165,6 +280,9 @@ vertex GaussianOutData vertexGaussianTBDRShader(
     const device EncodedGaussianSplat *splats     [[buffer(gaussianTBDRRenderSplatIndex)]],
     constant Uniforms                 &uniforms   [[buffer(gaussianTBDRRenderUniformIndex)]],
     constant float2                   &viewport   [[buffer(gaussianTBDRRenderViewPortIndex)]],
+    const device half                 *shCoefficients [[buffer(gaussianTBDRRenderSHIndex)]],
+    constant GaussianSHMetadata       &shMetadata [[buffer(gaussianTBDRRenderSHMetadataIndex)]],
+    constant float3                   &localCameraPosition [[buffer(gaussianTBDRRenderLocalCameraIndex)]],
     uint                               vid        [[vertex_id]],
     uint                               iid        [[instance_id]])
 {
@@ -174,6 +292,9 @@ vertex GaussianOutData vertexGaussianTBDRShader(
 
     uint64_t packed = packedKeys[iid];
     uint splatIndex = unpackIndex(packed);
+    if (splatIndex == 0xffffffffu) {
+        return out;
+    }
     const EncodedGaussianSplat splat = splats[splatIndex];
 
     float2 quad = getCurrentQuadVertex(vid);
@@ -216,7 +337,13 @@ vertex GaussianOutData vertexGaussianTBDRShader(
     float2 ndcOffset = quad * extent * 2.0f / viewport;
     out.position = centerClip;
     out.position.xy += ndcOffset * centerClip.w;
-    out.color = splat.color;
+    out.color = gaussianSRGBToLinear(evaluateGaussianSphericalHarmonics(
+        splat.color,
+        shCoefficients,
+        shMetadata,
+        splatIndex,
+        centerLocal - localCameraPosition
+    ));
     out.alpha = splat.opacity;
     out.valid = true;
 
diff --git a/Sources/UntoldEngine/Systems/GaussianSystem.swift b/Sources/UntoldEngine/Systems/GaussianSystem.swift
index 63a9cc2e..db534bb1 100644
--- a/Sources/UntoldEngine/Systems/GaussianSystem.swift
+++ b/Sources/UntoldEngine/Systems/GaussianSystem.swift
@@ -22,6 +22,16 @@ import simd
 
 let maxNumOfGaussians: UInt64 = 1024 * 1024 * 5
 
+private func activeGaussianSortCount(_ component: GaussianComponent) -> Int {
+    min(Int(component.visibleSplatCountForRendering), Int(component.splatCount))
+}
+
+private struct GaussianVisibleCountUpdate: @unchecked Sendable {
+    let component: GaussianComponent
+    let visibleCount: MTLBuffer
+    let splatCount: UInt
+}
+
 func initGuassianComputePipelines() {
     if renderInfo.device == nil {
         handleError(.metalDeviceNotFound)
@@ -33,6 +43,10 @@ func initGuassianComputePipelines() {
         return
     }
 
+    createComputePipeline(into: &gaussianResetVisibleCountPipeline, device: renderInfo.device, library: renderInfo.library, functionName: "gaussianResetVisibleCount", pipelineName: "Gaussian Reset Visible Count")
+
+    createComputePipeline(into: &gaussianFrustumCullPipeline, device: renderInfo.device, library: renderInfo.library, functionName: "gaussianFrustumCull", pipelineName: "Gaussian Frustum Cull")
+
     createComputePipeline(into: &gaussianDepthPipeline, device: renderInfo.device, library: renderInfo.library, functionName: "gaussianDepthKeys", pipelineName: "Gaussian Depth")
 
     createComputePipeline(into: &radixClearHistogramPipeline, device: renderInfo.device, library: renderInfo.library, functionName: "gaussianRadixClearHistogram", pipelineName: "Radix Clear")
@@ -46,7 +60,130 @@ func initGuassianComputePipelines() {
     createComputePipeline(into: &radixScatterPipeline, device: renderInfo.device, library: renderInfo.library, functionName: "gaussianRadixScatter", pipelineName: "Radix Scatter")
 }
 
+public func executeGaussianFrustumCulling(_ commandBuffer: MTLCommandBuffer) {
+    let profileStart = gaussianProfilingStartTime()
+    var profileTotals = GaussianProfileTotals()
+    var activeSplatTotal = 0
+
+    guard gaussianResetVisibleCountPipeline.success else {
+        handleError(.pipelineStateNulled, gaussianResetVisibleCountPipeline.name!); return
+    }
+    guard gaussianFrustumCullPipeline.success else {
+        handleError(.pipelineStateNulled, gaussianFrustumCullPipeline.name!); return
+    }
+    guard let camera = CameraSystem.shared.activeCamera,
+          let cameraComponent = scene.get(component: CameraComponent.self, for: camera)
+    else {
+        handleError(.noActiveCamera)
+        return
+    }
+    guard let resetPipelineState = gaussianResetVisibleCountPipeline.pipelineState,
+          let cullPipelineState = gaussianFrustumCullPipeline.pipelineState
+    else {
+        handleError(.pipelineStateNulled, "Gaussian culling pipeline state is nil")
+        return
+    }
+
+    let transformId = getComponentId(for: WorldTransformComponent.self)
+    let gaussianId = getComponentId(for: GaussianComponent.self)
+    let entities = queryEntitiesWithComponentIds([transformId, gaussianId], in: scene)
+    guard !entities.isEmpty else { return }
+
+    guard let computeEncoder = commandBuffer.makeComputeCommandEncoder() else { return }
+    computeEncoder.label = "Gaussian Frustum Culling"
+
+    var visibleCountUpdates: [GaussianVisibleCountUpdate] = []
+
+    for entityId in entities {
+        guard let gaussianComponent = scene.get(component: GaussianComponent.self, for: entityId) else {
+            handleError(.noGaussianComponent, entityId)
+            continue
+        }
+        guard let worldTransformComponent = scene.get(component: WorldTransformComponent.self, for: entityId) else {
+            handleError(.noWorldTransformComponent, entityId)
+            continue
+        }
+        guard let encodedSplatData = gaussianComponent.encodedSplatData,
+              let visibleIndices = gaussianComponent.gaussianVisibleIndices,
+              let visibleCount = gaussianComponent.gaussianVisibleCount
+        else {
+            handleError(.bufferAllocationFailed, "Gaussian culling buffers")
+            continue
+        }
+
+        profileTotals.include(component: gaussianComponent)
+        activeSplatTotal += activeGaussianSortCount(gaussianComponent)
+
+        computeEncoder.setComputePipelineState(resetPipelineState)
+        computeEncoder.setBuffer(visibleCount, offset: 0, index: Int(gaussianVisibleCountIndex.rawValue))
+        computeEncoder.dispatchThreadgroups(MTLSizeMake(1, 1, 1), threadsPerThreadgroup: MTLSizeMake(1, 1, 1))
+        profileTotals.dispatchCount += 1
+
+        let modelMatrix = simd_mul(worldTransformComponent.space, .identity)
+        let viewMatrix = cameraComponent.viewSpace
+        let modelViewMatrix = simd_mul(viewMatrix, modelMatrix)
+
+        var gaussianUniform = Uniforms()
+        gaussianUniform.modelViewMatrix = modelViewMatrix
+        gaussianUniform.viewMatrix = viewMatrix
+        gaussianUniform.modelMatrix = modelMatrix
+        gaussianUniform.cameraPosition = cameraComponent.localPosition
+        gaussianUniform.projectionMatrix = renderInfo.perspectiveSpace
+
+        var totalSplats = UInt32(gaussianComponent.splatCount)
+        var clipGuardBand: Float = 0.25
+
+        computeEncoder.setComputePipelineState(cullPipelineState)
+        computeEncoder.setBuffer(encodedSplatData, offset: 0, index: Int(gaussianEncodedSplatIndex.rawValue))
+        computeEncoder.setBytes(&gaussianUniform, length: MemoryLayout<Uniforms>.stride, index: Int(gaussianUniformIndex.rawValue))
+        computeEncoder.setBytes(&totalSplats, length: MemoryLayout<UInt32>.stride, index: Int(gaussianNumberOfSplatsIndex.rawValue))
+        computeEncoder.setBuffer(visibleIndices, offset: 0, index: Int(gaussianVisibleIndicesIndex.rawValue))
+        computeEncoder.setBuffer(visibleCount, offset: 0, index: Int(gaussianVisibleCountIndex.rawValue))
+        computeEncoder.setBytes(&clipGuardBand, length: MemoryLayout<Float>.stride, index: Int(gaussianIndicesIndex.rawValue))
+
+        let tew = cullPipelineState.threadExecutionWidth
+        let maxT = cullPipelineState.maxTotalThreadsPerThreadgroup
+        var block = min(256, maxT)
+        block = max((block / tew) * tew, tew)
+        let numThreadgroups = (Int(gaussianComponent.splatCount) + block - 1) / block
+        computeEncoder.dispatchThreadgroups(
+            MTLSizeMake(numThreadgroups, 1, 1),
+            threadsPerThreadgroup: MTLSizeMake(block, 1, 1)
+        )
+        profileTotals.dispatchCount += 1
+
+        visibleCountUpdates.append(
+            GaussianVisibleCountUpdate(
+                component: gaussianComponent,
+                visibleCount: visibleCount,
+                splatCount: gaussianComponent.splatCount
+            )
+        )
+    }
+
+    computeEncoder.endEncoding()
+
+    let updates = visibleCountUpdates
+    commandBuffer.addCompletedHandler { _ in
+        for update in updates {
+            let count = update.visibleCount.contents().load(as: UInt32.self)
+            update.component.visibleSplatCountForRendering = min(UInt(count), update.splatCount)
+        }
+    }
+
+    logGaussianProfile(
+        stage: "FrustumCull",
+        startTime: profileStart,
+        totals: profileTotals,
+        extra: "previousActiveSplats=\(activeSplatTotal)"
+    )
+}
+
 public func executeGaussianDepth(_ commandBuffer: MTLCommandBuffer) {
+    let profileStart = gaussianProfilingStartTime()
+    var profileTotals = GaussianProfileTotals()
+    var activeSplatTotal = 0
+
     if gaussianDepthPipeline.success == false {
         handleError(.pipelineStateNulled, gaussianDepthPipeline.name!)
         return
@@ -72,6 +209,11 @@ public func executeGaussianDepth(_ commandBuffer: MTLCommandBuffer) {
             handleError(.noGaussianComponent, entityId)
             continue
         }
+        profileTotals.include(component: gaussianComponent)
+
+        let activeCount = activeGaussianSortCount(gaussianComponent)
+        guard activeCount > 0 else { continue }
+        activeSplatTotal += activeCount
 
         guard let worldTransformComponent = scene.get(component: WorldTransformComponent.self, for: entityId) else {
             handleError(.noWorldTransformComponent, entityId)
@@ -84,7 +226,13 @@ public func executeGaussianDepth(_ commandBuffer: MTLCommandBuffer) {
         }
 
         computeEncoder.setBuffer(gaussianComponent.gaussianSortedIndices, offset: 0, index: Int(gaussianIndicesIndex.rawValue))
-        computeEncoder.setBuffer(gaussianComponent.splatData, offset: 0, index: Int(gaussianSplatIndex.rawValue))
+        computeEncoder.setBuffer(gaussianComponent.gaussianVisibleIndices, offset: 0, index: Int(gaussianVisibleIndicesIndex.rawValue))
+        computeEncoder.setBuffer(gaussianComponent.gaussianVisibleCount, offset: 0, index: Int(gaussianVisibleCountIndex.rawValue))
+        computeEncoder.setBuffer(
+            gaussianComponent.encodedSplatData,
+            offset: 0,
+            index: Int(gaussianEncodedSplatIndex.rawValue)
+        )
 
         // update uniforms
         var gaussianUniform = Uniforms()
@@ -132,7 +280,7 @@ public func executeGaussianDepth(_ commandBuffer: MTLCommandBuffer) {
             gaussianComponent.spaceUniform[uniformBufferIndex], offset: 0, index: Int(gaussianUniformIndex.rawValue)
         )
 
-        var localNumGaussians = UInt32(gaussianComponent.splatCount)
+        var localNumGaussians = UInt32(activeCount)
         computeEncoder.setBytes(&localNumGaussians, length: MemoryLayout<UInt32>.stride, index: Int(gaussianNumberOfSplatsIndex.rawValue))
 
         let tew = gaussianDepthPipeline.pipelineState?.threadExecutionWidth ?? 32
@@ -145,13 +293,20 @@ public func executeGaussianDepth(_ commandBuffer: MTLCommandBuffer) {
         let threadsPerThreadgroup: MTLSize = MTLSizeMake(block, 1, 1)
 
         // Use dispatchThreadgroups for broader device compatibility (including Vision Pro)
-        let numThreadgroups = (Int(gaussianComponent.splatCount) + block - 1) / block
+        let numThreadgroups = (activeCount + block - 1) / block
         let threadgroupsPerGrid: MTLSize = MTLSizeMake(numThreadgroups, 1, 1)
 
         computeEncoder.dispatchThreadgroups(threadgroupsPerGrid, threadsPerThreadgroup: threadsPerThreadgroup)
+        profileTotals.dispatchCount += 1
     }
 
     computeEncoder.endEncoding()
+    logGaussianProfile(
+        stage: "DepthKeys",
+        startTime: profileStart,
+        totals: profileTotals,
+        extra: "activeSplats=\(activeSplatTotal)"
+    )
 }
 
 // MARK: - Device Radix Sort
@@ -175,6 +330,10 @@ public func executeGaussianDepth(_ commandBuffer: MTLCommandBuffer) {
 // indexing (perTGStart[groupId * 256 + digit]).
 
 public func executeRadixSort(_ commandBuffer: MTLCommandBuffer) {
+    let profileStart = gaussianProfilingStartTime()
+    var profileTotals = GaussianProfileTotals()
+    var activeSplatTotal = 0
+
     guard radixClearHistogramPipeline.success else {
         handleError(.pipelineStateNulled, radixClearHistogramPipeline.name!); return
     }
@@ -212,8 +371,10 @@ public func executeRadixSort(_ commandBuffer: MTLCommandBuffer) {
         guard let gc = scene.get(component: GaussianComponent.self, for: entityId) else { continue }
         guard let sortedIndices = gc.gaussianSortedIndices else { continue }
 
-        let n = Int(gc.splatCount)
+        let n = activeGaussianSortCount(gc)
         guard n >= 2 else { continue }
+        profileTotals.include(component: gc)
+        activeSplatTotal += n
 
         // Ping-pong temp buffer (CPU alloc before encoding)
         let keyBufLen = n * MemoryLayout<UInt64>.stride
@@ -223,6 +384,7 @@ public func executeRadixSort(_ commandBuffer: MTLCommandBuffer) {
             )
         }
         guard let tempBuffer = radixSortTempBuffer else { continue }
+        profileTotals.scratchBytes = max(profileTotals.scratchBytes, tempBuffer.length)
 
         // Fixed block size: histogram and scatter MUST use the same value so
         // that histGroups == scatterGroups and perTGStart indexing is correct.
@@ -236,6 +398,10 @@ public func executeRadixSort(_ commandBuffer: MTLCommandBuffer) {
             )
         }
         guard let perTGBuf = radixPerTGHistBuffer else { continue }
+        profileTotals.scratchBytes = max(
+            profileTotals.scratchBytes,
+            tempBuffer.length + histBuffer.length + perTGBuf.length
+        )
 
         var numElems = UInt32(n)
         var numBuckets = UInt32(256)
@@ -251,6 +417,7 @@ public func executeRadixSort(_ commandBuffer: MTLCommandBuffer) {
             enc.setComputePipelineState(radixClearHistogramPipeline.pipelineState!)
             enc.setBuffer(histBuffer, offset: 0, index: Int(radixClearHistogramBuffer.rawValue))
             enc.dispatchThreadgroups(MTLSizeMake(1, 1, 1), threadsPerThreadgroup: MTLSizeMake(256, 1, 1))
+            profileTotals.dispatchCount += 1
 
             // ── 1. Histogram + per-TG counts ─────────────────────────────────
             enc.setComputePipelineState(radixHistogramPipeline.pipelineState!)
@@ -260,18 +427,21 @@ public func executeRadixSort(_ commandBuffer: MTLCommandBuffer) {
             enc.setBytes(&numElems, length: MemoryLayout<UInt32>.stride, index: Int(radixHistogramNumElems.rawValue))
             enc.setBytes(&passIdx, length: MemoryLayout<UInt32>.stride, index: Int(radixHistogramPassIndex.rawValue))
             enc.dispatchThreadgroups(MTLSizeMake(numGroups, 1, 1), threadsPerThreadgroup: MTLSizeMake(radixBlock, 1, 1))
+            profileTotals.dispatchCount += 1
 
             // ── 2. Per-TG column scan → per-TG starting offsets ─────────────
             enc.setComputePipelineState(radixScanPerTGPipeline.pipelineState!)
             enc.setBuffer(perTGBuf, offset: 0, index: Int(radixScanPerTGBuffer.rawValue))
             enc.setBytes(&numGroups32, length: MemoryLayout<UInt32>.stride, index: Int(radixScanPerTGNumGroups.rawValue))
             enc.dispatchThreadgroups(MTLSizeMake(1, 1, 1), threadsPerThreadgroup: MTLSizeMake(256, 1, 1))
+            profileTotals.dispatchCount += 1
 
             // ── 3. Global exclusive scan ─────────────────────────────────────
             enc.setComputePipelineState(radixScanPipeline.pipelineState!)
             enc.setBuffer(histBuffer, offset: 0, index: Int(radixScanHistogram.rawValue))
             enc.setBytes(&numBuckets, length: MemoryLayout<UInt32>.stride, index: Int(radixScanNumBuckets.rawValue))
             enc.dispatchThreadgroups(MTLSizeMake(1, 1, 1), threadsPerThreadgroup: MTLSizeMake(256, 1, 1))
+            profileTotals.dispatchCount += 1
 
             // ── 4. Stable scatter ────────────────────────────────────────────
             enc.setComputePipelineState(radixScatterPipeline.pipelineState!)
@@ -282,61 +452,16 @@ public func executeRadixSort(_ commandBuffer: MTLCommandBuffer) {
             enc.setBytes(&numElems, length: MemoryLayout<UInt32>.stride, index: Int(radixScatterNumElems.rawValue))
             enc.setBytes(&passIdx, length: MemoryLayout<UInt32>.stride, index: Int(radixScatterPassIdx.rawValue))
             enc.dispatchThreadgroups(MTLSizeMake(numGroups, 1, 1), threadsPerThreadgroup: MTLSizeMake(radixBlock, 1, 1))
+            profileTotals.dispatchCount += 1
+            profileTotals.radixPassCount += 1
         }
     }
 
     enc.endEncoding()
+    logGaussianProfile(
+        stage: "RadixSort",
+        startTime: profileStart,
+        totals: profileTotals,
+        extra: "activeSplats=\(activeSplatTotal)"
+    )
 }
-
-// MARK: - PLY Loading Helpers
-
-/*
- /// Loads Gaussian splats from a PLY file into the GPU buffer
- public func loadGaussianSplatsFromPLY(url: URL) throws -> Int {
-     // Load splats from PLY file
-     let splats = try PLYReader.readGaussianSplats(from: url)
-
-     // Check if we exceed the buffer capacity
-     guard splats.count <= Int(maxNumOfGaussians) else {
-         handleError(.bufferAllocationFailed, "Too many Gaussian splats: \(splats.count) exceeds maximum \(maxNumOfGaussians)")
-         throw PLYError.invalidData("Too many Gaussian splats: \(splats.count) exceeds maximum \(maxNumOfGaussians)")
-     }
-
-     // Copy to GPU buffer
-     guard let splatBuffer = bufferResources.splatData else {
-         handleError(.bufferAllocationFailed, "Gaussian splat buffer is nil")
-         throw PLYError.invalidData("Gaussian splat buffer not initialized")
-     }
-
-     let pointer = splatBuffer.contents().bindMemory(
-         to: GaussianSplat.self,
-         capacity: splats.count
-     )
-
-     for (index, splat) in splats.enumerated() {
-         pointer[index] = splat
-     }
-
-     // Update current count
-     currentNumOfGaussians = splats.count
-
-     print("✓ Loaded \(splats.count) Gaussian splats from \(url.lastPathComponent)")
-
-     return splats.count
- }
-
- /// Loads Gaussian splats from a PLY file path
- public func loadGaussianSplatsFromPLY(path: String) throws -> Int {
-     let url = URL(fileURLWithPath: path)
-     return try loadGaussianSplatsFromPLY(url: url)
- }
-
- /// Loads Gaussian splats from a PLY file in the app bundle
- public func loadGaussianSplatsFromBundle(filename: String, bundle: Bundle = .main) throws -> Int {
-     guard let url = bundle.url(forResource: filename, withExtension: "ply") else {
-         handleError(.assetDataMissing, "PLY file '\(filename).ply' not found in bundle")
-         throw PLYError.invalidFormat("PLY file '\(filename).ply' not found in bundle")
-     }
-     return try loadGaussianSplatsFromPLY(url: url)
- }
- */
diff --git a/Sources/UntoldEngine/Systems/RegistrationSystem.swift b/Sources/UntoldEngine/Systems/RegistrationSystem.swift
index 63e9df02..9c283560 100644
--- a/Sources/UntoldEngine/Systems/RegistrationSystem.swift
+++ b/Sources/UntoldEngine/Systems/RegistrationSystem.swift
@@ -2911,13 +2911,14 @@ public func setEntityGaussian(entityId: EntityID, filename: String, withExtensio
     }
 
     // Attempt to read Gaussian splats, handling errors internally
-    let splats: [GaussianSplat]
+    let asset: GaussianSplatAsset
     do {
-        splats = try PLYReader.readGaussianSplats(from: url)
+        asset = try PLYReader.readGaussianAsset(from: url)
     } catch {
         handleError(.assetDataMissing, "Failed to read Gaussian splats from \(filename): \(error.localizedDescription)")
         return
     }
+    let splats = asset.splats
 
     // Check if we exceed the buffer capacity
     guard splats.count <= Int(maxNumOfGaussians) else {
@@ -2926,36 +2927,76 @@ public func setEntityGaussian(entityId: EntityID, filename: String, withExtensio
     }
 
     let splatCount = UInt(splats.count)
-    var temSplatCount = splatCount
-    let tempPowerOfTwoSplatCount: UInt = nextPowerOf2(x: &temSplatCount)
+    guard splatCount > 0 else {
+        handleError(.assetDataMissing, "Gaussian splat file contains no vertices: \(filename)")
+        return
+    }
 
-    let gaussianSortedIndices = renderInfo.device.makeBuffer(length: MemoryLayout<UInt64>.stride * Int(tempPowerOfTwoSplatCount), options: .storageModeShared)
+    guard let gaussianSortedIndices = renderInfo.device.makeBuffer(
+        length: MemoryLayout<UInt64>.stride * Int(splatCount),
+        options: .storageModeShared
+    ) else {
+        handleError(.bufferAllocationFailed, "Gaussian sorted-index buffer is nil")
+        return
+    }
 
-    guard let splatBuffer = renderInfo.device.makeBuffer(length: MemoryLayout<GaussianSplat>.stride * Int(splatCount), options: .storageModeShared) else {
-        handleError(.bufferAllocationFailed, "Gaussian splat buffer is nil")
+    guard let gaussianVisibleIndices = renderInfo.device.makeBuffer(
+        length: MemoryLayout<UInt32>.stride * Int(splatCount),
+        options: .storageModeShared
+    ) else {
+        handleError(.bufferAllocationFailed, "Gaussian visible-index buffer is nil")
         return
     }
 
+    guard let gaussianVisibleCount = renderInfo.device.makeBuffer(
+        length: MemoryLayout<UInt32>.stride,
+        options: .storageModeShared
+    ) else {
+        handleError(.bufferAllocationFailed, "Gaussian visible-count buffer is nil")
+        return
+    }
+    gaussianVisibleCount.contents().storeBytes(of: UInt32(splatCount), as: UInt32.self)
+
     guard let encodedSplatBuffer = renderInfo.device.makeBuffer(length: MemoryLayout<EncodedGaussianSplat>.stride * Int(splatCount), options: .storageModeShared) else {
         handleError(.bufferAllocationFailed, "Encoded Gaussian splat buffer is nil")
         return
     }
 
-    let pointer = splatBuffer.contents().bindMemory(
-        to: GaussianSplat.self,
-        capacity: splats.count
-    )
-
     let encodedPointer = encodedSplatBuffer.contents().bindMemory(
         to: EncodedGaussianSplat.self,
         capacity: splats.count
     )
 
     for (index, splat) in splats.enumerated() {
-        pointer[index] = splat
         encodedPointer[index] = encodeGaussianSplatForTBDR(splat)
     }
 
+    let packedSphericalHarmonics: PackedGaussianSphericalHarmonics?
+    do {
+        packedSphericalHarmonics = try asset.sphericalHarmonics.map {
+            try packGaussianSphericalHarmonics($0, splatCount: splats.count)
+        }
+    } catch {
+        handleError(.assetDataMissing, "Failed to pack spherical harmonics from \(filename): \(error.localizedDescription)")
+        return
+    }
+
+    let sphericalHarmonicsBuffer: MTLBuffer?
+    if let packedSphericalHarmonics, !packedSphericalHarmonics.coefficients.isEmpty {
+        sphericalHarmonicsBuffer = renderInfo.device.makeBuffer(
+            bytes: packedSphericalHarmonics.coefficients,
+            length: packedSphericalHarmonics.coefficients.count * MemoryLayout<Float16>.stride,
+            options: .storageModeShared
+        )
+        guard sphericalHarmonicsBuffer != nil else {
+            handleError(.bufferAllocationFailed, "Gaussian spherical-harmonics buffer is nil")
+            return
+        }
+        sphericalHarmonicsBuffer?.label = "Gaussian Spherical Harmonics"
+    } else {
+        sphericalHarmonicsBuffer = nil
+    }
+
     let spaceUniform = (0 ..< totalPerMeshUniformBuffers()).compactMap { _ in
         renderInfo.device.makeBuffer(length: MemoryLayout<Uniforms>.stride,
                                      options: [MTLResourceOptions.storageModeShared])
@@ -2970,13 +3011,69 @@ public func setEntityGaussian(entityId: EntityID, filename: String, withExtensio
         }
 
         gaussianComponent.splatCount = splatCount
+        gaussianComponent.visibleSplatCountForRendering = splatCount
         gaussianComponent.gaussianSortedIndices = gaussianSortedIndices
-        gaussianComponent.splatData = splatBuffer
+        gaussianComponent.gaussianVisibleIndices = gaussianVisibleIndices
+        gaussianComponent.gaussianVisibleCount = gaussianVisibleCount
         gaussianComponent.encodedSplatData = encodedSplatBuffer
+        gaussianComponent.sphericalHarmonicsData = sphericalHarmonicsBuffer
+        gaussianComponent.sphericalHarmonicsMetadata = packedSphericalHarmonics?.metadata
         gaussianComponent.spaceUniform = spaceUniform
     }
 }
 
+struct PackedGaussianSphericalHarmonics {
+    let coefficients: [Float16]
+    let metadata: GaussianSHMetadata
+}
+
+/// Packs higher-order SH coefficients to the GPU contract while leaving DC
+/// color in `EncodedGaussianSplat`. Input and output are both channel-major.
+func packGaussianSphericalHarmonics(
+    _ sphericalHarmonics: GaussianSphericalHarmonics,
+    splatCount: Int
+) throws -> PackedGaussianSphericalHarmonics {
+    let coefficientsPerChannel = sphericalHarmonics.coefficientsPerChannel
+    let higherOrderPerChannel = coefficientsPerChannel - 1
+    let inputPerSplat = sphericalHarmonics.coefficientsPerSplat
+    let expectedInputCount = splatCount * inputPerSplat
+
+    guard (0 ... 3).contains(sphericalHarmonics.degree),
+          coefficientsPerChannel == (sphericalHarmonics.degree + 1) * (sphericalHarmonics.degree + 1),
+          sphericalHarmonics.coefficients.count == expectedInputCount
+    else {
+        throw PLYError.invalidData("Spherical-harmonic data does not match its degree or splat count")
+    }
+
+    let outputPerSplat = higherOrderPerChannel * 3
+    var packed: [Float16] = []
+    packed.reserveCapacity(splatCount * outputPerSplat)
+
+    for splatIndex in 0 ..< splatCount {
+        let splatBase = splatIndex * inputPerSplat
+        for channel in 0 ..< 3 {
+            let channelBase = splatBase + channel * coefficientsPerChannel
+            for coefficient in 1 ..< coefficientsPerChannel {
+                let value = Float16(sphericalHarmonics.coefficients[channelBase + coefficient])
+                guard value.isFinite else {
+                    throw PLYError.invalidData("Spherical-harmonic coefficient cannot be represented as Float16")
+                }
+                packed.append(value)
+            }
+        }
+    }
+
+    return PackedGaussianSphericalHarmonics(
+        coefficients: packed,
+        metadata: GaussianSHMetadata(
+            degree: UInt32(sphericalHarmonics.degree),
+            coefficientsPerChannel: UInt32(coefficientsPerChannel),
+            higherOrderCoefficientsPerSplat: UInt32(outputPerSplat),
+            _pad0: 0
+        )
+    )
+}
+
 private func encodeGaussianSplatForTBDR(_ splat: GaussianSplat) -> EncodedGaussianSplat {
     let scale = simd_float3(splat.scale.x, splat.scale.y, splat.scale.z)
     let rotation = simd_quatf(
@@ -3138,9 +3235,13 @@ func removeEntityLOD(entityId: EntityID) {
 func removeEntityGaussian(entityId: EntityID) {
     if let gaussianComponent = scene.get(component: GaussianComponent.self, for: entityId) {
         // Release Metal buffers
-        gaussianComponent.splatData = nil
         gaussianComponent.encodedSplatData = nil
+        gaussianComponent.sphericalHarmonicsData = nil
+        gaussianComponent.sphericalHarmonicsMetadata = nil
         gaussianComponent.gaussianSortedIndices = nil
+        gaussianComponent.gaussianVisibleIndices = nil
+        gaussianComponent.gaussianVisibleCount = nil
+        gaussianComponent.visibleSplatCountForRendering = 0
         gaussianComponent.spaceUniform.removeAll()
         scene.remove(component: GaussianComponent.self, from: entityId)
     }
diff --git a/Sources/UntoldEngine/Systems/RenderingSystem.swift b/Sources/UntoldEngine/Systems/RenderingSystem.swift
index 2114d7a8..86677bd2 100644
--- a/Sources/UntoldEngine/Systems/RenderingSystem.swift
+++ b/Sources/UntoldEngine/Systems/RenderingSystem.swift
@@ -70,6 +70,7 @@ func UpdateRenderingSystem(in view: MTKView) {
                 }
             #endif
 
+            executeGaussianFrustumCulling(commandBuffer)
             executeGaussianDepth(commandBuffer)
             executeRadixSort(commandBuffer)
             EngineProfiler.shared.endScope(.renderPrep)
diff --git a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-ios.air b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-ios.air
index e8eebb77..2ddcd3d6 100644
Binary files a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-ios.air and b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-ios.air differ
diff --git a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-ios.metallib b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-ios.metallib
index a7babd97..c469644c 100644
Binary files a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-ios.metallib and b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-ios.metallib differ
diff --git a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-iossim.metallib b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-iossim.metallib
index c8fd5a41..90fc676c 100644
Binary files a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-iossim.metallib and b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-iossim.metallib differ
diff --git a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-tvos.air b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-tvos.air
index 25da8a37..c2e480b9 100644
Binary files a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-tvos.air and b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-tvos.air differ
diff --git a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-tvos.metallib b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-tvos.metallib
index 220379b6..4ace51cd 100644
Binary files a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-tvos.metallib and b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-tvos.metallib differ
diff --git a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-tvossim.air b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-tvossim.air
index cf3dc327..4ebbff02 100644
Binary files a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-tvossim.air and b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-tvossim.air differ
diff --git a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-tvossim.metallib b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-tvossim.metallib
index f8c0b986..233d6122 100644
Binary files a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-tvossim.metallib and b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-tvossim.metallib differ
diff --git a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-xros.air b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-xros.air
index cd2c4a6e..dafe650b 100644
Binary files a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-xros.air and b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-xros.air differ
diff --git a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-xros.metallib b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-xros.metallib
index 96c62a16..b055cdf4 100644
Binary files a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-xros.metallib and b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-xros.metallib differ
diff --git a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-xrossim.air b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-xrossim.air
index f125450a..9bf3f0d8 100644
Binary files a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-xrossim.air and b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-xrossim.air differ
diff --git a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-xrossim.metallib b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-xrossim.metallib
index 277bc933..c9cc750d 100644
Binary files a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-xrossim.metallib and b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels-xrossim.metallib differ
diff --git a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels.metallib b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels.metallib
index 8f555696..3e5c00b2 100644
Binary files a/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels.metallib and b/Sources/UntoldEngine/UntoldEngineKernels/UntoldEngineKernels.metallib differ
diff --git a/Sources/UntoldEngine/Utils/GaussianSphericalHarmonics.swift b/Sources/UntoldEngine/Utils/GaussianSphericalHarmonics.swift
new file mode 100644
index 00000000..e493be7e
--- /dev/null
+++ b/Sources/UntoldEngine/Utils/GaussianSphericalHarmonics.swift
@@ -0,0 +1,127 @@
+//
+//  GaussianSphericalHarmonics.swift
+//  UntoldEngine
+//
+// Copyright (C) Untold Engine Studios
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at https://mozilla.org/MPL/2.0/.
+
+import simd
+
+private let gaussianSHC1: Float = 0.4886025119029199
+private let gaussianSHC2: [Float] = [
+    1.0925484305920792,
+    -1.0925484305920792,
+    0.31539156525252005,
+    -1.0925484305920792,
+    0.5462742152960396,
+]
+private let gaussianSHC3: [Float] = [
+    -0.5900435899266435,
+    2.890611442640554,
+    -0.4570457994644658,
+    0.3731763325901154,
+    -0.4570457994644658,
+    1.445305721320277,
+    -0.5900435899266435,
+]
+
+/// Decodes the display-referred RGB values represented by standard 3DGS SH
+/// coefficients before they enter the engine's linear working color space.
+func gaussianSRGBToLinear(_ color: simd_float3) -> simd_float3 {
+    let nonnegative = simd_max(color, .zero)
+    return simd_float3(
+        nonnegative.x <= 0.04045
+            ? nonnegative.x / 12.92
+            : pow((nonnegative.x + 0.055) / 1.055, 2.4),
+        nonnegative.y <= 0.04045
+            ? nonnegative.y / 12.92
+            : pow((nonnegative.y + 0.055) / 1.055, 2.4),
+        nonnegative.z <= 0.04045
+            ? nonnegative.z / 12.92
+            : pow((nonnegative.z + 0.055) / 1.055, 2.4)
+    )
+}
+
+/// Converts a world-space camera position to the Gaussian asset's local space.
+func gaussianLocalCameraPosition(
+    cameraWorldPosition: simd_float3,
+    modelMatrix: simd_float4x4
+) -> simd_float3 {
+    let local = modelMatrix.inverse * simd_float4(cameraWorldPosition, 1)
+    return simd_float3(local.x, local.y, local.z)
+}
+
+/// CPU reference for the degree 1–3 terms used by `Gaussians.metal`.
+/// The DC contribution and 0.5 bias must already be present in `baseColor`.
+func evaluateGaussianSphericalHarmonics(
+    baseColor: simd_float3,
+    higherOrderCoefficients: [Float],
+    degree: Int,
+    direction: simd_float3
+) -> simd_float3 {
+    guard degree > 0 else { return baseColor }
+
+    let coefficientsPerChannel = (degree + 1) * (degree + 1) - 1
+    guard (1 ... 3).contains(degree),
+          higherOrderCoefficients.count == coefficientsPerChannel * 3
+    else {
+        return baseColor
+    }
+
+    let directionLengthSquared = simd_length_squared(direction)
+    guard directionLengthSquared > Float.ulpOfOne else { return baseColor }
+    let d = direction * simd_rsqrt(directionLengthSquared)
+    let x = d.x
+    let y = d.y
+    let z = d.z
+
+    func coefficient(_ index: Int) -> simd_float3 {
+        let offset = index - 1
+        return simd_float3(
+            higherOrderCoefficients[offset],
+            higherOrderCoefficients[coefficientsPerChannel + offset],
+            higherOrderCoefficients[2 * coefficientsPerChannel + offset]
+        )
+    }
+
+    var result = baseColor
+    result += gaussianSHC1 * (-y * coefficient(1) + z * coefficient(2) - x * coefficient(3))
+
+    if degree > 1 {
+        let xx = x * x
+        let yy = y * y
+        let zz = z * z
+        let basis4 = gaussianSHC2[0] * x * y
+        let basis5 = gaussianSHC2[1] * y * z
+        let basis6 = gaussianSHC2[2] * (2 * zz - xx - yy)
+        let basis7 = gaussianSHC2[3] * x * z
+        let basis8 = gaussianSHC2[4] * (xx - yy)
+        result += basis4 * coefficient(4)
+        result += basis5 * coefficient(5)
+        result += basis6 * coefficient(6)
+        result += basis7 * coefficient(7)
+        result += basis8 * coefficient(8)
+
+        if degree > 2 {
+            let basis9 = gaussianSHC3[0] * y * (3 * xx - yy)
+            let basis10 = gaussianSHC3[1] * x * y * z
+            let basis11 = gaussianSHC3[2] * y * (4 * zz - xx - yy)
+            let basis12 = gaussianSHC3[3] * z * (2 * zz - 3 * xx - 3 * yy)
+            let basis13 = gaussianSHC3[4] * x * (4 * zz - xx - yy)
+            let basis14 = gaussianSHC3[5] * z * (xx - yy)
+            let basis15 = gaussianSHC3[6] * x * (xx - 3 * yy)
+            result += basis9 * coefficient(9)
+            result += basis10 * coefficient(10)
+            result += basis11 * coefficient(11)
+            result += basis12 * coefficient(12)
+            result += basis13 * coefficient(13)
+            result += basis14 * coefficient(14)
+            result += basis15 * coefficient(15)
+        }
+    }
+
+    return simd_max(result, .zero)
+}
diff --git a/Sources/UntoldEngine/Utils/Globals.swift b/Sources/UntoldEngine/Utils/Globals.swift
index e522dca9..5be9bd12 100644
--- a/Sources/UntoldEngine/Utils/Globals.swift
+++ b/Sources/UntoldEngine/Utils/Globals.swift
@@ -49,6 +49,8 @@ private final class CoreRuntimeGlobals: @unchecked Sendable {
     var reduceScanScatterCompactedPipeline = ComputePipeline()
     var hzbBuildPyramidPipeline = ComputePipeline()
     var hzbOcclusionCullingPipeline = ComputePipeline()
+    var gaussianResetVisibleCountPipeline = ComputePipeline()
+    var gaussianFrustumCullPipeline = ComputePipeline()
     var gaussianDepthPipeline = ComputePipeline()
     var radixClearHistogramPipeline = ComputePipeline()
     var radixHistogramPipeline = ComputePipeline()
@@ -489,6 +491,48 @@ var gaussianDepthPipeline: ComputePipeline {
     }
 }
 
+var gaussianResetVisibleCountPipeline: ComputePipeline {
+    get {
+        let state = CoreRuntimeGlobals.shared
+        state.lock.lock()
+        defer { state.lock.unlock() }
+        return state.gaussianResetVisibleCountPipeline
+    }
+    set {
+        let state = CoreRuntimeGlobals.shared
+        state.lock.lock()
+        state.gaussianResetVisibleCountPipeline = newValue
+        state.lock.unlock()
+    }
+    _modify {
+        let state = CoreRuntimeGlobals.shared
+        state.lock.lock()
+        defer { state.lock.unlock() }
+        yield &state.gaussianResetVisibleCountPipeline
+    }
+}
+
+var gaussianFrustumCullPipeline: ComputePipeline {
+    get {
+        let state = CoreRuntimeGlobals.shared
+        state.lock.lock()
+        defer { state.lock.unlock() }
+        return state.gaussianFrustumCullPipeline
+    }
+    set {
+        let state = CoreRuntimeGlobals.shared
+        state.lock.lock()
+        state.gaussianFrustumCullPipeline = newValue
+        state.lock.unlock()
+    }
+    _modify {
+        let state = CoreRuntimeGlobals.shared
+        state.lock.lock()
+        defer { state.lock.unlock() }
+        yield &state.gaussianFrustumCullPipeline
+    }
+}
+
 var radixClearHistogramPipeline: ComputePipeline {
     get {
         let state = CoreRuntimeGlobals.shared
diff --git a/Sources/UntoldEngine/Utils/Logger.swift b/Sources/UntoldEngine/Utils/Logger.swift
index 08068056..1aedfaaa 100644
--- a/Sources/UntoldEngine/Utils/Logger.swift
+++ b/Sources/UntoldEngine/Utils/Logger.swift
@@ -35,6 +35,7 @@ public enum LogCategory: String, CaseIterable, Sendable {
     case xrCamera = "XRCamera"
     case batching = "Batching"
     case lightPortal = "LightPortal"
+    case gaussian = "Gaussian"
 }
 
 public struct LogEvent: Identifiable, Sendable {
@@ -201,6 +202,7 @@ public enum Logger {
             LogCategory.xrCamera.rawValue,
             LogCategory.batching.rawValue,
             LogCategory.lightPortal.rawValue,
+            LogCategory.gaussian.rawValue,
         ]
         private var categoryOverrides: [String: Bool] = [:]
 
diff --git a/Sources/UntoldEngine/Utils/PLYReader.swift b/Sources/UntoldEngine/Utils/PLYReader.swift
index 5d8be9a1..4c98f754 100644
--- a/Sources/UntoldEngine/Utils/PLYReader.swift
+++ b/Sources/UntoldEngine/Utils/PLYReader.swift
@@ -44,6 +44,27 @@ public struct PLYHeader {
     var comments: [String]
 }
 
+/// Contiguous spherical-harmonic coefficients for a Gaussian asset.
+///
+/// Coefficients are stored per splat in channel-major order:
+/// `R[0 ... n], G[0 ... n], B[0 ... n]`, where coefficient zero is the DC
+/// term and `n + 1 == coefficientsPerChannel`.
+public struct GaussianSphericalHarmonics: Sendable {
+    public let degree: Int
+    public let coefficientsPerChannel: Int
+    public let coefficients: [Float]
+
+    public var coefficientsPerSplat: Int {
+        coefficientsPerChannel * 3
+    }
+}
+
+/// CPU import representation for a Gaussian PLY asset.
+public struct GaussianSplatAsset {
+    public let splats: [GaussianSplat]
+    public let sphericalHarmonics: GaussianSphericalHarmonics?
+}
+
 public class PLYReader {
     // MARK: - Public Methods
 
@@ -52,6 +73,11 @@ public class PLYReader {
     /// - Returns: Array of GaussianSplat structs
     /// - Throws: Error if file cannot be read or parsed
     public static func readGaussianSplats(from url: URL) throws -> [GaussianSplat] {
+        try readGaussianAsset(from: url).splats
+    }
+
+    /// Reads Gaussian geometry and preserves all spherical-harmonic coefficients.
+    public static func readGaussianAsset(from url: URL) throws -> GaussianSplatAsset {
         let data = try Data(contentsOf: url)
         let (header, bodyOffset) = try parseHeader(from: data)
 
@@ -60,15 +86,32 @@ public class PLYReader {
             throw PLYError.missingElement("vertex")
         }
 
+        let shSchema = try sphericalHarmonicSchema(for: vertexElement.properties)
+
         // Parse the body based on format
+        let parsed: ([GaussianSplat], [Float])
         switch header.format {
         case .ascii:
-            return try parseASCIIGaussians(data: data, bodyOffset: bodyOffset, element: vertexElement)
+            parsed = try parseASCIIGaussians(data: data, bodyOffset: bodyOffset, element: vertexElement, shSchema: shSchema)
         case .binaryLittleEndian:
-            return try parseBinaryGaussians(data: data, bodyOffset: bodyOffset, element: vertexElement, bigEndian: false)
+            parsed = try parseBinaryGaussians(data: data, bodyOffset: bodyOffset, element: vertexElement, bigEndian: false, shSchema: shSchema)
         case .binaryBigEndian:
-            return try parseBinaryGaussians(data: data, bodyOffset: bodyOffset, element: vertexElement, bigEndian: true)
+            parsed = try parseBinaryGaussians(data: data, bodyOffset: bodyOffset, element: vertexElement, bigEndian: true, shSchema: shSchema)
         }
+
+        let sphericalHarmonics = shSchema.map {
+            GaussianSphericalHarmonics(
+                degree: $0.degree,
+                coefficientsPerChannel: $0.coefficientsPerChannel,
+                coefficients: parsed.1
+            )
+        }
+        if let shSchema,
+           parsed.1.count != vertexElement.count * shSchema.coefficientsPerSplat
+        {
+            throw PLYError.invalidData("Spherical-harmonic coefficient data is incomplete")
+        }
+        return GaussianSplatAsset(splats: parsed.0, sphericalHarmonics: sphericalHarmonics)
     }
 
     // MARK: - Header Parsing
@@ -197,7 +240,12 @@ public class PLYReader {
 
     // MARK: - ASCII Parsing
 
-    private static func parseASCIIGaussians(data: Data, bodyOffset: Int, element: PLYElement) throws -> [GaussianSplat] {
+    private static func parseASCIIGaussians(
+        data: Data,
+        bodyOffset: Int,
+        element: PLYElement,
+        shSchema: SphericalHarmonicSchema?
+    ) throws -> ([GaussianSplat], [Float]) {
         // Get the body string
         let bodyData = data.suffix(from: bodyOffset)
         guard let bodyString = String(data: bodyData, encoding: .utf8) else {
@@ -207,6 +255,10 @@ public class PLYReader {
         let lines = bodyString.components(separatedBy: .newlines).filter { !$0.isEmpty }
         var splats: [GaussianSplat] = []
         splats.reserveCapacity(element.count)
+        var shCoefficients: [Float] = []
+        if let shSchema {
+            shCoefficients.reserveCapacity(element.count * shSchema.coefficientsPerSplat)
+        }
 
         // Create property index map
         let propertyMap = createPropertyIndexMap(properties: element.properties)
@@ -218,18 +270,36 @@ public class PLYReader {
 
             if values.isEmpty { continue }
 
-            let splat = try parseGaussianFromValues(values: values, propertyMap: propertyMap)
+            let splat = try parseGaussianFromValues(
+                values: values,
+                propertyMap: propertyMap,
+                shSchema: shSchema,
+                shCoefficients: &shCoefficients
+            )
             splats.append(splat)
         }
 
-        return splats
+        guard splats.count == element.count else {
+            throw PLYError.invalidData("Expected \(element.count) Gaussian vertices, found \(splats.count)")
+        }
+        return (splats, shCoefficients)
     }
 
     // MARK: - Binary Parsing
 
-    private static func parseBinaryGaussians(data: Data, bodyOffset: Int, element: PLYElement, bigEndian: Bool) throws -> [GaussianSplat] {
+    private static func parseBinaryGaussians(
+        data: Data,
+        bodyOffset: Int,
+        element: PLYElement,
+        bigEndian: Bool,
+        shSchema: SphericalHarmonicSchema?
+    ) throws -> ([GaussianSplat], [Float]) {
         var splats: [GaussianSplat] = []
         splats.reserveCapacity(element.count)
+        var shCoefficients: [Float] = []
+        if let shSchema {
+            shCoefficients.reserveCapacity(element.count * shSchema.coefficientsPerSplat)
+        }
 
         // Calculate stride for each vertex
         let stride = calculateStride(properties: element.properties)
@@ -249,14 +319,16 @@ public class PLYReader {
                 properties: element.properties,
                 propertyMap: propertyMap,
                 propertyOffsets: propertyOffsets,
-                bigEndian: bigEndian
+                bigEndian: bigEndian,
+                shSchema: shSchema,
+                shCoefficients: &shCoefficients
             )
             splats.append(splat)
 
             offset += stride
         }
 
-        return splats
+        return (splats, shCoefficients)
     }
 
     // MARK: - Gaussian Parsing Helpers
@@ -269,7 +341,59 @@ public class PLYReader {
         return map
     }
 
-    private static func parseGaussianFromValues(values: [String], propertyMap: [String: Int]) throws -> GaussianSplat {
+    private struct SphericalHarmonicSchema {
+        let degree: Int
+        let coefficientsPerChannel: Int
+        let restPropertyNames: [String]
+
+        var coefficientsPerSplat: Int {
+            coefficientsPerChannel * 3
+        }
+    }
+
+    private static func sphericalHarmonicSchema(for properties: [PLYProperty]) throws -> SphericalHarmonicSchema? {
+        let names = Set(properties.map(\.name))
+        let dcNames = (0 ..< 3).map { "f_dc_\($0)" }
+        let dcCount = dcNames.filter(names.contains).count
+        let restProperties = properties.filter { $0.name.hasPrefix("f_rest_") }
+        let parsedRestIndices = restProperties.map { Int($0.name.dropFirst("f_rest_".count)) }
+        guard parsedRestIndices.allSatisfy({ $0 != nil }) else {
+            throw PLYError.invalidData("f_rest_* properties must end in a numeric index")
+        }
+        let restIndices = parsedRestIndices.compactMap { $0 }.sorted()
+
+        guard dcCount == 0 || dcCount == 3 else {
+            throw PLYError.invalidData("Spherical harmonics require f_dc_0, f_dc_1, and f_dc_2")
+        }
+        guard dcCount == 3 || restIndices.isEmpty else {
+            throw PLYError.invalidData("f_rest_* properties require all three f_dc_* properties")
+        }
+        guard dcCount == 3 else { return nil }
+
+        guard restIndices == Array(0 ..< restIndices.count) else {
+            throw PLYError.invalidData("f_rest_* property indices must be contiguous starting at zero")
+        }
+
+        let supportedRestCounts = [0: 0, 9: 1, 24: 2, 45: 3]
+        guard let degree = supportedRestCounts[restIndices.count] else {
+            throw PLYError.invalidData(
+                "Unsupported spherical-harmonic coefficient count: \(restIndices.count) f_rest_* values"
+            )
+        }
+
+        return SphericalHarmonicSchema(
+            degree: degree,
+            coefficientsPerChannel: (degree + 1) * (degree + 1),
+            restPropertyNames: restIndices.map { "f_rest_\($0)" }
+        )
+    }
+
+    private static func parseGaussianFromValues(
+        values: [String],
+        propertyMap: [String: Int],
+        shSchema: SphericalHarmonicSchema?,
+        shCoefficients: inout [Float]
+    ) throws -> GaussianSplat {
         // Extract position
         let x = try getFloat(values: values, propertyMap: propertyMap, key: "x")
         let y = try getFloat(values: values, propertyMap: propertyMap, key: "y")
@@ -287,17 +411,36 @@ public class PLYReader {
 
         // Extract color (from spherical harmonics DC component or direct RGB)
         var r: Float, g: Float, b: Float
-        if let dcR = try? getFloat(values: values, propertyMap: propertyMap, key: "f_dc_0") {
+        if let shSchema {
             // Color from spherical harmonics
+            let dcR = try getFloat(values: values, propertyMap: propertyMap, key: "f_dc_0")
+            let dcG = try getFloat(values: values, propertyMap: propertyMap, key: "f_dc_1")
+            let dcB = try getFloat(values: values, propertyMap: propertyMap, key: "f_dc_2")
             r = dcR
-            g = try getFloat(values: values, propertyMap: propertyMap, key: "f_dc_1")
-            b = try getFloat(values: values, propertyMap: propertyMap, key: "f_dc_2")
+            g = dcG
+            b = dcB
 
             // Convert from SH to RGB (DC component of SH corresponds to RGB / C0 where C0 = 0.28209479177387814)
             let C0: Float = 0.28209479177387814
             r = (r * C0 + 0.5)
             g = (g * C0 + 0.5)
             b = (b * C0 + 0.5)
+
+            let dc = (dcR, dcG, dcB)
+            let restPerChannel = shSchema.coefficientsPerChannel - 1
+            for channel in 0 ..< 3 {
+                shCoefficients.append(channel == 0 ? dc.0 : (channel == 1 ? dc.1 : dc.2))
+                let start = channel * restPerChannel
+                for index in start ..< start + restPerChannel {
+                    try shCoefficients.append(
+                        getFloat(
+                            values: values,
+                            propertyMap: propertyMap,
+                            key: shSchema.restPropertyNames[index]
+                        )
+                    )
+                }
+            }
         } else {
             // Direct RGB
             r = try getFloat(values: values, propertyMap: propertyMap, key: "red", default: 1.0) / 255.0
@@ -387,7 +530,9 @@ public class PLYReader {
         properties: [PLYProperty],
         propertyMap: [String: Int],
         propertyOffsets: [Int],
-        bigEndian: Bool
+        bigEndian: Bool,
+        shSchema: SphericalHarmonicSchema?,
+        shCoefficients: inout [Float]
     ) throws -> GaussianSplat {
         func readFloat(propertyName: String, default defaultValue: Float? = nil) throws -> Float {
             guard let index = propertyMap[propertyName] else {
@@ -405,8 +550,12 @@ public class PLYReader {
                 throw PLYError.invalidData("Buffer overflow reading '\(propertyName)'")
             }
 
-            let bytes = data.subdata(in: offset ..< (offset + size))
-            return try convertToFloat(bytes: bytes, type: property.type, bigEndian: bigEndian)
+            return try convertToFloat(
+                data: data,
+                offset: offset,
+                type: property.type,
+                bigEndian: bigEndian
+            )
         }
 
         // Extract all properties
@@ -423,15 +572,30 @@ public class PLYReader {
         let scaleZ = exp(scale2)
 
         var r: Float, g: Float, b: Float
-        if let dcR = try? readFloat(propertyName: "f_dc_0") {
+        if let shSchema {
+            let dcR = try readFloat(propertyName: "f_dc_0")
+            let dcG = try readFloat(propertyName: "f_dc_1")
+            let dcB = try readFloat(propertyName: "f_dc_2")
             r = dcR
-            g = try readFloat(propertyName: "f_dc_1")
-            b = try readFloat(propertyName: "f_dc_2")
+            g = dcG
+            b = dcB
 
             let C0: Float = 0.28209479177387814
             r = (r * C0 + 0.5)
             g = (g * C0 + 0.5)
             b = (b * C0 + 0.5)
+
+            let dc = (dcR, dcG, dcB)
+            let restPerChannel = shSchema.coefficientsPerChannel - 1
+            for channel in 0 ..< 3 {
+                shCoefficients.append(channel == 0 ? dc.0 : (channel == 1 ? dc.1 : dc.2))
+                let start = channel * restPerChannel
+                for index in start ..< start + restPerChannel {
+                    try shCoefficients.append(
+                        readFloat(propertyName: shSchema.restPropertyNames[index])
+                    )
+                }
+            }
         } else {
             r = try readFloat(propertyName: "red", default: 1.0) / 255.0
             g = try readFloat(propertyName: "green", default: 1.0) / 255.0
@@ -457,64 +621,48 @@ public class PLYReader {
         )
     }
 
-    private static func convertToFloat(bytes: Data, type: String, bigEndian: Bool) throws -> Float {
-        switch type {
-        case "float", "float32":
-            var value: Float = 0
-            _ = withUnsafeMutableBytes(of: &value) { bytes.copyBytes(to: $0) }
-            if bigEndian {
-                value = Float(bitPattern: UInt32(bigEndian: value.bitPattern))
-            }
-            return value
-
-        case "double", "float64":
-            var value: Double = 0
-            _ = withUnsafeMutableBytes(of: &value) { bytes.copyBytes(to: $0) }
-            if bigEndian {
-                value = Double(bitPattern: UInt64(bigEndian: value.bitPattern))
-            }
-            return Float(value)
-
-        case "uchar", "uint8":
-            return Float(bytes[0])
-
-        case "char", "int8":
-            return Float(Int8(bitPattern: bytes[0]))
-
-        case "ushort", "uint16":
-            var value: UInt16 = 0
-            _ = withUnsafeMutableBytes(of: &value) { bytes.copyBytes(to: $0) }
-            if bigEndian {
-                value = UInt16(bigEndian: value)
-            }
-            return Float(value)
-
-        case "short", "int16":
-            var value: Int16 = 0
-            _ = withUnsafeMutableBytes(of: &value) { bytes.copyBytes(to: $0) }
-            if bigEndian {
-                value = Int16(bigEndian: value)
-            }
-            return Float(value)
-
-        case "uint", "uint32":
-            var value: UInt32 = 0
-            _ = withUnsafeMutableBytes(of: &value) { bytes.copyBytes(to: $0) }
-            if bigEndian {
-                value = UInt32(bigEndian: value)
-            }
-            return Float(value)
+    private static func convertToFloat(data: Data, offset: Int, type: String, bigEndian: Bool) throws -> Float {
+        try data.withUnsafeBytes { bytes in
+            switch type {
+            case "float", "float32":
+                var bits = bytes.loadUnaligned(fromByteOffset: offset, as: UInt32.self)
+                if bigEndian { bits = UInt32(bigEndian: bits) }
+                return Float(bitPattern: bits)
+
+            case "double", "float64":
+                var bits = bytes.loadUnaligned(fromByteOffset: offset, as: UInt64.self)
+                if bigEndian { bits = UInt64(bigEndian: bits) }
+                return Float(Double(bitPattern: bits))
+
+            case "uchar", "uint8":
+                return Float(bytes[offset])
+
+            case "char", "int8":
+                return Float(Int8(bitPattern: bytes[offset]))
+
+            case "ushort", "uint16":
+                var value = bytes.loadUnaligned(fromByteOffset: offset, as: UInt16.self)
+                if bigEndian { value = UInt16(bigEndian: value) }
+                return Float(value)
+
+            case "short", "int16":
+                var value = bytes.loadUnaligned(fromByteOffset: offset, as: Int16.self)
+                if bigEndian { value = Int16(bigEndian: value) }
+                return Float(value)
+
+            case "uint", "uint32":
+                var value = bytes.loadUnaligned(fromByteOffset: offset, as: UInt32.self)
+                if bigEndian { value = UInt32(bigEndian: value) }
+                return Float(value)
+
+            case "int", "int32":
+                var value = bytes.loadUnaligned(fromByteOffset: offset, as: Int32.self)
+                if bigEndian { value = Int32(bigEndian: value) }
+                return Float(value)
 
-        case "int", "int32":
-            var value: Int32 = 0
-            _ = withUnsafeMutableBytes(of: &value) { bytes.copyBytes(to: $0) }
-            if bigEndian {
-                value = Int32(bigEndian: value)
+            default:
+                throw PLYError.unsupportedType(type)
             }
-            return Float(value)
-
-        default:
-            throw PLYError.unsupportedType(type)
         }
     }
 }
diff --git a/Tests/UntoldEngineRenderTests/DeviceRadixSortTest.swift b/Tests/UntoldEngineRenderTests/DeviceRadixSortTest.swift
index 09c73017..21665587 100644
--- a/Tests/UntoldEngineRenderTests/DeviceRadixSortTest.swift
+++ b/Tests/UntoldEngineRenderTests/DeviceRadixSortTest.swift
@@ -708,6 +708,7 @@ final class DeviceRadixSortTest: BaseRenderSetup {
         guard let keyBuf = makeBuffer(inputKeys) else { XCTFail("Buffer allocation failed"); return }
         gc.gaussianSortedIndices = keyBuf
         gc.splatCount = UInt(n)
+        gc.visibleSplatCountForRendering = UInt(n)
 
         guard let queue = renderInfo.device.makeCommandQueue(),
               let cmd = queue.makeCommandBuffer()
@@ -731,7 +732,9 @@ final class DeviceRadixSortTest: BaseRenderSetup {
 
     /// Radix sort output must match Swift's sort for 256 random UInt64 values.
     func test_radixSort_matchesCPUSort_256elements() {
-        guard makePipeline(named: "gaussianRadixHistogram") != nil,
+        guard makePipeline(named: "gaussianRadixClearHistogram") != nil,
+              makePipeline(named: "gaussianRadixHistogram") != nil,
+              makePipeline(named: "gaussianRadixScanPerTG") != nil,
               makePipeline(named: "gaussianRadixScan") != nil,
               makePipeline(named: "gaussianRadixScatter") != nil
         else {
@@ -740,7 +743,11 @@ final class DeviceRadixSortTest: BaseRenderSetup {
         }
 
         let n = 256
-        let inputKeys: [UInt64] = (0 ..< n).map { _ in UInt64.random(in: 0 ..< UInt64.max) }
+        let inputKeys: [UInt64] = (0 ..< n).map { i in
+            let high = UInt64((n - 1 - i) * 2 + 1)
+            let low = UInt64((i &* 1_664_525 &+ 1_013_904_223) & 0xFFFF_FFFF)
+            return (high << 32) | low
+        }
         let expectedSorted = inputKeys.sorted()
 
         let entity = createEntity()
@@ -757,6 +764,7 @@ final class DeviceRadixSortTest: BaseRenderSetup {
         guard let keyBuf = makeBuffer(inputKeys) else { XCTFail("Buffer allocation failed"); return }
         gc.gaussianSortedIndices = keyBuf
         gc.splatCount = UInt(n)
+        gc.visibleSplatCountForRendering = UInt(n)
 
         guard let queue = renderInfo.device.makeCommandQueue(),
               let cmd = queue.makeCommandBuffer()
@@ -787,7 +795,11 @@ final class DeviceRadixSortTest: BaseRenderSetup {
         }
 
         let n = 256
-        let inputKeys: [UInt64] = (0 ..< n).map { _ in UInt64.random(in: 0 ..< UInt64.max) }
+        let inputKeys: [UInt64] = (0 ..< n).map { i in
+            let high = UInt64((n - 1 - i) * 2 + 1)
+            let low = UInt64((i &* 1_664_525 &+ 1_013_904_223) & 0xFFFF_FFFF)
+            return (high << 32) | low
+        }
         let expectedSorted = inputKeys.sorted()
 
         // First sort
@@ -804,6 +816,7 @@ final class DeviceRadixSortTest: BaseRenderSetup {
         wt1.space = matrix_identity_float4x4
         gc1.gaussianSortedIndices = buf1
         gc1.splatCount = UInt(n)
+        gc1.visibleSplatCountForRendering = UInt(n)
 
         guard let queue = renderInfo.device.makeCommandQueue(),
               let cmd1 = queue.makeCommandBuffer()
@@ -846,12 +859,15 @@ final class DeviceRadixSortTest: BaseRenderSetup {
         let numSplats = positions.count
 
         let splats = positions.map { pos in
-            GaussianSplat(
-                center: simd_float4(pos, 1.0),
-                scale: simd_float4(1, 1, 1, 1),
-                color: simd_float4(1, 0, 0, 1),
-                quat: simd_float4(0, 0, 0, 1),
-                opacity: 1.0
+            EncodedGaussianSplat(
+                position: pos,
+                opacity: 1.0,
+                color: simd_float3(1, 0, 0),
+                _pad0: 0.0,
+                covA: simd_float3(1, 0, 0),
+                _pad1: 0.0,
+                covB: simd_float3(1, 0, 1),
+                _pad2: 0.0
             )
         }
 
@@ -876,15 +892,28 @@ final class DeviceRadixSortTest: BaseRenderSetup {
 
         guard let splatBuf = renderInfo.device.makeBuffer(
             bytes: splats,
-            length: numSplats * MemoryLayout<GaussianSplat>.stride,
+            length: numSplats * MemoryLayout<EncodedGaussianSplat>.stride,
             options: .storageModeShared
         ),
-            let keyBuf = makeBuffer(count: numSplats, type: UInt64.self)
+            let keyBuf = makeBuffer(count: numSplats, type: UInt64.self),
+            let visibleIndexBuf = renderInfo.device.makeBuffer(
+                bytes: (0 ..< numSplats).map { UInt32($0) },
+                length: numSplats * MemoryLayout<UInt32>.stride,
+                options: .storageModeShared
+            ),
+            let visibleCountBuf = renderInfo.device.makeBuffer(
+                length: MemoryLayout<UInt32>.stride,
+                options: .storageModeShared
+            )
         else { XCTFail("Buffer allocation failed"); return }
 
-        gc.splatData = splatBuf
+        visibleCountBuf.contents().storeBytes(of: UInt32(numSplats), as: UInt32.self)
+        gc.encodedSplatData = splatBuf
         gc.gaussianSortedIndices = keyBuf
+        gc.gaussianVisibleIndices = visibleIndexBuf
+        gc.gaussianVisibleCount = visibleCountBuf
         gc.splatCount = UInt(numSplats)
+        gc.visibleSplatCountForRendering = UInt(numSplats)
         gc.spaceUniform = (0 ..< 2).compactMap { _ in
             renderInfo.device.makeBuffer(
                 length: MemoryLayout<Uniforms>.stride,
@@ -932,7 +961,11 @@ final class DeviceRadixSortTest: BaseRenderSetup {
         }
 
         let n = 1024
-        let inputKeys: [UInt64] = (0 ..< n).map { _ in UInt64.random(in: 0 ..< UInt64.max) }
+        let inputKeys: [UInt64] = (0 ..< n).map { i in
+            let high = UInt64((n - 1 - i) * 2 + 1)
+            let low = UInt64((i &* 1_664_525 &+ 1_013_904_223) & 0xFFFF_FFFF)
+            return (high << 32) | low
+        }
         let expectedSorted = inputKeys.sorted()
 
         let entity = createEntity()
@@ -948,6 +981,7 @@ final class DeviceRadixSortTest: BaseRenderSetup {
         wt.space = matrix_identity_float4x4
         gc.gaussianSortedIndices = keyBuf
         gc.splatCount = UInt(n)
+        gc.visibleSplatCountForRendering = UInt(n)
 
         guard let queue = renderInfo.device.makeCommandQueue(),
               let cmd = queue.makeCommandBuffer()
diff --git a/Tests/UntoldEngineRenderTests/GaussianRenderingTest.swift b/Tests/UntoldEngineRenderTests/GaussianRenderingTest.swift
index bda050b1..00f964c7 100644
--- a/Tests/UntoldEngineRenderTests/GaussianRenderingTest.swift
+++ b/Tests/UntoldEngineRenderTests/GaussianRenderingTest.swift
@@ -8,6 +8,7 @@
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
 
+import CShaderTypes
 import Foundation
 import Metal
 import simd
@@ -202,7 +203,183 @@ final class GaussianRenderingTest: BaseRenderSetup {
 
         XCTAssertNotNil(gaussianComponent.spaceUniform,
                         "Gaussian component should have space uniform array")
-        // Note: splatData and gaussianSortedIndices may be nil until loaded
+        // Note: encodedSplatData and gaussianSortedIndices may be nil until loaded
+    }
+
+    func testLoadedGaussianUsesExactSizeGPUBufferAllocations() {
+        let transformId = getComponentId(for: WorldTransformComponent.self)
+        let gaussianId = getComponentId(for: GaussianComponent.self)
+        let entities = queryEntitiesWithComponentIds([transformId, gaussianId], in: scene)
+
+        guard let entity = entities.first,
+              let component = scene.get(component: GaussianComponent.self, for: entity),
+              let encodedSplats = component.encodedSplatData,
+              let sortedIndices = component.gaussianSortedIndices
+        else {
+            XCTFail("Expected the Gaussian test asset to be loaded")
+            return
+        }
+
+        let count = Int(component.splatCount)
+        XCTAssertGreaterThan(count, 0)
+        XCTAssertEqual(
+            encodedSplats.length,
+            count * MemoryLayout<EncodedGaussianSplat>.stride
+        )
+        XCTAssertEqual(
+            sortedIndices.length,
+            count * MemoryLayout<UInt64>.stride
+        )
+
+        let metadata = component.sphericalHarmonicsMetadata
+        XCTAssertEqual(metadata?.degree, 0)
+        XCTAssertEqual(metadata?.coefficientsPerChannel, 1)
+        XCTAssertEqual(metadata?.higherOrderCoefficientsPerSplat, 0)
+        XCTAssertNil(component.sphericalHarmonicsData)
+    }
+
+    func testRemoveGaussianReleasesSphericalHarmonicsState() {
+        let entity = createEntity()
+        guard let component = scene.assign(to: entity, component: GaussianComponent.self) else {
+            XCTFail("Expected Gaussian component")
+            return
+        }
+        component.sphericalHarmonicsData = renderInfo.device.makeBuffer(
+            length: MemoryLayout<Float16>.stride,
+            options: .storageModeShared
+        )
+        component.sphericalHarmonicsMetadata = GaussianSHMetadata(
+            degree: 1,
+            coefficientsPerChannel: 4,
+            higherOrderCoefficientsPerSplat: 9,
+            _pad0: 0
+        )
+
+        removeEntityGaussian(entityId: entity)
+
+        XCTAssertNil(scene.get(component: GaussianComponent.self, for: entity))
+    }
+
+    func testPackedSphericalHarmonicsRoundTripsThroughMetalBuffer() throws {
+        let sphericalHarmonics = GaussianSphericalHarmonics(
+            degree: 1,
+            coefficientsPerChannel: 4,
+            coefficients: (0 ..< 12).map { Float($0) * 0.125 }
+        )
+        let packed = try packGaussianSphericalHarmonics(sphericalHarmonics, splatCount: 1)
+        guard let buffer = renderInfo.device.makeBuffer(
+            bytes: packed.coefficients,
+            length: packed.coefficients.count * MemoryLayout<Float16>.stride,
+            options: .storageModeShared
+        ) else {
+            XCTFail("Expected SH buffer allocation")
+            return
+        }
+
+        XCTAssertEqual(buffer.length, 9 * MemoryLayout<Float16>.stride)
+        let pointer = buffer.contents().bindMemory(to: Float16.self, capacity: 9)
+        let roundTripped = (0 ..< 9).map { pointer[$0] }
+        XCTAssertEqual(roundTripped, packed.coefficients)
+    }
+
+    func testPackedSphericalHarmonicsMatchesActiveMetalEvaluator() throws {
+        guard let library = renderInfo.library,
+              let function = library.makeFunction(name: "gaussianSphericalHarmonicsDiagnostic")
+        else {
+            XCTFail("Gaussian SH diagnostic kernel is missing from the active metallib")
+            return
+        }
+
+        let pipeline = try renderInfo.device.makeComputePipelineState(function: function)
+        let sphericalHarmonics = GaussianSphericalHarmonics(
+            degree: 3,
+            coefficientsPerChannel: 16,
+            coefficients: (0 ..< 48).map { Float($0 - 24) * 0.0078125 }
+        )
+        let packed = try packGaussianSphericalHarmonics(sphericalHarmonics, splatCount: 1)
+        let baseColor = simd_float4(0.35, 0.5, 0.65, 1)
+        let direction = simd_float4(0.25, -0.5, 0.75, 0)
+        let cpuResult = evaluateGaussianSphericalHarmonics(
+            baseColor: simd_float3(baseColor.x, baseColor.y, baseColor.z),
+            higherOrderCoefficients: packed.coefficients.map(Float.init),
+            degree: 3,
+            direction: simd_float3(direction.x, direction.y, direction.z)
+        )
+
+        guard let coefficients = renderInfo.device.makeBuffer(
+            bytes: packed.coefficients,
+            length: packed.coefficients.count * MemoryLayout<Float16>.stride,
+            options: .storageModeShared
+        ), let output = renderInfo.device.makeBuffer(
+            length: 2 * MemoryLayout<simd_float4>.stride,
+            options: .storageModeShared
+        ), let commandBuffer = renderInfo.commandQueue.makeCommandBuffer(),
+        let encoder = commandBuffer.makeComputeCommandEncoder()
+        else {
+            XCTFail("Failed to allocate Gaussian SH diagnostic resources")
+            return
+        }
+
+        var metadata = packed.metadata
+        var baseColorArgument = baseColor
+        var directionArgument = direction
+        encoder.setComputePipelineState(pipeline)
+        encoder.setBuffer(coefficients, offset: 0, index: 0)
+        encoder.setBytes(&metadata, length: MemoryLayout<GaussianSHMetadata>.stride, index: 1)
+        encoder.setBytes(&baseColorArgument, length: MemoryLayout<simd_float4>.stride, index: 2)
+        encoder.setBytes(&directionArgument, length: MemoryLayout<simd_float4>.stride, index: 3)
+        encoder.setBuffer(output, offset: 0, index: 4)
+        encoder.dispatchThreads(MTLSize(width: 1, height: 1, depth: 1),
+                                threadsPerThreadgroup: MTLSize(width: 1, height: 1, depth: 1))
+        encoder.endEncoding()
+        commandBuffer.commit()
+        commandBuffer.waitUntilCompleted()
+
+        XCTAssertEqual(commandBuffer.status, .completed)
+        let gpuResults = output.contents().bindMemory(to: simd_float4.self, capacity: 2)
+        let gpuResult = gpuResults[0]
+        XCTAssertEqual(gpuResult.x, cpuResult.x, accuracy: 1e-5)
+        XCTAssertEqual(gpuResult.y, cpuResult.y, accuracy: 1e-5)
+        XCTAssertEqual(gpuResult.z, cpuResult.z, accuracy: 1e-5)
+        XCTAssertNotEqual(simd_float3(gpuResult.x, gpuResult.y, gpuResult.z),
+                          simd_float3(baseColor.x, baseColor.y, baseColor.z))
+
+        let expectedLinear = gaussianSRGBToLinear(cpuResult)
+        let gpuLinear = gpuResults[1]
+        XCTAssertEqual(gpuLinear.x, expectedLinear.x, accuracy: 1e-5)
+        XCTAssertEqual(gpuLinear.y, expectedLinear.y, accuracy: 1e-5)
+        XCTAssertEqual(gpuLinear.z, expectedLinear.z, accuracy: 1e-5)
+    }
+
+    func testExternalGaussianDiagnosticAssetReachesRuntimeGPUBuffers() throws {
+        guard let path = ProcessInfo.processInfo.environment["UNTOLD_GAUSSIAN_DIAGNOSTIC_PLY"] else {
+            throw XCTSkip("Set UNTOLD_GAUSSIAN_DIAGNOSTIC_PLY to audit a production Gaussian asset")
+        }
+
+        let entity = createEntity()
+        let source = (path as NSString).deletingPathExtension
+        let fileExtension = (path as NSString).pathExtension
+        setEntityGaussian(entityId: entity, filename: source, withExtension: fileExtension)
+
+        let component = try XCTUnwrap(scene.get(component: GaussianComponent.self, for: entity))
+        let metadata = try XCTUnwrap(component.sphericalHarmonicsMetadata)
+        let coefficientBuffer = try XCTUnwrap(component.sphericalHarmonicsData)
+        let splatBuffer = try XCTUnwrap(component.encodedSplatData)
+        let splatCount = Int(component.splatCount)
+
+        XCTAssertEqual(splatCount, 1_732_378)
+        XCTAssertEqual(metadata.degree, 3)
+        XCTAssertEqual(metadata.coefficientsPerChannel, 16)
+        XCTAssertEqual(metadata.higherOrderCoefficientsPerSplat, 45)
+        XCTAssertEqual(coefficientBuffer.length, splatCount * 45 * MemoryLayout<Float16>.stride)
+        XCTAssertEqual(splatBuffer.length, splatCount * MemoryLayout<EncodedGaussianSplat>.stride)
+
+        let coefficients = coefficientBuffer.contents().bindMemory(
+            to: Float16.self,
+            capacity: splatCount * 45
+        )
+        let sampledValues = [0, splatCount * 45 / 2, splatCount * 45 - 1].map { coefficients[$0] }
+        XCTAssertTrue(sampledValues.contains { $0 != 0 })
     }
 
     func testGaussianExecution_HandlesEmptyScene() {
diff --git a/Tests/UntoldEngineTests/GaussianSphericalHarmonicsTests.swift b/Tests/UntoldEngineTests/GaussianSphericalHarmonicsTests.swift
new file mode 100644
index 00000000..504021b6
--- /dev/null
+++ b/Tests/UntoldEngineTests/GaussianSphericalHarmonicsTests.swift
@@ -0,0 +1,108 @@
+//
+//  GaussianSphericalHarmonicsTests.swift
+//  UntoldEngine
+//
+// Copyright (C) Untold Engine Studios
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at https://mozilla.org/MPL/2.0/.
+
+import simd
+@testable import UntoldEngine
+import XCTest
+
+final class GaussianSphericalHarmonicsTests: XCTestCase {
+    func testSRGBGaussianColorIsDecodedIntoLinearWorkingSpace() {
+        let result = gaussianSRGBToLinear(simd_float3(0, 0.04045, 0.5))
+
+        XCTAssertEqual(result.x, 0, accuracy: 1e-7)
+        XCTAssertEqual(result.y, 0.0031308, accuracy: 1e-7)
+        XCTAssertEqual(result.z, 0.21404114, accuracy: 1e-7)
+        XCTAssertEqual(gaussianSRGBToLinear(simd_float3(repeating: 1)), simd_float3(repeating: 1))
+    }
+
+    func testSRGBGaussianColorClampsNegativeSHResultsBeforeDecode() {
+        XCTAssertEqual(
+            gaussianSRGBToLinear(simd_float3(-1, -0.1, -.ulpOfOne)),
+            .zero
+        )
+    }
+
+    func testDegreeZeroReturnsBaseColorExactly() {
+        let base = simd_float3(0.2, 0.4, 0.6)
+        XCTAssertEqual(
+            evaluateGaussianSphericalHarmonics(
+                baseColor: base,
+                higherOrderCoefficients: [],
+                degree: 0,
+                direction: simd_float3(1, 0, 0)
+            ),
+            base
+        )
+    }
+
+    func testDegreeOneUsesStandardRealSHBasis() {
+        var coefficients = [Float](repeating: 0, count: 9)
+        coefficients[0] = 1 // Red coefficient 1: -C1 * y
+        let result = evaluateGaussianSphericalHarmonics(
+            baseColor: simd_float3(repeating: 1),
+            higherOrderCoefficients: coefficients,
+            degree: 1,
+            direction: simd_float3(0, 2, 0)
+        )
+
+        XCTAssertEqual(result.x, 1 - 0.4886025119, accuracy: 1e-6)
+        XCTAssertEqual(result.y, 1, accuracy: 1e-6)
+        XCTAssertEqual(result.z, 1, accuracy: 1e-6)
+    }
+
+    func testDegreeTwoUsesChannelMajorCoefficientSix() {
+        var coefficients = [Float](repeating: 0, count: 24)
+        coefficients[2 * 8 + 5] = 1 // Blue coefficient 6
+        let result = evaluateGaussianSphericalHarmonics(
+            baseColor: simd_float3(repeating: 0.5),
+            higherOrderCoefficients: coefficients,
+            degree: 2,
+            direction: simd_float3(0, 0, 1)
+        )
+
+        XCTAssertEqual(result.x, 0.5, accuracy: 1e-6)
+        XCTAssertEqual(result.y, 0.5, accuracy: 1e-6)
+        XCTAssertEqual(result.z, 0.5 + 2 * 0.3153915653, accuracy: 1e-6)
+    }
+
+    func testDegreeThreeUsesChannelMajorCoefficientTwelve() {
+        var coefficients = [Float](repeating: 0, count: 45)
+        coefficients[15 + 11] = 1 // Green coefficient 12
+        let result = evaluateGaussianSphericalHarmonics(
+            baseColor: simd_float3(repeating: 0.5),
+            higherOrderCoefficients: coefficients,
+            degree: 3,
+            direction: simd_float3(0, 0, 1)
+        )
+
+        XCTAssertEqual(result.x, 0.5, accuracy: 1e-6)
+        XCTAssertEqual(result.y, 0.5 + 2 * 0.3731763326, accuracy: 1e-6)
+        XCTAssertEqual(result.z, 0.5, accuracy: 1e-6)
+    }
+
+    func testLocalCameraPositionAccountsForEntityTransform() {
+        let rotation = simd_float4x4(simd_quatf(angle: .pi / 2, axis: simd_float3(0, 0, 1)))
+        let scale = simd_float4x4(diagonal: simd_float4(2, 3, 4, 1))
+        var translation = matrix_identity_float4x4
+        translation.columns.3 = simd_float4(10, -5, 7, 1)
+        let model = translation * rotation * scale
+        let expectedLocal = simd_float3(1, 2, 3)
+        let cameraWorld4 = model * simd_float4(expectedLocal, 1)
+
+        let actual = gaussianLocalCameraPosition(
+            cameraWorldPosition: simd_float3(cameraWorld4.x, cameraWorld4.y, cameraWorld4.z),
+            modelMatrix: model
+        )
+
+        XCTAssertEqual(actual.x, expectedLocal.x, accuracy: 1e-5)
+        XCTAssertEqual(actual.y, expectedLocal.y, accuracy: 1e-5)
+        XCTAssertEqual(actual.z, expectedLocal.z, accuracy: 1e-5)
+    }
+}
diff --git a/Tests/UntoldEngineTests/LoggerCategoryTests.swift b/Tests/UntoldEngineTests/LoggerCategoryTests.swift
index e87a7f61..b76d69e3 100644
--- a/Tests/UntoldEngineTests/LoggerCategoryTests.swift
+++ b/Tests/UntoldEngineTests/LoggerCategoryTests.swift
@@ -35,6 +35,7 @@ final class LoggerCategoryTests: XCTestCase {
         XCTAssertFalse(Logger.isEnabled(category: .textureStreaming))
         XCTAssertFalse(Logger.isEnabled(category: .textureLoading))
         XCTAssertFalse(Logger.isEnabled(category: .lightPortal))
+        XCTAssertFalse(Logger.isEnabled(category: .gaussian))
     }
 
     func testStreamingCategoriesCanBeEnabledIndividually() {
@@ -47,6 +48,17 @@ final class LoggerCategoryTests: XCTestCase {
         XCTAssertFalse(Logger.isEnabled(category: .textureLoading))
         XCTAssertFalse(Logger.isEnabled(category: .streamingHeartbeat))
         XCTAssertFalse(Logger.isEnabled(category: .lightPortal))
+        XCTAssertFalse(Logger.isEnabled(category: .gaussian))
+    }
+
+    func testGaussianCategoryCanBeEnabledIndividually() {
+        Logger.resetCategoryToggles()
+
+        setLogger(.category(.gaussian, true))
+
+        XCTAssertTrue(Logger.isEnabled(category: .gaussian))
+        XCTAssertFalse(Logger.isEnabled(category: .tileStreaming))
+        XCTAssertFalse(Logger.isEnabled(category: .batching))
     }
 
     func testWarningsRespectCategoryToggles() {
diff --git a/Tests/UntoldEngineTests/PLYReaderTest.swift b/Tests/UntoldEngineTests/PLYReaderTest.swift
index 507bb02d..801efaf7 100644
--- a/Tests/UntoldEngineTests/PLYReaderTest.swift
+++ b/Tests/UntoldEngineTests/PLYReaderTest.swift
@@ -200,6 +200,219 @@ final class PLYReaderTest: XCTestCase {
         XCTAssertEqual(quatMagnitude, 1.0, accuracy: 0.001)
     }
 
+    // MARK: - Spherical Harmonics
+
+    func test_readGaussianAsset_preservesDegreeZeroThroughThreeSphericalHarmonics() throws {
+        for degree in 0 ... 3 {
+            let coefficientsPerChannel = (degree + 1) * (degree + 1)
+            let restCount = (coefficientsPerChannel - 1) * 3
+            let restValues: [Float] = (0 ..< restCount).map { Float($0) }
+            let restProperties = (0 ..< restCount).map { "property float f_rest_\($0)" }
+
+            var lines: [String] = [
+                "ply",
+                "format ascii 1.0",
+                "element vertex 1",
+                "property float x",
+                "property float y",
+                "property float z",
+                "property float f_dc_0",
+                "property float f_dc_1",
+                "property float f_dc_2",
+            ]
+            lines.append(contentsOf: restProperties)
+            lines.append("end_header")
+            let valueStrings = ["1", "2", "3", "100", "200", "300"] + restValues.map { String($0) }
+            lines.append(valueStrings.joined(separator: " "))
+            let plyContent = lines.joined(separator: "\n")
+
+            let tempURL = createTempFile(content: plyContent)
+            tempFileURL = tempURL
+            let asset = try PLYReader.readGaussianAsset(from: tempURL)
+            let sh = try XCTUnwrap(asset.sphericalHarmonics)
+
+            XCTAssertEqual(sh.degree, degree)
+            XCTAssertEqual(sh.coefficientsPerChannel, coefficientsPerChannel)
+            XCTAssertEqual(sh.coefficientsPerSplat, coefficientsPerChannel * 3)
+
+            var expected: [Float] = []
+            let restPerChannel = coefficientsPerChannel - 1
+            for channel in 0 ..< 3 {
+                expected.append(Float((channel + 1) * 100))
+                let start = channel * restPerChannel
+                expected.append(contentsOf: restValues[start ..< start + restPerChannel])
+            }
+            XCTAssertEqual(sh.coefficients, expected)
+        }
+    }
+
+    func test_externalGaussianDiagnosticAssetPreservesAndPacksDegreeThreeSH() throws {
+        guard let path = ProcessInfo.processInfo.environment["UNTOLD_GAUSSIAN_DIAGNOSTIC_PLY"] else {
+            throw XCTSkip("Set UNTOLD_GAUSSIAN_DIAGNOSTIC_PLY to audit a production Gaussian asset")
+        }
+
+        let asset = try PLYReader.readGaussianAsset(from: URL(fileURLWithPath: path))
+        let sphericalHarmonics = try XCTUnwrap(asset.sphericalHarmonics)
+        XCTAssertEqual(sphericalHarmonics.degree, 3)
+        XCTAssertEqual(sphericalHarmonics.coefficientsPerChannel, 16)
+        XCTAssertEqual(sphericalHarmonics.coefficients.count, asset.splats.count * 48)
+
+        let packed = try packGaussianSphericalHarmonics(
+            sphericalHarmonics,
+            splatCount: asset.splats.count
+        )
+        XCTAssertEqual(packed.metadata.degree, 3)
+        XCTAssertEqual(packed.metadata.coefficientsPerChannel, 16)
+        XCTAssertEqual(packed.metadata.higherOrderCoefficientsPerSplat, 45)
+        XCTAssertEqual(packed.coefficients.count, asset.splats.count * 45)
+
+        let sampleIndices = [0, asset.splats.count / 2, asset.splats.count - 1]
+        var observedDirectionalChange = false
+        for splatIndex in sampleIndices {
+            let sourceBase = splatIndex * 48
+            let packedBase = splatIndex * 45
+            let higherOrder = packed.coefficients[packedBase ..< packedBase + 45].map(Float.init)
+            for channel in 0 ..< 3 {
+                for coefficient in 0 ..< 15 {
+                    XCTAssertEqual(
+                        higherOrder[channel * 15 + coefficient],
+                        Float(Float16(sphericalHarmonics.coefficients[sourceBase + channel * 16 + coefficient + 1]))
+                    )
+                }
+            }
+
+            let splat = asset.splats[splatIndex]
+            let baseColor = simd_float3(splat.color.x, splat.color.y, splat.color.z)
+            let evaluated = evaluateGaussianSphericalHarmonics(
+                baseColor: baseColor,
+                higherOrderCoefficients: higherOrder,
+                degree: 3,
+                direction: simd_float3(0.25, -0.5, 0.75)
+            )
+            observedDirectionalChange = observedDirectionalChange || simd_distance(evaluated, baseColor) > 1e-5
+        }
+        XCTAssertTrue(observedDirectionalChange, "Production SH data must affect evaluated color")
+
+        print("Gaussian diagnostic: splats=\(asset.splats.count), degree=\(sphericalHarmonics.degree), "
+            + "sourceCoefficients=\(sphericalHarmonics.coefficients.count), "
+            + "packedCoefficients=\(packed.coefficients.count), packedBytes=\(packed.coefficients.count * MemoryLayout<Float16>.stride)")
+    }
+
+    func test_readGaussianAsset_binaryAndASCIIHaveIdenticalSphericalHarmonics() throws {
+        let propertyLines = [
+            "property float x", "property float y", "property float z",
+            "property float f_dc_0", "property float f_dc_1", "property float f_dc_2",
+        ] + (0 ..< 9).map { "property float f_rest_\($0)" }
+        let values: [Float] = (1 ... 15).map { Float($0) }
+
+        var asciiLines = ["ply", "format ascii 1.0", "element vertex 1"]
+        asciiLines.append(contentsOf: propertyLines)
+        asciiLines.append("end_header")
+        asciiLines.append(values.map { String($0) }.joined(separator: " "))
+        let ascii = asciiLines.joined(separator: "\n")
+        let asciiURL = createTempFile(content: ascii)
+
+        var binaryHeaderLines = ["ply", "format binary_little_endian 1.0", "element vertex 1"]
+        binaryHeaderLines.append(contentsOf: propertyLines)
+        binaryHeaderLines.append("end_header")
+        let binaryHeader = binaryHeaderLines.joined(separator: "\n") + "\n"
+        let binaryURL = createTempBinaryFile(header: binaryHeader, values: values)
+        defer {
+            try? FileManager.default.removeItem(at: asciiURL)
+            try? FileManager.default.removeItem(at: binaryURL)
+        }
+
+        let asciiAsset = try PLYReader.readGaussianAsset(from: asciiURL)
+        let binaryAsset = try PLYReader.readGaussianAsset(from: binaryURL)
+        XCTAssertEqual(asciiAsset.sphericalHarmonics?.degree, 1)
+        XCTAssertEqual(
+            asciiAsset.sphericalHarmonics?.coefficients,
+            binaryAsset.sphericalHarmonics?.coefficients
+        )
+    }
+
+    func test_readGaussianAsset_rejectsMalformedSphericalHarmonicSchemas() {
+        let malformedPropertySets = [
+            ["property float f_dc_0", "property float f_dc_1"],
+            ["property float f_rest_0"],
+            ["property float f_dc_0", "property float f_dc_1", "property float f_dc_2", "property float f_rest_0", "property float f_rest_2"],
+            ["property float f_dc_0", "property float f_dc_1", "property float f_dc_2", "property float f_rest_bad"],
+        ]
+
+        for properties in malformedPropertySets {
+            let values = Array(repeating: "0", count: 3 + properties.count).joined(separator: " ")
+            let plyContent = ([
+                "ply", "format ascii 1.0", "element vertex 1",
+                "property float x", "property float y", "property float z",
+            ] + properties + ["end_header", values]).joined(separator: "\n")
+            let tempURL = createTempFile(content: plyContent)
+            defer { try? FileManager.default.removeItem(at: tempURL) }
+
+            XCTAssertThrowsError(try PLYReader.readGaussianAsset(from: tempURL)) { error in
+                guard case PLYError.invalidData = error else {
+                    XCTFail("Expected invalidData, got \(error)")
+                    return
+                }
+            }
+        }
+    }
+
+    func test_packGaussianSphericalHarmonics_preservesHigherOrderChannelMajorLayout() throws {
+        XCTAssertEqual(MemoryLayout<Float16>.stride, 2)
+
+        let first = (0 ..< 12).map { Float($0) }
+        let second = (100 ..< 112).map { Float($0) }
+        let sphericalHarmonics = GaussianSphericalHarmonics(
+            degree: 1,
+            coefficientsPerChannel: 4,
+            coefficients: first + second
+        )
+
+        let packed = try packGaussianSphericalHarmonics(sphericalHarmonics, splatCount: 2)
+        let expected = [
+            1, 2, 3, 5, 6, 7, 9, 10, 11,
+            101, 102, 103, 105, 106, 107, 109, 110, 111,
+        ].map(Float16.init)
+
+        XCTAssertEqual(packed.coefficients, expected)
+        XCTAssertEqual(packed.metadata.degree, 1)
+        XCTAssertEqual(packed.metadata.coefficientsPerChannel, 4)
+        XCTAssertEqual(packed.metadata.higherOrderCoefficientsPerSplat, 9)
+    }
+
+    func test_packGaussianSphericalHarmonics_degreeZeroRequiresNoCoefficientBuffer() throws {
+        let sphericalHarmonics = GaussianSphericalHarmonics(
+            degree: 0,
+            coefficientsPerChannel: 1,
+            coefficients: [0.1, 0.2, 0.3]
+        )
+
+        let packed = try packGaussianSphericalHarmonics(sphericalHarmonics, splatCount: 1)
+
+        XCTAssertTrue(packed.coefficients.isEmpty)
+        XCTAssertEqual(packed.metadata.degree, 0)
+        XCTAssertEqual(packed.metadata.coefficientsPerChannel, 1)
+        XCTAssertEqual(packed.metadata.higherOrderCoefficientsPerSplat, 0)
+    }
+
+    func test_packGaussianSphericalHarmonics_rejectsInvalidCountAndFloat16Overflow() {
+        let invalidCount = GaussianSphericalHarmonics(
+            degree: 1,
+            coefficientsPerChannel: 4,
+            coefficients: [0]
+        )
+        XCTAssertThrowsError(try packGaussianSphericalHarmonics(invalidCount, splatCount: 1))
+
+        var overflowCoefficients = [Float](repeating: 0, count: 12)
+        overflowCoefficients[1] = Float.greatestFiniteMagnitude
+        let overflow = GaussianSphericalHarmonics(
+            degree: 1,
+            coefficientsPerChannel: 4,
+            coefficients: overflowCoefficients
+        )
+        XCTAssertThrowsError(try packGaussianSphericalHarmonics(overflow, splatCount: 1))
+    }
+
     // MARK: - Test Error Handling
 
     func test_invalidPLY_missingMagicNumber() {
@@ -346,4 +559,20 @@ final class PLYReaderTest: XCTestCase {
 
         return fileURL
     }
+
+    private func createTempBinaryFile(header: String, values: [Float]) -> URL {
+        let fileURL = FileManager.default.temporaryDirectory
+            .appendingPathComponent("test_\(UUID().uuidString).ply")
+        var data = Data(header.utf8)
+        for value in values {
+            var bits = value.bitPattern.littleEndian
+            withUnsafeBytes(of: &bits) { data.append(contentsOf: $0) }
+        }
+        do {
+            try data.write(to: fileURL)
+        } catch {
+            XCTFail("Failed to create binary temporary file: \(error)")
+        }
+        return fileURL
+    }
 }
diff --git a/Tools/UntoldEngineCLI/Sources/UntoldEngineCLI/ExportCommand.swift b/Tools/UntoldEngineCLI/Sources/UntoldEngineCLI/ExportCommand.swift
index 33af6648..77629442 100644
--- a/Tools/UntoldEngineCLI/Sources/UntoldEngineCLI/ExportCommand.swift
+++ b/Tools/UntoldEngineCLI/Sources/UntoldEngineCLI/ExportCommand.swift
@@ -219,7 +219,10 @@ enum ExportError: LocalizedError {
         case let .inputNotFound(path):
             return "Input asset does not exist: \(path)"
         case .blenderNotFound:
-            return "Blender was not found. Install Blender, use --blender, or set BLENDER_BIN."
+            return """
+            Blender was not found. Download it from https://www.blender.org/download/ \
+            (tested with Blender 5.1.0), then use --blender /path/to/Blender or set BLENDER_BIN.
+            """
         case let .blenderNotExecutable(path):
             return "Blender is not executable at: \(path)"
         case let .exporterNotInstalled(path):