From 880517d91a29e3cc6d71a6b1787f126ac939d6a1 Mon Sep 17 00:00:00 2001 From: snakamura <88454117+snkmcb@users.noreply.github.com> Date: Sun, 27 Sep 2026 21:21:26 +0900 Subject: [PATCH] Share Gaussian preparation and radix sort with Metal compute --- CHANGELOG.md | 12 + .../shaders/gaussian-prepare-metal.slang | 27 ++ .../shaders/gaussian-sort-metal.slang | 37 ++ .../merlin-metal/src/gaussian_compute_abi.hpp | 91 +++++ .../shaders/gaussian-prepare.slang | 369 +---------------- .../merlin-vulkan/shaders/gaussian-sort.slang | 320 +-------------- cmake/MerlinMetalShaders.cmake | 33 +- .../shaders/gaussian-prepare-abi.slang | 52 +++ .../shaders/gaussian-prepare-common.slang | 304 ++++++++++++++ .../shaders/gaussian-prepared-record.slang | 18 + .../shaders/gaussian-sort-abi.slang | 58 +++ .../shaders/gaussian-sort-common.slang | 257 ++++++++++++ docs/design/metal-gaussian-execution.md | 27 +- docs/reference/support-matrix.md | 2 +- docs/roadmap/current.md | 6 +- tests/CMakeLists.txt | 16 +- tests/compile-gaussian-spirv.cmake | 20 + tests/metal_gaussian_compute_test.mm | 378 ++++++++++++++++++ tests/run_install_consumer.cmake | 7 + tests/verify-metal-gaussian-abi.cmake | 81 ++++ 20 files changed, 1420 insertions(+), 695 deletions(-) create mode 100644 backend/merlin-metal/shaders/gaussian-prepare-metal.slang create mode 100644 backend/merlin-metal/shaders/gaussian-sort-metal.slang create mode 100644 backend/merlin-metal/src/gaussian_compute_abi.hpp create mode 100644 core/merlin-render-backend/shaders/gaussian-prepare-abi.slang create mode 100644 core/merlin-render-backend/shaders/gaussian-prepare-common.slang create mode 100644 core/merlin-render-backend/shaders/gaussian-prepared-record.slang create mode 100644 core/merlin-render-backend/shaders/gaussian-sort-abi.slang create mode 100644 core/merlin-render-backend/shaders/gaussian-sort-common.slang create mode 100644 tests/metal_gaussian_compute_test.mm diff --git a/CHANGELOG.md b/CHANGELOG.md index 947a44d..af9a57a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,13 @@ after its public API and release process are established. ### Added +- Shared Slang Gaussian projection, covariance, SH, culling/compaction and + deterministic radix-sort kernels now compile into the Metal Gaussian library. + Native Metal buffer bindings and a scalar 64-byte prepared-record ABI have + reflection, install-package and Apple GPU comparisons against the CPU + reference. These kernels are not yet connected to Metal renderer scheduling; + interactive rendering still uses CPU preparation and sorting. + - Metal Gaussian reference rasterization now uses shared Slang ellipse/alpha math with Vulkan. Slang 2026.8.x and Xcode compile an embedded metallib during the build; installed consumers do not compile Gaussian shaders at runtime. @@ -283,6 +290,11 @@ after its public API and release process are established. ### Fixed +- Shared GPU Gaussian preparation now applies the inverse local-to-camera + matrix when evaluating directional SH. The previous inverse-transpose + calculation changed radiance under rotation and nonuniform scale; a Metal + compute regression compares transformed degree-three SH with the CPU path. + - Zooming into a Gaussian scene no longer floods the view with a single color in usdview or the development viewport. CPU and GPU preparation cull a kernel whose center lies in front of the camera's near plane, as diff --git a/backend/merlin-metal/shaders/gaussian-prepare-metal.slang b/backend/merlin-metal/shaders/gaussian-prepare-metal.slang new file mode 100644 index 0000000..cdfdb0a --- /dev/null +++ b/backend/merlin-metal/shaders/gaussian-prepare-metal.slang @@ -0,0 +1,27 @@ +#include "../../../core/merlin-render-backend/shaders/gaussian-prepare-abi.slang" + +// Metal uses one explicit buffer namespace for constants and storage. +[[vk::binding(7, 0)]] +ConstantBuffer gaussian_prepare_constants : register(b7); + +// Byte-addressed inputs match the existing tightly packed arena payloads: +// float3 positions, six-float covariance, float opacity, and float3 SH terms. +[[vk::binding(0, 0)]] ByteAddressBuffer gaussian_positions : register(t0); +[[vk::binding(1, 0)]] ByteAddressBuffer gaussian_covariances : register(t1); +[[vk::binding(2, 0)]] ByteAddressBuffer gaussian_opacities : register(t2); +[[vk::binding(3, 0)]] ByteAddressBuffer gaussian_radiance : register(t3); +[[vk::binding(4, 0)]] RWStructuredBuffer gaussian_candidate_results : register(u4); +[[vk::binding(5, 0)]] RWByteAddressBuffer gaussian_prepared_records : register(u5); +[[vk::binding(6, 0)]] RWStructuredBuffer gaussian_prepare_counters : register(u6); + +// Raw records keep the Vulkan 64-byte ABI despite Metal float3 alignment. +void StorePreparedRecord(uint index, GaussianPreparedRecord record) +{ + uint offset = index * 64u; + gaussian_prepared_records.Store4(offset, asuint(float4(record.center_pixels, record.radius_pixels, record.depth))); + gaussian_prepared_records.Store4(offset + 16u, asuint(float4(record.inverse_conic, record.opacity))); + gaussian_prepared_records.Store4(offset + 32u, asuint(float4(record.radiance, record.sort_key))); + gaussian_prepared_records.Store4(offset + 48u, uint4(record.resource_id_low, record.resource_id_high, record.particle_id, record.padding)); +} + +#include "../../../core/merlin-render-backend/shaders/gaussian-prepare-common.slang" diff --git a/backend/merlin-metal/shaders/gaussian-sort-metal.slang b/backend/merlin-metal/shaders/gaussian-sort-metal.slang new file mode 100644 index 0000000..28f89c8 --- /dev/null +++ b/backend/merlin-metal/shaders/gaussian-sort-metal.slang @@ -0,0 +1,37 @@ +#include "../../../core/merlin-render-backend/shaders/gaussian-sort-abi.slang" + +[[vk::binding(4, 0)]] +ConstantBuffer gaussian_sort_constants : register(b4); + +[[vk::binding(0, 0)]] RWStructuredBuffer gaussian_sort_source : register(u0); +[[vk::binding(1, 0)]] RWStructuredBuffer gaussian_sort_destination : register(u1); +// Verification words, per-resource visible counts, then the digit-major +// histogram and its scan levels. +[[vk::binding(2, 0)]] RWStructuredBuffer gaussian_sort_scan : register(u2); +[[vk::binding(3, 0)]] ByteAddressBuffer gaussian_sort_prepared_records : register(t3); + +// Match the raw 64-byte stream written by Metal preparation. +[ForceInline] +GaussianPreparedRecord LoadSortPreparedRecord(uint index) +{ + uint offset = index * 64u; + float4 geometry = asfloat(gaussian_sort_prepared_records.Load4(offset)); + float4 conic = asfloat(gaussian_sort_prepared_records.Load4(offset + 16u)); + float4 radiance = asfloat(gaussian_sort_prepared_records.Load4(offset + 32u)); + uint4 identity = gaussian_sort_prepared_records.Load4(offset + 48u); + GaussianPreparedRecord record; + record.center_pixels = geometry.xy; + record.radius_pixels = geometry.z; + record.depth = geometry.w; + record.inverse_conic = conic.xyz; + record.opacity = conic.w; + record.radiance = radiance.xyz; + record.sort_key = radiance.w; + record.resource_id_low = identity.x; + record.resource_id_high = identity.y; + record.particle_id = identity.z; + record.padding = identity.w; + return record; +} + +#include "../../../core/merlin-render-backend/shaders/gaussian-sort-common.slang" diff --git a/backend/merlin-metal/src/gaussian_compute_abi.hpp b/backend/merlin-metal/src/gaussian_compute_abi.hpp new file mode 100644 index 0000000..00a51f1 --- /dev/null +++ b/backend/merlin-metal/src/gaussian_compute_abi.hpp @@ -0,0 +1,91 @@ +#pragma once + +#include +#include + +#include + +// Private host records for the shared Gaussian compute kernels. Metal accesses +// prepared records through byte-addressed buffers to preserve the scalar ABI. +namespace merlin::metal::gaussian_compute { + +struct alignas(16) PrepareConstants { + Mat4 local_to_camera; + Mat4 projection; + Vec2 viewport_size; + float sigma_extent{3.0F}; + float minimum_variance_pixels{0.25F}; + std::uint32_t resource_id_low{}; + std::uint32_t resource_id_high{}; + std::uint32_t particle_count{}; + std::uint32_t coefficients_per_particle{}; + std::uint32_t spherical_harmonics_degree{}; + std::uint32_t projection_mode{}; + std::uint32_t sorting_mode{}; + std::uint32_t padding{}; +}; + +struct alignas(16) PreparedRecord { + Vec2 center_pixels; + float radius_pixels{}; + float depth{}; + Vec3 inverse_conic; + float opacity{}; + Vec3 radiance; + float sort_key{}; + std::uint32_t resource_id_low{}; + std::uint32_t resource_id_high{}; + std::uint32_t particle_id{}; + std::uint32_t padding{}; +}; + +struct PrepareCounters { + std::uint32_t candidate_count{}; + std::uint32_t visible_count{}; + std::uint32_t opacity_culled_count{}; + std::uint32_t frustum_culled_count{}; + std::uint32_t invalid_culled_count{}; + std::uint32_t padding[3]{}; +}; + +struct SortElement { + std::uint32_t key_low{}; + std::uint32_t key_high{}; + std::uint32_t value{}; +}; + +struct SortConstants { + std::uint32_t element_count{}; + std::uint32_t block_count{}; + std::uint32_t digit_shift{}; + std::uint32_t digit_word{}; + std::uint32_t scan_offset{}; + std::uint32_t scan_count{}; + std::uint32_t scan_sums_offset{}; + std::uint32_t candidate_base{}; + std::uint32_t prepared_base{}; + std::uint32_t visible_count_offset{}; + std::uint32_t count_word{}; + std::uint32_t flags{}; +}; + +static_assert(sizeof(PrepareConstants) == 176); +static_assert(offsetof(PrepareConstants, projection) == 64); +static_assert(offsetof(PrepareConstants, viewport_size) == 128); +static_assert(offsetof(PrepareConstants, particle_count) == 152); +static_assert(offsetof(PrepareConstants, sorting_mode) == 168); +static_assert(sizeof(PreparedRecord) == 64); +static_assert(offsetof(PreparedRecord, inverse_conic) == 16); +static_assert(offsetof(PreparedRecord, opacity) == 28); +static_assert(offsetof(PreparedRecord, radiance) == 32); +static_assert(offsetof(PreparedRecord, sort_key) == 44); +static_assert(offsetof(PreparedRecord, resource_id_low) == 48); +static_assert(offsetof(PreparedRecord, particle_id) == 56); +static_assert(sizeof(PrepareCounters) == 32); +static_assert(sizeof(SortElement) == 12); +static_assert(sizeof(SortConstants) == 48); +static_assert(offsetof(SortConstants, scan_offset) == 16); +static_assert(offsetof(SortConstants, visible_count_offset) == 36); +static_assert(offsetof(SortConstants, flags) == 44); + +} // namespace merlin::metal::gaussian_compute diff --git a/backend/merlin-vulkan/shaders/gaussian-prepare.slang b/backend/merlin-vulkan/shaders/gaussian-prepare.slang index 2006faf..9e8f5b8 100644 --- a/backend/merlin-vulkan/shaders/gaussian-prepare.slang +++ b/backend/merlin-vulkan/shaders/gaussian-prepare.slang @@ -1,68 +1,4 @@ -// First GPU-driven Gaussian preparation contract. One dispatch processes one -// resident Gaussian resource directly from its four device-local attribute -// ranges, classifies every particle, and atomically compacts visible records. -// A later radix-sort/tile slice consumes the stable sort key and identity. - -static const uint kGaussianProjectionPerspective = 0u; -static const uint kGaussianProjectionTangential = 1u; -static const uint kGaussianSortingZDepth = 0u; -static const uint kGaussianSortingCameraDistance = 1u; - -static const uint kGaussianCandidateVisible = 0u; -static const uint kGaussianCandidateOpacityCulled = 1u; -static const uint kGaussianCandidateFrustumCulled = 2u; -static const uint kGaussianCandidateInvalidCulled = 3u; - -static const float kProjectionEpsilon = 1.0e-6f; -// Matches the CPU reference: the 3DGS guard band bounds the perspective -// Jacobian of close off-axis kernels. -static const float kJacobianGuardBand = 1.3f; - -struct GaussianPrepareConstants -{ - float4x4 local_to_camera; - float4x4 projection; - float2 viewport_size; - float sigma_extent; - float minimum_variance_pixels; - uint resource_id_low; - uint resource_id_high; - uint particle_count; - uint coefficients_per_particle; - uint spherical_harmonics_degree; - uint projection_mode; - uint sorting_mode; - uint padding; -}; - -// Kept at 64 bytes so it can become the direct input of sorted/tiled raster -// without repacking. Resource + particle identity is a deterministic tie break. -struct GaussianPreparedRecord -{ - float2 center_pixels; - float radius_pixels; - float depth; - float3 inverse_conic; - float opacity; - float3 radiance; - float sort_key; - uint resource_id_low; - uint resource_id_high; - uint particle_id; - uint padding; -}; - -struct GaussianPrepareDispatchCounters -{ - uint candidate_count; - uint visible_count; - uint opacity_culled_count; - uint frustum_culled_count; - uint invalid_culled_count; - uint padding0; - uint padding1; - uint padding2; -}; +#include "../../../core/merlin-render-backend/shaders/gaussian-prepare-abi.slang" // This block intentionally uses a uniform descriptor: its two matrices make // it 176 bytes, larger than Vulkan's guaranteed 128-byte push-constant limit. @@ -79,306 +15,9 @@ ConstantBuffer gaussian_prepare_constants; [[vk::binding(5, 3)]] RWStructuredBuffer gaussian_prepared_records; [[vk::binding(6, 3)]] RWStructuredBuffer gaussian_prepare_counters; -bool IsFinite1(float value) -{ - return !isnan(value) && !isinf(value); -} - -bool IsFinite2(float2 value) -{ - return all(!isnan(value)) && all(!isinf(value)); -} - -bool IsFinite3(float3 value) +void StorePreparedRecord(uint index, GaussianPreparedRecord record) { - return all(!isnan(value)) && all(!isinf(value)); + gaussian_prepared_records[index] = record; } -float3 SafeNormalize(float3 value) -{ - float length_squared = dot(value, value); - if (!IsFinite1(length_squared) || length_squared <= kProjectionEpsilon * kProjectionEpsilon) - { - return float3(0.0f); - } - return value * rsqrt(length_squared); -} - -float3 LoadPosition(uint particle) -{ - return gaussian_positions.Load(particle * 12u); -} - -float3x3 LoadCovariance(uint particle) -{ - uint offset = particle * 24u; - float3 first = gaussian_covariances.Load(offset); - float3 second = gaussian_covariances.Load(offset + 12u); - return float3x3( - first.x, first.y, first.z, - first.y, second.x, second.y, - first.z, second.y, second.z); -} - -float3 LoadCoefficient(uint index) -{ - return gaussian_radiance.Load(index * 12u); -} - -void AddScaled(inout float3 result, uint coefficient, float basis) -{ - result += LoadCoefficient(coefficient) * basis; -} - -float3 EvaluateRadiance(uint particle, float3 direction) -{ - static const float c0 = 0.28209479177387814f; - static const float c1 = 0.4886025119029199f; - uint first = particle * gaussian_prepare_constants.coefficients_per_particle; - float3 result = float3(0.5f) + LoadCoefficient(first) * c0; - uint degree = min(gaussian_prepare_constants.spherical_harmonics_degree, 3u); - if (degree >= 1u) - { - AddScaled(result, first + 1u, -c1 * direction.y); - AddScaled(result, first + 2u, c1 * direction.z); - AddScaled(result, first + 3u, -c1 * direction.x); - } - float xx = direction.x * direction.x; - float yy = direction.y * direction.y; - float zz = direction.z * direction.z; - if (degree >= 2u) - { - AddScaled(result, first + 4u, 1.0925484305920792f * direction.x * direction.y); - AddScaled(result, first + 5u, 1.0925484305920792f * direction.y * direction.z); - AddScaled(result, first + 6u, 0.31539156525252005f * (2.0f * zz - xx - yy)); - AddScaled(result, first + 7u, 1.0925484305920792f * direction.x * direction.z); - AddScaled(result, first + 8u, 0.5462742152960396f * (xx - yy)); - } - if (degree >= 3u) - { - AddScaled(result, first + 9u, -0.5900435899266435f * direction.y * (3.0f * xx - yy)); - AddScaled(result, first + 10u, 2.890611442640554f * direction.x * direction.y * direction.z); - AddScaled(result, first + 11u, -0.4570457994644658f * direction.y * (4.0f * zz - xx - yy)); - AddScaled(result, first + 12u, 0.3731763325901154f * direction.z * (2.0f * zz - 3.0f * xx - 3.0f * yy)); - AddScaled(result, first + 13u, -0.4570457994644658f * direction.x * (4.0f * zz - xx - yy)); - AddScaled(result, first + 14u, 1.445305721320277f * direction.z * (xx - yy)); - AddScaled(result, first + 15u, -0.5900435899266435f * direction.x * (xx - 3.0f * yy)); - } - return max(result, float3(0.0f)); -} - -struct GaussianProjectionJacobian -{ - float3 dx; - float3 dy; -}; - -GaussianProjectionJacobian ProjectionJacobian(float3 camera, float4 clip) -{ - float3 row_x = gaussian_prepare_constants.projection[0].xyz; - float3 row_y = gaussian_prepare_constants.projection[1].xyz; - float3 row_w = gaussian_prepare_constants.projection[3].xyz; - float inverse_w_squared = 1.0f / (clip.w * clip.w); - float guard = kJacobianGuardBand * clip.w; - float3 dx = (row_x * clip.w - clamp(clip.x, -guard, guard) * row_w) * inverse_w_squared; - float3 dy = (row_y * clip.w - clamp(clip.y, -guard, guard) * row_w) * inverse_w_squared; - - if (gaussian_prepare_constants.projection_mode == kGaussianProjectionTangential && - dot(row_w, row_w) > kProjectionEpsilon * kProjectionEpsilon) - { - float distance = length(camera); - float3 direction = SafeNormalize(camera); - float3 tangent_x = SafeNormalize(float3( - 1.0f - direction.x * direction.x, - -direction.x * direction.y, - -direction.x * direction.z)); - if (dot(tangent_x, tangent_x) <= kProjectionEpsilon * kProjectionEpsilon) - { - tangent_x = SafeNormalize(float3( - -direction.y * direction.x, - 1.0f - direction.y * direction.y, - -direction.y * direction.z)); - } - float3 tangent_y = SafeNormalize(cross(tangent_x, direction)); - if (tangent_y.y < 0.0f) - { - tangent_y = -tangent_y; - } - dx = tangent_x * (max(abs(gaussian_prepare_constants.projection[0][0]), kProjectionEpsilon) / distance); - dy = tangent_y * (max(abs(gaussian_prepare_constants.projection[1][1]), kProjectionEpsilon) / distance); - } - GaussianProjectionJacobian result; - result.dx = dx; - result.dy = dy; - return result; -} - -bool ProjectCovariance(float3x3 covariance, GaussianProjectionJacobian jacobian, - out float3 inverse_conic, out float radius_pixels) -{ - inverse_conic = float3(0.0f); - radius_pixels = 0.0f; - float3 dx = jacobian.dx; - float3 dy = jacobian.dy; - float width_scale = gaussian_prepare_constants.viewport_size.x * 0.5f; - float height_scale = gaussian_prepare_constants.viewport_size.y * 0.5f; - float xx = dot(dx, mul(covariance, dx)) * width_scale * width_scale + - gaussian_prepare_constants.minimum_variance_pixels; - float xy = dot(dx, mul(covariance, dy)) * width_scale * height_scale; - float yy = dot(dy, mul(covariance, dy)) * height_scale * height_scale + - gaussian_prepare_constants.minimum_variance_pixels; - float determinant = xx * yy - xy * xy; - if (!IsFinite1(xx) || !IsFinite1(xy) || !IsFinite1(yy) || - !IsFinite1(determinant) || determinant <= kProjectionEpsilon) - { - return false; - } - float discriminant = sqrt(max(0.0f, (xx - yy) * (xx - yy) + 4.0f * xy * xy)); - float largest_eigenvalue = 0.5f * (xx + yy + discriminant); - radius_pixels = gaussian_prepare_constants.sigma_extent * sqrt(max(largest_eigenvalue, 0.0f)); - inverse_conic = float3(yy / determinant, -xy / determinant, xx / determinant); - return IsFinite1(radius_pixels) && radius_pixels > 0.0f && IsFinite3(inverse_conic); -} - -bool ProjectDepthExtent(float3x3 covariance, float3 camera, float4 clip, - out float radius) -{ - radius = 0.0f; - float3 row_z = gaussian_prepare_constants.projection[2].xyz; - float3 row_w = gaussian_prepare_constants.projection[3].xyz; - float inverse_w_squared = 1.0f / (clip.w * clip.w); - float3 dz = (row_z * clip.w - clip.z * row_w) * inverse_w_squared; - float variance = dot(dz, mul(covariance, dz)); - if (!IsFinite1(variance) || variance < -kProjectionEpsilon) - { - return false; - } - radius = gaussian_prepare_constants.sigma_extent * sqrt(max(variance, 0.0f)); - return IsFinite1(radius); -} - -float3 LocalCameraDirection(float3x3 local_to_camera, float3 camera) -{ - float3 row0 = local_to_camera[0]; - float3 row1 = local_to_camera[1]; - float3 row2 = local_to_camera[2]; - float3 cofactor0 = cross(row1, row2); - float3 cofactor1 = cross(row2, row0); - float3 cofactor2 = cross(row0, row1); - float determinant_value = dot(row0, cofactor0); - if (!IsFinite1(determinant_value) || abs(determinant_value) <= kProjectionEpsilon) - { - return SafeNormalize(camera); - } - float inverse_determinant = 1.0f / determinant_value; - return SafeNormalize(float3( - dot(cofactor0, camera), - dot(cofactor1, camera), - dot(cofactor2, camera)) * inverse_determinant); -} - -void RecordRejection(uint particle, uint result, inout uint counter) -{ - gaussian_candidate_results[particle] = result; - uint ignored; - InterlockedAdd(counter, 1u, ignored); -} - -[shader("compute")] -[numthreads(64, 1, 1)] -void gaussian_prepare_compact(uint3 dispatch_thread_id : SV_DispatchThreadID) -{ - uint particle = dispatch_thread_id.x; - if (particle >= gaussian_prepare_constants.particle_count) - { - return; - } - if (particle == 0u) - { - gaussian_prepare_counters[0].candidate_count = gaussian_prepare_constants.particle_count; - } - - float opacity = gaussian_opacities.Load(particle * 4u); - if (!IsFinite1(opacity) || opacity <= 0.0f) - { - RecordRejection(particle, kGaussianCandidateOpacityCulled, - gaussian_prepare_counters[0].opacity_culled_count); - return; - } - - float3 local_position = LoadPosition(particle); - float4 camera4 = mul(gaussian_prepare_constants.local_to_camera, float4(local_position, 1.0f)); - float3 camera = camera4.xyz; - float4 clip = mul(gaussian_prepare_constants.projection, float4(camera, 1.0f)); - if (!IsFinite3(camera) || any(isnan(clip)) || any(isinf(clip)) || clip.w <= kProjectionEpsilon) - { - RecordRejection(particle, kGaussianCandidateFrustumCulled, - gaussian_prepare_counters[0].frustum_culled_count); - return; - } - float3 ndc = clip.xyz / clip.w; - // Kernels centered in front of the near plane are clipped like Mesh - // geometry; their center-evaluated footprint would cover the view. - if (ndc.z < 0.0f) - { - RecordRejection(particle, kGaussianCandidateFrustumCulled, - gaussian_prepare_counters[0].frustum_culled_count); - return; - } - float3x3 local_linear = (float3x3)gaussian_prepare_constants.local_to_camera; - float3x3 covariance = LoadCovariance(particle); - float3x3 camera_covariance = mul(mul(local_linear, covariance), transpose(local_linear)); - float depth_radius; - if (!ProjectDepthExtent(camera_covariance, camera, clip, depth_radius)) - { - RecordRejection(particle, kGaussianCandidateInvalidCulled, - gaussian_prepare_counters[0].invalid_culled_count); - return; - } - if (ndc.z - depth_radius > 1.0f) - { - RecordRejection(particle, kGaussianCandidateFrustumCulled, - gaussian_prepare_counters[0].frustum_culled_count); - return; - } - float3 inverse_conic; - float radius_pixels; - if (!ProjectCovariance(camera_covariance, ProjectionJacobian(camera, clip), - inverse_conic, radius_pixels)) - { - RecordRejection(particle, kGaussianCandidateInvalidCulled, - gaussian_prepare_counters[0].invalid_culled_count); - return; - } - - float2 center_pixels = (ndc.xy * 0.5f + 0.5f) * gaussian_prepare_constants.viewport_size; - if (!IsFinite2(center_pixels) || - center_pixels.x + radius_pixels < 0.0f || - center_pixels.x - radius_pixels > gaussian_prepare_constants.viewport_size.x || - center_pixels.y + radius_pixels < 0.0f || - center_pixels.y - radius_pixels > gaussian_prepare_constants.viewport_size.y) - { - RecordRejection(particle, kGaussianCandidateFrustumCulled, - gaussian_prepare_counters[0].frustum_culled_count); - return; - } - - uint compact_index; - InterlockedAdd(gaussian_prepare_counters[0].visible_count, 1u, compact_index); - GaussianPreparedRecord prepared; - prepared.center_pixels = center_pixels; - prepared.radius_pixels = radius_pixels; - prepared.depth = clamp(ndc.z, 0.0f, 1.0f); - prepared.inverse_conic = inverse_conic; - prepared.opacity = opacity; - prepared.radiance = EvaluateRadiance( - particle, LocalCameraDirection(local_linear, camera)); - prepared.sort_key = gaussian_prepare_constants.sorting_mode == kGaussianSortingCameraDistance - ? dot(camera, camera) : ndc.z; - prepared.resource_id_low = gaussian_prepare_constants.resource_id_low; - prepared.resource_id_high = gaussian_prepare_constants.resource_id_high; - prepared.particle_id = particle; - prepared.padding = 0u; - gaussian_prepared_records[compact_index] = prepared; - gaussian_candidate_results[particle] = kGaussianCandidateVisible; -} +#include "../../../core/merlin-render-backend/shaders/gaussian-prepare-common.slang" diff --git a/backend/merlin-vulkan/shaders/gaussian-sort.slang b/backend/merlin-vulkan/shaders/gaussian-sort.slang index a07a074..244471b 100644 --- a/backend/merlin-vulkan/shaders/gaussian-sort.slang +++ b/backend/merlin-vulkan/shaders/gaussian-sort.slang @@ -1,73 +1,4 @@ -// Deterministic global sort of GPU-prepared Gaussian records. Each resident -// resource contributes one key segment; an LSD radix sort over 8-bit digits -// then orders the whole frame back to front, and a verification pass checks -// the result against the prepared records it indexes. -// -// The 64-bit key reproduces the CPU reference order exactly: the high word -// maps a larger sort key to a smaller unsigned value, and the low word is the -// frame-global candidate index, assigned by ascending resource identity and -// then particle index. Atomic compaction order therefore never reaches the -// sorted result, and the low word is unique for every real record. -// -// Every kernel uses one portable reduce-then-scan formulation without -// subgroup operations or forward-progress assumptions, and binds at most four -// storage buffers, the Vulkan guaranteed minimum. - -static const uint kGaussianSortWorkgroupSize = 256u; -static const uint kGaussianSortRadixBins = 256u; -static const uint kGaussianSortScanElementsPerThread = 4u; -static const uint kGaussianSortScanWorkgroupElements = 1024u; -static const uint kGaussianSortInvalidValue = 0xFFFFFFFFu; -static const uint kGaussianSortDynamicCount = 1u; - -// Verification words at the start of the scan/control buffer. -static const uint kGaussianSortSortedCount = 0u; -static const uint kGaussianSortOrderViolationCount = 1u; -static const uint kGaussianSortKeyMismatchCount = 2u; -static const uint kGaussianSortIdentityChecksum = 3u; - -struct GaussianSortConstants -{ - // Keys: segment length. Histogram, scatter, and verify: padded key count, - // or the buffer capacity when flags selects a dynamic count. - uint element_count; - uint block_count; - uint digit_shift; - uint digit_word; - uint scan_offset; - uint scan_count; - uint scan_sums_offset; - uint candidate_base; - uint prepared_base; - uint visible_count_offset; - // Histogram and scatter: with kGaussianSortDynamicCount, only the first - // min(element_count, scan[count_word]) elements participate. - uint count_word; - uint flags; -}; - -struct GaussianSortElement -{ - uint key_low; - uint key_high; - uint value; -}; - -// Must match gaussian-prepare.slang; reflection checks both layouts. -struct GaussianPreparedRecord -{ - float2 center_pixels; - float radius_pixels; - float depth; - float3 inverse_conic; - float opacity; - float3 radiance; - float sort_key; - uint resource_id_low; - uint resource_id_high; - uint particle_id; - uint padding; -}; +#include "../../../core/merlin-render-backend/shaders/gaussian-sort-abi.slang" [[vk::push_constant]] ConstantBuffer gaussian_sort_constants; @@ -79,252 +10,9 @@ ConstantBuffer gaussian_sort_constants; [[vk::binding(2, 0)]] RWStructuredBuffer gaussian_sort_scan; [[vk::binding(3, 0)]] StructuredBuffer gaussian_sort_prepared_records; -groupshared uint gs_sort_bins[kGaussianSortRadixBins]; -groupshared uint gs_sort_digits[kGaussianSortWorkgroupSize]; -groupshared uint gs_sort_scan[kGaussianSortWorkgroupSize]; - -// Back to front: a larger sort key must produce a smaller unsigned key. The -// CPU reference compares floats, so -0 and +0 have to share one key. -uint GaussianSortKeyHigh(float sort_key) -{ - uint bits = sort_key == 0.0f ? 0u : asuint(sort_key); - uint ordered = (bits & 0x80000000u) != 0u ? ~bits : (bits | 0x80000000u); - return ~ordered; -} - -// Order-sensitive identity mix shared with the CPU reference checksum. -uint GaussianSortMix(uint value) -{ - value ^= value >> 16u; - value *= 0x85EBCA6Bu; - value ^= value >> 13u; - value *= 0xC2B2AE35u; - value ^= value >> 16u; - return value; -} - -uint GaussianSortIdentityHash(uint resource_low, uint resource_high, uint particle) -{ - uint hash = GaussianSortMix(particle * 0x9E3779B9u ^ resource_low); - return GaussianSortMix(hash ^ resource_high * 0x85EBCA6Bu); -} - -GaussianSortElement SentinelElement() -{ - GaussianSortElement element; - element.key_low = 0xFFFFFFFFu; - element.key_high = 0xFFFFFFFFu; - element.value = kGaussianSortInvalidValue; - return element; -} - -uint Digit(GaussianSortElement element) -{ - uint word = gaussian_sort_constants.digit_word == 0u - ? element.key_low : element.key_high; - return (word >> gaussian_sort_constants.digit_shift) & (kGaussianSortRadixBins - 1u); -} - -uint ActiveElementCount() -{ - uint count = gaussian_sort_constants.element_count; - if ((gaussian_sort_constants.flags & kGaussianSortDynamicCount) != 0u) - { - count = min(count, gaussian_sort_scan[gaussian_sort_constants.count_word]); - } - return count; -} - -bool KeyLess(GaussianSortElement lhs, GaussianSortElement rhs) -{ - return lhs.key_high < rhs.key_high || - (lhs.key_high == rhs.key_high && lhs.key_low < rhs.key_low); -} - -// One dispatch per resource. Threads past the resource's visible count, and -// the padding appended to the last segment, write sentinels that sort last. -[shader("compute")] -[numthreads(256, 1, 1)] -void gaussian_sort_keys(uint3 dispatch_thread_id : SV_DispatchThreadID) -{ - uint index = dispatch_thread_id.x; - if (index >= gaussian_sort_constants.element_count) - { - return; - } - GaussianSortElement element = SentinelElement(); - if (index < gaussian_sort_scan[gaussian_sort_constants.visible_count_offset]) - { - uint record_index = gaussian_sort_constants.prepared_base + index; - GaussianPreparedRecord record = gaussian_sort_prepared_records[record_index]; - element.key_low = gaussian_sort_constants.candidate_base + record.particle_id; - element.key_high = GaussianSortKeyHigh(record.sort_key); - element.value = record_index; - } - gaussian_sort_destination[gaussian_sort_constants.candidate_base + index] = element; -} - -// Per-workgroup digit counts, stored digit-major so one exclusive scan over -// the whole array yields every workgroup's global scatter base. -[shader("compute")] -[numthreads(256, 1, 1)] -void gaussian_sort_histogram(uint3 group_id : SV_GroupID, - uint3 group_thread_id : SV_GroupThreadID, - uint3 dispatch_thread_id : SV_DispatchThreadID) +GaussianPreparedRecord LoadSortPreparedRecord(uint index) { - uint lane = group_thread_id.x; - gs_sort_bins[lane] = 0u; - GroupMemoryBarrierWithGroupSync(); - if (dispatch_thread_id.x < ActiveElementCount()) - { - uint ignored; - InterlockedAdd(gs_sort_bins[Digit(gaussian_sort_source[dispatch_thread_id.x])], - 1u, ignored); - } - GroupMemoryBarrierWithGroupSync(); - gaussian_sort_scan[gaussian_sort_constants.scan_offset + - lane * gaussian_sort_constants.block_count + group_id.x] = - gs_sort_bins[lane]; + return gaussian_sort_prepared_records[index]; } -// Exclusive scan of one 1024-element block per workgroup; each block total is -// written to the next level so a later scan and add make the scan global. -[shader("compute")] -[numthreads(256, 1, 1)] -void gaussian_sort_scan_blocks(uint3 group_id : SV_GroupID, - uint3 group_thread_id : SV_GroupThreadID) -{ - uint lane = group_thread_id.x; - uint first = group_id.x * kGaussianSortScanWorkgroupElements + - lane * kGaussianSortScanElementsPerThread; - uint values[kGaussianSortScanElementsPerThread]; - uint total = 0u; - for (uint i = 0u; i < kGaussianSortScanElementsPerThread; ++i) - { - uint index = first + i; - values[i] = index < gaussian_sort_constants.scan_count - ? gaussian_sort_scan[gaussian_sort_constants.scan_offset + index] : 0u; - total += values[i]; - } - gs_sort_scan[lane] = total; - GroupMemoryBarrierWithGroupSync(); - for (uint offset = 1u; offset < kGaussianSortWorkgroupSize; offset <<= 1u) - { - uint addend = lane >= offset ? gs_sort_scan[lane - offset] : 0u; - GroupMemoryBarrierWithGroupSync(); - gs_sort_scan[lane] += addend; - GroupMemoryBarrierWithGroupSync(); - } - uint running = lane == 0u ? 0u : gs_sort_scan[lane - 1u]; - for (uint i = 0u; i < kGaussianSortScanElementsPerThread; ++i) - { - uint index = first + i; - if (index < gaussian_sort_constants.scan_count) - { - gaussian_sort_scan[gaussian_sort_constants.scan_offset + index] = running; - } - running += values[i]; - } - if (lane == kGaussianSortWorkgroupSize - 1u) - { - gaussian_sort_scan[gaussian_sort_constants.scan_sums_offset + group_id.x] = - gs_sort_scan[lane]; - } -} - -[shader("compute")] -[numthreads(256, 1, 1)] -void gaussian_sort_scan_add(uint3 group_id : SV_GroupID, - uint3 group_thread_id : SV_GroupThreadID) -{ - uint addend = gaussian_sort_scan[gaussian_sort_constants.scan_sums_offset + group_id.x]; - for (uint i = 0u; i < kGaussianSortScanElementsPerThread; ++i) - { - uint index = group_id.x * kGaussianSortScanWorkgroupElements + - i * kGaussianSortWorkgroupSize + group_thread_id.x; - if (index < gaussian_sort_constants.scan_count) - { - gaussian_sort_scan[gaussian_sort_constants.scan_offset + index] += addend; - } - } -} - -// Stable scatter: the destination is the workgroup's global digit base plus -// the number of earlier elements in the same workgroup with the same digit. -[shader("compute")] -[numthreads(256, 1, 1)] -void gaussian_sort_scatter(uint3 group_id : SV_GroupID, - uint3 group_thread_id : SV_GroupThreadID, - uint3 dispatch_thread_id : SV_DispatchThreadID) -{ - uint lane = group_thread_id.x; - bool active = dispatch_thread_id.x < ActiveElementCount(); - GaussianSortElement element = SentinelElement(); - uint digit = kGaussianSortRadixBins; - if (active) - { - element = gaussian_sort_source[dispatch_thread_id.x]; - digit = Digit(element); - } - gs_sort_digits[lane] = digit; - gs_sort_bins[lane] = gaussian_sort_scan[ - gaussian_sort_constants.scan_offset + - lane * gaussian_sort_constants.block_count + group_id.x]; - GroupMemoryBarrierWithGroupSync(); - if (!active) - { - return; - } - uint rank = 0u; - for (uint i = 0u; i < lane; ++i) - { - rank += gs_sort_digits[i] == digit ? 1u : 0u; - } - gaussian_sort_destination[gs_sort_bins[digit] + rank] = element; -} - -// Checks strict key order, sentinel placement, and each key against the -// prepared record it indexes, and accumulates an order-sensitive identity -// checksum the CPU reference can reproduce. -[shader("compute")] -[numthreads(256, 1, 1)] -void gaussian_sort_verify(uint3 dispatch_thread_id : SV_DispatchThreadID) -{ - uint index = dispatch_thread_id.x; - if (index >= gaussian_sort_constants.element_count) - { - return; - } - GaussianSortElement element = gaussian_sort_source[index]; - bool valid = element.value != kGaussianSortInvalidValue; - uint ignored; - if (index + 1u < gaussian_sort_constants.element_count) - { - GaussianSortElement next = gaussian_sort_source[index + 1u]; - bool ordered = valid ? KeyLess(element, next) : !KeyLess(next, element); - if (!ordered || (!valid && next.value != kGaussianSortInvalidValue)) - { - InterlockedAdd(gaussian_sort_scan[kGaussianSortOrderViolationCount], 1u, ignored); - } - } - if (!valid) - { - if (element.key_low != 0xFFFFFFFFu || element.key_high != 0xFFFFFFFFu) - { - InterlockedAdd(gaussian_sort_scan[kGaussianSortKeyMismatchCount], 1u, ignored); - } - return; - } - GaussianPreparedRecord record = gaussian_sort_prepared_records[element.value]; - if (element.key_high != GaussianSortKeyHigh(record.sort_key) || - element.key_low < record.particle_id) - { - InterlockedAdd(gaussian_sort_scan[kGaussianSortKeyMismatchCount], 1u, ignored); - } - InterlockedAdd(gaussian_sort_scan[kGaussianSortSortedCount], 1u, ignored); - InterlockedAdd(gaussian_sort_scan[kGaussianSortIdentityChecksum], - (index + 1u) * GaussianSortIdentityHash( - record.resource_id_low, record.resource_id_high, - record.particle_id), - ignored); -} +#include "../../../core/merlin-render-backend/shaders/gaussian-sort-common.slang" diff --git a/cmake/MerlinMetalShaders.cmake b/cmake/MerlinMetalShaders.cmake index 5fbcc87..1cfe6ab 100644 --- a/cmake/MerlinMetalShaders.cmake +++ b/cmake/MerlinMetalShaders.cmake @@ -13,19 +13,39 @@ macro(merlin_add_metal_shaders) set(_metal_artifacts "") set(_metal_air "") set(_metal_records "") - foreach(_stage vertex fragment) - set(_entry "gaussian_metal_${_stage}") - set(_msl "${_metal_output}/gaussian.${_stage}.metal") + set(_metal_entries + "gaussian-metal.slang|gaussian_metal_vertex|vertex|gaussian.vertex|gaussian-reference|prepared_stream+ellipse+color+ids" + "gaussian-metal.slang|gaussian_metal_fragment|fragment|gaussian.fragment|gaussian-reference|prepared_stream+ellipse+color+ids" + "gaussian-prepare-metal.slang|gaussian_prepare_compact|compute|gaussian.prepare|gaussian-prepare|projection+sh+compaction") + foreach(_sort keys histogram scan_blocks scan_add scatter verify) + list(APPEND _metal_entries + "gaussian-sort-metal.slang|gaussian_sort_${_sort}|compute|gaussian.sort-${_sort}|gaussian-sort-${_sort}|deterministic_radix_sort") + endforeach() + foreach(_record IN LISTS _metal_entries) + string(REPLACE "|" ";" _fields "${_record}") + list(GET _fields 0 _source) + list(GET _fields 1 _entry) + list(GET _fields 2 _stage) + list(GET _fields 3 _name) + list(GET _fields 4 _permutation) + list(GET _fields 5 _features) + set(_msl "${_metal_output}/${_name}.metal") set(_reflection "${_msl}.reflection.json") set(_depfile "${_msl}.d") + set(_metal_compile_options "") + if(_stage STREQUAL "compute") + # Preparation classifies non-finite inputs. Preserve those checks. + # Slang emits unused temporaries for the return value of atomics. + set(_metal_compile_options -fno-fast-math -Wno-unused-variable -Werror) + endif() add_custom_command( OUTPUT "${_msl}" "${_reflection}" COMMAND ${CMAKE_COMMAND} -E make_directory "${_metal_output}" - COMMAND "${MERLIN_SLANGC_EXECUTABLE}" gaussian-metal.slang + COMMAND "${MERLIN_SLANGC_EXECUTABLE}" "${_source}" -entry "${_entry}" -stage "${_stage}" -target metal -profile metallib_2_4 -matrix-layout-column-major -O2 -warnings-as-errors all -reflection-json "${_reflection}" -depfile "${_depfile}" -o "${_msl}" - DEPENDS "${_metal_source}/gaussian-metal.slang" + DEPENDS "${_metal_source}/${_source}" "${_metal_shared}/gaussian-raster-common.slang" "${_metal_shared}/gaussian-raster-abi.slang" DEPFILE "${_depfile}" @@ -35,6 +55,7 @@ macro(merlin_add_metal_shaders) OUTPUT "${_msl}.air" COMMAND "${MERLIN_XCRUN_EXECUTABLE}" -sdk macosx metal -std=macos-metal2.4 + ${_metal_compile_options} "-mmacosx-version-min=${_metal_deployment_target}" "-fmodules-cache-path=${CMAKE_CURRENT_BINARY_DIR}/metal-module-cache" -c "${_msl}" -o "${_msl}.air" @@ -43,7 +64,7 @@ macro(merlin_add_metal_shaders) list(APPEND _metal_artifacts "${_msl}" "${_reflection}") list(APPEND _metal_air "${_msl}.air") string(APPEND _metal_records - "${_msl}|${_reflection}|${_depfile}|gaussian-metal.slang|${_entry}|${_stage}|metal|metallib_2_4|none|gaussian-reference|prepared_stream+ellipse+color+ids\n") + "${_msl}|${_reflection}|${_depfile}|${_source}|${_entry}|${_stage}|metal|metallib_2_4|none|${_permutation}|${_features}\n") endforeach() set(_metal_library "${_metal_output}/gaussian.metallib") add_custom_command( diff --git a/core/merlin-render-backend/shaders/gaussian-prepare-abi.slang b/core/merlin-render-backend/shaders/gaussian-prepare-abi.slang new file mode 100644 index 0000000..bc05469 --- /dev/null +++ b/core/merlin-render-backend/shaders/gaussian-prepare-abi.slang @@ -0,0 +1,52 @@ +#pragma once + +// First GPU-driven Gaussian preparation contract. One dispatch processes one +// resident Gaussian resource directly from its four device-local attribute +// ranges, classifies every particle, and atomically compacts visible records. +// A later radix-sort/tile slice consumes the stable sort key and identity. + +static const uint kGaussianProjectionPerspective = 0u; +static const uint kGaussianProjectionTangential = 1u; +static const uint kGaussianSortingZDepth = 0u; +static const uint kGaussianSortingCameraDistance = 1u; + +static const uint kGaussianCandidateVisible = 0u; +static const uint kGaussianCandidateOpacityCulled = 1u; +static const uint kGaussianCandidateFrustumCulled = 2u; +static const uint kGaussianCandidateInvalidCulled = 3u; + +static const float kProjectionEpsilon = 1.0e-6f; +// Matches the CPU reference: the 3DGS guard band bounds the perspective +// Jacobian of close off-axis kernels. +static const float kJacobianGuardBand = 1.3f; + +struct GaussianPrepareConstants +{ + float4x4 local_to_camera; + float4x4 projection; + float2 viewport_size; + float sigma_extent; + float minimum_variance_pixels; + uint resource_id_low; + uint resource_id_high; + uint particle_count; + uint coefficients_per_particle; + uint spherical_harmonics_degree; + uint projection_mode; + uint sorting_mode; + uint padding; +}; + +#include "gaussian-prepared-record.slang" + +struct GaussianPrepareDispatchCounters +{ + uint candidate_count; + uint visible_count; + uint opacity_culled_count; + uint frustum_culled_count; + uint invalid_culled_count; + uint padding0; + uint padding1; + uint padding2; +}; diff --git a/core/merlin-render-backend/shaders/gaussian-prepare-common.slang b/core/merlin-render-backend/shaders/gaussian-prepare-common.slang new file mode 100644 index 0000000..c8022c7 --- /dev/null +++ b/core/merlin-render-backend/shaders/gaussian-prepare-common.slang @@ -0,0 +1,304 @@ +// Resource declarations belong to the backend wrapper. + +bool IsFinite1(float value) +{ + return !isnan(value) && !isinf(value); +} + +bool IsFinite2(float2 value) +{ + return all(!isnan(value)) && all(!isinf(value)); +} + +bool IsFinite3(float3 value) +{ + return all(!isnan(value)) && all(!isinf(value)); +} + +float3 SafeNormalize(float3 value) +{ + float length_squared = dot(value, value); + if (!IsFinite1(length_squared) || length_squared <= kProjectionEpsilon * kProjectionEpsilon) + { + return float3(0.0f); + } + return value * rsqrt(length_squared); +} + +float3 LoadPosition(uint particle) +{ + return gaussian_positions.Load(particle * 12u); +} + +float3x3 LoadCovariance(uint particle) +{ + uint offset = particle * 24u; + float3 first = gaussian_covariances.Load(offset); + float3 second = gaussian_covariances.Load(offset + 12u); + return float3x3( + first.x, first.y, first.z, + first.y, second.x, second.y, + first.z, second.y, second.z); +} + +float3 LoadCoefficient(uint index) +{ + return gaussian_radiance.Load(index * 12u); +} + +void AddScaled(inout float3 result, uint coefficient, float basis) +{ + result += LoadCoefficient(coefficient) * basis; +} + +float3 EvaluateRadiance(uint particle, float3 direction) +{ + static const float c0 = 0.28209479177387814f; + static const float c1 = 0.4886025119029199f; + uint first = particle * gaussian_prepare_constants.coefficients_per_particle; + float3 result = float3(0.5f) + LoadCoefficient(first) * c0; + uint degree = min(gaussian_prepare_constants.spherical_harmonics_degree, 3u); + if (degree >= 1u) + { + AddScaled(result, first + 1u, -c1 * direction.y); + AddScaled(result, first + 2u, c1 * direction.z); + AddScaled(result, first + 3u, -c1 * direction.x); + } + float xx = direction.x * direction.x; + float yy = direction.y * direction.y; + float zz = direction.z * direction.z; + if (degree >= 2u) + { + AddScaled(result, first + 4u, 1.0925484305920792f * direction.x * direction.y); + AddScaled(result, first + 5u, 1.0925484305920792f * direction.y * direction.z); + AddScaled(result, first + 6u, 0.31539156525252005f * (2.0f * zz - xx - yy)); + AddScaled(result, first + 7u, 1.0925484305920792f * direction.x * direction.z); + AddScaled(result, first + 8u, 0.5462742152960396f * (xx - yy)); + } + if (degree >= 3u) + { + AddScaled(result, first + 9u, -0.5900435899266435f * direction.y * (3.0f * xx - yy)); + AddScaled(result, first + 10u, 2.890611442640554f * direction.x * direction.y * direction.z); + AddScaled(result, first + 11u, -0.4570457994644658f * direction.y * (4.0f * zz - xx - yy)); + AddScaled(result, first + 12u, 0.3731763325901154f * direction.z * (2.0f * zz - 3.0f * xx - 3.0f * yy)); + AddScaled(result, first + 13u, -0.4570457994644658f * direction.x * (4.0f * zz - xx - yy)); + AddScaled(result, first + 14u, 1.445305721320277f * direction.z * (xx - yy)); + AddScaled(result, first + 15u, -0.5900435899266435f * direction.x * (xx - 3.0f * yy)); + } + return max(result, float3(0.0f)); +} + +struct GaussianProjectionJacobian +{ + float3 dx; + float3 dy; +}; + +GaussianProjectionJacobian ProjectionJacobian(float3 camera, float4 clip) +{ + float3 row_x = gaussian_prepare_constants.projection[0].xyz; + float3 row_y = gaussian_prepare_constants.projection[1].xyz; + float3 row_w = gaussian_prepare_constants.projection[3].xyz; + float inverse_w_squared = 1.0f / (clip.w * clip.w); + float guard = kJacobianGuardBand * clip.w; + float3 dx = (row_x * clip.w - clamp(clip.x, -guard, guard) * row_w) * inverse_w_squared; + float3 dy = (row_y * clip.w - clamp(clip.y, -guard, guard) * row_w) * inverse_w_squared; + + if (gaussian_prepare_constants.projection_mode == kGaussianProjectionTangential && + dot(row_w, row_w) > kProjectionEpsilon * kProjectionEpsilon) + { + float distance = length(camera); + float3 direction = SafeNormalize(camera); + float3 tangent_x = SafeNormalize(float3( + 1.0f - direction.x * direction.x, + -direction.x * direction.y, + -direction.x * direction.z)); + if (dot(tangent_x, tangent_x) <= kProjectionEpsilon * kProjectionEpsilon) + { + tangent_x = SafeNormalize(float3( + -direction.y * direction.x, + 1.0f - direction.y * direction.y, + -direction.y * direction.z)); + } + float3 tangent_y = SafeNormalize(cross(tangent_x, direction)); + if (tangent_y.y < 0.0f) + { + tangent_y = -tangent_y; + } + dx = tangent_x * (max(abs(gaussian_prepare_constants.projection[0][0]), kProjectionEpsilon) / distance); + dy = tangent_y * (max(abs(gaussian_prepare_constants.projection[1][1]), kProjectionEpsilon) / distance); + } + GaussianProjectionJacobian result; + result.dx = dx; + result.dy = dy; + return result; +} + +bool ProjectCovariance(float3x3 covariance, GaussianProjectionJacobian jacobian, + out float3 inverse_conic, out float radius_pixels) +{ + inverse_conic = float3(0.0f); + radius_pixels = 0.0f; + float3 dx = jacobian.dx; + float3 dy = jacobian.dy; + float width_scale = gaussian_prepare_constants.viewport_size.x * 0.5f; + float height_scale = gaussian_prepare_constants.viewport_size.y * 0.5f; + float xx = dot(dx, mul(covariance, dx)) * width_scale * width_scale + + gaussian_prepare_constants.minimum_variance_pixels; + float xy = dot(dx, mul(covariance, dy)) * width_scale * height_scale; + float yy = dot(dy, mul(covariance, dy)) * height_scale * height_scale + + gaussian_prepare_constants.minimum_variance_pixels; + float determinant = xx * yy - xy * xy; + if (!IsFinite1(xx) || !IsFinite1(xy) || !IsFinite1(yy) || + !IsFinite1(determinant) || determinant <= kProjectionEpsilon) + { + return false; + } + float discriminant = sqrt(max(0.0f, (xx - yy) * (xx - yy) + 4.0f * xy * xy)); + float largest_eigenvalue = 0.5f * (xx + yy + discriminant); + radius_pixels = gaussian_prepare_constants.sigma_extent * sqrt(max(largest_eigenvalue, 0.0f)); + inverse_conic = float3(yy / determinant, -xy / determinant, xx / determinant); + return IsFinite1(radius_pixels) && radius_pixels > 0.0f && IsFinite3(inverse_conic); +} + +bool ProjectDepthExtent(float3x3 covariance, float3 camera, float4 clip, + out float radius) +{ + radius = 0.0f; + float3 row_z = gaussian_prepare_constants.projection[2].xyz; + float3 row_w = gaussian_prepare_constants.projection[3].xyz; + float inverse_w_squared = 1.0f / (clip.w * clip.w); + float3 dz = (row_z * clip.w - clip.z * row_w) * inverse_w_squared; + float variance = dot(dz, mul(covariance, dz)); + if (!IsFinite1(variance) || variance < -kProjectionEpsilon) + { + return false; + } + radius = gaussian_prepare_constants.sigma_extent * sqrt(max(variance, 0.0f)); + return IsFinite1(radius); +} + +float3 LocalCameraDirection(float3x3 local_to_camera, float3 camera) +{ + float3 row0 = local_to_camera[0]; + float3 row1 = local_to_camera[1]; + float3 row2 = local_to_camera[2]; + float3 cofactor0 = cross(row1, row2); + float3 cofactor1 = cross(row2, row0); + float3 cofactor2 = cross(row0, row1); + float determinant_value = dot(row0, cofactor0); + if (!IsFinite1(determinant_value) || abs(determinant_value) <= kProjectionEpsilon) + { + return SafeNormalize(camera); + } + float inverse_determinant = 1.0f / determinant_value; + // Cofactors of the rows are columns of the inverse. Dotting them with + // camera instead applies the inverse transpose and rotates SH incorrectly. + return SafeNormalize((cofactor0 * camera.x + cofactor1 * camera.y + + cofactor2 * camera.z) * inverse_determinant); +} + +void RecordRejection(uint particle, uint result, inout uint counter) +{ + gaussian_candidate_results[particle] = result; + InterlockedAdd(counter, 1u); +} + +[shader("compute")] +[numthreads(64, 1, 1)] +void gaussian_prepare_compact(uint3 dispatch_thread_id : SV_DispatchThreadID) +{ + uint particle = dispatch_thread_id.x; + if (particle >= gaussian_prepare_constants.particle_count) + { + return; + } + if (particle == 0u) + { + gaussian_prepare_counters[0].candidate_count = gaussian_prepare_constants.particle_count; + } + + float opacity = gaussian_opacities.Load(particle * 4u); + if (!IsFinite1(opacity) || opacity <= 0.0f) + { + RecordRejection(particle, kGaussianCandidateOpacityCulled, + gaussian_prepare_counters[0].opacity_culled_count); + return; + } + + float3 local_position = LoadPosition(particle); + float4 camera4 = mul(gaussian_prepare_constants.local_to_camera, float4(local_position, 1.0f)); + float3 camera = camera4.xyz; + float4 clip = mul(gaussian_prepare_constants.projection, float4(camera, 1.0f)); + if (!IsFinite3(camera) || any(isnan(clip)) || any(isinf(clip)) || clip.w <= kProjectionEpsilon) + { + RecordRejection(particle, kGaussianCandidateFrustumCulled, + gaussian_prepare_counters[0].frustum_culled_count); + return; + } + float3 ndc = clip.xyz / clip.w; + // Kernels centered in front of the near plane are clipped like Mesh + // geometry; their center-evaluated footprint would cover the view. + if (ndc.z < 0.0f) + { + RecordRejection(particle, kGaussianCandidateFrustumCulled, + gaussian_prepare_counters[0].frustum_culled_count); + return; + } + float3x3 local_linear = (float3x3)gaussian_prepare_constants.local_to_camera; + float3x3 covariance = LoadCovariance(particle); + float3x3 camera_covariance = mul(mul(local_linear, covariance), transpose(local_linear)); + float depth_radius; + if (!ProjectDepthExtent(camera_covariance, camera, clip, depth_radius)) + { + RecordRejection(particle, kGaussianCandidateInvalidCulled, + gaussian_prepare_counters[0].invalid_culled_count); + return; + } + if (ndc.z - depth_radius > 1.0f) + { + RecordRejection(particle, kGaussianCandidateFrustumCulled, + gaussian_prepare_counters[0].frustum_culled_count); + return; + } + float3 inverse_conic; + float radius_pixels; + if (!ProjectCovariance(camera_covariance, ProjectionJacobian(camera, clip), + inverse_conic, radius_pixels)) + { + RecordRejection(particle, kGaussianCandidateInvalidCulled, + gaussian_prepare_counters[0].invalid_culled_count); + return; + } + + float2 center_pixels = (ndc.xy * 0.5f + 0.5f) * gaussian_prepare_constants.viewport_size; + if (!IsFinite2(center_pixels) || + center_pixels.x + radius_pixels < 0.0f || + center_pixels.x - radius_pixels > gaussian_prepare_constants.viewport_size.x || + center_pixels.y + radius_pixels < 0.0f || + center_pixels.y - radius_pixels > gaussian_prepare_constants.viewport_size.y) + { + RecordRejection(particle, kGaussianCandidateFrustumCulled, + gaussian_prepare_counters[0].frustum_culled_count); + return; + } + + uint compact_index; + InterlockedAdd(gaussian_prepare_counters[0].visible_count, 1u, compact_index); + GaussianPreparedRecord prepared; + prepared.center_pixels = center_pixels; + prepared.radius_pixels = radius_pixels; + prepared.depth = clamp(ndc.z, 0.0f, 1.0f); + prepared.inverse_conic = inverse_conic; + prepared.opacity = opacity; + prepared.radiance = EvaluateRadiance( + particle, LocalCameraDirection(local_linear, camera)); + prepared.sort_key = gaussian_prepare_constants.sorting_mode == kGaussianSortingCameraDistance + ? dot(camera, camera) : ndc.z; + prepared.resource_id_low = gaussian_prepare_constants.resource_id_low; + prepared.resource_id_high = gaussian_prepare_constants.resource_id_high; + prepared.particle_id = particle; + prepared.padding = 0u; + StorePreparedRecord(compact_index, prepared); + gaussian_candidate_results[particle] = kGaussianCandidateVisible; +} diff --git a/core/merlin-render-backend/shaders/gaussian-prepared-record.slang b/core/merlin-render-backend/shaders/gaussian-prepared-record.slang new file mode 100644 index 0000000..f79fb1e --- /dev/null +++ b/core/merlin-render-backend/shaders/gaussian-prepared-record.slang @@ -0,0 +1,18 @@ +#pragma once + +// Kept at 64 bytes so it can become the direct input of sorted/tiled raster +// without repacking. Resource + particle identity is a deterministic tie break. +struct GaussianPreparedRecord +{ + float2 center_pixels; + float radius_pixels; + float depth; + float3 inverse_conic; + float opacity; + float3 radiance; + float sort_key; + uint resource_id_low; + uint resource_id_high; + uint particle_id; + uint padding; +}; diff --git a/core/merlin-render-backend/shaders/gaussian-sort-abi.slang b/core/merlin-render-backend/shaders/gaussian-sort-abi.slang new file mode 100644 index 0000000..cb785c3 --- /dev/null +++ b/core/merlin-render-backend/shaders/gaussian-sort-abi.slang @@ -0,0 +1,58 @@ +#pragma once + +// Deterministic global sort of GPU-prepared Gaussian records. Each resident +// resource contributes one key segment; an LSD radix sort over 8-bit digits +// then orders the whole frame back to front, and a verification pass checks +// the result against the prepared records it indexes. +// +// The 64-bit key reproduces the CPU reference order exactly: the high word +// maps a larger sort key to a smaller unsigned value, and the low word is the +// frame-global candidate index, assigned by ascending resource identity and +// then particle index. Atomic compaction order therefore never reaches the +// sorted result, and the low word is unique for every real record. +// +// Every kernel uses one portable reduce-then-scan formulation without +// subgroup operations or forward-progress assumptions, and binds at most four +// storage buffers, the Vulkan guaranteed minimum. + +static const uint kGaussianSortWorkgroupSize = 256u; +static const uint kGaussianSortRadixBins = 256u; +static const uint kGaussianSortScanElementsPerThread = 4u; +static const uint kGaussianSortScanWorkgroupElements = 1024u; +static const uint kGaussianSortInvalidValue = 0xFFFFFFFFu; +static const uint kGaussianSortDynamicCount = 1u; + +// Verification words at the start of the scan/control buffer. +static const uint kGaussianSortSortedCount = 0u; +static const uint kGaussianSortOrderViolationCount = 1u; +static const uint kGaussianSortKeyMismatchCount = 2u; +static const uint kGaussianSortIdentityChecksum = 3u; + +struct GaussianSortConstants +{ + // Keys: segment length. Histogram, scatter, and verify: padded key count, + // or the buffer capacity when flags selects a dynamic count. + uint element_count; + uint block_count; + uint digit_shift; + uint digit_word; + uint scan_offset; + uint scan_count; + uint scan_sums_offset; + uint candidate_base; + uint prepared_base; + uint visible_count_offset; + // Histogram and scatter: with kGaussianSortDynamicCount, only the first + // min(element_count, scan[count_word]) elements participate. + uint count_word; + uint flags; +}; + +struct GaussianSortElement +{ + uint key_low; + uint key_high; + uint value; +}; + +#include "gaussian-prepared-record.slang" diff --git a/core/merlin-render-backend/shaders/gaussian-sort-common.slang b/core/merlin-render-backend/shaders/gaussian-sort-common.slang new file mode 100644 index 0000000..1e57709 --- /dev/null +++ b/core/merlin-render-backend/shaders/gaussian-sort-common.slang @@ -0,0 +1,257 @@ +// Resource declarations belong to the backend wrapper. +// Inline helpers so independently compiled Metal entry points can link into +// one metallib without duplicate externally visible helper symbols. + +groupshared uint gs_sort_bins[kGaussianSortRadixBins]; +groupshared uint gs_sort_digits[kGaussianSortWorkgroupSize]; +groupshared uint gs_sort_scan[kGaussianSortWorkgroupSize]; + +// Back to front: a larger sort key must produce a smaller unsigned key. The +// CPU reference compares floats, so -0 and +0 have to share one key. +[ForceInline] +uint GaussianSortKeyHigh(float sort_key) +{ + uint bits = sort_key == 0.0f ? 0u : asuint(sort_key); + uint ordered = (bits & 0x80000000u) != 0u ? ~bits : (bits | 0x80000000u); + return ~ordered; +} + +// Order-sensitive identity mix shared with the CPU reference checksum. +[ForceInline] +uint GaussianSortMix(uint value) +{ + value ^= value >> 16u; + value *= 0x85EBCA6Bu; + value ^= value >> 13u; + value *= 0xC2B2AE35u; + value ^= value >> 16u; + return value; +} + +[ForceInline] +uint GaussianSortIdentityHash(uint resource_low, uint resource_high, uint particle) +{ + uint hash = GaussianSortMix(particle * 0x9E3779B9u ^ resource_low); + return GaussianSortMix(hash ^ resource_high * 0x85EBCA6Bu); +} + +[ForceInline] +GaussianSortElement SentinelElement() +{ + GaussianSortElement element; + element.key_low = 0xFFFFFFFFu; + element.key_high = 0xFFFFFFFFu; + element.value = kGaussianSortInvalidValue; + return element; +} + +[ForceInline] +uint Digit(GaussianSortElement element) +{ + uint word = gaussian_sort_constants.digit_word == 0u + ? element.key_low : element.key_high; + return (word >> gaussian_sort_constants.digit_shift) & (kGaussianSortRadixBins - 1u); +} + +[ForceInline] +uint ActiveElementCount() +{ + uint count = gaussian_sort_constants.element_count; + if ((gaussian_sort_constants.flags & kGaussianSortDynamicCount) != 0u) + { + count = min(count, gaussian_sort_scan[gaussian_sort_constants.count_word]); + } + return count; +} + +[ForceInline] +bool KeyLess(GaussianSortElement lhs, GaussianSortElement rhs) +{ + return lhs.key_high < rhs.key_high || + (lhs.key_high == rhs.key_high && lhs.key_low < rhs.key_low); +} + +// One dispatch per resource. Threads past the resource's visible count, and +// the padding appended to the last segment, write sentinels that sort last. +[shader("compute")] +[numthreads(256, 1, 1)] +void gaussian_sort_keys(uint3 dispatch_thread_id : SV_DispatchThreadID) +{ + uint index = dispatch_thread_id.x; + if (index >= gaussian_sort_constants.element_count) + { + return; + } + GaussianSortElement element = SentinelElement(); + if (index < gaussian_sort_scan[gaussian_sort_constants.visible_count_offset]) + { + uint record_index = gaussian_sort_constants.prepared_base + index; + GaussianPreparedRecord record = LoadSortPreparedRecord(record_index); + element.key_low = gaussian_sort_constants.candidate_base + record.particle_id; + element.key_high = GaussianSortKeyHigh(record.sort_key); + element.value = record_index; + } + gaussian_sort_destination[gaussian_sort_constants.candidate_base + index] = element; +} + +// Per-workgroup digit counts, stored digit-major so one exclusive scan over +// the whole array yields every workgroup's global scatter base. +[shader("compute")] +[numthreads(256, 1, 1)] +void gaussian_sort_histogram(uint3 group_id : SV_GroupID, + uint3 group_thread_id : SV_GroupThreadID, + uint3 dispatch_thread_id : SV_DispatchThreadID) +{ + uint lane = group_thread_id.x; + gs_sort_bins[lane] = 0u; + GroupMemoryBarrierWithGroupSync(); + if (dispatch_thread_id.x < ActiveElementCount()) + { + InterlockedAdd(gs_sort_bins[Digit(gaussian_sort_source[dispatch_thread_id.x])], + 1u); + } + GroupMemoryBarrierWithGroupSync(); + gaussian_sort_scan[gaussian_sort_constants.scan_offset + + lane * gaussian_sort_constants.block_count + group_id.x] = + gs_sort_bins[lane]; +} + +// Exclusive scan of one 1024-element block per workgroup; each block total is +// written to the next level so a later scan and add make the scan global. +[shader("compute")] +[numthreads(256, 1, 1)] +void gaussian_sort_scan_blocks(uint3 group_id : SV_GroupID, + uint3 group_thread_id : SV_GroupThreadID) +{ + uint lane = group_thread_id.x; + uint first = group_id.x * kGaussianSortScanWorkgroupElements + + lane * kGaussianSortScanElementsPerThread; + uint values[kGaussianSortScanElementsPerThread]; + uint total = 0u; + for (uint i = 0u; i < kGaussianSortScanElementsPerThread; ++i) + { + uint index = first + i; + values[i] = index < gaussian_sort_constants.scan_count + ? gaussian_sort_scan[gaussian_sort_constants.scan_offset + index] : 0u; + total += values[i]; + } + gs_sort_scan[lane] = total; + GroupMemoryBarrierWithGroupSync(); + for (uint offset = 1u; offset < kGaussianSortWorkgroupSize; offset <<= 1u) + { + uint addend = lane >= offset ? gs_sort_scan[lane - offset] : 0u; + GroupMemoryBarrierWithGroupSync(); + gs_sort_scan[lane] += addend; + GroupMemoryBarrierWithGroupSync(); + } + uint running = lane == 0u ? 0u : gs_sort_scan[lane - 1u]; + for (uint i = 0u; i < kGaussianSortScanElementsPerThread; ++i) + { + uint index = first + i; + if (index < gaussian_sort_constants.scan_count) + { + gaussian_sort_scan[gaussian_sort_constants.scan_offset + index] = running; + } + running += values[i]; + } + if (lane == kGaussianSortWorkgroupSize - 1u) + { + gaussian_sort_scan[gaussian_sort_constants.scan_sums_offset + group_id.x] = + gs_sort_scan[lane]; + } +} + +[shader("compute")] +[numthreads(256, 1, 1)] +void gaussian_sort_scan_add(uint3 group_id : SV_GroupID, + uint3 group_thread_id : SV_GroupThreadID) +{ + uint addend = gaussian_sort_scan[gaussian_sort_constants.scan_sums_offset + group_id.x]; + for (uint i = 0u; i < kGaussianSortScanElementsPerThread; ++i) + { + uint index = group_id.x * kGaussianSortScanWorkgroupElements + + i * kGaussianSortWorkgroupSize + group_thread_id.x; + if (index < gaussian_sort_constants.scan_count) + { + gaussian_sort_scan[gaussian_sort_constants.scan_offset + index] += addend; + } + } +} + +// Stable scatter: the destination is the workgroup's global digit base plus +// the number of earlier elements in the same workgroup with the same digit. +[shader("compute")] +[numthreads(256, 1, 1)] +void gaussian_sort_scatter(uint3 group_id : SV_GroupID, + uint3 group_thread_id : SV_GroupThreadID, + uint3 dispatch_thread_id : SV_DispatchThreadID) +{ + uint lane = group_thread_id.x; + bool active = dispatch_thread_id.x < ActiveElementCount(); + GaussianSortElement element = SentinelElement(); + uint digit = kGaussianSortRadixBins; + if (active) + { + element = gaussian_sort_source[dispatch_thread_id.x]; + digit = Digit(element); + } + gs_sort_digits[lane] = digit; + gs_sort_bins[lane] = gaussian_sort_scan[ + gaussian_sort_constants.scan_offset + + lane * gaussian_sort_constants.block_count + group_id.x]; + GroupMemoryBarrierWithGroupSync(); + if (!active) + { + return; + } + uint rank = 0u; + for (uint i = 0u; i < lane; ++i) + { + rank += gs_sort_digits[i] == digit ? 1u : 0u; + } + gaussian_sort_destination[gs_sort_bins[digit] + rank] = element; +} + +// Checks strict key order, sentinel placement, and each key against the +// prepared record it indexes, and accumulates an order-sensitive identity +// checksum the CPU reference can reproduce. +[shader("compute")] +[numthreads(256, 1, 1)] +void gaussian_sort_verify(uint3 dispatch_thread_id : SV_DispatchThreadID) +{ + uint index = dispatch_thread_id.x; + if (index >= gaussian_sort_constants.element_count) + { + return; + } + GaussianSortElement element = gaussian_sort_source[index]; + bool valid = element.value != kGaussianSortInvalidValue; + if (index + 1u < gaussian_sort_constants.element_count) + { + GaussianSortElement next = gaussian_sort_source[index + 1u]; + bool ordered = valid ? KeyLess(element, next) : !KeyLess(next, element); + if (!ordered || (!valid && next.value != kGaussianSortInvalidValue)) + { + InterlockedAdd(gaussian_sort_scan[kGaussianSortOrderViolationCount], 1u); + } + } + if (!valid) + { + if (element.key_low != 0xFFFFFFFFu || element.key_high != 0xFFFFFFFFu) + { + InterlockedAdd(gaussian_sort_scan[kGaussianSortKeyMismatchCount], 1u); + } + return; + } + GaussianPreparedRecord record = LoadSortPreparedRecord(element.value); + if (element.key_high != GaussianSortKeyHigh(record.sort_key) || + element.key_low < record.particle_id) + { + InterlockedAdd(gaussian_sort_scan[kGaussianSortKeyMismatchCount], 1u); + } + InterlockedAdd(gaussian_sort_scan[kGaussianSortSortedCount], 1u); + InterlockedAdd(gaussian_sort_scan[kGaussianSortIdentityChecksum], + (index + 1u) * GaussianSortIdentityHash( + record.resource_id_low, record.resource_id_high, + record.particle_id)); +} diff --git a/docs/design/metal-gaussian-execution.md b/docs/design/metal-gaussian-execution.md index 1a3393f..467d96f 100644 --- a/docs/design/metal-gaussian-execution.md +++ b/docs/design/metal-gaussian-execution.md @@ -51,10 +51,29 @@ with Vulkan through Slang. Metal-specific entry points preserve the scalar 52-byte stream, Y convention and combined color/ID attachments. Slang 2026.8.x emits MSL/reflection, Xcode produces an embedded metallib, and the versioned install package retains those artifacts and their identity/checksum evidence. -Host layout assertions, generated binding checks and local Apple GPU image tests -cover this first step. Projection/SH/sort shader sharing, controlled performance -captures and updated native viewport/HgiMetal comparisons remain unfinished; -this is not completion of the phase gate below. +Projection, covariance, SH evaluation, culling/compaction and portable radix-sort +kernels now also live in shared Slang sources, with backend-owned bindings and +record access. Metal uses explicit buffer slots and byte-addressed 64-byte +prepared records because native MSL `float3` storage has different padding. +The embedded/installed library includes these compute entry points. Preparation +preserves non-finite classification with native Metal fast math disabled. + +Host layout assertions, generated binding checks and local Apple GPU tests +compare preparation and sorting with the CPU reference, including SH degrees +0–3, projection/sorting policies, nonidentity transforms, culling boundaries, +multiple resources, hierarchical scans, device-written counts and empty input. +The kernel harness allows absolute error of `2e-4 * max(1, abs(reference))` +per floating-point field; classifications, counters, identities and sentinel +placement must match exactly. Intentional distance ties use binary-exact input +coordinates so rounding of nearly equal CPU/GPU keys is not mistaken for a +radix-sort error. These checks do not establish renderer image parity or a +general exact ordering guarantee for numerically near-equal keys. +The compute test submits preparation through sorting without an intermediate +CPU readback. It is a kernel correctness harness, not renderer GPU execution: +persistent attributes, completion-safe frame scheduling, gather/indirect raster +and telemetry integration remain Phase 2 work. Controlled performance captures +and updated native viewport/HgiMetal comparisons remain unfinished; this is +not completion of the phase gate below. Establish repeatable static, camera-motion and particle-edit captures before changing execution. Use the existing public Gaussian corpus and deterministic diff --git a/docs/reference/support-matrix.md b/docs/reference/support-matrix.md index a0fe52d..52f5324 100644 --- a/docs/reference/support-matrix.md +++ b/docs/reference/support-matrix.md @@ -74,7 +74,7 @@ between separately produced OpenUSD SDKs remain the operator's responsibility. | MaterialXGenSlang material-function prototype | Available: optional `Merlin::MaterialX` emits deterministic graph-only Slang functions and renderer-owned minimum Standard Surface results for constants, image/UV0/world-normal, add/multiply/mix, `base`, `base_color`, `metalness`, `specular_roughness`, and `normal`. Portable library/include fingerprints feed topology-only module keys separated from parameter/resource state. Registered parameter-only and texture/sampler artifacts execute in renderer-owned Vulkan Forward after ABI/reflection checks, reuse pipelines across value and texture-content edits, and report structured fallback/capability telemetry. The same sources retain installed SPIR-V, Metal-target, and reflection evidence. General MaterialX documents, tangent-space normal mapping, production IBL, and Hydra MaterialX ingestion are not claimed | | Material ABI agreement | `merlin.material-abi/v1` is available in Core. A consumer declares the result fields it reads and the geometry inputs it can build, and a module is checked against that rather than in isolation. A compiled artifact's reflected interface is checked back against the module's logical one by name, type, and array size, so SPIR-V and Metal describe one material through their own native bindings and agree with each other by agreeing with the module. The same contract owns the pass-neutrality rule, which `Merlin::MaterialX` applies to its own output; both generated modules and all four of their SPIR-V/Metal artifacts are checked in the test suite, including the dropped, retyped, and undeclared-parameter cases | | Native Metal backend and residency | Available for offscreen Mesh Forward: native device/queue, runtime MSL, buffers/textures/samplers, heap residency, generation-checked argument-buffer tables with conventional fallback, frames-in-flight retirement, basic material/opacity mask, color/depth/primId/instanceId AOVs, CPU readback, capacity diagnostics, and Metal-specific telemetry | -| Metal Gaussian rendering | Available: shared CPU projection/covariance/SH evaluation and sorting, shared Slang raster math compiled into an embedded Metal library, alpha composition with opaque Mesh depth, resource/particle IDs, and immutable prepared-stream reuse. Image tests cover edits, camera/resize changes and in-flight lifetime; the prior embedded-MSL path was checked in the OpenUSD 26.08 native viewport and usdview/HgiMetal GPU-copy display with a 5.8-million-particle stage. The Slang path was also manually checked in usdview/HgiMetal with that stage; the native viewport comparison remains open. GPU preparation/sorting/tiling are not implemented: `prefer` falls back and `require` rejects. | +| Metal Gaussian rendering | Available: shared CPU projection/covariance/SH evaluation and sorting, shared Slang raster math compiled into an embedded Metal library, alpha composition with opaque Mesh depth, resource/particle IDs, and immutable prepared-stream reuse. Image tests cover edits, camera/resize changes and in-flight lifetime; the prior embedded-MSL path was checked in the OpenUSD 26.08 native viewport and usdview/HgiMetal GPU-copy display with a 5.8-million-particle stage. The Slang path was also manually checked in usdview/HgiMetal with that stage; the native viewport comparison remains open. Shared GPU preparation/sort kernels have Metal ABI and CPU-reference compute tests, but are not connected to renderer execution; persistent residency, indirect raster and tiling remain open. `prefer` falls back and `require` rejects. | | Native Metal viewport presentation | Available: adapter-owned `CAMetalLayer`, renderer-owned drawable encoding, GPU-only offscreen-to-drawable presentation, resize/frames-in-flight safety, sRGB/Display P3 SDR policy with an explicit future HDR boundary, vsync/drawable-count pacing, Dear ImGui integration, presentation telemetry, and exact offscreen reference parity | | HgiVulkan host presentation bridge | Available: public-driver discovery, Hgi-owned color targets, Tier 0 upload fallback, and selected color GPU copy on validated OpenUSD 26.05/26.08 packages when their imported `hgiVulkan` target is present. Merlin borrows the Hgi Vulkan 1.3 device/graphics queue for its conventional renderer path, exports one color image, copies it with explicit barriers, and releases its lease from Hgi command-buffer completion; depth and id AOVs stay on CPU readback. Runtime comparison evidence covers image parity and no color Map/upload or coarse wait through resize. Merlin-owned Vulkan remains 1.4. Direct sharing reports `public-texture-import-unavailable` because public Hgi cannot import a Merlin-owned `VkImage`, so GPU copy remains selected. | | HgiMetal host presentation bridge | Available on validated OpenUSD 26.05/26.08 packages: Hgi-owned color targets receive same-device Metal GPU copies from leased renderer AOVs with completion-safe resize retirement and Tier 0 fallback. Direct sharing remains rejected because the public host texture-import contract is unavailable. | diff --git a/docs/roadmap/current.md b/docs/roadmap/current.md index e664614..c0cfe0e 100644 --- a/docs/roadmap/current.md +++ b/docs/roadmap/current.md @@ -20,8 +20,10 @@ ownership, dependency, and fallback contracts. - 🚧 Advance the [Metal Gaussian execution plan](../design/metal-gaussian-execution.md). Reference raster math is shared through Slang and an embedded metallib, with ABI, artifact identity, install-consumer and local Apple GPU image checks. - Controlled static/motion/edit captures, shared GPU preparation/sort shaders, - native viewport/HgiMetal rechecks, and complete GPU execution remain open. + GPU preparation/SH/radix-sort kernels are also shared and checked against the + CPU reference on Apple GPU. Controlled static/motion/edit captures, native + viewport/HgiMetal rechecks, persistent attribute residency and complete GPU + execution through raster remain open. - 🚧 Recheck the Hydra CPU-readback floor in usdview after preferring host-cached coherent Vulkan AOV buffers. Local RTX A5000 headless evidence diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 90c1bf3..d6f7440 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -210,7 +210,7 @@ if(TARGET Merlin::Metal) add_dependencies(merlin-metal-shader-artifact-key-tests merlin-metal-shaders) add_test(NAME merlin-metal-shader-artifact-key COMMAND merlin-metal-shader-artifact-key-tests - "${MERLIN_METAL_SHADER_OUTPUT_DIR}/manifest.json" 2 1) + "${MERLIN_METAL_SHADER_OUTPUT_DIR}/manifest.json" 9 3) set_tests_properties(merlin-metal-shader-artifact-key PROPERTIES LABELS "metal;shader;cache") add_executable(merlin-metal-gaussian-tests metal_gaussian_test.cpp) @@ -221,6 +221,20 @@ if(TARGET Merlin::Metal) LABELS "metal;gaussian;gpu;image;aov" SKIP_RETURN_CODE 77 TIMEOUT 60) + add_executable(merlin-metal-gaussian-compute-tests metal_gaussian_compute_test.mm) + merlin_target_defaults(merlin-metal-gaussian-compute-tests) + target_link_libraries(merlin-metal-gaussian-compute-tests PRIVATE Merlin::Metal) + target_compile_options(merlin-metal-gaussian-compute-tests PRIVATE -fobjc-arc) + set_target_properties(merlin-metal-gaussian-compute-tests PROPERTIES + OBJCXX_STANDARD 20 OBJCXX_STANDARD_REQUIRED ON OBJCXX_EXTENSIONS OFF) + add_dependencies(merlin-metal-gaussian-compute-tests merlin-metal-shaders) + add_test(NAME merlin-metal-gaussian-compute + COMMAND merlin-metal-gaussian-compute-tests + "${MERLIN_METAL_SHADER_OUTPUT_DIR}/gaussian.metallib") + set_tests_properties(merlin-metal-gaussian-compute PROPERTIES + LABELS "metal;gaussian;gpu;compute;cpu-reference" + SKIP_RETURN_CODE 77 TIMEOUT 60) + add_executable(merlin-metal-resource-table-tests metal_resource_table_test.cpp ) diff --git a/tests/compile-gaussian-spirv.cmake b/tests/compile-gaussian-spirv.cmake index 9249911..30a1809 100644 --- a/tests/compile-gaussian-spirv.cmake +++ b/tests/compile-gaussian-spirv.cmake @@ -15,3 +15,23 @@ foreach(_entry gaussian_vertex gaussian_id_vertex gaussian_fragment gaussian_id_ message(FATAL_ERROR "Vulkan Gaussian compilation failed: ${_entry}") endif() endforeach() + +# Metal-only builders also compile the production Vulkan wrappers of shared +# preparation and sorting, so changes cannot silently break the other target. +get_filename_component(_source_dir "${MERLIN_GAUSSIAN_SOURCE}" DIRECTORY) +foreach(_kernel prepare_compact sort_keys sort_histogram sort_scan_blocks sort_scan_add sort_scatter sort_verify) + if(_kernel MATCHES "^prepare") + set(_source gaussian-prepare.slang) + else() + set(_source gaussian-sort.slang) + endif() + execute_process( + COMMAND "${MERLIN_SLANGC_EXECUTABLE}" "${_source_dir}/${_source}" + -entry "gaussian_${_kernel}" -stage compute -target spirv -profile sm_6_6 + -capability spirv_1_5 -matrix-layout-column-major -O2 -warnings-as-errors all + -o "${MERLIN_TEST_OUTPUT_DIR}/gaussian_${_kernel}.spv" + RESULT_VARIABLE _result) + if(NOT _result EQUAL 0) + message(FATAL_ERROR "Vulkan Gaussian compilation failed: gaussian_${_kernel}") + endif() +endforeach() diff --git a/tests/metal_gaussian_compute_test.mm b/tests/metal_gaussian_compute_test.mm new file mode 100644 index 0000000..b1ed6ec --- /dev/null +++ b/tests/metal_gaussian_compute_test.mm @@ -0,0 +1,378 @@ +#include "../backend/merlin-metal/src/gaussian_compute_abi.hpp" +#include + +#import + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { +using namespace merlin::metal::gaussian_compute; +using merlin::extraction::FrameSnapshot; + +void Require(bool condition, const char* message) { + if (!condition) throw std::runtime_error(message); +} + +void Near(float actual, float expected, const char* field) { + // Field tolerance is documented in the Metal Gaussian execution design. + if (!std::isfinite(actual) || + std::abs(actual - expected) > 2.0e-4F * std::max(1.0F, std::abs(expected))) { + std::cerr << field << " actual=" << actual << " expected=" << expected << '\n'; + throw std::runtime_error("Metal preparation differs from the CPU reference"); + } +} + +std::uint32_t Groups(std::uint32_t count, std::uint32_t size) { + return (count + size - 1) / size; +} + +id Buffer(id device, std::size_t bytes, const void* data = nullptr) { + auto buffer = [device newBufferWithLength:std::max(bytes, std::size_t{16}) + options:MTLResourceStorageModeShared]; + Require(buffer != nil, "Metal buffer allocation failed"); + std::memset(buffer.contents, 0, buffer.length); + if (data && bytes) std::memcpy(buffer.contents, data, bytes); + return buffer; +} + +template +id Upload(id device, const std::vector& values) { + return Buffer(device, values.size() * sizeof(T), values.data()); +} + +struct Kernels { + id prepare; + id keys; + id histogram; + id scan; + id add; + id scatter; + id verify; +}; + +id Pipeline(id device, id library, + NSString* name) { + auto function = [library newFunctionWithName:name]; + Require(function != nil, "Missing Metal compute entry point"); + NSError* error = nil; + auto pipeline = [device newComputePipelineStateWithFunction:function error:&error]; + if (!pipeline) throw std::runtime_error(error.localizedDescription.UTF8String); + Require(pipeline.maxTotalThreadsPerThreadgroup >= 256, + "Metal device cannot run the portable radix workgroup"); + return pipeline; +} + +// A complete preparation + global sort in one submission. Counts are copied +// device-to-device; no intermediate readback schedules subsequent kernels. +void Compare(id device, id queue, const Kernels& kernels, + FrameSnapshot snapshot, bool dynamic_count = false) { + std::vector ordered(snapshot.gaussians.begin(), snapshot.gaussians.end()); + std::sort(ordered.begin(), ordered.end(), + [](const auto& a, const auto& b) { return a.gaussian < b.gaussian; }); + snapshot.gaussians.assign(std::move(ordered)); + const auto reference = merlin::extraction::PrepareGaussianFrame(snapshot, {320, 192}); + std::uint32_t total = 0; + for (const auto& record : snapshot.gaussians) + total += static_cast(record.positions->size()); + const std::uint32_t padded = std::max(256U, Groups(total, 256) * 256); + auto prepared = Buffer(device, padded * sizeof(PreparedRecord)); + auto source = Buffer(device, padded * sizeof(SortElement)); + auto destination = Buffer(device, padded * sizeof(SortElement)); + // Four verification words, then per-resource counts, then hierarchical scan. + const auto histogram_offset = 4U + static_cast(snapshot.gaussians.size()); + auto control = Buffer(device, (histogram_offset + padded * 2U + 16U) * sizeof(std::uint32_t)); + auto command = [queue commandBuffer]; + Require(command != nil, "Metal command allocation failed"); + std::vector> counters; + std::vector> classifications; + const auto policy = merlin::extraction::SelectGaussianSortingPolicy(snapshot); + auto sort_dispatch = [&](id pipeline, + const SortConstants& constants, std::uint32_t groups) { + auto encoder = [command computeCommandEncoder]; + [encoder setComputePipelineState:pipeline]; + [encoder setBuffer:source offset:0 atIndex:0]; + [encoder setBuffer:destination offset:0 atIndex:1]; + [encoder setBuffer:control offset:0 atIndex:2]; + [encoder setBuffer:prepared offset:0 atIndex:3]; + [encoder setBytes:&constants length:sizeof(constants) atIndex:4]; + [encoder dispatchThreadgroups:MTLSizeMake(groups, 1, 1) + threadsPerThreadgroup:MTLSizeMake(256, 1, 1)]; + [encoder endEncoding]; + }; + std::uint32_t base = 0; + for (std::size_t i = 0; i < snapshot.gaussians.size(); ++i) { + const auto& record = snapshot.gaussians[i]; + const auto count = static_cast(record.positions->size()); + PrepareConstants constants; + constants.local_to_camera = record.transform; // Fixtures use identity view. + constants.projection = snapshot.projection; + constants.viewport_size = {320, 192}; + constants.resource_id_low = static_cast(record.gaussian); + constants.resource_id_high = static_cast(record.gaussian >> 32U); + constants.particle_count = count; + constants.spherical_harmonics_degree = record.spherical_harmonics_degree; + constants.coefficients_per_particle = (record.spherical_harmonics_degree + 1U) * + (record.spherical_harmonics_degree + 1U); + constants.projection_mode = static_cast(record.projection_mode); + constants.sorting_mode = static_cast(policy.mode); + counters.push_back(Buffer(device, sizeof(PrepareCounters))); + classifications.push_back(Buffer(device, count * sizeof(std::uint32_t))); + if (count && record.visible) { + auto encoder = [command computeCommandEncoder]; + [encoder setComputePipelineState:kernels.prepare]; + [encoder setBuffer:Upload(device, *record.positions) offset:0 atIndex:0]; + [encoder setBuffer:Upload(device, *record.covariances) offset:0 atIndex:1]; + [encoder setBuffer:Upload(device, *record.opacities) offset:0 atIndex:2]; + [encoder setBuffer:Upload(device, *record.spherical_harmonics_coefficients) offset:0 atIndex:3]; + [encoder setBuffer:classifications.back() offset:0 atIndex:4]; + [encoder setBuffer:prepared offset:base * sizeof(PreparedRecord) atIndex:5]; + [encoder setBuffer:counters.back() offset:0 atIndex:6]; + [encoder setBytes:&constants length:sizeof(constants) atIndex:7]; + [encoder dispatchThreadgroups:MTLSizeMake(Groups(count, 64), 1, 1) + threadsPerThreadgroup:MTLSizeMake(64, 1, 1)]; + [encoder endEncoding]; + } + auto blit = [command blitCommandEncoder]; + [blit copyFromBuffer:counters.back() sourceOffset:offsetof(PrepareCounters, visible_count) + toBuffer:control destinationOffset:(4 + i) * sizeof(std::uint32_t) + size:sizeof(std::uint32_t)]; + [blit endEncoding]; + SortConstants keys; + keys.element_count = i + 1 == snapshot.gaussians.size() ? padded - base : count; + keys.candidate_base = base; + keys.prepared_base = base; + keys.visible_count_offset = 4U + static_cast(i); + if (keys.element_count) sort_dispatch(kernels.keys, keys, Groups(keys.element_count, 256)); + base += count; + } + // A resource-free frame still initializes the sentinel stream on the GPU. + if (snapshot.gaussians.empty()) { + SortConstants keys; + keys.element_count = padded; + sort_dispatch(kernels.keys, keys, Groups(padded, 256)); + } + std::swap(source, destination); + for (std::uint32_t digit = 0; digit < 8; ++digit) { + SortConstants constants; + constants.element_count = padded; + constants.block_count = Groups(padded, 256); + constants.digit_word = digit / 4; + constants.digit_shift = (digit % 4) * 8; + constants.scan_offset = histogram_offset; + if (dynamic_count) { + // Single-resource fixtures use the GPU-written visible count. Padding + // remains in the buffers but does not participate in this sort. + Require(snapshot.gaussians.size() == 1, "Dynamic-count fixture must use one resource"); + constants.flags = 1; + constants.count_word = 4; + } + sort_dispatch(kernels.histogram, constants, constants.block_count); + std::vector levels; + auto level = constants; + level.scan_count = padded; + for (;;) { + level.scan_sums_offset = level.scan_offset + level.scan_count; + const auto groups = Groups(level.scan_count, 1024); + sort_dispatch(kernels.scan, level, groups); + levels.push_back(level); + if (groups == 1) break; + level.scan_offset = level.scan_sums_offset; + level.scan_count = groups; + } + for (std::size_t i = levels.size() - 1; i > 0; --i) + sort_dispatch(kernels.add, levels[i - 1], Groups(levels[i - 1].scan_count, 1024)); + sort_dispatch(kernels.scatter, constants, constants.block_count); + std::swap(source, destination); + } + SortConstants verify; + verify.element_count = dynamic_count ? static_cast(reference.gaussians.size()) : padded; + if (verify.element_count) sort_dispatch(kernels.verify, verify, Groups(verify.element_count, 256)); + [command commit]; + [command waitUntilCompleted]; + if (command.status != MTLCommandBufferStatusCompleted) + throw std::runtime_error(command.error.localizedDescription.UTF8String); + const auto* verification = static_cast(control.contents); + Require(verification[0] == reference.gaussians.size(), "Wrong sorted count"); + Require(verification[1] == 0 && verification[2] == 0, "GPU sort verification failed"); + const auto* sorted = static_cast(source.contents); + const auto* records = static_cast(prepared.contents); + for (std::size_t i = 0; i < reference.gaussians.size(); ++i) { + Require(sorted[i].value < padded, "Sorted record index exceeds capacity"); + const auto& actual = records[sorted[i].value]; + const auto& expected = reference.gaussians[i]; + if (actual.particle_id != expected.particle) { + std::cerr << "identity at " << i << " particles=" << total << " degree=" + << snapshot.gaussians[0].spherical_harmonics_degree + << " projection=" << int(snapshot.gaussians[0].projection_mode) + << " sorting=" << int(policy.mode) << " dynamic=" << dynamic_count + << " actual=" << actual.particle_id << " key=" << actual.sort_key + << " expected=" << expected.particle << " key=" << expected.sort_key << '\n'; + } + Require((std::uint64_t{actual.resource_id_high} << 32U | actual.resource_id_low) == expected.resource && + actual.particle_id == expected.particle, "Non-deterministic sorted identity"); + Near(actual.center_pixels.x, expected.center_pixels.x, "center_pixels.x"); + Near(actual.center_pixels.y, expected.center_pixels.y, "center_pixels.y"); + Near(actual.radius_pixels, expected.radius_pixels, "radius_pixels"); + Near(actual.depth, expected.depth, "depth"); + Near(actual.inverse_conic.x, expected.inverse_conic.x, "inverse_conic.x"); + Near(actual.inverse_conic.y, expected.inverse_conic.y, "inverse_conic.y"); + Near(actual.inverse_conic.z, expected.inverse_conic.z, "inverse_conic.z"); + Near(actual.opacity, expected.opacity, "opacity"); + Near(actual.radiance.x, expected.radiance.x, "radiance.x"); + Near(actual.radiance.y, expected.radiance.y, "radiance.y"); + Near(actual.radiance.z, expected.radiance.z, "radiance.z"); + Near(actual.sort_key, expected.sort_key, "sort_key"); + } + if (!dynamic_count) + for (std::size_t i = reference.gaussians.size(); i < padded; ++i) + Require(sorted[i].value == UINT32_MAX && sorted[i].key_low == UINT32_MAX && + sorted[i].key_high == UINT32_MAX, "Invalid sentinel tail"); + std::array sums{}; + for (std::size_t i = 0; i < counters.size(); ++i) { + const auto& record = snapshot.gaussians[i]; + const auto& actual = *static_cast(counters[i].contents); + Require(actual.candidate_count == (record.visible ? record.positions->size() : 0), + "Wrong candidate count"); + const std::array counts{actual.visible_count, actual.opacity_culled_count, + actual.frustum_culled_count, actual.invalid_culled_count}; + std::array classified{}; + if (record.visible) { + const auto* values = static_cast(classifications[i].contents); + for (std::size_t j = 0; j < record.positions->size(); ++j) { + Require(values[j] < 4, "Invalid candidate classification"); + ++classified[values[j]]; + } + } + Require(counts == classified, "Candidate classifications disagree with counters"); + for (std::size_t j = 0; j < counts.size(); ++j) sums[j] += counts[j]; + } + Require(sums[0] == reference.counters.visible_count && + sums[1] == reference.counters.opacity_culled_count && + sums[2] == reference.counters.frustum_culled_count && + sums[3] == reference.counters.invalid_culled_count, "GPU culling differs from CPU"); +} + +FrameSnapshot Fixture(std::uint32_t count, std::uint32_t degree, bool perspective) { + FrameSnapshot snapshot; + if (perspective) snapshot.projection.values = { + 1.2F, 0, 0, 0, 0, 1.5F, 0, 0, 0, 0, -1.001001F, -1, + 0, 0, -0.1001001F, 0}; + merlin::extraction::GaussianRecord record; + record.gaussian = 0x200000003ULL; + record.spherical_harmonics_degree = degree; + std::vector positions(count); + std::vector covariances(count, {0.0004F, 0.0001F, 0, 0.0002F, 0, 0.0001F}); + std::vector opacities(count, 0.65F); + const auto coefficients = (degree + 1) * (degree + 1); + std::vector radiance(count * coefficients); + for (std::uint32_t i = 0; i < count; ++i) { + // Binary-exact coordinates keep intentional distance ties exact on CPU + // and GPU; unrelated floating-point near-ties are not sort failures. + positions[i] = {float(int(i % 17) - 8) * 0.0625F, + float(int(i % 13) - 6) * 0.0625F, + perspective ? -1.0F - float(i % 7) * 0.25F : 0.25F + float(i % 7) * 0.0625F}; + for (std::uint32_t j = 0; j < coefficients; ++j) + radiance[i * coefficients + j] = {float(j + 1) * 0.02F, -0.03F, 0.04F}; + // Duplicate depths exercise deterministic tie-breaking across workgroups. + if (i % 31 == 0) opacities[i] = 0; + if (i % 37 == 0) positions[i].x = 1000; + } + if (count >= 8) { + opacities[1] = std::numeric_limits::quiet_NaN(); + positions[2].z = perspective ? -0.05F : -0.01F; // Near plane. + positions[3].z = perspective ? -101.0F : 1.01F; // Conservative far bound. + positions[4].z = perspective ? -200.0F : 2.0F; + positions[5].x = std::numeric_limits::infinity(); + covariances[6].xx = std::numeric_limits::quiet_NaN(); + covariances[7].zz = -1.0F; + } + record.positions = std::make_shared>(positions); + record.covariances = std::make_shared>(covariances); + record.opacities = std::make_shared>(opacities); + record.spherical_harmonics_coefficients = std::make_shared>(radiance); + snapshot.gaussians.push_back(record); + return snapshot; +} +} // namespace + +int main(int argc, char** argv) { + @autoreleasepool { + try { + Require(argc == 2, "Expected compiled Metal library path"); + auto device = MTLCreateSystemDefaultDevice(); + if (!device) { std::cerr << "skip: no Metal device\n"; return 77; } + NSError* error = nil; + auto library = [device newLibraryWithURL:[NSURL fileURLWithPath:@(argv[1])] error:&error]; + if (!library) throw std::runtime_error(error.localizedDescription.UTF8String); + auto queue = [device newCommandQueue]; + Require(queue != nil, "Metal queue creation failed"); + const Kernels kernels{ + Pipeline(device, library, @"gaussian_prepare_compact"), + Pipeline(device, library, @"gaussian_sort_keys"), + Pipeline(device, library, @"gaussian_sort_histogram"), + Pipeline(device, library, @"gaussian_sort_scan_blocks"), + Pipeline(device, library, @"gaussian_sort_scan_add"), + Pipeline(device, library, @"gaussian_sort_scatter"), + Pipeline(device, library, @"gaussian_sort_verify")}; + Compare(device, queue, kernels, {}); + for (auto count : {0U, 1U, 63U, 64U, 65U, 255U, 256U, 257U, 1301U}) + Compare(device, queue, kernels, Fixture(count, 0, false)); + for (std::uint32_t degree = 0; degree <= 3; ++degree) { + for (bool perspective : {false, true}) { + for (auto projection : {merlin::GaussianProjectionMode::Perspective, + merlin::GaussianProjectionMode::Tangential}) { + for (auto sorting : {merlin::GaussianSortingMode::ZDepth, + merlin::GaussianSortingMode::CameraDistance}) { + auto snapshot = Fixture(1301, degree, perspective); + auto first = snapshot.gaussians[0]; + first.projection_mode = projection; + first.sorting_mode = sorting; + snapshot.gaussians.assign({first}); + Compare(device, queue, kernels, snapshot, true); + auto second = first; + second.gaussian = 0x100000004ULL; + snapshot.gaussians.push_back(second); + Compare(device, queue, kernels, snapshot); + } + } + } + } + auto transformed = Fixture(257, 3, true); + auto transformed_record = transformed.gaussians[0]; + transformed_record.transform.values = { + 0.8F, 0.2F, 0, 0, -0.3F, 1.1F, 0, 0, 0, 0, 1.2F, 0, + 0.1F, -0.1F, -0.2F, 1}; + transformed.gaussians.assign({transformed_record}); + Compare(device, queue, kernels, transformed); + auto mixed = Fixture(257, 1, true); + auto second = mixed.gaussians[0]; + second.gaussian = 1; + second.sorting_mode = merlin::GaussianSortingMode::CameraDistance; + mixed.gaussians.push_back(second); + Compare(device, queue, kernels, mixed); + auto hidden = mixed.gaussians[0]; + hidden.visible = false; + mixed.gaussians.assign({hidden, second}); + Compare(device, queue, kernels, mixed); + auto culled = Fixture(257, 0, false); + auto culled_record = culled.gaussians[0]; + culled_record.opacities = std::make_shared>(257, 0); + culled.gaussians.assign({culled_record}); + Compare(device, queue, kernels, culled); + Compare(device, queue, kernels, culled, true); + std::cout << "Metal preparation and radix sort match CPU: " << device.name.UTF8String << '\n'; + return 0; + } catch (const std::exception& error) { + std::cerr << error.what() << '\n'; + return 1; + } + } +} diff --git a/tests/run_install_consumer.cmake b/tests/run_install_consumer.cmake index e450c0f..426cd71 100644 --- a/tests/run_install_consumer.cmake +++ b/tests/run_install_consumer.cmake @@ -168,6 +168,13 @@ if(MERLIN_EXPECTED_METAL) message(FATAL_ERROR "Metal install is missing ${_artifact}") endif() endforeach() + foreach(_kernel prepare sort-keys sort-histogram sort-scan_blocks sort-scan_add sort-scatter sort-verify) + foreach(_suffix metal metal.reflection.json) + if(NOT EXISTS "${_metal_shader_dir}/gaussian.${_kernel}.${_suffix}") + message(FATAL_ERROR "Metal install is missing gaussian.${_kernel}.${_suffix}") + endif() + endforeach() + endforeach() file(SHA256 "${_metal_shader_dir}/gaussian.metallib" _metal_library_hash) file(READ "${_metal_shader_dir}/gaussian.metallib.sha256" _metal_library_checksum) if(NOT _metal_library_checksum STREQUAL "${_metal_library_hash} gaussian.metallib\n") diff --git a/tests/verify-metal-gaussian-abi.cmake b/tests/verify-metal-gaussian-abi.cmake index 384361e..8202a21 100644 --- a/tests/verify-metal-gaussian-abi.cmake +++ b/tests/verify-metal-gaussian-abi.cmake @@ -39,3 +39,84 @@ foreach(_binding "gaussian_instances_[0-9]+ \\[\\[buffer\\(0\\)\\]\\]" message(FATAL_ERROR "Generated Metal binding is missing: ${_binding}") endif() endforeach() + +# Compute wrappers bind one Metal buffer namespace. Check both reflection and +# emitted attributes, including kernels where dead resources are optimized out. +function(check_compute _file _entry _threads) + file(READ "${MERLIN_METAL_SHADER_DIR}/${_file}.metal.reflection.json" _json) + file(READ "${MERLIN_METAL_SHADER_DIR}/${_file}.metal" _source) + require_json("${_json}" "${_entry}" entryPoints 0 name) + require_json("${_json}" compute entryPoints 0 stage) + require_json("${_json}" "${_threads}" entryPoints 0 threadGroupSize 0) + require_json("${_json}" 1 entryPoints 0 threadGroupSize 1) + require_json("${_json}" 1 entryPoints 0 threadGroupSize 2) + set(_index 0) + foreach(_binding IN LISTS ARGN) + string(REPLACE ":" ";" _parts "${_binding}") + list(GET _parts 0 _name) + list(GET _parts 1 _slot) + require_json("${_json}" "${_name}" parameters ${_index} name) + if(_name MATCHES "_constants$") + require_json("${_json}" constantBuffer parameters ${_index} binding kind) + require_json("${_json}" "${_slot}" parameters ${_index} binding index) + else() + require_json("${_json}" constantBuffer parameters ${_index} bindings 0 kind) + require_json("${_json}" "${_slot}" parameters ${_index} bindings 0 index) + endif() + # Only bindings actually used by this entry point appear in the MSL. + if(_source MATCHES "${_name}_[0-9]+ \\[\\[buffer") + if(NOT _source MATCHES "${_name}_[0-9]+ \\[\\[buffer\\(${_slot}\\)\\]\\]") + message(FATAL_ERROR "Generated compute binding differs: ${_name}") + endif() + endif() + math(EXPR _index "${_index} + 1") + endforeach() + set(_compute_json "${_json}" PARENT_SCOPE) +endfunction() + +function(check_fields _json _parameter _names _offsets) + set(_index 0) + foreach(_name IN LISTS _names) + list(GET _offsets ${_index} _offset) + require_json("${_json}" "${_name}" parameters ${_parameter} type elementType fields ${_index} name) + require_json("${_json}" "${_offset}" parameters ${_parameter} type elementType fields ${_index} binding offset) + math(EXPR _index "${_index} + 1") + endforeach() +endfunction() + +check_compute(gaussian.prepare gaussian_prepare_compact 64 + gaussian_prepare_constants:7 gaussian_positions:0 gaussian_covariances:1 + gaussian_opacities:2 gaussian_radiance:3 gaussian_candidate_results:4 + gaussian_prepared_records:5 gaussian_prepare_counters:6) +require_json("${_compute_json}" 176 parameters 0 type elementVarLayout binding size) +check_fields("${_compute_json}" 0 + "local_to_camera;projection;viewport_size;sigma_extent;minimum_variance_pixels;resource_id_low;resource_id_high;particle_count;coefficients_per_particle;spherical_harmonics_degree;projection_mode;sorting_mode;padding" + "0;64;128;136;140;144;148;152;156;160;164;168;172") +foreach(_matrix 0 1) + require_json("${_compute_json}" matrix parameters 0 type elementType fields ${_matrix} type kind) + require_json("${_compute_json}" 4 parameters 0 type elementType fields ${_matrix} type rowCount) + require_json("${_compute_json}" 4 parameters 0 type elementType fields ${_matrix} type columnCount) +endforeach() +require_json("${_compute_json}" byteAddressBuffer parameters 6 type baseShape) +# The native counter struct consists entirely of 32-bit words. +foreach(_field RANGE 0 7) + math(EXPR _offset "${_field} * 4") + require_json("${_compute_json}" "${_offset}" parameters 7 type resultType fields ${_field} binding offset) +endforeach() +foreach(_sort keys histogram scan_blocks scan_add scatter verify) + check_compute("gaussian.sort-${_sort}" "gaussian_sort_${_sort}" 256 + gaussian_sort_constants:4 gaussian_sort_source:0 gaussian_sort_destination:1 + gaussian_sort_scan:2 gaussian_sort_prepared_records:3) + require_json("${_compute_json}" 48 parameters 0 type elementVarLayout binding size) + check_fields("${_compute_json}" 0 + "element_count;block_count;digit_shift;digit_word;scan_offset;scan_count;scan_sums_offset;candidate_base;prepared_base;visible_count_offset;count_word;flags" + "0;4;8;12;16;20;24;28;32;36;40;44") + require_json("${_compute_json}" byteAddressBuffer parameters 4 type baseShape) + foreach(_parameter 1 2) + foreach(_field RANGE 0 2) + math(EXPR _offset "${_field} * 4") + require_json("${_compute_json}" "${_offset}" parameters ${_parameter} type resultType fields ${_field} binding offset) + require_json("${_compute_json}" 4 parameters ${_parameter} type resultType fields ${_field} binding size) + endforeach() + endforeach() +endforeach()