diff --git a/CMakeLists.txt b/CMakeLists.txt index b55d145..9d05def 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -9,8 +9,6 @@ set(OPENSPLAT_BUILD_VISUALIZER OFF CACHE BOOL "Build visualizer application") set(OPENSPLAT_USE_FAST_MATH OFF CACHE BOOL "Enable fast math optimizations for GPU kernels (-use_fast_math / -ffast-math)") set(OPENSPLAT_USE_PCH ON CACHE BOOL "Use precompiled headers to speed up compilation") -set(FETCH_DEPENDENCIES ON CACHE BOOL "Fetch additional dependencies from the Internet during configuration") - set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}) # Read version @@ -53,41 +51,59 @@ include(FetchContent) set(NANOFLANN_BUILD_EXAMPLES OFF) set(NANOFLANN_BUILD_TESTS OFF) -set(SPZ_BUILD_PYTHON_BINDINGS OFF) -set(SPZ_BUILD_TOOLS OFF) -set(SPZ_BUILD_EXTENSIONS OFF) -set(SPZ_BUILD_WASM OFF) -set(SPZ_SHOULD_INSTALL OFF) +set(ZLIB_BUILD_TESTING OFF) +set(ZLIB_BUILD_SHARED OFF) +set(ZLIB_INSTALL OFF) -if(FETCH_DEPENDENCIES) +find_package(nlohmann_json QUIET) +if(NOT nlohmann_json_FOUND) + message(STATUS "nlohmann_json not found, fetching it") FetchContent_Declare(nlohmann_json URL https://github.com/nlohmann/json/archive/refs/tags/v3.11.3.zip ) + FetchContent_MakeAvailable(nlohmann_json) +endif() + +find_package(nanoflann QUIET) +if(NOT nanoflann_FOUND) + message(STATUS "nanoflann not found, fetching it") FetchContent_Declare(nanoflann URL https://github.com/jlblancoc/nanoflann/archive/refs/tags/v1.5.5.zip ) + FetchContent_MakeAvailable(nanoflann) +endif() + +find_package(cxxopts QUIET) +if(NOT cxxopts_FOUND) + message(STATUS "cxxopts not found, fetching it") FetchContent_Declare(cxxopts URL https://github.com/jarro2783/cxxopts/archive/refs/tags/v3.2.0.zip ) - FetchContent_Declare(spz - URL https://github.com/nianticlabs/spz/archive/affd0ecea7fbb4c265ee119475af7ee5b2997482.zip - ) - FetchContent_MakeAvailable(nlohmann_json nanoflann cxxopts spz) + FetchContent_MakeAvailable(cxxopts) +endif() + +find_package(ZLIB QUIET) +if(ZLIB_FOUND) + set(ZLIB_LIB ZLIB::ZLIB) else() - find_package(nlohmann_json) - find_package(nanoflann) - find_package(cxxopts) - find_package(spz) + message(STATUS "ZLIB not found, fetching it") + FetchContent_Declare(zlib + URL https://github.com/pierotofy/OpenSplat/releases/download/v1.1.4/zlib-1.3.2.tar.gz + ) + FetchContent_MakeAvailable(zlib) + set(ZLIB_LIB zlibstatic) endif() +add_subdirectory(vendor/spz) + if((GPU_RUNTIME STREQUAL "CUDA") OR (GPU_RUNTIME STREQUAL "HIP")) - if(FETCH_DEPENDENCIES) + find_package(glm QUIET) + if(NOT glm_FOUND) + message(STATUS "glm not found, fetching it") FetchContent_Declare(glm URL https://github.com/g-truc/glm/archive/refs/tags/1.0.1.zip ) FetchContent_MakeAvailable(glm) - else() - find_package(glm) endif() endif() @@ -265,7 +281,7 @@ target_include_directories(gsplat_cpu PRIVATE ${TORCH_INCLUDE_DIRS}) set(OPENSPLAT_SRC_FILES opensplat.cpp point_io.cpp nerfstudio.cpp model.cpp kdtree_tensor.cpp spherical_harmonics.cpp cv_utils.cpp utils.cpp project_gaussians.cpp rasterize_gaussians.cpp ssim.cpp optim_scheduler.cpp colmap.cpp opensfm.cpp openmvg.cpp input_data.cpp -tensor_math.cpp) +tensor_math.cpp rad.cpp) if (OPENSPLAT_BUILD_VISUALIZER) if (Pangolin_FOUND) @@ -302,6 +318,7 @@ target_link_libraries(opensplat PRIVATE nlohmann_json::nlohmann_json cxxopts::cxxopts nanoflann::nanoflann + ${ZLIB_LIB} spz::spz ) if (NOT WIN32) diff --git a/README.md b/README.md index e0fb674..46ed690 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ A free and open source implementation of 3D [gaussian splatting](https://www.you -OpenSplat takes camera poses + sparse points in [COLMAP](https://colmap.github.io/), [OpenSfM](https://github.com/mapillary/OpenSfM), [ODX](https://github.com/WebODM/ODX), [OpenMVG](https://github.com/OpenMVG/OpenMVG) or [nerfstudio](https://docs.nerf.studio/quickstart/custom_dataset.html) project format and computes a [scene file](https://drive.google.com/file/d/12lmvVWpFlFPL6nxl2e2d-4u4a31RCSKT/view?usp=sharing) (.ply, .splat, or .spz) that can be later imported for [viewing](https://antimatter15.com/splat/?url=https://splat.uav4geo.com/banana.splat), editing and rendering in other [software](https://github.com/MrNeRF/awesome-3D-gaussian-splatting?tab=readme-ov-file#open-source-implementations). +OpenSplat takes camera poses + sparse points in [COLMAP](https://colmap.github.io/), [OpenSfM](https://github.com/mapillary/OpenSfM), [ODX](https://github.com/WebODM/ODX), [OpenMVG](https://github.com/OpenMVG/OpenMVG) or [nerfstudio](https://docs.nerf.studio/quickstart/custom_dataset.html) project format and computes a [scene file](https://drive.google.com/file/d/12lmvVWpFlFPL6nxl2e2d-4u4a31RCSKT/view?usp=sharing) (.ply, .splat, .spz, or .rad) that can be later imported for [viewing](https://antimatter15.com/splat/?url=https://splat.uav4geo.com/banana.splat), editing and rendering in other [software](https://github.com/MrNeRF/awesome-3D-gaussian-splatting?tab=readme-ov-file#open-source-implementations). Graphics card recommended, but not required! OpenSplat runs the fastest on NVIDIA, AMD and Apple (Metal) GPUs, but can also run entirely on the CPU (~100x slower). diff --git a/model.cpp b/model.cpp index ca880f9..71b44ab 100644 --- a/model.cpp +++ b/model.cpp @@ -1,4 +1,5 @@ #include +#include #include "model.hpp" #include "constants.hpp" #include "splat-types.h" @@ -8,7 +9,7 @@ #include "tensor_math.hpp" #include "gsplat.hpp" #include "utils.hpp" -#include +#include "rad.hpp" #ifdef USE_MPS #include @@ -510,6 +511,14 @@ void Model::save(const std::string &filename, int step){ savePly(filename, step); std::cout << "Wrote " << filename << std::endl; } + else if (extension == ".rad") { + if (saveRad(filename)) { + std::cout << "Wrote " << filename << std::endl; + } + else { + std::cerr << "Failed to write " << filename << ", aborting save." << std::endl; + } + } else { bool success = saveSpz(filename); if (success) { @@ -621,14 +630,14 @@ bool Model::saveSpz(const std::string &filename){ torch::Tensor meansCpu = keepCrs ? (means.cpu() / scale) + translation : means.cpu(); - torch::Tensor scalesCpu = keepCrs ? (scales.cpu() / scale) : scales.cpu(); - + torch::Tensor scalesCpu = keepCrs ? torch::log(torch::exp(scales.cpu()) / scale) : scales.cpu(); + torch::Tensor meansFlat = meansCpu.flatten(); torch::Tensor scalesFlat = scalesCpu.flatten(); torch::Tensor colorsFlat = featuresDc.cpu().flatten(); // raw DC coefficients torch::Tensor opacFlat = opacities.flatten().cpu(); - torch::Tensor quatsFlat = quats.flatten().cpu(); - torch::Tensor shRestFlat = featuresRest.cpu().transpose(1, 2).flatten(); + torch::Tensor quatsFlat = torch::roll(quats.cpu(), -1, 1).flatten(); + torch::Tensor shRestFlat = featuresRest.cpu().flatten(); spz::GaussianCloud gaussians; gaussians.numPoints = meansCpu.size(0); @@ -641,16 +650,41 @@ bool Model::saveSpz(const std::string &filename){ gaussians.colors = tensor_to_vector(colorsFlat); gaussians.sh = tensor_to_vector(shRestFlat); - auto options = spz::PackOptions{ - .version = 3, // V4 available but not handled by many viewers - .from = spz::CoordinateSystem::RUB, - .sh1Bits = 6, - .shRestBits = 5 - }; + spz::PackOptions options; + options.version = 3; // V4 available but not handled by many viewers + options.from = spz::CoordinateSystem::RUB; + options.sh1Bits = 6; + options.shRestBits = 5; bool success = spz::saveSpz(gaussians, options, filename); return success; } +bool Model::saveRad(const std::string &filename){ + size_t numPoints = means.size(0); + + // Same value preparation as savePly; rad.cpp expects exactly the values + // a saved PLY would contain + torch::Tensor meansCpu = (keepCrs ? (means.cpu() / scale) + translation : means.cpu()).contiguous(); + torch::Tensor featuresDcCpu = featuresDc.cpu().contiguous(); + // featuresRest [N, K, 3] coefficient-major matches rad's SH layout; no transpose + torch::Tensor featuresRestCpu = featuresRest.cpu().contiguous(); + torch::Tensor opacitiesCpu = opacities.cpu().contiguous(); + torch::Tensor scalesCpu = (keepCrs ? torch::log(torch::exp(scales.cpu()) / scale) : scales.cpu()).contiguous(); + torch::Tensor quatsCpu = quats.cpu().contiguous(); + + rad::SplatData data; + data.numPoints = numPoints; + data.numRestCoeffs = featuresRest.size(1); + data.means = tensor_to_vector(meansCpu); + data.featuresDc = tensor_to_vector(featuresDcCpu); + data.featuresRest = tensor_to_vector(featuresRestCpu); + data.opacities = tensor_to_vector(opacitiesCpu); + data.scales = tensor_to_vector(scalesCpu); + data.quats = tensor_to_vector(quatsCpu); + + return rad::saveRad(filename, data); +} + void Model::saveDebugPly(const std::string &filename, int step){ // A standard PLY std::ofstream o(filename, std::ios::binary); diff --git a/model.hpp b/model.hpp index bdec258..ecb33e2 100644 --- a/model.hpp +++ b/model.hpp @@ -73,6 +73,7 @@ struct Model{ void savePly(const std::string &filename, int step); void saveSplat(const std::string &filename); bool saveSpz(const std::string &filename); + bool saveRad(const std::string &filename); void saveDebugPly(const std::string &filename, int step); int loadPly(const std::string &filename); torch::Tensor mainLoss(torch::Tensor &rgb, torch::Tensor >, float ssimWeight); diff --git a/opensplat.cpp b/opensplat.cpp index 197d101..bcc505a 100644 --- a/opensplat.cpp +++ b/opensplat.cpp @@ -160,7 +160,7 @@ int main(int argc, char *argv[]){ if (step % displayStep == 0) { const float percentage = static_cast(step) / numIters; - std::cout << "Step " << step << ": " << mainLoss.item() << " (" << floor(percentage * 100) << "%)" << std::endl; + std::cout << "Step " << step << ": " << mainLoss.item() << " [" << floor(percentage * 100) << "%]" << std::endl; } model.optimizersStep(); diff --git a/rad.cpp b/rad.cpp new file mode 100644 index 0000000..8a6a7a6 --- /dev/null +++ b/rad.cpp @@ -0,0 +1,2314 @@ +// Writer for Spark's .RAD LoD gaussian splat format. +// +// C++ port of the Rust build-lod tool (https://github.com/sparkjsdev/spark/tree/main/rust/build-lod) +// +// Splat attributes are stored as IEEE 754 half-precision values during +// processing; this quantization is part of the format design and shapes the +// LoD construction. + +#include "rad.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace rad { +namespace { + +// insertion-ordered to control metadata key order +using json = nlohmann::ordered_json; + +// IEEE 754 half-precision conversions (round-to-nearest-even) + +inline uint32_t f32Bits(float f){ + uint32_t x; + std::memcpy(&x, &f, 4); + return x; +} + +inline float f32FromBits(uint32_t x){ + float f; + std::memcpy(&f, &x, 4); + return f; +} + +// half-2.6.0 f32_to_f16_fallback +uint16_t f16FromF32(float value){ + uint32_t x = f32Bits(value); + + uint32_t sign = x & 0x80000000u; + uint32_t exp = x & 0x7F800000u; + uint32_t man = x & 0x007FFFFFu; + + if (exp == 0x7F800000u){ + uint32_t nanBit = (man == 0) ? 0 : 0x0200u; + return static_cast((sign >> 16) | 0x7C00u | nanBit | (man >> 13)); + } + + uint32_t halfSign = sign >> 16; + int32_t unbiasedExp = static_cast(exp >> 23) - 127; + int32_t halfExp = unbiasedExp + 15; + + if (halfExp >= 0x1F){ + return static_cast(halfSign | 0x7C00u); + } + + if (halfExp <= 0){ + if (14 - halfExp > 24){ + return static_cast(halfSign); + } + uint32_t man2 = man | 0x00800000u; + uint32_t halfMan = man2 >> (14 - halfExp); + uint32_t roundBit = 1u << (13 - halfExp); + if ((man2 & roundBit) != 0 && (man2 & (3 * roundBit - 1)) != 0){ + halfMan += 1; + } + return static_cast(halfSign | halfMan); + } + + uint32_t halfExpBits = static_cast(halfExp) << 10; + uint32_t halfMan = man >> 13; + uint32_t roundBit = 0x00001000u; + if ((man & roundBit) != 0 && (man & (3 * roundBit - 1)) != 0){ + return static_cast((halfSign | halfExpBits | halfMan) + 1); + }else{ + return static_cast(halfSign | halfExpBits | halfMan); + } +} + +// half-2.6.0 f16_to_f32_fallback +float f16ToF32(uint16_t i){ + if ((i & 0x7FFFu) == 0){ + return f32FromBits(static_cast(i) << 16); + } + + uint32_t halfSign = i & 0x8000u; + uint32_t halfExp = i & 0x7C00u; + uint32_t halfMan = i & 0x03FFu; + + if (halfExp == 0x7C00u){ + if (halfMan == 0){ + return f32FromBits((halfSign << 16) | 0x7F800000u); + }else{ + return f32FromBits((halfSign << 16) | 0x7FC00000u | (halfMan << 13)); + } + } + + uint32_t sign = halfSign << 16; + int32_t unbiasedExp = (static_cast(halfExp) >> 10) - 15; + + if (halfExp == 0){ + // leading_zeros_u16(halfMan) - 6; halfMan != 0 here + int32_t lz = 0; + uint16_t m = static_cast(halfMan); + while (!(m & 0x8000u)){ lz++; m <<= 1; } + int32_t e = lz - 6; + uint32_t exp = static_cast(127 - 15 - e) << 23; + uint32_t man = (halfMan << (14 + e)) & 0x7FFFFFu; + return f32FromBits(sign | exp | man); + } + + uint32_t exp = static_cast(unbiasedExp + 127) << 23; + uint32_t man = (halfMan & 0x03FFu) << 13; + return f32FromBits(sign | exp | man); +} + +// A stored f16 value (bit pattern). Mirrors half::f16. +struct F16 { + uint16_t bits = 0; + + F16() = default; + static F16 fromF32(float v){ F16 h; h.bits = f16FromF32(v); return h; } + static F16 fromBits(uint16_t b){ F16 h; h.bits = b; return h; } + float toF32() const { return f16ToF32(bits); } + bool isNan() const { return (bits & 0x7C00u) == 0x7C00u && (bits & 0x03FFu) != 0; } +}; + +// half-2.6.0 f16::max: if other > self && !other.is_nan() { other } else { self } +// (non-NaN f16 ordering matches f32 ordering of the converted values, incl. -0 == 0) +inline F16 f16Max(F16 self, F16 other){ + if (!self.isNan() && !other.isNan() && other.toF32() > self.toF32()) return other; + return self; +} + +// Rust semantics helpers + +// Rust `as` cast f32 -> integer: truncate toward zero, saturate, NaN -> 0. +inline int64_t rustCastI64(float v){ + if (std::isnan(v)) return 0; + if (v <= -9223372036854775808.0f) return INT64_MIN; + if (v >= 9223372036854775807.0f) return INT64_MAX; + return static_cast(v); +} + +inline int16_t rustCastI16(float v){ + if (std::isnan(v)) return 0; + if (v <= -32768.0f) return INT16_MIN; + if (v >= 32767.0f) return INT16_MAX; + return static_cast(v); +} + +inline uint8_t rustCastU8(float v){ + if (std::isnan(v)) return 0; + if (v <= 0.0f) return 0; + if (v >= 255.0f) return 255; + return static_cast(v); +} + +inline int8_t rustCastI8(float v){ + if (std::isnan(v)) return 0; + if (v <= -128.0f) return INT8_MIN; + if (v >= 127.0f) return INT8_MAX; + return static_cast(v); +} + +inline uint32_t rustCastU32(uint64_t v){ return static_cast(v); } + +// Vector / quaternion / matrix math (scalar port of glam) + +struct Vec3 { + float x = 0.0f, y = 0.0f, z = 0.0f; + + Vec3() = default; + Vec3(float x_, float y_, float z_) : x(x_), y(y_), z(z_) {} + static Vec3 splat(float v){ return Vec3(v, v, v); } + + float operator[](int i) const { return i == 0 ? x : (i == 1 ? y : z); } + float &operator[](int i){ return i == 0 ? x : (i == 1 ? y : z); } + + Vec3 operator+(const Vec3 &o) const { return Vec3(x + o.x, y + o.y, z + o.z); } + Vec3 operator-(const Vec3 &o) const { return Vec3(x - o.x, y - o.y, z - o.z); } + Vec3 operator*(const Vec3 &o) const { return Vec3(x * o.x, y * o.y, z * o.z); } + Vec3 operator*(float s) const { return Vec3(x * s, y * s, z * s); } + Vec3 operator/(float s) const { return Vec3(x / s, y / s, z / s); } + + // glam mul_add without the fma target feature calls f32::mul_add, which is + // a correctly-rounded fused multiply-add (fmaf), NOT mul-then-add. + Vec3 mulAdd(const Vec3 &a, const Vec3 &b) const { + return Vec3(std::fmaf(x, a.x, b.x), std::fmaf(y, a.y, b.y), std::fmaf(z, a.z, b.z)); + } + + // glam sse2 max/min use _mm_max_ps/_mm_min_ps: a > b ? a : b per element + Vec3 max(const Vec3 &o) const { + return Vec3(x > o.x ? x : o.x, y > o.y ? y : o.y, z > o.z ? z : o.z); + } + + Vec3 floorv() const { return Vec3(std::floor(x), std::floor(y), std::floor(z)); } + + float maxElement() const { + // glam sse2 max_element: max(max(x, y), z) via _mm_max_ss (a > b ? a : b) + float m = x > y ? x : y; + return m > z ? m : z; + } + + // glam sse2 dot3: (x*x' + y*y') + z*z' + float dot(const Vec3 &o) const { return (x * o.x + y * o.y) + z * o.z; } + float lengthSquared() const { return dot(*this); } + float length() const { return std::sqrt(dot(*this)); } + float distance(const Vec3 &o) const { return (*this - o).length(); } + + bool isFinite() const { return std::isfinite(x) && std::isfinite(y) && std::isfinite(z); } +}; + +struct I64Vec3 { + int64_t x = 0, y = 0, z = 0; + + I64Vec3() = default; + I64Vec3(int64_t x_, int64_t y_, int64_t z_) : x(x_), y(y_), z(z_) {} + + int64_t operator[](int i) const { return i == 0 ? x : (i == 1 ? y : z); } + bool operator==(const I64Vec3 &o) const { return x == o.x && y == o.y && z == o.z; } + + I64Vec3 minv(const I64Vec3 &o) const { + return I64Vec3(std::min(x, o.x), std::min(y, o.y), std::min(z, o.z)); + } + I64Vec3 maxv(const I64Vec3 &o) const { + return I64Vec3(std::max(x, o.x), std::max(y, o.y), std::max(z, o.z)); + } +}; + +struct I64Vec3Hash { + size_t operator()(const I64Vec3 &v) const { + // Iteration order of the cell map is never observed, so any hash works. + uint64_t h = 0x9E3779B97F4A7C15ull; + for (uint64_t c : {static_cast(v.x), static_cast(v.y), static_cast(v.z)}){ + h ^= c + 0x9E3779B97F4A7C15ull + (h << 6) + (h >> 2); + } + return static_cast(h); + } +}; + +struct Quat { + float x = 0.0f, y = 0.0f, z = 0.0f, w = 1.0f; + + Quat() = default; + Quat(float x_, float y_, float z_, float w_) : x(x_), y(y_), z(z_), w(w_) {} + + // glam sse2 dot4: (x*x' + z*z') + (y*y' + w*w') + float dot(const Quat &o) const { return (x * o.x + z * o.z) + (y * o.y + w * o.w); } + float length() const { return std::sqrt(dot(*this)); } + + bool isFinite() const { + return std::isfinite(x) && std::isfinite(y) && std::isfinite(z) && std::isfinite(w); + } +}; + +// Column-major 3x3 matrix: m[c] is column c. col(c)[r] = m[c][r]. +struct Mat3 { + float m[3][3]; + + static Mat3 identity(){ + Mat3 r{}; + r.m[0][0] = 1.0f; r.m[1][1] = 1.0f; r.m[2][2] = 1.0f; + return r; + } + static Mat3 fromCols(const Vec3 &c0, const Vec3 &c1, const Vec3 &c2){ + Mat3 r; + for (int i = 0; i < 3; i++){ r.m[0][i] = c0[i]; r.m[1][i] = c1[i]; r.m[2][i] = c2[i]; } + return r; + } + Vec3 col(int c) const { return Vec3(m[c][0], m[c][1], m[c][2]); } + + // glam Mat3A::from_quat + static Mat3 fromQuat(const Quat &q){ + float x2 = q.x + q.x; + float y2 = q.y + q.y; + float z2 = q.z + q.z; + float xx = q.x * x2; + float xy = q.x * y2; + float xz = q.x * z2; + float yy = q.y * y2; + float yz = q.y * z2; + float zz = q.z * z2; + float wx = q.w * x2; + float wy = q.w * y2; + float wz = q.w * z2; + return fromCols( + Vec3(1.0f - (yy + zz), xy + wz, xz - wy), + Vec3(xy - wz, 1.0f - (xx + zz), yz + wx), + Vec3(xz + wy, yz - wx, 1.0f - (xx + yy))); + } +}; + +// glam Quat::from_rotation_axes (Quat::from_mat3a passes the matrix columns) +Quat quatFromMat3(const Mat3 &mat){ + float m00 = mat.m[0][0], m01 = mat.m[0][1], m02 = mat.m[0][2]; + float m10 = mat.m[1][0], m11 = mat.m[1][1], m12 = mat.m[1][2]; + float m20 = mat.m[2][0], m21 = mat.m[2][1], m22 = mat.m[2][2]; + if (m22 <= 0.0f){ + float dif10 = m11 - m00; + float omm22 = 1.0f - m22; + if (dif10 <= 0.0f){ + float fourXsq = omm22 - dif10; + float inv4x = 0.5f / std::sqrt(fourXsq); + return Quat(fourXsq * inv4x, (m01 + m10) * inv4x, (m02 + m20) * inv4x, (m12 - m21) * inv4x); + }else{ + float fourYsq = omm22 + dif10; + float inv4y = 0.5f / std::sqrt(fourYsq); + return Quat((m01 + m10) * inv4y, fourYsq * inv4y, (m12 + m21) * inv4y, (m20 - m02) * inv4y); + } + }else{ + float sum10 = m11 + m00; + float opm22 = 1.0f + m22; + if (sum10 <= 0.0f){ + float fourZsq = opm22 - sum10; + float inv4z = 0.5f / std::sqrt(fourZsq); + return Quat((m02 + m20) * inv4z, (m12 + m21) * inv4z, fourZsq * inv4z, (m01 - m10) * inv4z); + }else{ + float fourWsq = opm22 + sum10; + float inv4w = 0.5f / std::sqrt(fourWsq); + return Quat((m12 - m21) * inv4w, (m20 - m02) * inv4w, (m01 - m10) * inv4w, fourWsq * inv4w); + } + } +} + +// SymMat3 (port of spark-lib/src/symmat3.rs) +// Storage [xx, yy, zz, xy] + [xz, yz] + +struct SymMat3 { + float v0[4] = {0, 0, 0, 0}; // xx, yy, zz, xy + float v1[2] = {0, 0}; // xz, yz + + static SymMat3 make(float xx, float yy, float zz, float xy, float xz, float yz){ + SymMat3 r; + r.v0[0] = xx; r.v0[1] = yy; r.v0[2] = zz; r.v0[3] = xy; + r.v1[0] = xz; r.v1[1] = yz; + return r; + } + + float xx() const { return v0[0]; } + float yy() const { return v0[1]; } + float zz() const { return v0[2]; } + float xy() const { return v0[3]; } + float xz() const { return v1[0]; } + float yz() const { return v1[1]; } + + static SymMat3 newScaleQuaternion(const Vec3 &scale, const Quat &quat){ + Mat3 rot = Mat3::fromQuat(quat); + Vec3 sx = rot.col(0) * scale.x; + Vec3 sy = rot.col(1) * scale.y; + Vec3 sz = rot.col(2) * scale.z; + float xx = sx.x * sx.x + sy.x * sy.x + sz.x * sz.x; + float yy = sx.y * sx.y + sy.y * sy.y + sz.y * sz.y; + float zz = sx.z * sx.z + sy.z * sy.z + sz.z * sz.z; + float xy = sx.x * sx.y + sy.x * sy.y + sz.x * sz.y; + float xz = sx.x * sx.z + sy.x * sy.z + sz.x * sz.z; + float yz = sx.y * sx.z + sy.y * sy.z + sz.y * sz.z; + return make(xx, yy, zz, xy, xz, yz); + } + + // Vec4/Vec2 mul_add = correctly-rounded fused multiply-add (see Vec3::mulAdd) + void addWeighted(const SymMat3 &other, float weight){ + for (int i = 0; i < 4; i++) v0[i] = std::fmaf(other.v0[i], weight, v0[i]); + for (int i = 0; i < 2; i++) v1[i] = std::fmaf(other.v1[i], weight, v1[i]); + } + + static SymMat3 newAverage(const SymMat3 &a, const SymMat3 &b){ + SymMat3 r; + for (int i = 0; i < 4; i++) r.v0[i] = std::fmaf(a.v0[i], 0.5f, b.v0[i] * 0.5f); + for (int i = 0; i < 2; i++) r.v1[i] = std::fmaf(a.v1[i], 0.5f, b.v1[i] * 0.5f); + return r; + } + + float determinant() const { + float m00 = v0[0], m11 = v0[1], m22 = v0[2]; + float m01 = v0[3], m02 = v1[0], m12 = v1[1]; + return m00 * (m11 * m22 - m12 * m12) - + m01 * (m01 * m22 - m12 * m02) + + m02 * (m01 * m12 - m11 * m02); + } + + bool inverse(SymMat3 &out) const { + float m00 = v0[0], m11 = v0[1], m22 = v0[2]; + float m01 = v0[3], m02 = v1[0], m12 = v1[1]; + + float det = determinant(); + float diagMax = std::fmaxf(std::fmaxf(std::fabs(v0[0]), std::fabs(v0[1])), std::fabs(v0[2])); + float relTol = 1e-9f * std::fmaxf(diagMax * diagMax * diagMax, 1e-30f); + if (std::fabs(det) < relTol){ + return false; + } + float invDet = 1.0f / det; + + out = make((m11 * m22 - m12 * m12) * invDet, + (m00 * m22 - m02 * m02) * invDet, + (m00 * m11 - m01 * m01) * invDet, + (m02 * m12 - m01 * m22) * invDet, + (m01 * m12 - m02 * m11) * invDet, + (m01 * m02 - m00 * m12) * invDet); + return true; + } + + void eigens(float vals[3], Vec3 vecs[3]) const { + const int MAX_ITERS = 32; + float eps; + { + float s = std::fabs(v0[0]) + std::fabs(v0[1]) + std::fabs(v0[2]); + eps = 1e-6f * std::fmaxf(s, 1.0f); + } + + Mat3 current = Mat3::fromCols( + Vec3(v0[0], v0[3], v1[0]), + Vec3(v0[3], v0[1], v1[1]), + Vec3(v1[0], v1[1], v0[2])); + Mat3 eigs = Mat3::identity(); + + auto offDiagNorm2 = [](const Mat3 &a) -> float { + float a01 = a.m[0][1]; + float a02 = a.m[0][2]; + float a12 = a.m[1][2]; + return a01 * a01 + a02 * a02 + a12 * a12; + }; + + int k = 0; + while (k < MAX_ITERS && offDiagNorm2(current) > (eps * eps)){ + int p = 0, q = 1; + float maxVal = std::fabs(current.m[0][1]); + const struct { int i, j; } cand[2] = {{0, 2}, {1, 2}}; + const float candVal[2] = {std::fabs(current.m[0][2]), std::fabs(current.m[1][2])}; + for (int c = 0; c < 2; c++){ + if (candVal[c] > maxVal){ + maxVal = candVal[c]; + p = cand[c].i; + q = cand[c].j; + } + } + + float apq = current.m[p][q]; + if (std::fabs(apq) > eps){ + float app = current.m[p][p]; + float aqq = current.m[q][q]; + float tau = aqq - app; + float phi = 0.5f * std::atan2(2.0f * apq, tau); + float c = std::cos(phi); + float s = std::sin(phi); + + for (int r = 0; r < 3; r++){ + float arp = current.m[r][p]; + float arq = current.m[r][q]; + current.m[r][p] = c * arp - s * arq; + current.m[r][q] = s * arp + c * arq; + } + for (int r = 0; r < 3; r++){ + float apr = current.m[p][r]; + float aqr = current.m[q][r]; + current.m[p][r] = c * apr - s * aqr; + current.m[q][r] = s * apr + c * aqr; + } + current.m[p][q] = 0.0f; + current.m[q][p] = 0.0f; + + for (int r = 0; r < 3; r++){ + float vrp = eigs.m[r][p]; + float vrq = eigs.m[r][q]; + eigs.m[r][p] = c * vrp - s * vrq; + eigs.m[r][q] = s * vrp + c * vrq; + } + } + k += 1; + } + + float rawVals[3] = {current.m[0][0], current.m[1][1], current.m[2][2]}; + Vec3 rawVecs[3] = { + Vec3(eigs.m[0][0], eigs.m[1][0], eigs.m[2][0]), + Vec3(eigs.m[0][1], eigs.m[1][1], eigs.m[2][1]), + Vec3(eigs.m[0][2], eigs.m[1][2], eigs.m[2][2]), + }; + + for (int j = 0; j < 3; j++){ + float n = std::sqrt(rawVecs[j].x * rawVecs[j].x + rawVecs[j].y * rawVecs[j].y + rawVecs[j].z * rawVecs[j].z); + if (n > 0.0f){ + rawVecs[j].x /= n; + rawVecs[j].y /= n; + rawVecs[j].z /= n; + } + } + + // Stable sort of [0,1,2] by descending eigenvalue + int idx[3] = {0, 1, 2}; + std::stable_sort(idx, idx + 3, [&](int a, int b){ + return rawVals[b] < rawVals[a]; + }); + for (int j = 0; j < 3; j++){ + vals[j] = rawVals[idx[j]]; + vecs[j] = rawVecs[idx[j]]; + } + } + + void positiveEigens(float vals[3], Vec3 vecs[3]) const { + eigens(vals, vecs); + float det = + vecs[0][0] * (vecs[1][1] * vecs[2][2] - vecs[1][2] * vecs[2][1]) - + vecs[0][1] * (vecs[1][0] * vecs[2][2] - vecs[1][2] * vecs[2][0]) + + vecs[0][2] * (vecs[1][0] * vecs[2][1] - vecs[1][1] * vecs[2][0]); + if (det < 0.0f){ + vecs[2] = Vec3(-vecs[2][0], -vecs[2][1], -vecs[2][2]); + } + } +}; + +// Max-heap keyed by (float, index), compared lexicographically. Keys are +// unique, so pop order is fully deterministic; the internal array order only +// affects tie-breaking in neighbor scans. + +struct HeapKey { + float key; + size_t index; + + bool operator<(const HeapKey &o) const { + if (key != o.key) return key < o.key; + return index < o.index; + } +}; + +struct MaxHeap { + std::vector data; + + size_t len() const { return data.size(); } + bool isEmpty() const { return data.empty(); } + + void push(const HeapKey &item){ + data.push_back(item); + std::push_heap(data.begin(), data.end()); + } + + bool pop(HeapKey &out){ + if (data.empty()) return false; + std::pop_heap(data.begin(), data.end()); + out = data.back(); + data.pop_back(); + return true; + } + + void extend(const std::vector &items){ + data.insert(data.end(), items.begin(), items.end()); + std::make_heap(data.begin(), data.end()); + } +}; + +// Gsplat / GsplatArray (port of spark-lib/src/gsplat.rs + tsplat.rs) + +// tsplat.rs ellipsoid_area (Knud Thomsen approximation) +float ellipsoidArea(const Vec3 &scales){ + const float P = 1.6075f; + float numerator = std::pow(scales.x * scales.y, P) + std::pow(scales.x * scales.z, P) + + std::pow(scales.y * scales.z, P); + return 4.0f * 3.14159265358979323846264338327950288f * std::pow(numerator / 3.0f, 1.0f / P); +} + +struct Gsplat { + Vec3 center; + F16 opacity; + F16 rgb[3]; + F16 lnScales[3]; + F16 quaternion[4]; // x, y, z, w + + static Gsplat make(const Vec3 ¢er, float opacity, const Vec3 &rgb, const Vec3 &scales, + const Quat &quaternion){ + Gsplat s; + s.center = center; + s.opacity = F16::fromF32(opacity); + s.rgb[0] = F16::fromF32(rgb.x); + s.rgb[1] = F16::fromF32(rgb.y); + s.rgb[2] = F16::fromF32(rgb.z); + s.lnScales[0] = F16::fromF32(std::log(scales.x)); + s.lnScales[1] = F16::fromF32(std::log(scales.y)); + s.lnScales[2] = F16::fromF32(std::log(scales.z)); + s.quaternion[0] = F16::fromF32(quaternion.x); + s.quaternion[1] = F16::fromF32(quaternion.y); + s.quaternion[2] = F16::fromF32(quaternion.z); + s.quaternion[3] = F16::fromF32(quaternion.w); + return s; + } + + Vec3 getCenter() const { return center; } + float getOpacity() const { return opacity.toF32(); } + Vec3 getRgb() const { return Vec3(rgb[0].toF32(), rgb[1].toF32(), rgb[2].toF32()); } + Vec3 getScales() const { + return Vec3(std::exp(lnScales[0].toF32()), std::exp(lnScales[1].toF32()), + std::exp(lnScales[2].toF32())); + } + Quat getQuaternion() const { + return Quat(quaternion[0].toF32(), quaternion[1].toF32(), quaternion[2].toF32(), + quaternion[3].toF32()); + } + // gsplat.rs: max over f16 ln_scales, then exp + float maxScale() const { + return std::exp(f16Max(f16Max(lnScales[0], lnScales[1]), lnScales[2]).toF32()); + } + + void setCenter(const Vec3 &c){ center = c; } + void setOpacity(float v){ opacity = F16::fromF32(v); } + void setRgb(const Vec3 &v){ + rgb[0] = F16::fromF32(v.x); + rgb[1] = F16::fromF32(v.y); + rgb[2] = F16::fromF32(v.z); + } + void setScales(const Vec3 &scales){ + lnScales[0] = F16::fromF32(std::log(scales.x)); + lnScales[1] = F16::fromF32(std::log(scales.y)); + lnScales[2] = F16::fromF32(std::log(scales.z)); + } + void setQuaternion(const Quat &q){ + quaternion[0] = F16::fromF32(q.x); + quaternion[1] = F16::fromF32(q.y); + quaternion[2] = F16::fromF32(q.z); + quaternion[3] = F16::fromF32(q.w); + } + + float area() const { return ellipsoidArea(getScales()); } + + // tsplat.rs lod_opacity + float lodOpacity() const { + float op = getOpacity(); + if (op > 1.0f){ + return std::sqrt(1.0f + 2.71828182845904523536028747135266250f * std::log(op)); + } + return 1.0f; + } + + float featureSize() const { return 2.0f * maxScale() * lodOpacity(); } + + // tsplat.rs grid: (center / step).floor().as_i64vec3() + I64Vec3 grid(float stepSize) const { + Vec3 g = (center / stepSize).floorv(); + return I64Vec3(rustCastI64(g.x), rustCastI64(g.y), rustCastI64(g.z)); + } +}; + +// tsplat.rs bhattacharyya_distance +float bhattacharyyaDistance(const Gsplat &a, const Gsplat &b){ + SymMat3 covA = SymMat3::newScaleQuaternion(a.getScales(), a.getQuaternion()); + SymMat3 covB = SymMat3::newScaleQuaternion(b.getScales(), b.getQuaternion()); + SymMat3 sigma = SymMat3::newAverage(covA, covB); + SymMat3 inv; + if (!sigma.inverse(inv)){ + return 0.0f; + } + + Vec3 delta = b.getCenter() - a.getCenter(); + float quad = inv.xx() * delta.x * delta.x + + inv.yy() * delta.y * delta.y + + inv.zz() * delta.z * delta.z + + 2.0f * inv.xy() * delta.x * delta.y + + 2.0f * inv.xz() * delta.x * delta.z + + 2.0f * inv.yz() * delta.y * delta.z; + float term1 = 0.125f * quad; + + float detSigma = sigma.determinant(); + float detA = covA.determinant(); + float detB = covB.determinant(); + float term2 = 0.5f * std::log(detSigma / std::sqrt(detA * detB)); + + return term1 + term2; +} + +// tsplat.rs similarity_metric +float similarityMetric(const Gsplat &a, const Gsplat &b){ + float spatial = std::exp(-bhattacharyyaDistance(a, b)); + + Vec3 colorDelta = a.getRgb() - b.getRgb(); + float colorDelta2 = colorDelta.lengthSquared(); + + float metric = spatial * std::exp(-colorDelta2); + if (std::isnan(metric)){ + return 0.0f; + } + return metric; +} + +// tsplat.rs compute_swaps +std::vector> computeSwaps(const std::vector &indexMap){ + size_t n = indexMap.size(); + std::vector destOfSrc(n, 0); + for (size_t newI = 0; newI < n; newI++){ + destOfSrc[indexMap[newI]] = newI; + } + + std::vector> swaps; + for (size_t i = 0; i < n; i++){ + while (destOfSrc[i] != i){ + size_t j = destOfSrc[i]; + swaps.push_back({i, j}); + std::swap(destOfSrc[i], destOfSrc[j]); + } + } + return swaps; +} + +template +void applySwaps(std::vector &data, const std::vector> &swaps){ + for (const auto &s : swaps){ + std::swap(data[s.first], data[s.second]); + } +} + +typedef std::array GsplatSH1; // 3 coeffs x rgb +typedef std::array GsplatSH2; // 5 coeffs x rgb +typedef std::array GsplatSH3; // 7 coeffs x rgb + +struct GsplatArray { + size_t maxShDegree = 0; + std::vector splats; + std::vector> children; + std::vector sh1; + std::vector sh2; + std::vector sh3; + + size_t len() const { return splats.size(); } + + void prepareChildren(){ children.resize(len()); } + bool hasChildren() const { return !children.empty(); } + bool hasLodTree() const { return !children.empty(); } + + // gsplat.rs new_merged (step is always 0.0 from bhatt_lod) + size_t newMerged(const size_t *indices, size_t numIndices, float step){ + size_t newIndex = splats.size(); + + std::vector weights(numIndices); + for (size_t i = 0; i < numIndices; i++){ + const Gsplat &splat = splats[indices[i]]; + weights[i] = splat.area() * splat.getOpacity(); + } + float sum = 0.0f; + for (size_t i = 0; i < numIndices; i++) sum += weights[i]; + float totalWeight = std::fmaxf(sum, 1.0e-30f); + for (size_t i = 0; i < numIndices; i++) weights[i] /= totalWeight; + + Vec3 center = Vec3(0, 0, 0); + Vec3 rgb = Vec3(0, 0, 0); + + for (size_t i = 0; i < numIndices; i++){ + const Gsplat &splat = splats[indices[i]]; + float weight = weights[i]; + center = splat.getCenter().mulAdd(Vec3::splat(weight), center); + rgb = splat.getRgb().mulAdd(Vec3::splat(weight), rgb); + } + + SymMat3 totalCov; + float filter2 = (0.5f * step) * (0.5f * step); // powi(2) + + for (size_t i = 0; i < numIndices; i++){ + const Gsplat &splat = splats[indices[i]]; + float weight = weights[i]; + Vec3 delta = splat.getCenter() - center; + SymMat3 cov = SymMat3::newScaleQuaternion(splat.getScales(), splat.getQuaternion()); + float xx = delta.x * delta.x + cov.xx() + filter2; + float yy = delta.y * delta.y + cov.yy() + filter2; + float zz = delta.z * delta.z + cov.zz() + filter2; + float xy = delta.x * delta.y + cov.xy(); + float xz = delta.x * delta.z + cov.xz(); + float yz = delta.y * delta.z + cov.yz(); + totalCov.addWeighted(SymMat3::make(xx, yy, zz, xy, xz, yz), weight); + } + + float vals[3]; + Vec3 vecs[3]; + totalCov.positiveEigens(vals, vecs); + Vec3 scales = Vec3(std::sqrt(std::fmaxf(vals[0], 0.0f)), std::sqrt(std::fmaxf(vals[1], 0.0f)), + std::sqrt(std::fmaxf(vals[2], 0.0f))); + scales = scales.max(Vec3::splat(1.0e-30f)); + + Mat3 basis = Mat3::fromCols(vecs[0], vecs[1], vecs[2]); + Quat quaternion = quatFromMat3(basis); + float opacity = totalWeight / ellipsoidArea(scales); + opacity = std::clamp(opacity, 0.000001f, 1000.0f); + + splats.push_back(Gsplat::make(center, opacity, rgb, scales, quaternion)); + children.push_back(std::vector(indices, indices + numIndices)); + + if (maxShDegree >= 1){ + Vec3 total[3] = {Vec3(0, 0, 0), Vec3(0, 0, 0), Vec3(0, 0, 0)}; + for (size_t i = 0; i < numIndices; i++){ + float weight = weights[i]; + const GsplatSH1 &s = sh1[indices[i]]; + for (int c = 0; c < 3; c++){ + Vec3 v(s[c * 3 + 0].toF32(), s[c * 3 + 1].toF32(), s[c * 3 + 2].toF32()); + total[c] = v.mulAdd(Vec3::splat(weight), total[c]); + } + } + GsplatSH1 out; + for (int c = 0; c < 3; c++){ + out[c * 3 + 0] = F16::fromF32(total[c].x); + out[c * 3 + 1] = F16::fromF32(total[c].y); + out[c * 3 + 2] = F16::fromF32(total[c].z); + } + sh1.push_back(out); + } + + if (maxShDegree >= 2){ + Vec3 total[5] = {Vec3(0, 0, 0), Vec3(0, 0, 0), Vec3(0, 0, 0), Vec3(0, 0, 0), Vec3(0, 0, 0)}; + for (size_t i = 0; i < numIndices; i++){ + float weight = weights[i]; + const GsplatSH2 &s = sh2[indices[i]]; + for (int c = 0; c < 5; c++){ + Vec3 v(s[c * 3 + 0].toF32(), s[c * 3 + 1].toF32(), s[c * 3 + 2].toF32()); + total[c] = v.mulAdd(Vec3::splat(weight), total[c]); + } + } + GsplatSH2 out; + for (int c = 0; c < 5; c++){ + out[c * 3 + 0] = F16::fromF32(total[c].x); + out[c * 3 + 1] = F16::fromF32(total[c].y); + out[c * 3 + 2] = F16::fromF32(total[c].z); + } + sh2.push_back(out); + } + + if (maxShDegree >= 3){ + Vec3 total[7]; + for (size_t i = 0; i < numIndices; i++){ + float weight = weights[i]; + const GsplatSH3 &s = sh3[indices[i]]; + for (int c = 0; c < 7; c++){ + Vec3 v(s[c * 3 + 0].toF32(), s[c * 3 + 1].toF32(), s[c * 3 + 2].toF32()); + total[c] = v.mulAdd(Vec3::splat(weight), total[c]); + } + } + GsplatSH3 out; + for (int c = 0; c < 7; c++){ + out[c * 3 + 0] = F16::fromF32(total[c].x); + out[c * 3 + 1] = F16::fromF32(total[c].y); + out[c * 3 + 2] = F16::fromF32(total[c].z); + } + sh3.push_back(out); + } + + return newIndex; + } + + void setChildren(size_t parent, const std::vector &c){ children[parent] = c; } + + std::vector getChildren(size_t parent) const { return children[parent]; } + + float similarity(size_t a, size_t b) const { return similarityMetric(splats[a], splats[b]); } + + template + void retain(F f){ + std::vector keep(splats.size()); + for (size_t i = 0; i < splats.size(); i++){ + keep[i] = f(splats[i]); + } + retainByMask(splats, keep); + if (!children.empty()) retainByMask(children, keep); + if (!sh1.empty()) retainByMask(sh1, keep); + if (!sh2.empty()) retainByMask(sh2, keep); + if (!sh3.empty()) retainByMask(sh3, keep); + } + + void permute(const std::vector &indexMap){ + assert(indexMap.size() == splats.size()); + auto swaps = computeSwaps(indexMap); + applySwaps(splats, swaps); + if (!children.empty()) applySwaps(children, swaps); + if (!sh1.empty()) applySwaps(sh1, swaps); + if (!sh2.empty()) applySwaps(sh2, swaps); + if (!sh3.empty()) applySwaps(sh3, swaps); + } + + void truncate(size_t count){ + if (splats.size() > count) splats.resize(count); + if (!children.empty() && children.size() > count) children.resize(count); + if (!sh1.empty() && sh1.size() > count) sh1.resize(count); + if (!sh2.empty() && sh2.size() > count) sh2.resize(count); + if (!sh3.empty() && sh3.size() > count) sh3.resize(count); + } + + // tsplat.rs sort_by: stable sort of the index map by key, then permute + void sortByFeatureSize(){ + std::vector indexMap(len()); + for (size_t i = 0; i < indexMap.size(); i++) indexMap[i] = i; + std::vector keys(len()); + for (size_t i = 0; i < keys.size(); i++) keys[i] = splats[i].featureSize(); + // OrderedFloat ordering; keys are finite here so plain < works, + // and stable_sort preserves equal-key order like Rust's stable sort. + std::stable_sort(indexMap.begin(), indexMap.end(), + [&](size_t a, size_t b){ return keys[a] < keys[b]; }); + permute(indexMap); + } + + // tsplat.rs encode_lod_opacity + void encodeLodOpacity(){ + for (size_t i = 0; i < len(); i++){ + Gsplat &splat = splats[i]; + if (splat.getOpacity() > 1.0f){ + float d = splat.lodOpacity(); + splat.setOpacity(std::clamp(0.25f * (d - 1.0f) + 1.0f, 1.0f, 2.0f)); + } + } + } + +private: + template + static void retainByMask(std::vector &v, const std::vector &keep){ + size_t out = 0; + for (size_t i = 0; i < v.size(); i++){ + if (keep[i]){ + if (out != i) v[out] = std::move(v[i]); + out++; + } + } + v.resize(out); + } +}; + +// bhatt_lod (port of spark-lib/src/bhatt_lod.rs) + +const float MERGE_BASE = 2.0f; + +void bhattRecurseToOutput(GsplatArray &splats, size_t index, std::vector &toOutput, + float lodBase, float &featureSizeOut, std::vector &childrenOut){ + float featureSize; + { + const Gsplat &splat = splats.splats[index]; + featureSize = splat.area() * splat.getOpacity(); + } + + std::vector children = splats.getChildren(index); + if (children.empty()){ + featureSizeOut = featureSize; + childrenOut.assign(1, index); + return; + } + + std::vector newChildren; + float maxChildFeatureSize = -std::numeric_limits::infinity(); + + for (size_t child : children){ + float childFeatureSize; + std::vector childChildren; + bhattRecurseToOutput(splats, child, toOutput, lodBase, childFeatureSize, childChildren); + maxChildFeatureSize = std::fmaxf(maxChildFeatureSize, childFeatureSize); + newChildren.insert(newChildren.end(), childChildren.begin(), childChildren.end()); + } + + if (featureSize >= (maxChildFeatureSize * lodBase)){ + toOutput[index] = true; + } + + if (toOutput[index]){ + assert(newChildren.size() <= 65535); + splats.setChildren(index, newChildren); + featureSizeOut = featureSize; + childrenOut.assign(1, index); + }else{ + splats.setChildren(index, std::vector()); + featureSizeOut = maxChildFeatureSize; + childrenOut = std::move(newChildren); + } +} + +void bhattRecurseIndices(GsplatArray &splats, size_t index, std::vector &indices, + float limitSize, std::vector &frontier){ + if (splats.splats[index].featureSize() < limitSize){ + frontier.push_back(index); + return; + } + + std::vector children = splats.getChildren(index); + if (children.empty()){ + return; + } + + std::vector newChildren(children.size()); + for (size_t i = 0; i < children.size(); i++) newChildren[i] = indices.size() + i; + splats.setChildren(index, newChildren); + + std::sort(children.begin(), children.end()); + for (size_t child : children){ + indices.push_back(child); + } + + for (size_t child : children){ + bhattRecurseIndices(splats, child, indices, limitSize, frontier); + } +} + +void bhattComputeLodTree(GsplatArray &splats, float lodBase){ + size_t initialLen = splats.len(); + if (initialLen == 0){ + return; + } + + splats.sortByFeatureSize(); + splats.prepareChildren(); + + std::vector isActive(splats.len(), true); + + float minFeatureSize = std::fmaxf(splats.splats[0].featureSize(), 0.000001f); + // Rust: min_feature_size.log(MERGE_BASE).ceil() as i16 + int16_t levelMin = rustCastI16(std::ceil(std::log(minFeatureSize) / std::log(MERGE_BASE))); + + int32_t level = levelMin; + size_t frontier = 0; + MaxHeap active; + std::unordered_map, I64Vec3Hash> cells; + + for (;;){ + float step = std::pow(MERGE_BASE, static_cast(level)); + + size_t frontierStart = frontier; + while (frontier < initialLen){ + if (splats.splats[frontier].featureSize() > step){ + break; + } + frontier += 1; + } + + if (frontier > frontierStart){ + std::vector newSplats; + newSplats.reserve(frontier - frontierStart); + for (size_t i = frontierStart; i < frontier; i++){ + newSplats.push_back(HeapKey{-splats.splats[i].featureSize(), i}); + } + active.extend(newSplats); + } + + cells.clear(); + + // Iterate the heap's internal array order, like Rust's active.iter() + for (const HeapKey &hk : active.data){ + size_t index = hk.index; + I64Vec3 grid = splats.splats[index].grid(step); + cells[grid].push_back(index); + } + + std::vector nextActive; + + HeapKey top; + while (active.pop(top)){ + float negSize = top.key; + size_t index = top.index; + if (!isActive[index]){ + continue; + } + + I64Vec3 grid = splats.splats[index].grid(step); + + size_t bestIndex = SIZE_MAX; + float bestMetric = -std::numeric_limits::infinity(); + I64Vec3 bestGrid(INT64_MAX, INT64_MAX, INT64_MAX); + + for (int64_t z = grid.z - 1; z <= grid.z + 1; z++){ + for (int64_t y = grid.y - 1; y <= grid.y + 1; y++){ + for (int64_t x = grid.x - 1; x <= grid.x + 1; x++){ + I64Vec3 g(x, y, z); + auto it = cells.find(g); + if (it != cells.end()){ + for (size_t neighbor : it->second){ + if (isActive[neighbor] && neighbor != index){ + float metric = splats.similarity(index, neighbor); + if (metric > bestMetric){ + bestIndex = neighbor; + bestMetric = metric; + bestGrid = g; + } + } + } + } + } + } + } + + if (bestIndex != SIZE_MAX){ + size_t bestNeighbor = bestIndex; + size_t mergeIndices[2] = {index, bestNeighbor}; + size_t merged = splats.newMerged(mergeIndices, 2, 0.0f); + + isActive[index] = false; + { + std::vector &cellIndex = cells[grid]; + cellIndex.erase(std::remove(cellIndex.begin(), cellIndex.end(), index), + cellIndex.end()); + } + + isActive[bestNeighbor] = false; + { + std::vector &cellBest = cells[bestGrid]; + cellBest.erase(std::remove(cellBest.begin(), cellBest.end(), bestNeighbor), + cellBest.end()); + } + + isActive.push_back(true); + + float featureSize = splats.splats[merged].featureSize(); + if (featureSize > step){ + nextActive.push_back(HeapKey{-featureSize, merged}); + }else{ + I64Vec3 mergedGrid = splats.splats[merged].grid(step); + cells[mergedGrid].push_back(merged); + + active.push(HeapKey{-featureSize, merged}); + } + }else{ + // Can't find a neighbor to merge, so kick to next level + nextActive.push_back(HeapKey{negSize, index}); + } + } + + level += 1; + active.extend(nextActive); + + if (frontier < initialLen){ + continue; + } + + if (active.len() <= 1){ + break; + } + } + + size_t rootIndex = splats.len() - 1; + + std::vector toOutput; + toOutput.resize(initialLen, true); + toOutput.resize(splats.len(), false); + toOutput[rootIndex] = true; + + { + float rootFeatureSize; + std::vector rootChildren; + bhattRecurseToOutput(splats, rootIndex, toOutput, lodBase, rootFeatureSize, rootChildren); + } + + size_t outputCount = 0; + for (bool b : toOutput){ + if (b) outputCount++; + } + + std::vector indices; + indices.push_back(rootIndex); + float limitSize = splats.splats[rootIndex].featureSize(); + std::vector frontierIdx{rootIndex}; + + for (;;){ + std::vector nextFrontier; + for (size_t index : frontierIdx){ + bhattRecurseIndices(splats, index, indices, limitSize, nextFrontier); + } + frontierIdx.clear(); + + if (nextFrontier.empty()){ + break; + } + limitSize = limitSize / 4.0f; + frontierIdx = std::move(nextFrontier); + } + + assert(indices.size() == outputCount); + + for (size_t index = 0; index < toOutput.size(); index++){ + if (!toOutput[index]){ + indices.push_back(index); + } + } + + splats.permute(indices); + splats.truncate(outputCount); +} + +// chunk_tree (port of spark-lib/src/chunk_tree.rs chunk_tree_size) + +const size_t BATCH_SIZE = 64 * 1024; +const size_t MIN_BATCH_SIZE = 8 * 1024; +const float STD_DEVS = 1.5f; + +struct Aabb { + Vec3 min; + Vec3 max; + + static Aabb empty(){ + Aabb a; + a.min = Vec3::splat(std::numeric_limits::infinity()); + a.max = Vec3::splat(-std::numeric_limits::infinity()); + return a; + } + + // glam Vec3A::min/max are _mm_min_ps/_mm_max_ps: a < b ? a : b / a > b ? a : b + Aabb extend(const Aabb &other) const { + Aabb r; + r.min = Vec3(min.x < other.min.x ? min.x : other.min.x, + min.y < other.min.y ? min.y : other.min.y, + min.z < other.min.z ? min.z : other.min.z); + r.max = max.max(other.max); + return r; + } + + Vec3 center() const { return (min + max) * 0.5f; } + Vec3 extent() const { return max - min; } + + float minElement() const { + Vec3 e = extent(); + float m = e.x < e.y ? e.x : e.y; + return m < e.z ? m : e.z; + } + + // longest_axis: x if ex >= ey && ex >= ez, else y if ey >= ez, else z + int longestAxis() const { + Vec3 e = extent(); + if (e.x >= e.y && e.x >= e.z) return 0; + if (e.y >= e.z) return 1; + return 2; + } + + static Aabb fromSplat(const Gsplat &s, float stdDevs){ + Vec3 clampedScales = s.getScales().max(Vec3::splat(1.0e-3f)); + Vec3 r = clampedScales * stdDevs; + + Mat3 rmat = Mat3::fromQuat(s.getQuaternion()); + // half = |R| * r: glam Mat3A * Vec3A = x_axis*v.x + y_axis*v.y + z_axis*v.z + Vec3 ax(std::fabs(rmat.m[0][0]), std::fabs(rmat.m[0][1]), std::fabs(rmat.m[0][2])); + Vec3 ay(std::fabs(rmat.m[1][0]), std::fabs(rmat.m[1][1]), std::fabs(rmat.m[1][2])); + Vec3 az(std::fabs(rmat.m[2][0]), std::fabs(rmat.m[2][1]), std::fabs(rmat.m[2][2])); + Vec3 half = ax * r.x + ay * r.y + az * r.z; + + Vec3 c = s.getCenter(); + Aabb a; + a.min = c - half; + a.max = c + half; + return a; + } +}; + +void chunkTreeSize(GsplatArray &splats, size_t root){ + std::vector indices; + indices.push_back(root); + + std::vector> batches; + size_t batchesHead = 0; // VecDeque via head index + batches.push_back({root}); + + while (batchesHead < batches.size()){ + std::vector batch = std::move(batches[batchesHead]); + batchesHead++; + + MaxHeap priority; + for (size_t parent : batch){ + priority.push(HeapKey{splats.splats[parent].featureSize(), parent}); + } + + size_t startIndex = indices.size(); + size_t endIndex = ((startIndex + MIN_BATCH_SIZE + BATCH_SIZE - 1) / BATCH_SIZE) * BATCH_SIZE; + + HeapKey top; + while (priority.pop(top)){ + std::vector children = splats.getChildren(top.index); + if ((indices.size() + children.size()) > endIndex){ + priority.push(top); + break; + } + std::vector newChildren(children.size()); + for (size_t i = 0; i < children.size(); i++) newChildren[i] = indices.size() + i; + splats.setChildren(top.index, newChildren); + + for (size_t child : children){ + indices.push_back(child); + priority.push(HeapKey{splats.splats[child].featureSize(), child}); + } + } + + if (!priority.isEmpty()){ + Aabb aabb = Aabb::empty(); + for (const HeapKey &hk : priority.data){ + aabb = aabb.extend(Aabb::fromSplat(splats.splats[hk.index], STD_DEVS)); + } + + if (aabb.extent().maxElement() >= (3.0f * aabb.minElement())){ + int axis = aabb.longestAxis(); + float split = aabb.center()[axis]; + // partition keeps relative (internal array) order; pop order in the + // next round is key-determined, so the batch order within a/b is + // not observable — but keep it faithful anyway. + std::vector a, b; + for (const HeapKey &hk : priority.data){ + if (splats.splats[hk.index].getCenter()[axis] < split) a.push_back(hk.index); + else b.push_back(hk.index); + } + + // sort_by_key(-(len)): stable descending by length, a first on tie + if (b.size() > a.size()) std::swap(a, b); + batches.push_back(std::move(a)); + batches.push_back(std::move(b)); + continue; + } + + std::vector octants[8]; + Vec3 split = aabb.center(); + + for (const HeapKey &hk : priority.data){ + Vec3 center = splats.splats[hk.index].getCenter(); + int octant = (center.x < split.x ? 0 : 1) + (center.y < split.y ? 0 : 2) + + (center.z < split.z ? 0 : 4); + octants[octant].push_back(hk.index); + } + + // Hilbert order + const int hilbert[8] = {0, 1, 3, 2, 6, 7, 5, 4}; + for (int i = 0; i < 8; i++){ + batches.push_back(std::move(octants[hilbert[i]])); + } + } + } + + assert(indices.size() == splats.len()); + splats.permute(indices); +} + +// Ingestion: SplatData -> GsplatArray, applying the same transforms +// spark-lib/src/ply.rs applies when reading a 3DGS PLY + +const float SH_C0 = 0.28209479177387814f; + +size_t maxShDegreeFromRestCoeffs(size_t numRestCoeffs){ + // f_rest count = numRestCoeffs * 3: 0 -> 0, 9 -> 1, 24 -> 2, 45 -> 3 + switch (numRestCoeffs){ + case 0: return 0; + case 3: return 1; + case 8: return 2; + case 15: return 3; + default: return SIZE_MAX; + } +} + +bool ingestSplats(const SplatData &data, GsplatArray &splats){ + size_t maxSh = maxShDegreeFromRestCoeffs(data.numRestCoeffs); + if (maxSh == SIZE_MAX){ + std::cerr << "saveRad: invalid number of SH rest coefficients: " << data.numRestCoeffs + << std::endl; + return false; + } + + size_t n = data.numPoints; + splats.maxShDegree = maxSh; + splats.splats.resize(n); + if (maxSh >= 1) splats.sh1.resize(n); + if (maxSh >= 2) splats.sh2.resize(n); + if (maxSh >= 3) splats.sh3.resize(n); + + size_t K = data.numRestCoeffs; + + for (size_t i = 0; i < n; i++){ + size_t i3 = i * 3; + size_t i4 = i * 4; + + Vec3 center(data.means[i3], data.means[i3 + 1], data.means[i3 + 2]); + + float opLogistic = data.opacities[i]; + float opacity = 1.0f / (1.0f + std::exp(-opLogistic)); + + Vec3 rgb(0.5f + data.featuresDc[i3] * SH_C0, + 0.5f + data.featuresDc[i3 + 1] * SH_C0, + 0.5f + data.featuresDc[i3 + 2] * SH_C0); + + Vec3 scale(std::exp(data.scales[i3]), std::exp(data.scales[i3 + 1]), + std::exp(data.scales[i3 + 2])); + + // stored order w,x,y,z (rot_0..rot_3); spark reads [rot_1, rot_2, rot_3, rot_0] + float quatArr[4] = {data.quats[i4 + 1], data.quats[i4 + 2], data.quats[i4 + 3], + data.quats[i4]}; + float quatMagnitude = std::sqrt(((quatArr[0] * quatArr[0] + quatArr[1] * quatArr[1]) + + quatArr[2] * quatArr[2]) + quatArr[3] * quatArr[3]); + Quat quat(quatArr[0] / quatMagnitude, quatArr[1] / quatMagnitude, + quatArr[2] / quatMagnitude, quatArr[3] / quatMagnitude); + + // set_batch applies each setter (f16 quantization) individually. + Gsplat &s = splats.splats[i]; + s.setCenter(center); + s.setOpacity(opacity); + s.setRgb(rgb); + s.setScales(scale); // f16(ln(exp(stored))) + s.setQuaternion(quat); + + // featuresRest [i][k][d] linear order equals spark's interleaved SH slots: + // sh1 = coeffs 0..2, sh2 = coeffs 3..7, sh3 = coeffs 8..14 + const float *fr = data.featuresRest.data() + i * K * 3; + if (maxSh >= 1){ + for (int d = 0; d < 9; d++) splats.sh1[i][d] = F16::fromF32(fr[d]); + } + if (maxSh >= 2){ + for (int d = 0; d < 15; d++) splats.sh2[i][d] = F16::fromF32(fr[9 + d]); + } + if (maxSh >= 3){ + for (int d = 0; d < 21; d++) splats.sh3[i][d] = F16::fromF32(fr[24 + d]); + } + } + + return true; +} + +// Pipeline (port of build-lod/src/main.rs process_file_lod_tsplat, default +// options) minus the encoder, which is invoked by saveRad + +struct PipelineResult { + size_t inputSplatCount = 0; + size_t inputShDegree = 0; + bool removedAny = false; + size_t emptySplatCount = 0; + size_t initialSplatCount = 0; + double lodDuration = 0.0; + size_t finalSplatCount = 0; + double chunkDuration = 0.0; + size_t maxShDegree = 0; +}; + +bool runPipeline(const SplatData &data, GsplatArray &splats, PipelineResult &result){ + if (!ingestSplats(data, splats)){ + return false; + } + + result.inputSplatCount = splats.len(); + result.inputShDegree = splats.maxShDegree; + + // Validation (main.rs, no --skip-validate) + { + size_t invalidCount = 0; + for (size_t i = 0; i < splats.len(); i++){ + const Gsplat &s = splats.splats[i]; + if (!s.getCenter().isFinite() || !s.getScales().isFinite() || + !s.getQuaternion().isFinite() || !std::isfinite(s.getOpacity()) || + !s.getRgb().isFinite()){ + invalidCount += 1; + } + } + if (invalidCount > 0){ + std::cerr << "saveRad: found " << invalidCount << " invalid (non-finite) splats" + << std::endl; + return false; + } + } + + splats.retain([](const Gsplat &s){ + return (s.getOpacity() > 0.0f) && (s.maxScale() > 0.0f) && + (s.getQuaternion().isFinite() && s.getQuaternion().length() > 0.0f); + }); + + if (result.inputSplatCount != splats.len()){ + result.removedAny = true; + result.emptySplatCount = result.inputSplatCount - splats.len(); + result.initialSplatCount = splats.len(); + } + + // LoD: default method Quality -> BhattLod { lod_base: 1.75 } + auto lodStart = std::chrono::steady_clock::now(); + bhattComputeLodTree(splats, 1.75f); + result.lodDuration = std::chrono::duration(std::chrono::steady_clock::now() - lodStart).count(); + + result.finalSplatCount = splats.len(); + + auto chunkStart = std::chrono::steady_clock::now(); + chunkTreeSize(splats, 0); + result.chunkDuration = std::chrono::duration(std::chrono::steady_clock::now() - chunkStart).count(); + + result.maxShDegree = splats.maxShDegree; + + splats.encodeLodOpacity(); + + return true; +} + +// Deflate via zlib: raw stream (no zlib header, windowBits -15), level 6 + +std::vector compressToVec(const std::vector &data){ + z_stream strm; + std::memset(&strm, 0, sizeof(strm)); + if (deflateInit2(&strm, 6, Z_DEFLATED, -15, 8, Z_DEFAULT_STRATEGY) != Z_OK){ + throw std::runtime_error("saveRad: deflateInit2 failed"); + } + uLong bound = deflateBound(&strm, static_cast(data.size())); + std::vector out(bound); + strm.next_in = const_cast(data.data()); + strm.avail_in = static_cast(data.size()); + strm.next_out = out.data(); + strm.avail_out = static_cast(bound); + int ret = deflate(&strm, Z_FINISH); + if (ret != Z_STREAM_END){ + deflateEnd(&strm); + throw std::runtime_error("saveRad: deflate failed"); + } + out.resize(strm.total_out); + deflateEnd(&strm); + return out; +} + +// Property encoders (port of rad.rs encode_*). All planar (dimension-major) +// except oct88r8, which is 3 bytes per splat interleaved. + +std::vector encodeF32(const std::vector &data, size_t dims, size_t count){ + std::vector result; + result.reserve(4 * dims * count); + for (size_t d = 0; d < dims; d++){ + size_t index = d; + for (size_t i = 0; i < count; i++){ + uint32_t bits = f32Bits(data[index]); + result.push_back(static_cast(bits)); + result.push_back(static_cast(bits >> 8)); + result.push_back(static_cast(bits >> 16)); + result.push_back(static_cast(bits >> 24)); + index += dims; + } + } + return result; +} + +std::vector encodeF16Prop(const std::vector &data, size_t dims, size_t count){ + std::vector result; + result.reserve(2 * dims * count); + for (size_t d = 0; d < dims; d++){ + size_t index = d; + for (size_t i = 0; i < count; i++){ + uint16_t bits = f16FromF32(data[index]); + result.push_back(static_cast(bits)); + result.push_back(static_cast(bits >> 8)); + index += dims; + } + } + return result; +} + +std::vector encodeF32LeBytes(const std::vector &data, size_t dims, size_t count){ + std::vector result; + result.reserve(4 * dims * count); + for (int b = 0; b < 4; b++){ + for (size_t d = 0; d < dims; d++){ + size_t index = d; + for (size_t i = 0; i < count; i++){ + uint32_t bits = f32Bits(data[index]); + result.push_back(static_cast(bits >> (8 * b))); + index += dims; + } + } + } + return result; +} + +std::vector encodeR8(const std::vector &data, size_t dims, size_t count, float min, + float max){ + std::vector result; + result.reserve(dims * count); + for (size_t d = 0; d < dims; d++){ + size_t index = d; + for (size_t i = 0; i < count; i++){ + float value = (data[index] - min) / (max - min) * 255.0f; + result.push_back(rustCastU8(std::round(std::clamp(value, 0.0f, 255.0f)))); + index += dims; + } + } + return result; +} + +std::vector encodeS8(const std::vector &data, size_t dims, size_t count, float max){ + std::vector result; + result.reserve(dims * count); + for (size_t d = 0; d < dims; d++){ + size_t index = d; + for (size_t i = 0; i < count; i++){ + float value = data[index] / max * 127.0f; + result.push_back(static_cast( + rustCastI8(std::round(std::clamp(value, -127.0f, 127.0f))))); + index += dims; + } + } + return result; +} + +std::vector encodeR8Delta(const std::vector &data, size_t dims, size_t count, + float min, float max){ + std::vector result; + result.reserve(dims * count); + for (size_t d = 0; d < dims; d++){ + size_t index = d; + uint8_t last = 0; + for (size_t i = 0; i < count; i++){ + float v = (data[index] - min) / (max - min) * 255.0f; + uint8_t value = rustCastU8(std::round(std::clamp(v, 0.0f, 255.0f))); + result.push_back(static_cast(value - last)); // wrapping_sub + last = value; + index += dims; + } + } + return result; +} + +// splat_encode.rs encode_scale8_zero +uint8_t encodeScale8Zero(float scale, float lnScaleZero, float lnScaleMin, float lnScaleMax){ + if (scale <= 0.0f) return 0; + float lnScale = std::log(scale); + if (lnScale <= lnScaleZero) return 0; + float value = (lnScale - lnScaleMin) / (lnScaleMax - lnScaleMin) * 254.0f; + return static_cast(rustCastU8(std::round(std::clamp(value, 0.0f, 254.0f))) + 1); +} + +std::vector encodeLn0R8(const std::vector &data, size_t dims, size_t count, + float zero, float min, float max){ + std::vector result; + result.reserve(dims * count); + for (size_t d = 0; d < dims; d++){ + size_t index = d; + for (size_t i = 0; i < count; i++){ + result.push_back(encodeScale8Zero(data[index], zero, min, max)); + index += dims; + } + } + return result; +} + +std::vector encodeLnF16(const std::vector &data, size_t dims, size_t count){ + std::vector result; + result.reserve(2 * dims * count); + for (size_t d = 0; d < dims; d++){ + size_t index = d; + for (size_t i = 0; i < count; i++){ + uint16_t bits = f16FromF32(std::log(data[index])); + result.push_back(static_cast(bits)); + result.push_back(static_cast(bits >> 8)); + index += dims; + } + } + return result; +} + +// splat_encode.rs float_to_u8 +uint8_t floatToU8(float value, float min, float max){ + return rustCastU8(std::round(std::clamp((value - min) / (max - min) * 255.0f, 0.0f, 255.0f))); +} + +// splat_encode.rs encode_quat_oct888 (quat in xyzw order) +void encodeQuatOct888(const float quatXyzw[4], uint8_t out[3]){ + float quat[4]; + if (quatXyzw[3] < 0.0f){ + for (int i = 0; i < 4; i++) quat[i] = -quatXyzw[i]; + }else{ + for (int i = 0; i < 4; i++) quat[i] = quatXyzw[i]; + } + float theta = 2.0f * std::acos(std::clamp(quat[3], 0.0f, 1.0f)); + float s = std::sin(theta * 0.5f); + + float axis[3]; + if (std::fabs(s) < 1e-6f){ + axis[0] = 1.0f; + axis[1] = 0.0f; + axis[2] = 0.0f; + }else{ + for (int i = 0; i < 3; i++) axis[i] = quat[i] / s; + } + float sum = std::fabs(axis[0]) + std::fabs(axis[1]) + std::fabs(axis[2]); + float p[2] = {axis[0] / sum, axis[1] / sum}; + if (axis[2] < 0.0f){ + float p0 = (1.0f - std::fabs(p[1])) * (p[0] >= 0.0f ? 1.0f : -1.0f); + float p1 = (1.0f - std::fabs(p[0])) * (p[1] >= 0.0f ? 1.0f : -1.0f); + p[0] = p0; + p[1] = p1; + } + out[0] = floatToU8(p[0], -1.0f, 1.0f); + out[1] = floatToU8(p[1], -1.0f, 1.0f); + out[2] = floatToU8(theta, 0.0f, 3.14159265358979323846264338327950288f); +} + +std::vector encodeQuatOct88R8(const std::vector &data, size_t count){ + std::vector result; + result.reserve(3 * count); + for (size_t i = 0; i < count; i++){ + uint8_t enc[3]; + encodeQuatOct888(&data[i * 4], enc); + result.push_back(enc[0]); + result.push_back(enc[1]); + result.push_back(enc[2]); + } + return result; +} + +std::vector encodeU16Prop(const std::vector &data, size_t count){ + std::vector result; + result.reserve(2 * count); + for (size_t i = 0; i < count; i++){ + result.push_back(static_cast(data[i])); + result.push_back(static_cast(data[i] >> 8)); + } + return result; +} + +std::vector encodeUsizeAsU32(const std::vector &data, size_t count){ + std::vector result; + result.reserve(4 * count); + for (size_t i = 0; i < count; i++){ + uint32_t v = rustCastU32(data[i]); + result.push_back(static_cast(v)); + result.push_back(static_cast(v >> 8)); + result.push_back(static_cast(v >> 16)); + result.push_back(static_cast(v >> 24)); + } + return result; +} + +// RadEncoder (port of rad.rs RadEncoder for the default configuration: +// no SH clusters, resolve everything from Auto) + +const uint32_t RAD_MAGIC = 0x30444152; // 'RAD0' +const uint32_t RAD_CHUNK_MAGIC = 0x43444152; // 'RADC' + +inline size_t roundup8(size_t size){ return (size + 7) & ~static_cast(7); } +inline size_t pad8(size_t size){ return (8 - (size & 7)) & 7; } + +struct SplatEncodingParams { + float rgbMin = 0.0f; + float rgbMax = 1.0f; + float lnScaleMin = -12.0f; + float lnScaleMax = 9.0f; + float sh1Max = 1.0f; + float sh2Max = 1.0f; + float sh3Max = 1.0f; + bool lodOpacity = false; +}; + +enum class PropEncoding { F32LeBytes, F16, R8, R8Delta, Ln0R8, LnF16, Oct88R8, S8, U16, U32 }; + +const char *propEncodingName(PropEncoding e){ + switch (e){ + case PropEncoding::F32LeBytes: return "f32_lebytes"; + case PropEncoding::F16: return "f16"; + case PropEncoding::R8: return "r8"; + case PropEncoding::R8Delta: return "r8_delta"; + case PropEncoding::Ln0R8: return "ln_0r8"; + case PropEncoding::LnF16: return "ln_f16"; + case PropEncoding::Oct88R8: return "oct88r8"; + case PropEncoding::S8: return "s8"; + case PropEncoding::U16: return "u16"; + case PropEncoding::U32: return "u32"; + } + return ""; +} + +struct EncodedProp { + const char *name; + PropEncoding encoding; + bool hasMinMax = false; + float min = 0.0f; + float max = 0.0f; + uint64_t offset = 0; + uint64_t bytes = 0; + std::vector data; +}; + +struct RadEncoder { + GsplatArray &splats; + size_t maxSh; // min(getter max_sh_degree, 3) + bool hasEncoding = false; + SplatEncodingParams encoding; + + // Rust-enum-variant names of the current encoding state (for the comment) + std::string centerEncoding = "Auto"; + std::string alphaEncoding = "Auto"; + std::string rgbEncoding = "Auto"; + std::string scalesEncoding = "Auto"; + std::string orientationEncoding = "Auto"; + std::string shEncoding = "Auto"; + std::string shLabelEncoding = "Auto"; + + std::string comment; + + explicit RadEncoder(GsplatArray &s) : splats(s), maxSh(std::min(s.maxShDegree, size_t(3))) {} + + // Percentile helper: n-th order statistic with OrderedFloat comparator + static float selectNth(std::vector &v, size_t n){ + std::nth_element(v.begin(), v.begin() + n, v.end()); + return v[n]; + } + + void resolveEncoding(){ + // resolve_center_encoding + centerEncoding = "F32LeBytes"; + + // resolve_alpha_encoding + { + float maxAlpha = -std::numeric_limits::infinity(); + for (const Gsplat &s : splats.splats){ + maxAlpha = std::fmaxf(maxAlpha, s.getOpacity()); + } + alphaEncoding = (maxAlpha > 1.0f) ? "F16" : "R8"; + } + + // resolve_rgb_encoding + { + size_t n = splats.len(); + std::vector allRgb(n * 3); + for (size_t i = 0; i < n; i++){ + Vec3 rgb = splats.splats[i].getRgb(); + allRgb[i * 3] = rgb.x; + allRgb[i * 3 + 1] = rgb.y; + allRgb[i * 3 + 2] = rgb.z; + } + size_t len = allRgb.size(); + size_t n1 = std::min(static_cast(std::round(static_cast(len) * 0.01f)), len - 1); + size_t n99 = std::min(static_cast(std::round(static_cast(len) * 0.99f)), len - 1); + float rgb1 = selectNth(allRgb, n1); + float rgb99 = selectNth(allRgb, n99); + float rgbMin = std::fminf(rgb1, 0.0f); + float rgbMax = std::fmaxf(rgb99, 1.0f); + + if (rgbMin < -1.0f || rgbMax > 2.0f){ + rgbEncoding = "F16"; + }else{ + rgbEncoding = "R8Delta"; + hasEncoding = true; + encoding.rgbMin = rgbMin; + encoding.rgbMax = rgbMax; + } + } + + // resolve_scales_encoding + { + std::vector scales; + scales.reserve(splats.len() * 2); + for (const Gsplat &s : splats.splats){ + Vec3 sc = s.getScales(); + float splatScales[3] = {sc.x, sc.y, sc.z}; + // 3-element stable ascending sort (Rust sort_by_key) + std::stable_sort(splatScales, splatScales + 3); + scales.push_back(splatScales[1]); + scales.push_back(splatScales[2]); + } + size_t len = scales.size(); + size_t n1 = std::min(static_cast(std::round(static_cast(len) * 0.01f)), len - 1); + size_t n99 = std::min(static_cast(std::round(static_cast(len) * 0.99f)), len - 1); + float scale1 = selectNth(scales, n1); + float scale99 = selectNth(scales, n99); + float lnScaleMin = std::fminf(std::log(std::fmaxf(scale1, 1.0e-30f)), -12.0f); + float lnScaleMax = std::fmaxf(std::log(std::fmaxf(scale99, 1.0e-30f)), 9.0f); + + if ((lnScaleMax - lnScaleMin) > 25.0f){ + scalesEncoding = "LnF16"; + }else{ + scalesEncoding = "Ln0R8"; + hasEncoding = true; + encoding.lnScaleMin = lnScaleMin; + encoding.lnScaleMax = lnScaleMax; + } + } + + // resolve_orientation_encoding + orientationEncoding = "Oct88R8"; + + // resolve_sh_encoding + { + size_t numSh = std::min(maxSh, splats.maxShDegree); + if (numSh == 0){ + shEncoding = "S8"; + }else{ + hasEncoding = true; + + size_t n = splats.len(); + std::vector allRgb(n * 9); + for (size_t i = 0; i < n; i++){ + for (int d = 0; d < 9; d++) allRgb[i * 9 + d] = splats.sh1[i][d].toF32(); + } + { + size_t len = allRgb.size(); + size_t n5 = std::min(static_cast(std::round(static_cast(len) * 0.05f)), len - 1); + size_t n95 = std::min(static_cast(std::round(static_cast(len) * 0.95f)), len - 1); + float sh5 = selectNth(allRgb, n5); + float sh95 = selectNth(allRgb, n95); + encoding.sh1Max = std::fmaxf(std::fmaxf(std::fabs(sh5), std::fabs(sh95)), 1.0f); + } + + if (numSh >= 2){ + allRgb.assign(n * 15, 0.0f); + for (size_t i = 0; i < n; i++){ + for (int d = 0; d < 15; d++) allRgb[i * 15 + d] = splats.sh2[i][d].toF32(); + } + size_t len = allRgb.size(); + size_t n5 = std::min(static_cast(std::round(static_cast(len) * 0.05f)), len - 1); + size_t n95 = std::min(static_cast(std::round(static_cast(len) * 0.95f)), len - 1); + float sh5 = selectNth(allRgb, n5); + float sh95 = selectNth(allRgb, n95); + encoding.sh2Max = std::fmaxf(std::fmaxf(std::fabs(sh5), std::fabs(sh95)), 1.0f); + } + + if (numSh >= 3){ + allRgb.assign(n * 21, 0.0f); + for (size_t i = 0; i < n; i++){ + for (int d = 0; d < 21; d++) allRgb[i * 21 + d] = splats.sh3[i][d].toF32(); + } + size_t len = allRgb.size(); + size_t n5 = std::min(static_cast(std::round(static_cast(len) * 0.05f)), len - 1); + size_t n95 = std::min(static_cast(std::round(static_cast(len) * 0.95f)), len - 1); + float sh5 = selectNth(allRgb, n5); + float sh95 = selectNth(allRgb, n95); + encoding.sh3Max = std::fmaxf(std::fmaxf(std::fabs(sh5), std::fabs(sh95)), 1.0f); + } + + shEncoding = "S8"; + } + } + + // resolve_sh_label_encoding: no clusters -> stays Auto + } + + json encodingStateJson() const { + json obj = json::object(); + obj["alpha"] = alphaEncoding; + obj["center"] = centerEncoding; + if (!hasEncoding){ + obj["encoding"] = nullptr; + }else{ + json enc = json::object(); + enc["lnScaleMax"] = encoding.lnScaleMax; + enc["lnScaleMin"] = encoding.lnScaleMin; + enc["lodOpacity"] = encoding.lodOpacity; + enc["rgbMax"] = encoding.rgbMax; + enc["rgbMin"] = encoding.rgbMin; + enc["sh1Max"] = encoding.sh1Max; + enc["sh2Max"] = encoding.sh2Max; + enc["sh3Max"] = encoding.sh3Max; + obj["encoding"] = std::move(enc); + } + obj["orientation"] = orientationEncoding; + obj["rgb"] = rgbEncoding; + obj["scales"] = scalesEncoding; + obj["sh"] = shEncoding; + obj["sh_label"] = shLabelEncoding; + return obj; + } + + json setSplatEncodingJson(bool lodOpacity) const { + json enc = json::object(); + enc["rgbMin"] = encoding.rgbMin; + enc["rgbMax"] = encoding.rgbMax; + enc["lnScaleMin"] = encoding.lnScaleMin; + enc["lnScaleMax"] = encoding.lnScaleMax; + enc["sh1Max"] = encoding.sh1Max; + enc["sh2Max"] = encoding.sh2Max; + enc["sh3Max"] = encoding.sh3Max; + enc["lodOpacity"] = lodOpacity; + return enc; + } + + EncodedProp encodeChunkCenter(size_t base, size_t count, std::vector &buffer){ + if (buffer.size() < count * 3) buffer.resize(count * 3, 0.0f); + for (size_t i = 0; i < count; i++){ + const Vec3 &c = splats.splats[base + i].center; + buffer[i * 3] = c.x; + buffer[i * 3 + 1] = c.y; + buffer[i * 3 + 2] = c.z; + } + EncodedProp p; + p.name = "center"; + p.encoding = PropEncoding::F32LeBytes; + p.data = compressToVec(encodeF32LeBytes(buffer, 3, count)); + return p; + } + + EncodedProp encodeChunkAlpha(size_t base, size_t count, std::vector &buffer){ + if (buffer.size() < count) buffer.resize(count, 0.0f); + for (size_t i = 0; i < count; i++){ + buffer[i] = splats.splats[base + i].getOpacity(); + } + float maxAlpha = splats.hasLodTree() ? 2.0f : 1.0f; + EncodedProp p; + p.name = "alpha"; + if (alphaEncoding == "R8"){ + p.encoding = PropEncoding::R8; + p.hasMinMax = true; + p.min = 0.0f; + p.max = maxAlpha; + p.data = compressToVec(encodeR8(buffer, 1, count, 0.0f, maxAlpha)); + }else{ + p.encoding = PropEncoding::F16; + p.data = compressToVec(encodeF16Prop(buffer, 1, count)); + } + return p; + } + + EncodedProp encodeChunkRgb(size_t base, size_t count, std::vector &buffer){ + if (buffer.size() < count * 3) buffer.resize(count * 3, 0.0f); + for (size_t i = 0; i < count; i++){ + Vec3 rgb = splats.splats[base + i].getRgb(); + buffer[i * 3] = rgb.x; + buffer[i * 3 + 1] = rgb.y; + buffer[i * 3 + 2] = rgb.z; + } + EncodedProp p; + p.name = "rgb"; + if (rgbEncoding == "F16"){ + p.encoding = PropEncoding::F16; + p.data = compressToVec(encodeF16Prop(buffer, 3, count)); + }else{ + p.encoding = PropEncoding::R8Delta; + p.hasMinMax = true; + p.min = encoding.rgbMin; + p.max = encoding.rgbMax; + p.data = compressToVec(encodeR8Delta(buffer, 3, count, encoding.rgbMin, encoding.rgbMax)); + } + return p; + } + + EncodedProp encodeChunkScales(size_t base, size_t count, std::vector &buffer){ + if (buffer.size() < count * 3) buffer.resize(count * 3, 0.0f); + for (size_t i = 0; i < count; i++){ + Vec3 sc = splats.splats[base + i].getScales(); + buffer[i * 3] = sc.x; + buffer[i * 3 + 1] = sc.y; + buffer[i * 3 + 2] = sc.z; + } + EncodedProp p; + p.name = "scales"; + if (scalesEncoding == "LnF16"){ + p.encoding = PropEncoding::LnF16; + p.data = compressToVec(encodeLnF16(buffer, 3, count)); + }else{ + p.encoding = PropEncoding::Ln0R8; + p.hasMinMax = true; + p.min = encoding.lnScaleMin; + p.max = encoding.lnScaleMax; + p.data = compressToVec( + encodeLn0R8(buffer, 3, count, -30.0f, encoding.lnScaleMin, encoding.lnScaleMax)); + } + return p; + } + + EncodedProp encodeChunkOrientation(size_t base, size_t count, std::vector &buffer){ + if (buffer.size() < count * 4) buffer.resize(count * 4, 0.0f); + for (size_t i = 0; i < count; i++){ + Quat q = splats.splats[base + i].getQuaternion(); + buffer[i * 4] = q.x; + buffer[i * 4 + 1] = q.y; + buffer[i * 4 + 2] = q.z; + buffer[i * 4 + 3] = q.w; + } + EncodedProp p; + p.name = "orientation"; + p.encoding = PropEncoding::Oct88R8; + p.data = compressToVec(encodeQuatOct88R8(buffer, count)); + return p; + } + + EncodedProp encodeChunkSh(size_t base, size_t count, std::vector &buffer, int degree){ + size_t elements = degree == 1 ? 9 : (degree == 2 ? 15 : 21); + float shMax = degree == 1 ? encoding.sh1Max : (degree == 2 ? encoding.sh2Max : encoding.sh3Max); + if (buffer.size() < count * elements) buffer.resize(count * elements, 0.0f); + for (size_t i = 0; i < count; i++){ + if (degree == 1){ + for (size_t d = 0; d < 9; d++) buffer[i * 9 + d] = splats.sh1[base + i][d].toF32(); + }else if (degree == 2){ + for (size_t d = 0; d < 15; d++) buffer[i * 15 + d] = splats.sh2[base + i][d].toF32(); + }else{ + for (size_t d = 0; d < 21; d++) buffer[i * 21 + d] = splats.sh3[base + i][d].toF32(); + } + } + EncodedProp p; + p.name = degree == 1 ? "sh1" : (degree == 2 ? "sh2" : "sh3"); + p.encoding = PropEncoding::S8; + p.hasMinMax = true; + p.min = -shMax; + p.max = shMax; + p.data = compressToVec(encodeS8(buffer, elements, count, shMax)); + return p; + } + + EncodedProp encodeChunkChildCount(size_t base, size_t count, std::vector &buffer){ + if (buffer.size() < count) buffer.resize(count, 0); + for (size_t i = 0; i < count; i++){ + buffer[i] = static_cast(splats.children[base + i].size()); + } + EncodedProp p; + p.name = "child_count"; + p.encoding = PropEncoding::U16; + p.data = compressToVec(encodeU16Prop(buffer, count)); + return p; + } + + EncodedProp encodeChunkChildStart(size_t base, size_t count, std::vector &buffer){ + if (buffer.size() < count) buffer.resize(count, 0); + for (size_t i = 0; i < count; i++){ + const std::vector &c = splats.children[base + i]; + buffer[i] = c.empty() ? 0 : c[0]; + } + EncodedProp p; + p.name = "child_start"; + p.encoding = PropEncoding::U32; + p.data = compressToVec(encodeUsizeAsU32(buffer, count)); + return p; + } + + std::vector encodeChunk(size_t base, size_t count, std::vector &buffer, + std::vector &bufferU16, + std::vector &bufferUsize){ + size_t chunkMaxSh = std::min(splats.maxShDegree, maxSh); + + std::vector props; + props.push_back(encodeChunkCenter(base, count, buffer)); + props.push_back(encodeChunkAlpha(base, count, buffer)); + props.push_back(encodeChunkRgb(base, count, buffer)); + props.push_back(encodeChunkScales(base, count, buffer)); + props.push_back(encodeChunkOrientation(base, count, buffer)); + + if (chunkMaxSh >= 1) props.push_back(encodeChunkSh(base, count, buffer, 1)); + if (chunkMaxSh >= 2) props.push_back(encodeChunkSh(base, count, buffer, 2)); + if (chunkMaxSh >= 3) props.push_back(encodeChunkSh(base, count, buffer, 3)); + + if (splats.hasLodTree()){ + props.push_back(encodeChunkChildCount(base, count, bufferU16)); + props.push_back(encodeChunkChildStart(base, count, bufferUsize)); + } + + uint64_t offset = 0; + for (EncodedProp &p : props){ + p.offset = offset; + p.bytes = p.data.size(); + offset += roundup8(p.data.size()); + } + uint64_t payloadBytes = offset; + + json meta = json::object(); + meta["version"] = 1; + meta["base"] = base; + meta["count"] = count; + meta["payloadBytes"] = payloadBytes; + meta["maxSh"] = chunkMaxSh; + if (splats.hasLodTree()) meta["lodTree"] = true; + if (hasEncoding){ + meta["splatEncoding"] = setSplatEncodingJson(splats.hasLodTree()); + } + json propsJson = json::array(); + for (const EncodedProp &p : props){ + json pj = json::object(); + pj["offset"] = p.offset; + pj["bytes"] = p.bytes; + pj["property"] = p.name; + pj["encoding"] = propEncodingName(p.encoding); + pj["compression"] = "gz"; + if (p.hasMinMax){ + pj["min"] = p.min; + pj["max"] = p.max; + } + propsJson.push_back(std::move(pj)); + } + meta["properties"] = std::move(propsJson); + + std::string metaStr = meta.dump(); + + std::vector encoded; + encoded.reserve(8 + roundup8(metaStr.size()) + 8 + static_cast(payloadBytes)); + auto pushU32 = [&](uint32_t v){ + for (int b = 0; b < 4; b++) encoded.push_back(static_cast(v >> (8 * b))); + }; + auto pushU64 = [&](uint64_t v){ + for (int b = 0; b < 8; b++) encoded.push_back(static_cast(v >> (8 * b))); + }; + pushU32(RAD_CHUNK_MAGIC); + pushU32(static_cast(metaStr.size())); + encoded.insert(encoded.end(), metaStr.begin(), metaStr.end()); + encoded.insert(encoded.end(), pad8(metaStr.size()), 0); + pushU64(payloadBytes); + for (const EncodedProp &p : props){ + encoded.insert(encoded.end(), p.data.begin(), p.data.end()); + encoded.insert(encoded.end(), pad8(p.data.size()), 0); + } + return encoded; + } + + bool encode(std::ofstream &out){ + const size_t CHUNK_SIZE = 65536; + + size_t numSplats = splats.len(); + size_t fileMaxSh = std::min(splats.maxShDegree, maxSh); + + std::vector buffer; + size_t bufferDim = fileMaxSh == 0 ? 4 : (fileMaxSh == 1 ? 9 : (fileMaxSh == 2 ? 15 : 21)); + buffer.resize(CHUNK_SIZE * bufferDim, 0.0f); + + std::vector bufferU16; + std::vector bufferUsize; + if (splats.hasLodTree()){ + bufferU16.resize(CHUNK_SIZE, 0); + bufferUsize.resize(CHUNK_SIZE, 0); + } + + size_t numChunks = (numSplats + CHUNK_SIZE - 1) / CHUNK_SIZE; + std::vector> chunks; + chunks.reserve(numChunks); + std::vector> chunkRanges; // offset, bytes + chunkRanges.reserve(numChunks); + uint64_t offset = 0; + + for (size_t chunkIndex = 0; chunkIndex < numChunks; chunkIndex++){ + size_t base = chunkIndex * CHUNK_SIZE; + size_t count = std::min(numSplats - base, CHUNK_SIZE); + std::vector chunk = encodeChunk(base, count, buffer, bufferU16, bufferUsize); + chunkRanges.push_back({offset, chunk.size()}); + offset += chunk.size(); + chunks.push_back(std::move(chunk)); + } + uint64_t allChunkBytes = offset; + + json meta = json::object(); + meta["version"] = 1; + meta["type"] = "gsplat"; + meta["count"] = numSplats; + meta["maxSh"] = fileMaxSh; + if (splats.hasLodTree()) meta["lodTree"] = true; + meta["chunkSize"] = 65536; + meta["allChunkBytes"] = allChunkBytes; + json chunksJson = json::array(); + for (const auto &cr : chunkRanges){ + json cj = json::object(); + cj["offset"] = cr.first; + cj["bytes"] = cr.second; + chunksJson.push_back(std::move(cj)); + } + meta["chunks"] = std::move(chunksJson); + if (hasEncoding){ + meta["splatEncoding"] = setSplatEncodingJson(splats.hasLodTree()); + } + if (!comment.empty()){ + meta["comment"] = comment; + } + + std::string metaStr = meta.dump(2); + metaStr.push_back('\n'); + size_t metaBytesSize = metaStr.size(); + + auto writeU32 = [&](uint32_t v){ + char b[4]; + for (int i = 0; i < 4; i++) b[i] = static_cast(v >> (8 * i)); + out.write(b, 4); + }; + writeU32(RAD_MAGIC); + writeU32(static_cast(metaBytesSize)); + out.write(metaStr.data(), static_cast(metaStr.size())); + { + size_t pad = pad8(metaBytesSize); + if (pad != 0){ + char zeros[8] = {0}; + out.write(zeros, static_cast(pad)); + } + } + for (const std::vector &chunk : chunks){ + assert((chunk.size() & 7) == 0); + out.write(reinterpret_cast(chunk.data()), + static_cast(chunk.size())); + } + return out.good(); + } +}; + +} + +bool saveRad(const std::string &filename, const SplatData &data){ + GsplatArray splats; + PipelineResult result; + + if (!runPipeline(data, splats, result)){ + return false; + } + + RadEncoder encoder(splats); + + // nlohmann::json objects keep keys sorted alphabetically + nlohmann::json description; + description["input_splat_count"] = result.inputSplatCount; + description["input_sh_degree"] = result.inputShDegree; + if (result.removedAny){ + description["empty_splat_count"] = result.emptySplatCount; + description["initial_splat_count"] = result.initialSplatCount; + } + description["method"] = "BhattLod { lod_base: 1.75 }"; + description["lod_duration"] = result.lodDuration; + description["final_splat_count"] = result.finalSplatCount; + description["chunk_duration"] = result.chunkDuration; + description["max_sh_degree"] = result.maxShDegree; + + description["input_encoding"] = encoder.encodingStateJson(); + + encoder.resolveEncoding(); + description["resolved_encoding"] = encoder.encodingStateJson(); + + encoder.comment = description.dump(2); + + std::ofstream out(filename, std::ios::binary); + if (!out.is_open()){ + std::cerr << "saveRad: cannot open " << filename << " for writing" << std::endl; + return false; + } + try { + return encoder.encode(out); + } catch (const std::exception &e){ + std::cerr << "saveRad: " << e.what() << std::endl; + return false; + } +} + +} + diff --git a/rad.hpp b/rad.hpp new file mode 100644 index 0000000..3454f43 --- /dev/null +++ b/rad.hpp @@ -0,0 +1,35 @@ +#ifndef RAD_H +#define RAD_H + +#include +#include +#include + +namespace rad { + +// Input splats, laid out exactly as Model::savePly writes them (post keepCrs): +// means n*3 world positions +// featuresDc n*3 SH degree-0 coefficients (not RGB) +// featuresRest n*K*3, coefficient-major (K in {0, 3, 8, 15}) +// opacities n logits +// scales n*3 log-space +// quats n*4 unnormalized, PLY order rot_0..rot_3 (w, x, y, z) +struct SplatData { + size_t numPoints = 0; + size_t numRestCoeffs = 0; // K: 0, 3, 8 or 15 (SH degree 0..3) + std::vector means; + std::vector featuresDc; + std::vector featuresRest; + std::vector opacities; + std::vector scales; + std::vector quats; +}; + +// Runs the full build-lod pipeline (LoD tree, chunk reordering, encoding) +// and writes the .rad file. Returns false on failure (non-finite input +// values or I/O error), matching build-lod's abort behavior. +bool saveRad(const std::string &filename, const SplatData &data); + +} + +#endif diff --git a/vendor/spz/CMakeLists.txt b/vendor/spz/CMakeLists.txt new file mode 100644 index 0000000..34c5765 --- /dev/null +++ b/vendor/spz/CMakeLists.txt @@ -0,0 +1,31 @@ +find_package(zstd QUIET) +if(zstd_FOUND) + if(TARGET zstd::libzstd_static) + set(SPZ_ZSTD_LIB zstd::libzstd_static) + else() + set(SPZ_ZSTD_LIB zstd::libzstd_shared) + endif() +else() + include(FetchContent) + FetchContent_Declare(zstd + URL https://github.com/facebook/zstd/releases/download/v1.5.6/zstd-1.5.6.tar.gz + URL_HASH SHA256=8c29e06cf42aacc1eafc4077ae2ec6c6fcb96a626157e0593d5e82a34fd403c1 + SOURCE_SUBDIR build/cmake + ) + set(ZSTD_BUILD_SHARED OFF) + set(ZSTD_BUILD_PROGRAMS OFF) + set(ZSTD_BUILD_TESTS OFF) + set(ZSTD_LEGACY_SUPPORT OFF) + FetchContent_MakeAvailable(zstd) + set(SPZ_ZSTD_LIB libzstd_static) +endif() + +add_library(spz STATIC load-spz.cc splat-types.cc) +add_library(spz::spz ALIAS spz) +set_target_properties(spz PROPERTIES + CXX_STANDARD 17 + CXX_STANDARD_REQUIRED ON + POSITION_INDEPENDENT_CODE ON +) +target_include_directories(spz PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) +target_link_libraries(spz PRIVATE ${ZLIB_LIB} ${SPZ_ZSTD_LIB}) diff --git a/vendor/spz/LICENSE b/vendor/spz/LICENSE new file mode 100644 index 0000000..0a161eb --- /dev/null +++ b/vendor/spz/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2024 Niantic Labs + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/vendor/spz/load-spz.cc b/vendor/spz/load-spz.cc new file mode 100644 index 0000000..59930fb --- /dev/null +++ b/vendor/spz/load-spz.cc @@ -0,0 +1,1484 @@ +/* +MIT License + +Copyright (c) 2025 Niantic Labs +Copyright (c) 2025 Adobe Inc. + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +*/ + +#include "load-spz.h" +#include "splat-utils.h" +#include "splat-types.h" +#ifdef SPZ_BUILD_EXTENSIONS +#include "splat-extensions.h" +#endif + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace spz { + +namespace { + +uint8_t toUint8(float x) { return static_cast(std::clamp(std::round(x), 0.0f, 255.0f)); } + +// Quantizes to 8 bits, then rounds to nearest bucket center. 0 always maps to a bucket center. +uint8_t quantizeSH(float x, int32_t bucketSize) { + int32_t q = static_cast(std::round(x * 128.0f) + 128.0f); + q = (q + bucketSize / 2) / bucketSize * bucketSize; + return static_cast(std::clamp(q, 0, 255)); +} + +float sigmoid(float x) { return 1 / (1 + std::exp(-x)); } + +template +size_t countBytes(std::vector vec) { + return vec.size() * sizeof(vec[0]); +} + +#define CHECK(x) \ + { \ + if (!(x)) { \ + SpzLog("[SPZ: ERROR] Check failed: %s:%d: %s", __FILE__, __LINE__, #x); \ + return false; \ + } \ + } + +#define CHECK_GE(x, y) CHECK((x) >= (y)) +#define CHECK_LE(x, y) CHECK((x) <= (y)) +#define CHECK_EQ(x, y) CHECK((x) == (y)) + +bool checkSizes(const GaussianCloud &g) { + CHECK_GE(g.numPoints, 0); + CHECK_GE(g.shDegree, 0); + CHECK_LE(g.shDegree, SH_MAX_DEGREE); + CHECK_EQ(g.positions.size(), g.numPoints * 3); + CHECK_EQ(g.scales.size(), g.numPoints * 3); + CHECK_EQ(g.rotations.size(), g.numPoints * 4); + CHECK_EQ(g.alphas.size(), g.numPoints); + CHECK_EQ(g.colors.size(), g.numPoints * 3); + CHECK_EQ(g.sh.size(), g.numPoints * dimForDegree(g.shDegree) * 3); + return true; +} + +bool checkSizes(const PackedGaussians &packed, int32_t numPoints, int32_t shDim, bool usesFloat16) { + CHECK_EQ(packed.positions.size(), numPoints * 3 * (usesFloat16 ? 2 : 3)); + CHECK_EQ(packed.scales.size(), numPoints * 3); + CHECK_EQ(packed.rotations.size(), numPoints * (packed.usesQuaternionSmallestThree ? 4 : 3)); + CHECK_EQ(packed.alphas.size(), numPoints); + CHECK_EQ(packed.colors.size(), numPoints * 3); + CHECK_EQ(packed.sh.size(), numPoints * shDim * 3); + return true; +} + +constexpr uint8_t FlagAntialiased = 0x1; +constexpr uint8_t FlagHasExtensions = 0x2; + +// Generous upper bound on the zstd compression ratio (uncompressed / compressed) used to +// sanity-check numPoints against the actual file size. Real packed-gaussian streams compress +// well below this; the bound only exists to reject headers that claim vastly more points than +// any legitimate file of the given size could contain. Effective ceiling is +// maxPoints = (fileSize * kMaxCompressionRatio) / kMinBytesPerPoint +// Examples (kMaxCompressionRatio=1024, kMinBytesPerPoint=9): +// file size max numPoints accepted +// 1 KB ~116 K +// 1 MB ~119 M +// 100 MB ~11.9 B (far above any realistic scene) +// 1 GB ~119 B (effectively unlimited) +// A 32-byte file is capped at ~3.6 K points, so a header claiming billions of points in a +// tiny file is rejected immediately. +// 64-bit so the `size * kMaxCompressionRatio` product below can't overflow on 32-bit targets +// (e.g. wasm32, where size_t is 32-bit and the multiply wraps for files >= 4 MiB). +constexpr uint64_t kMaxCompressionRatio = 1024; +constexpr uint64_t kMinBytesPerPoint = 9; // positions stream alone: 3 components * 3 bytes + +struct NgspFileHeader { + uint32_t magic = NGSP_MAGIC; + uint32_t version = LATEST_SPZ_HEADER_VERSION; + uint32_t numPoints = 0; + uint8_t shDegree = 0; + uint8_t fractionalBits = 0; + uint8_t flags = 0; + uint8_t numStreams = 0; // number of ZSTD-compressed attribute streams + uint32_t tocByteOffset = 0; // byte offset from file start to the TOC (table of contents) + uint8_t reserved[12] = {}; +}; +static_assert(sizeof(NgspFileHeader) == 32, "NgspFileHeader must be 32 bytes"); + +// TODO: After v4 is released, move legacy logic to separate files. +// Legacy 16-byte header used in gzip single-stream files (pre-v4). Read-only path. +struct LegacyPackedGaussiansHeader { + uint32_t magic = NGSP_MAGIC; + uint32_t version = 0; + uint32_t numPoints = 0; + uint8_t shDegree = 0; + uint8_t fractionalBits = 0; + uint8_t flags = 0; + uint8_t reserved = 0; +}; + +bool decompressGzippedImpl( + const uint8_t *compressed, size_t size, int32_t windowSize, std::vector *out) { + std::vector buffer(8192); + z_stream stream = {}; + stream.next_in = const_cast(compressed); + stream.avail_in = size; + if (inflateInit2(&stream, windowSize) != Z_OK) { + return false; + } + out->clear(); + out->reserve(size * 3); + bool success = false; + while (true) { + stream.next_out = buffer.data(); + stream.avail_out = buffer.size(); + int32_t res = inflate(&stream, Z_NO_FLUSH); + if (res != Z_OK && res != Z_STREAM_END) { + break; + } + out->insert(out->end(), buffer.data(), buffer.data() + buffer.size() - stream.avail_out); + if (res == Z_STREAM_END) { + success = true; + break; + } + } + inflateEnd(&stream); + return success; +} + +bool decompressGzipped(const uint8_t *compressed, size_t size, std::vector *out) { + // Here 16 means enable automatic gzip header detection; consider switching this to 32 to enable + // both automated gzip and zlib header detection. + return decompressGzippedImpl(compressed, size, 16 | MAX_WBITS, out); +} + +// A read-only streambuf over a contiguous byte range, avoiding any copy. +struct membuf : std::streambuf { + membuf(const uint8_t *data, size_t size) { + auto *p = reinterpret_cast(const_cast(data)); + setg(p, p, p + size); + } +}; + +} // namespace + +bool compressGzipped(const uint8_t *data, size_t size, std::vector *out) { + std::vector buffer(8192); + z_stream stream = {}; + if (deflateInit2(&stream, Z_DEFAULT_COMPRESSION, Z_DEFLATED, 16 + MAX_WBITS, 9, Z_DEFAULT_STRATEGY) + != Z_OK) { + return false; + } + out->clear(); + out->reserve(size / 4); + stream.next_in = const_cast(reinterpret_cast(data)); + stream.avail_in = static_cast(size); + bool success = false; + while (true) { + stream.next_out = buffer.data(); + stream.avail_out = static_cast(buffer.size()); + int32_t res = deflate(&stream, Z_FINISH); + if (res != Z_OK && res != Z_STREAM_END) { + break; + } + out->insert(out->end(), buffer.data(), buffer.data() + buffer.size() - stream.avail_out); + if (res == Z_STREAM_END) { + success = true; + break; + } + } + deflateEnd(&stream); + return success; +} + +bool compressZstd(const uint8_t *data, size_t size, std::vector *out, + int compressionLevel = 12) { + size_t const bound = ZSTD_compressBound(size); + out->resize(bound); + + ZSTD_CCtx *cctx = ZSTD_createCCtx(); + if (!cctx) return false; + ZSTD_CCtx_setParameter(cctx, ZSTD_c_compressionLevel, compressionLevel); + + size_t const compressedSize = ZSTD_compress2(cctx, out->data(), bound, data, size); + ZSTD_freeCCtx(cctx); + + if (ZSTD_isError(compressedSize)) return false; + out->resize(compressedSize); + return true; +} + +// Backward compatibility function for version 2. In version 2, rotations are represented as the +// (x, y, z) components of the normalized rotation quaternion, with each component encoded as an +// 8-bit signed integer. Version 3+ uses packQuaternionSmallestThree for better accuracy. +void packQuaternionFirstThree(uint8_t r[3], const float rotation[4], const CoordinateConverter& c) { + // Normalize the quaternion, make w positive, then store xyz. w can be derived from xyz. + // NOTE: These are already in xyzw order. + Quat4f q = normalized(quat4f(rotation)); + if (c.rotFlipQFunc) { + c.rotFlipQFunc(q.data()); + } else { + q[0] *= c.flipQ[0]; + q[1] *= c.flipQ[1]; + q[2] *= c.flipQ[2]; + } + q = times(q, (q[3] < 0 ? -127.5f : 127.5f)); + q = plus(q, Quat4f{127.5f, 127.5f, 127.5f, 127.5f}); + r[0] = toUint8(q[0]); + r[1] = toUint8(q[1]); + r[2] = toUint8(q[2]); +} + +void packQuaternionSmallestThree(uint8_t r[4], const float rotation[4], const CoordinateConverter& c) { + // Normalize the quaternion + Quat4f q = normalized(quat4f(&rotation[0])); + if (c.rotFlipQFunc) { c.rotFlipQFunc(q.data()); } + else { q[0] *= c.flipQ[0]; q[1] *= c.flipQ[1]; q[2] *= c.flipQ[2]; } + // Find largest component + unsigned iLargest = 0; + for (unsigned i = 1; i < 4; ++i) + { + if (std::abs(q[i]) > std::abs(q[iLargest])) + { + iLargest = i; + } + } + + // since -q represents the same rotation as q, transform the quaternion so the largest element + // is positive. This avoids having to send its sign bit. + unsigned negate = q[iLargest] < 0; + + // Do compression using sign bit and 9-bit precision per element. + uint32_t comp = iLargest; + for (unsigned i = 0; i < 4; ++i) + { + if (i != iLargest) + { + uint32_t negbit = (q[i] < 0) ^ negate; + uint32_t mag = + (uint32_t)(float((1u << 9u) - 1u) * (std::fabs(q[i]) / sqrt1_2) + 0.5f); + comp = (comp << 10u) | (negbit << 9u) | mag; + } + } + + // Ensure little-endianness on all platforms + r[0] = comp & 0xff; + r[1] = (comp >> 8) & 0xff; + r[2] = (comp >> 16) & 0xff; + r[3] = (comp >> 24) & 0xff; +} + +PackedGaussians packGaussians(const GaussianCloud &g, const PackOptions &o) { + if (!checkSizes(g)) { + return {}; + } + + // Validate SH quantization bit parameters + if (o.sh1Bits > 8 || o.shRestBits > 8 || o.sh1Bits < 1 || o.shRestBits < 1) { + SpzLog("[SPZ ERROR] SH quantization bits cannot exceed 8 or be less than 1 (sh1Bits=%d, shRestBits=%d)", + o.sh1Bits, o.shRestBits); + return {}; + } + + const int32_t numPoints = g.numPoints; + const int32_t shDim = dimForDegree(g.shDegree); + if (o.version < 3 && g.shDegree > 3) { + SpzLog("[SPZ WARNING] SPZ with SH degrees %d will not be loadable in a legacy loader of version %d", + g.shDegree, o.version); + } +#ifdef SPZ_BUILD_EXTENSIONS + // If the cloud carries a coordinate-system extension, use its value as the storage target + // instead of RUB so that callers can persist data in any coordinate system they choose. + CoordinateSystem packToCoord = getPackedCoordinateSystem(g.extensions); + CoordinateConverter c = coordinateConverter(o.from, packToCoord, g.shDegree); +#else + CoordinateConverter c = coordinateConverter(o.from, CoordinateSystem::RUB, g.shDegree); +#endif + + // Use 12 bits for the fractional part of coordinates (~0.25 millimeter resolution). In the future + // we can use different values on a per-splat basis and still be compatible with the decoder. + PackedGaussians packed; + packed.version = o.version; + packed.numPoints = g.numPoints; + packed.shDegree = g.shDegree; + packed.fractionalBits = 12; + packed.antialiased = g.antialiased; + // Turn off quaternion-smallest-three for backward compatibility, since version 2 does not + // support it. + packed.usesQuaternionSmallestThree = o.version >= MIN_SMALLEST_THREE_QUATERNIONS_VERSION; + + packed.rotations.resize(numPoints * (packed.usesQuaternionSmallestThree ? 4 : 3)); + packed.positions.resize(numPoints * 3 * 3); + packed.scales.resize(numPoints * 3); + packed.alphas.resize(numPoints); + packed.colors.resize(numPoints * 3); + packed.sh.resize(numPoints * shDim * 3); + +#ifdef SPZ_BUILD_EXTENSIONS + packed.extensions = g.extensions; +#endif + + // Store coordinates as 24-bit fixed point values. + const float scale = (1 << packed.fractionalBits); + std::array bufPos = {}; + for (int32_t pi = 0; pi < numPoints; ++pi) { + const size_t base = static_cast(pi) * 3; + bufPos[0] = g.positions[base + 0] * scale; + bufPos[1] = g.positions[base + 1] * scale; + bufPos[2] = g.positions[base + 2] * scale; + if (c.rotFlipPFunc) { c.rotFlipPFunc(bufPos.data()); } + else { bufPos[0] *= c.flipP[0]; bufPos[1] *= c.flipP[1]; bufPos[2] *= c.flipP[2]; } + for (size_t j = 0; j < 3; ++j) { + const int32_t fixed32 = + static_cast(std::round(bufPos[j])); + packed.positions[(base + j) * 3 + 0] = fixed32 & 0xff; + packed.positions[(base + j) * 3 + 1] = (fixed32 >> 8) & 0xff; + packed.positions[(base + j) * 3 + 2] = (fixed32 >> 16) & 0xff; + } + } + + for (size_t i = 0; i < numPoints * 3; i++) { + packed.scales[i] = toUint8((g.scales[i] + 10.0f) * 16.0f); + } + + if (packed.usesQuaternionSmallestThree) { + for (size_t i = 0; i < numPoints; i++) + { + packQuaternionSmallestThree(&packed.rotations[4 * i], &g.rotations[4 * i], c); + } + } else { + for (size_t i = 0; i < numPoints; i++) + { + packQuaternionFirstThree(&packed.rotations[3 * i], &g.rotations[4 * i], c); + } + } + + for (size_t i = 0; i < numPoints; i++) { + // Apply sigmoid activation to alpha + packed.alphas[i] = toUint8(sigmoid(g.alphas[i]) * 255.0f); + } + + for (size_t i = 0; i < numPoints * 3; i++) { + // Convert SH DC component to wide RGB (allowing values that are a bit above 1 and below 0). + packed.colors[i] = toUint8(g.colors[i] * (colorScale * 255.0f) + (0.5f * 255.0f)); + } + + if (g.shDegree > 0) { + // Use configurable spherical harmonics quantization parameters from PackOptions. + // Quantization reduces information entropy for better g-zipping compression. + // Note: Unpacking doesn't need these bits since g-unzipping fills zero bits automatically. + const uint8_t sh1Bits = o.sh1Bits; + const uint8_t shRestBits = o.shRestBits; + const int32_t shPerPoint = dimForDegree(g.shDegree) * 3; + std::array bufSh = {}; + for (size_t i = 0; i < numPoints * shPerPoint; i += shPerPoint) { + for (size_t channel = 0; channel < 3; channel++) { + for (size_t k = 0; k < static_cast(shDim); ++k) { + bufSh[k] = g.sh[i + k * 3 + channel]; + } + for (size_t band = 0; band < static_cast(g.shDegree) && band < static_cast(SH_MAX_DEGREE); ++band) { + if (c.rotFlipShFuncs[band]) { c.rotFlipShFuncs[band](bufSh.data() + band * (band + 2)); } + } + for (size_t k = 0; k < static_cast(shDim); ++k) { bufSh[k] *= c.flipSh[k]; } + size_t j = 0, k = 0; + for (; j < 9; j += 3, k++) { // degree-1: 3 coefficients × 3 RGB channels = 9 slots + packed.sh[i + j + channel] = quantizeSH(bufSh[k], 1 << (8 - sh1Bits)); + } + for (; j < shPerPoint; j += 3, k++) { + packed.sh[i + j + channel] = quantizeSH(bufSh[k], 1 << (8 - shRestBits)); + } + } + } + } + + return packed; +} + +UnpackedGaussian PackedGaussian::unpack( + bool usesFloat16, bool usesQuaternionSmallestThree, int32_t fractionalBits, const CoordinateConverter &c) const { + UnpackedGaussian result; + if (usesFloat16) { + // Decode legacy float16 format. We can remove this at some point as it was never released. + const auto *halfData = reinterpret_cast(position.data()); + for (size_t i = 0; i < 3; i++) { + result.position[i] = halfToFloat(halfData[i]); + } + } else { + // Decode 24-bit fixed point coordinates + float scale = 1.0 / (1 << fractionalBits); + for (size_t i = 0; i < 3; i++) { + int32_t fixed32 = position[i * 3 + 0]; + fixed32 |= position[i * 3 + 1] << 8; + fixed32 |= position[i * 3 + 2] << 16; + fixed32 |= (fixed32 & 0x800000) ? 0xff000000 : 0; // sign extension + result.position[i] = static_cast(fixed32) * scale; + } + } + if (c.rotFlipPFunc) { c.rotFlipPFunc(result.position.data()); } + else { for (size_t i = 0; i < 3; i++) result.position[i] *= c.flipP[i]; } + + for (size_t i = 0; i < 3; i++) { + result.scale[i] = (scale[i] / 16.0f - 10.0f); + } + + if (usesQuaternionSmallestThree) + { + unpackQuaternionSmallestThree(&result.rotation[0], &rotation[0], c); + } + else + { + unpackQuaternionFirstThree(&result.rotation[0], &rotation[0], c); + } + + result.alpha = invSigmoid(alpha / 255.0f); + + for (size_t i = 0; i < 3; i++) { + result.color[i] = ((color[i] / 255.0f) - 0.5f) / colorScale; + } + + for (size_t i = 0; i < SH_MAX_COEFFS; i++) { + result.shR[i] = unquantizeSH(shR[i]); + result.shG[i] = unquantizeSH(shG[i]); + result.shB[i] = unquantizeSH(shB[i]); + } + + if (c.rotFlipShFuncs[0]) { + for (size_t i = 0; i < SH_MAX_DEGREE; i++) { + if (!c.rotFlipShFuncs[i]) { break; } + const size_t baseIndex = i * (i + 2); + c.rotFlipShFuncs[i](result.shR.data() + baseIndex); + c.rotFlipShFuncs[i](result.shG.data() + baseIndex); + c.rotFlipShFuncs[i](result.shB.data() + baseIndex); + } + } else { + for (size_t i = 0; i < SH_MAX_COEFFS; i++) { + result.shR[i] *= c.flipSh[i]; + result.shG[i] *= c.flipSh[i]; + result.shB[i] *= c.flipSh[i]; + } + } + + return result; +} + +PackedGaussian PackedGaussians::at(int32_t i) const { + PackedGaussian result; + int32_t positionBytes = usesFloat16() ? 6 : 9; + int32_t start3 = i * 3; + const auto *p = &positions[i * positionBytes]; + std::copy(p, p + positionBytes, result.position.data()); + std::copy(&scales[start3], &scales[start3] + 3, result.scale.data()); + int32_t rotationBytes = usesQuaternionSmallestThree ? 4 : 3; + const auto& r = &rotations[i * rotationBytes]; + std::copy(r, r + rotationBytes, result.rotation.data()); + std::copy(&colors[start3], &colors[start3] + 3, result.color.data()); + result.alpha = alphas[i]; + + int32_t shDim = dimForDegree(shDegree); + const auto *sh = &this->sh[i * shDim * 3]; + for (int32_t j = 0; j < shDim; ++j, sh += 3) { + result.shR[j] = sh[0]; + result.shG[j] = sh[1]; + result.shB[j] = sh[2]; + } + for (int32_t j = shDim; j < SH_MAX_COEFFS; ++j) { + result.shR[j] = 128; + result.shG[j] = 128; + result.shB[j] = 128; + } + + return result; +} + +UnpackedGaussian PackedGaussians::unpack(int32_t i, const CoordinateConverter &c) const { + return at(i).unpack(usesFloat16(), usesQuaternionSmallestThree, fractionalBits, c); +} + +bool PackedGaussians::usesFloat16() const { return positions.size() == numPoints * 3 * 2; } + +GaussianCloud unpackGaussians(const PackedGaussians &packed, const UnpackOptions &o) { + const int32_t numPoints = packed.numPoints; + const int32_t shDim = dimForDegree(packed.shDegree); + const bool usesFloat16 = packed.usesFloat16(); + const bool usesQuaternionSmallestThree = packed.usesQuaternionSmallestThree; + if (!checkSizes(packed, numPoints, shDim, usesFloat16)) { + return {}; + } + + GaussianCloud result; + result.numPoints = packed.numPoints; + result.shDegree = packed.shDegree; + result.antialiased = packed.antialiased; + +#ifdef SPZ_BUILD_EXTENSIONS + // Copy all extensions from PackedGaussians to GaussianCloud. + // Note: Some extensions (like SH quantization) are only used during packing and may not + // be needed in the unpacked cloud, but we preserve them for metadata completeness + // and future extensibility. + result.extensions = packed.extensions; +#endif + + result.positions.resize(numPoints * 3); + result.scales.resize(numPoints * 3); + result.rotations.resize(numPoints * 4); + result.alphas.resize(numPoints); + result.colors.resize(numPoints * 3); + result.sh.resize(numPoints * shDim * 3); + + if (usesFloat16) { + // Decode legacy float16 format. We can remove this at some point as it was never released. + const auto *halfData = reinterpret_cast(packed.positions.data()); + for (size_t i = 0; i < numPoints * 3; i++) { + result.positions[i] = halfToFloat(halfData[i]); + } + } else { + // Decode 24-bit fixed point coordinates + float scale = 1.0 / (1 << packed.fractionalBits); + for (size_t i = 0; i < numPoints * 3; i++) { + int32_t fixed32 = packed.positions[i * 3 + 0]; + fixed32 |= packed.positions[i * 3 + 1] << 8; + fixed32 |= packed.positions[i * 3 + 2] << 16; + fixed32 |= (fixed32 & 0x800000) ? 0xff000000 : 0; // sign extension + result.positions[i] = static_cast(fixed32) * scale; + } + } + + for (size_t i = 0; i < numPoints * 3; i++) { + result.scales[i] = packed.scales[i] / 16.0f - 10.0f; + } + + for (size_t i = 0; i < numPoints; i++) { + if (usesQuaternionSmallestThree) { + unpackQuaternionSmallestThree(&result.rotations[4 * i], &packed.rotations[4 * i]); + } else { + unpackQuaternionFirstThree(&result.rotations[4 * i], &packed.rotations[3 * i]); + } + } + + for (size_t i = 0; i < numPoints; i++) { + result.alphas[i] = invSigmoid(packed.alphas[i] / 255.0f); + } + + for (size_t i = 0; i < numPoints * 3; i++) { + result.colors[i] = ((packed.colors[i] / 255.0f) - 0.5f) / colorScale; + } + + for (size_t i = 0; i < packed.sh.size(); i++) { + result.sh[i] = unquantizeSH(packed.sh[i]); + } + +#ifdef SPZ_BUILD_EXTENSIONS + { + CoordinateSystem fromCoord = getPackedCoordinateSystem(result.extensions); + result.convertCoordinates(fromCoord, o.to); + } +#else + if (packed.hadSkippedExtensions) { + SpzLog("[SPZ WARNING] unpackGaussians: extensions were skipped at load time — " + "unpacked data may be incorrect due to unknown packing behavior; " + "build with SPZ_BUILD_EXTENSIONS to ensure correct results"); + } + result.convertCoordinates(CoordinateSystem::RUB, o.to); +#endif + return result; +} + +void serializePackedGaussians(const PackedGaussians &packed, std::ostream *out) { + LegacyPackedGaussiansHeader header; + header.version = packed.version; + header.numPoints = static_cast(packed.numPoints); + header.shDegree = static_cast(packed.shDegree); + header.fractionalBits = static_cast(packed.fractionalBits); + header.flags = static_cast(packed.antialiased ? FlagAntialiased : 0) +#ifdef SPZ_BUILD_EXTENSIONS + | static_cast(packed.extensions.empty() ? 0 : FlagHasExtensions) +#endif + ; + out->write(reinterpret_cast(&header), sizeof(header)); + + out->write(reinterpret_cast(packed.positions.data()), countBytes(packed.positions)); + out->write(reinterpret_cast(packed.alphas.data()), countBytes(packed.alphas)); + out->write(reinterpret_cast(packed.colors.data()), countBytes(packed.colors)); + out->write(reinterpret_cast(packed.scales.data()), countBytes(packed.scales)); + out->write(reinterpret_cast(packed.rotations.data()), countBytes(packed.rotations)); + out->write(reinterpret_cast(packed.sh.data()), countBytes(packed.sh)); + + // Write extensions at the end +#ifdef SPZ_BUILD_EXTENSIONS + writeAllExtensions(packed.extensions, *out); +#endif +} + +// Decompresses the NGSP attribute streams directly into the caller-provided destination buffers, +// avoiding any intermediate combined-buffer allocation. dests must have exactly header.numStreams +// entries, each pre-sized to the expected uncompressed byte count for that stream. +bool decompressNgspStreams(const uint8_t *data, size_t size, + const NgspFileHeader &header, + const std::vector> &dests) { + if (header.tocByteOffset < sizeof(NgspFileHeader)) { + SpzLog("[SPZ ERROR] decompressNgspStreams: TOC byte offset is less than the size of the header"); + return false; + } + const size_t tocSize = header.numStreams * 2 * sizeof(uint64_t); + const size_t tocEnd = header.tocByteOffset + tocSize; + if (tocEnd > size) { + SpzLog("[SPZ ERROR] decompressNgspStreams: TOC end is greater than the size of the data"); + return false; + } + if (dests.size() != header.numStreams) { + SpzLog("[SPZ ERROR] decompressNgspStreams: stream count mismatch"); + return false; + } + + struct StreamInfo { + uint64_t compressedSize; + uint64_t uncompressedSize; + size_t compressedOffset; + }; + std::vector infos(header.numStreams); + size_t compressedOffset = tocEnd; + for (uint8_t i = 0; i < header.numStreams; i++) { + const size_t e = header.tocByteOffset + i * 16; + std::memcpy(&infos[i].compressedSize, data + e, sizeof(uint64_t)); + std::memcpy(&infos[i].uncompressedSize, data + e + sizeof(uint64_t), sizeof(uint64_t)); + infos[i].compressedOffset = compressedOffset; + if (infos[i].compressedSize > size - compressedOffset) { + SpzLog("[SPZ ERROR] decompressNgspStreams: stream extends past end of data"); + return false; + } + compressedOffset += infos[i].compressedSize; + if (infos[i].uncompressedSize != dests[i].second) { + SpzLog("[SPZ ERROR] decompressNgspStreams: stream size mismatch"); + return false; + } + } + if (compressedOffset != size) { + SpzLog("[SPZ ERROR] decompressNgspStreams: compressed data size mismatch"); + return false; + } + +#if defined(__EMSCRIPTEN__) + // TODO: Add support for parallel decompression on WASM. + for (size_t i = 0; i < infos.size(); i++) { + const size_t ret = ZSTD_decompress( + dests[i].first, dests[i].second, + data + infos[i].compressedOffset, infos[i].compressedSize); + if (ZSTD_isError(ret) || ret != dests[i].second) { + SpzLog("[SPZ ERROR] decompressNgspStreams: ZSTD decompression failed"); + return false; + } + } +#else + std::vector> futures; + for (size_t i = 0; i < infos.size(); i++) { + const auto info = infos[i]; + const auto dest = dests[i]; + futures.push_back(std::async(std::launch::async, [data, info, dest]() -> bool { + const size_t ret = ZSTD_decompress( + dest.first, dest.second, + data + info.compressedOffset, info.compressedSize); + if (ZSTD_isError(ret) || ret != dest.second) { + SpzLog("[SPZ ERROR] decompressNgspStreams: ZSTD decompression failed"); + return false; + } + return true; + })); + } + for (auto &f : futures) { + if (!f.get()) return false; + } +#endif + return true; +} + +bool compressNgspStreams(const std::vector> &srcs, + std::vector> *chunks, + std::vector *uncompressedSizes) { +#if defined(__EMSCRIPTEN__) + // TODO: Add support for parallel compression on WASM. + for (const auto &s : srcs) { + if (s.second == 0) continue; + uncompressedSizes->push_back(s.second); + std::vector chunk; + if (!compressZstd(s.first, s.second, &chunk)) { + SpzLog("[SPZ ERROR] compressNgspStreams: ZSTD compression failed"); + return false; + } + chunks->push_back(std::move(chunk)); + } +#else + std::vector>> futures; + for (const auto &s : srcs) { + if (s.second == 0) continue; + uncompressedSizes->push_back(s.second); + futures.push_back(std::async(std::launch::async, [s]() -> std::vector { + std::vector chunk; + if (!compressZstd(s.first, s.second, &chunk)) { + SpzLog("[SPZ ERROR] compressNgspStreams: ZSTD compression failed"); + return {}; + } + return chunk; + })); + } + for (auto &f : futures) { + chunks->push_back(f.get()); + if (chunks->back().empty()) { + SpzLog("[SPZ ERROR] compressNgspStreams: compression failed"); + return false; + } + } +#endif + return true; +} + +PackedGaussians loadPackedGaussiansFromNgsp(const uint8_t *data, size_t size, + const NgspFileHeader &header) { + const int32_t numPoints = header.numPoints; + const int32_t shDim = dimForDegree(header.shDegree); + const bool usesQuaternionSmallestThree = header.version >= MIN_SMALLEST_THREE_QUATERNIONS_VERSION; + + PackedGaussians result; + result.version = header.version; + result.numPoints = numPoints; + result.shDegree = header.shDegree; + result.fractionalBits = header.fractionalBits; + result.antialiased = (header.flags & FlagAntialiased) != 0; + result.usesQuaternionSmallestThree = usesQuaternionSmallestThree; + + // Pre-size attribute vectors so decompressNgspStreams can write directly into them, + // avoiding both the intermediate combined buffer and the readChunk copies. + result.positions.resize(static_cast(numPoints) * 9); + result.alphas.resize(static_cast(numPoints)); + result.colors.resize(static_cast(numPoints) * 3); + result.scales.resize(static_cast(numPoints) * 3); + result.rotations.resize(static_cast(numPoints) * (usesQuaternionSmallestThree ? 4u : 3u)); + result.sh.resize(static_cast(numPoints) * shDim * 3); + + // Build destination list in the same order saveSpz writes streams, skipping zero-size buffers. + std::vector> dests; + for (auto attr : kAllSplatAttributes) { + auto &v = packedBuffer(result, attr); + if (!v.empty()) dests.push_back({v.data(), v.size()}); + } + + if (!decompressNgspStreams(data, size, header, dests)) { + SpzLog("[SPZ ERROR] loadSpzPacked: NGSP stream decompression failed"); + return {}; + } + + if ((header.flags & FlagHasExtensions) != 0) { + const size_t extStart = sizeof(NgspFileHeader); + const size_t extEnd = header.tocByteOffset; + if (extStart < extEnd && extEnd <= size) { +#ifdef SPZ_BUILD_EXTENSIONS + std::string extStr(reinterpret_cast(data + extStart), extEnd - extStart); + std::istringstream extStream(std::move(extStr)); + readAllExtensions(extStream, result.extensions); +#else + SpzLog("[SPZ WARNING] loadSpzPacked: file has extensions but extension support is disabled — " + "skipped extensions may affect how data was packed or will be unpacked; " + "build with SPZ_BUILD_EXTENSIONS to ensure correct results"); + result.hadSkippedExtensions = true; +#endif + } + } + + return result; +} + +PackedGaussians deserializePackedGaussians(std::istream &in) { + LegacyPackedGaussiansHeader header; + in.read(reinterpret_cast(&header), sizeof(header)); + if (!in || header.magic != NGSP_MAGIC) { + SpzLog("[SPZ ERROR] deserializePackedGaussians: header not found"); + return {}; + } + if (header.version < 1 || header.version > LATEST_SPZ_HEADER_VERSION) { + SpzLog("[SPZ ERROR] deserializePackedGaussians: version not supported: %d", header.version); + return {}; + } + // Bound numPoints against the bytes actually remaining in the decompressed legacy stream. + // Legacy files are already gzip-decompressed at this point, so the ratio of header-claimed + // points to remaining bytes can't exceed ~1 point per kMinBytesPerPoint bytes. + size_t remaining = 0; + { + const std::streampos cur = in.tellg(); + if (cur != std::streampos(-1)) { + in.seekg(0, std::ios::end); + const std::streampos end = in.tellg(); + in.seekg(cur); + if (end != std::streampos(-1)) remaining = static_cast(end - cur); + } + } + // Same INT32_MAX cap as the NGSP path: numPoints is consumed as int32_t below, so values + // above INT32_MAX would wrap negative. This also bounds numPoints when the size probe above + // fails (remaining == 0), in which case the ratio check is skipped. + if (header.numPoints == 0 || + header.numPoints > static_cast(std::numeric_limits::max()) || + (remaining > 0 && + static_cast(header.numPoints) > remaining / kMinBytesPerPoint)) { + SpzLog("[SPZ ERROR] deserializePackedGaussians: invalid point count: %u", header.numPoints); + return {}; + } + if (header.shDegree > SH_MAX_DEGREE) { + SpzLog("[SPZ ERROR] deserializePackedGaussians: Unsupported SH degree: %d", header.shDegree); + return {}; + } + SpzLog( + "[SPZ] deserializePackedGaussians: version=%d, numPoints=%d, shDegree=%d, fractionalBits=%d, antialiased=%d, hasExtensions=%d", + header.version, + header.numPoints, + header.shDegree, + header.fractionalBits, + int((header.flags & FlagAntialiased) != 0), + int((header.flags & FlagHasExtensions) != 0) + ); + + const int32_t numPoints = header.numPoints; + const int32_t shDim = dimForDegree(header.shDegree); + const bool usesFloat16 = header.version == 1; + const bool usesQuaternionSmallestThree = header.version >= MIN_SMALLEST_THREE_QUATERNIONS_VERSION; + const bool hasExtensions = (header.flags & FlagHasExtensions) != 0; + + PackedGaussians result; + result.version = header.version; + result.numPoints = numPoints; + result.shDegree = header.shDegree; + result.fractionalBits = header.fractionalBits; + result.antialiased = (header.flags & FlagAntialiased) != 0; + result.positions.resize(numPoints * 3 * (usesFloat16 ? 2 : 3)); + result.scales.resize(numPoints * 3); + result.usesQuaternionSmallestThree = usesQuaternionSmallestThree; + result.rotations.resize(numPoints * (usesQuaternionSmallestThree ? 4 : 3)); + result.alphas.resize(numPoints); + result.colors.resize(numPoints * 3); + result.sh.resize(numPoints * shDim * 3); + in.read(reinterpret_cast(result.positions.data()), countBytes(result.positions)); + in.read(reinterpret_cast(result.alphas.data()), countBytes(result.alphas)); + in.read(reinterpret_cast(result.colors.data()), countBytes(result.colors)); + in.read(reinterpret_cast(result.scales.data()), countBytes(result.scales)); + in.read(reinterpret_cast(result.rotations.data()), countBytes(result.rotations)); + in.read(reinterpret_cast(result.sh.data()), countBytes(result.sh)); + + // Read extensions at the end + if (hasExtensions) { +#ifdef SPZ_BUILD_EXTENSIONS + readAllExtensions(in, result.extensions); +#else + SpzLog("[SPZ WARNING] deserializePackedGaussians: stream has extensions but extension support is disabled — " + "skipped extensions may affect how data was packed or will be unpacked; " + "build with SPZ_BUILD_EXTENSIONS to ensure correct results"); + result.hadSkippedExtensions = true; +#endif + } + + if (!in) { + SpzLog("[SPZ ERROR] deserializePackedGaussians: read error"); + return {}; + } + + return result; +} + +bool saveSpz(const GaussianCloud &g, const PackOptions &o, std::vector *out) { + PackedGaussians packed = packGaussians(g, o); + + if (g.numPoints > 0 && packed.numPoints == 0) { + return false; + } + + if (o.version < MIN_ZSTD_SPZ_HEADER_VERSION) { + // Legacy gzip path for versions 1–3. + std::stringstream ss; + serializePackedGaussians(packed, &ss); + const std::string data = ss.str(); + return compressGzipped(reinterpret_cast(data.data()), data.size(), out); + } + + std::vector extensionData; +#ifdef SPZ_BUILD_EXTENSIONS + if (!packed.extensions.empty()) { + std::ostringstream extStream; + writeAllExtensions(packed.extensions, extStream); + const std::string &s = extStream.str(); + extensionData.assign(s.begin(), s.end()); + } +#endif + + const std::vector> srcs = { + {packed.positions.data(), packed.positions.size()}, + {packed.alphas.data(), packed.alphas.size()}, + {packed.colors.data(), packed.colors.size()}, + {packed.scales.data(), packed.scales.size()}, + {packed.rotations.data(), packed.rotations.size()}, + {packed.sh.data(), packed.sh.size()}, + }; + + std::vector> chunks; + std::vector uncompressedSizes; + if (!compressNgspStreams(srcs, &chunks, &uncompressedSizes)) { + return false; + } + + const uint8_t numStreams = static_cast(chunks.size()); + const uint32_t tocByteOffset = static_cast(sizeof(NgspFileHeader) + extensionData.size()); + NgspFileHeader header; + header.version = o.version; + header.numPoints = packed.numPoints; + header.shDegree = packed.shDegree; + header.fractionalBits = packed.fractionalBits; + header.flags = static_cast(packed.antialiased ? FlagAntialiased : 0) +#ifdef SPZ_BUILD_EXTENSIONS + | static_cast(extensionData.empty() ? 0 : FlagHasExtensions) +#endif + ; + header.numStreams = numStreams; + header.tocByteOffset = tocByteOffset; + + // Write plaintext zone: [header][extensions][TOC] + out->resize(tocByteOffset + numStreams * 2 * sizeof(uint64_t)); + uint8_t *buf = out->data(); + std::memcpy(buf, &header, sizeof(header)); + if (!extensionData.empty()) { + std::memcpy(buf + sizeof(NgspFileHeader), extensionData.data(), extensionData.size()); + } + for (uint8_t i = 0; i < numStreams; i++) { + const uint64_t cs = static_cast(chunks[i].size()); + const uint64_t us = uncompressedSizes[i]; + const size_t e = tocByteOffset + i * 2 * sizeof(uint64_t); + std::memcpy(buf + e, &cs, sizeof(cs)); + std::memcpy(buf + e + 8, &us, sizeof(us)); + } + for (const auto &chunk : chunks) { + out->insert(out->end(), chunk.begin(), chunk.end()); + } + return true; +} + +PackedGaussians loadSpzPacked(const uint8_t *data, size_t size) { + uint32_t magic = 0; + if (size >= sizeof(magic)) std::memcpy(&magic, data, sizeof(magic)); + + if (magic == NGSP_MAGIC) { + if (size < sizeof(NgspFileHeader)) { + SpzLog("[SPZ ERROR] loadSpzPacked: NGSP file too short"); + return {}; + } + + NgspFileHeader header; + std::memcpy(&header, data, sizeof(header)); + + if (header.version < MIN_ZSTD_SPZ_HEADER_VERSION || header.version > LATEST_SPZ_HEADER_VERSION) { + SpzLog("[SPZ ERROR] loadSpzPacked: unsupported version: %d", header.version); + return {}; + } + // numPoints is stored unsigned but consumed as int32_t downstream (loadPackedGaussiansFromNgsp), + // so reject anything above INT32_MAX up front — otherwise the cast wraps negative. The + // file-size ratio check below is evaluated in 64-bit to avoid overflowing the multiply on + // 32-bit targets. + if (header.numPoints == 0 || + header.numPoints > static_cast(std::numeric_limits::max()) || + static_cast(header.numPoints) > + (static_cast(size) * kMaxCompressionRatio) / kMinBytesPerPoint) { + SpzLog("[SPZ ERROR] loadSpzPacked: invalid point count: %u (file size %zu)", + header.numPoints, size); + return {}; + } + SpzLog( + "[SPZ] loadSpzPacked (NGSP): version=%d, numPoints=%d, shDegree=%d, numStreams=%d, tocByteOffset=%d", + header.version, header.numPoints, header.shDegree, header.numStreams, header.tocByteOffset); + + return loadPackedGaussiansFromNgsp(data, size, header); + + } else if (size >= 2 && data[0] == 0x1f && data[1] == 0x8b) { + // Legacy single-stream GZip format (pre-v4). + std::vector decompressed; + if (!decompressGzipped(data, size, &decompressed)) return {}; + membuf buf(decompressed.data(), decompressed.size()); + std::istream stream(&buf); + return deserializePackedGaussians(stream); + + } else { + SpzLog("[SPZ ERROR] loadSpzPacked: unrecognized format"); + return {}; + } +} + +PackedGaussians loadSpzPacked(const std::vector &data) { + return loadSpzPacked(data.data(), data.size()); +} + +PackedGaussians loadSpzPacked(const std::string &filename) { + std::ifstream in(filename, std::ios::binary | std::ios::ate); + if (!in.good()) + return {}; + std::vector data(in.tellg()); + in.seekg(0, std::ios::beg); + in.read(reinterpret_cast(data.data()), data.size()); + if (!in.good()) { + return {}; + } + return loadSpzPacked(data); +} + +GaussianCloud loadSpz(const std::vector &data, const UnpackOptions &o) { + return unpackGaussians(loadSpzPacked(data), o); +} + +GaussianCloud loadSpz(const uint8_t *data, size_t size, const UnpackOptions &o) { + return unpackGaussians(loadSpzPacked(data, size), o); +} + +bool saveSpz(const GaussianCloud &g, const PackOptions &o, const std::string &filename) { + std::vector data; + if (!saveSpz(g, o, &data)) { + return false; + } + std::ofstream out(filename, std::ios::binary | std::ios::out); + out.write(reinterpret_cast(data.data()), data.size()); + out.close(); + return out.good(); +} + +GaussianCloud loadSpz(const std::string &filename, const UnpackOptions &o) { + std::ifstream in(filename, std::ios::binary | std::ios::ate); + if (!in.good()) { + SpzLog("[SPZ ERROR] Unable to open: %s", filename.c_str()); + return {}; + } + std::vector data(in.tellg()); + in.seekg(0, std::ios::beg); + in.read(reinterpret_cast(data.data()), data.size()); + in.close(); + if (!in.good()) { + SpzLog("[SPZ ERROR] Unable to load data from: %s", filename.c_str()); + return {}; + } + return loadSpz(data, o); +} + +bool getNextHeaderLine(std::ifstream &in, std::string &line) { + while (std::getline(in, line)) { + // Find the first non-whitespace character + size_t start = line.find_first_not_of(" \t\n\r\f\v"); + // If line is empty or whitespace-only, skip it and continue reading. + if (std::string::npos == start) { + continue; + } + // Trim leading whitespace and check for 'comment' + std::string trimmed_line = line.substr(start); + if (trimmed_line.rfind("comment", 0) == 0) { + continue; // Skip comment line + } + // Found a valid non-comment, non-empty line + line = trimmed_line; // Update the reference string with the trimmed line + return true; + } + // Failed to read a line (EOF or error) + return false; +} + +GaussianCloud loadSplatFromPly(const std::string &filename, const UnpackOptions &o) { + SpzLog("[SPZ] Loading: %s", filename.c_str()); + std::ifstream in(filename, std::ios::binary); + if (!in.good()) { + SpzLog("[SPZ ERROR] Unable to open: %s", filename.c_str()); + in.close(); + return {}; + } + std::string line; + std::getline(in, line); + if (line != "ply") { + SpzLog("[SPZ ERROR] %s: not a .ply file", filename.c_str()); + in.close(); + return {}; + } + if (!getNextHeaderLine(in, line) || line != "format binary_little_endian 1.0") { + SpzLog("[SPZ ERROR] %s: unsupported .ply format", filename.c_str()); + in.close(); + return {}; + } + if (!getNextHeaderLine(in, line) || line.find("element vertex ") != 0) { + SpzLog("[SPZ ERROR] %s: missing vertex count", filename.c_str()); + in.close(); + return {}; + } + int32_t numPoints = std::stoi(line.substr(std::strlen("element vertex "))); + if (numPoints <= 0) { + SpzLog("[SPZ ERROR] %s: invalid vertex count: %d", filename.c_str(), numPoints); + in.close(); + return {}; + } + + SpzLog("[SPZ] Loading %d points", numPoints); + std::unordered_map fields; // name -> index + + // Helper function to get property size from PLY type string + auto getPropertySize = [](const std::string& line) -> size_t { + if (line.find("property float ") == 0 || line.find("property int ") == 0 || + line.find("property uint ") == 0) { + return 4; + } else if (line.find("property double ") == 0) { + return 8; + } else if (line.find("property char ") == 0 || line.find("property uchar ") == 0) { + return 1; + } else if (line.find("property short ") == 0 || line.find("property ushort ") == 0) { + return 2; + } + return 4; // Default assumption + }; + + // Track extra elements (non-vertex) to handle their data + std::vector extraElements; + + // State machine for parsing header + enum class ParseState { IN_VERTEX, IN_EXTRA_ELEMENT }; + ParseState state = ParseState::IN_VERTEX; + std::string currentElementName; + int32_t currentElementCount = 0; + size_t currentElementBytes = 0; + bool currentElementIsKnown = false; + + for (int32_t i = 0;; i++) { + if (!getNextHeaderLine(in, line)) { + SpzLog("[SPZ ERROR] %s: unexpected EOF while reading header properties.", filename.c_str()); + in.close(); + return {}; + } + if (line == "end_header") { + // Finalize any pending extra element + if (state == ParseState::IN_EXTRA_ELEMENT && currentElementCount > 0) { + extraElements.push_back({currentElementName, currentElementCount, currentElementBytes, currentElementIsKnown}); + } + break; + } + + // Check for new element definitions (non-vertex) + if (line.find("element ") == 0 && line.find("element vertex ") != 0) { + // Finalize previous extra element if any + if (state == ParseState::IN_EXTRA_ELEMENT && currentElementCount > 0) { + extraElements.push_back({currentElementName, currentElementCount, currentElementBytes, currentElementIsKnown}); + } + + // Parse element name and count + size_t spacePos = line.find(' ', 8); // After "element " + if (spacePos != std::string::npos) { + currentElementName = line.substr(8, spacePos - 8); + currentElementCount = std::stoi(line.substr(spacePos + 1)); + currentElementBytes = 0; + + // Check if this is a known element we handle specially (via extensions) +#ifdef SPZ_BUILD_EXTENSIONS + currentElementIsKnown = isKnownPlyExtensionElement(currentElementName); +#else + currentElementIsKnown = false; +#endif + + state = ParseState::IN_EXTRA_ELEMENT; + if (!currentElementIsKnown) { + SpzLog("[SPZ] Found extra element: %s (%d items)", currentElementName.c_str(), currentElementCount); + } + } + continue; + } + + // Handle properties based on current state + if (state == ParseState::IN_EXTRA_ELEMENT) { + if (line.find("property ") == 0) { + currentElementBytes += getPropertySize(line); + } + continue; + } + + // We're in vertex element - only accept float properties + if (line.find("property float ") != 0) { + SpzLog("[SPZ ERROR] %s: unsupported vertex property type: %s", filename.c_str(), line.c_str()); + in.close(); + return {}; + } + std::string name = line.substr(std::strlen("property float ")); + fields[name] = i; + } + + // Returns the index for a given field name, ensuring the name exists. + const auto index = [&fields](const std::string &name) { + const auto &itr = fields.find(name); + if (itr == fields.end()) { + SpzLog("[SPZ ERROR] Missing field: %s", name.c_str()); + return -1; + } + return itr->second; + }; + + const std::vector positionIdx = {index("x"), index("y"), index("z")}; + const std::vector scaleIdx = {index("scale_0"), index("scale_1"), index("scale_2")}; + const std::vector rotIdx = {index("rot_1"), index("rot_2"), index("rot_3"), index("rot_0")}; + const std::vector alphaIdx = {index("opacity")}; + const std::vector colorIdx = {index("f_dc_0"), index("f_dc_1"), index("f_dc_2")}; + + // Check that only valid indices were returned. + for (auto idx : positionIdx) { + if (idx < 0) { + in.close(); + return {}; + } + } + for (auto idx : scaleIdx) { + if (idx < 0) { + in.close(); + return {}; + } + } + for (auto idx : rotIdx) { + if (idx < 0) { + in.close(); + return {}; + } + } + for (auto idx : alphaIdx) { + if (idx < 0) { + in.close(); + return {}; + } + } + for (auto idx : colorIdx) { + if (idx < 0) { + in.close(); + return {}; + } + } + + // Spherical harmonics are optional and variable in size (depending on degree) + std::vector shIdx; + const int32_t shMaxCoeffsRGB = SH_MAX_COEFFS * 3; + for (int32_t i = 0; i < shMaxCoeffsRGB; i++) { + const auto &itr = fields.find("f_rest_" + std::to_string(i)); + if (itr == fields.end()) + break; + shIdx.push_back(itr->second); + } + const int32_t shDim = static_cast(shIdx.size() / 3); + + std::vector values(numPoints * fields.size()); + in.read(reinterpret_cast(values.data()), values.size() * sizeof(float)); + if (!in.good()) { + SpzLog("[SPZ ERROR] Unable to load data from: %s", filename.c_str()); + in.close(); + return {}; + } + + GaussianCloud result; +#ifdef SPZ_BUILD_EXTENSIONS + readExtensionsFromPly(in, extraElements, result.extensions); +#endif + + // Skip data for extra elements (they appear after vertex and safe orbit data in the file) + for (const auto& elem : extraElements) { + if (elem.isKnown) continue; // Already handled above + size_t bytesToSkip = elem.count * elem.bytesPerElement; + if (bytesToSkip > 0) { + in.seekg(bytesToSkip, std::ios::cur); + SpzLog("[SPZ] Skipped %zu bytes for element '%s'", bytesToSkip, elem.name.c_str()); + } + } + + in.close(); + + result.numPoints = numPoints; + result.shDegree = degreeForDim(shDim); + result.positions.reserve(numPoints * 3); + result.scales.reserve(numPoints * 3); + result.rotations.reserve(numPoints * 4); + result.alphas.reserve(numPoints * 1); + result.colors.reserve(numPoints * 3); + for (size_t i = 0; i < values.size(); i += fields.size()) { + for (int32_t j = 0; j < positionIdx.size(); j++) { + result.positions.push_back(values[i + positionIdx[j]]); + } + for (int32_t j = 0; j < scaleIdx.size(); j++) { + result.scales.push_back(values[i + scaleIdx[j]]); + } + for (int32_t j = 0; j < rotIdx.size(); j++) { + result.rotations.push_back(values[i + rotIdx[j]]); + } + for (int32_t j = 0; j < alphaIdx.size(); j++) { + result.alphas.push_back(values[i + alphaIdx[j]]); + } + for (int32_t j = 0; j < colorIdx.size(); j++) { + result.colors.push_back(values[i + colorIdx[j]]); + } + // Convert from [N,C,S] to [N,S,C] (where C is color channel, S is SH coeff). + for (int32_t j = 0; j < shDim; j++) { + result.sh.push_back(values[i + shIdx[j]]); + result.sh.push_back(values[i + shIdx[j + shDim]]); + result.sh.push_back(values[i + shIdx[j + 2 * shDim]]); + } + } + + result.convertCoordinates(CoordinateSystem::RDF, o.to); + return result; +} + +bool saveSplatToPly(const GaussianCloud &data, const PackOptions &o, const std::string &filename) { + // Use int64_t for N so that N*D never overflows int32_t for clouds with >35M gaussians. + const int64_t N = data.numPoints; + CHECK_EQ(data.positions.size(), N * 3); + CHECK_EQ(data.scales.size(), N * 3); + CHECK_EQ(data.rotations.size(), N * 4); + CHECK_EQ(data.alphas.size(), N); + CHECK_EQ(data.colors.size(), N * 3); + const int64_t shDim = (N > 0) ? static_cast(data.sh.size() / N / 3) : 0; + const int64_t D = 17 + shDim * 3; + + CoordinateConverter c = coordinateConverter(o.from, CoordinateSystem::RDF, data.shDegree); + + std::ofstream out(filename, std::ios::binary); + if (!out.good()) { + SpzLog("[SPZ ERROR] Unable to open for writing: %s", filename.c_str()); + return false; + } + out << "ply\n"; + out << "format binary_little_endian 1.0\n"; + out << "element vertex " << N << "\n"; + out << "property float x\n"; + out << "property float y\n"; + out << "property float z\n"; + out << "property float nx\n"; + out << "property float ny\n"; + out << "property float nz\n"; + out << "property float f_dc_0\n"; + out << "property float f_dc_1\n"; + out << "property float f_dc_2\n"; + for (int64_t i = 0; i < shDim * 3; i++) { + out << "property float f_rest_" << i << "\n"; + } + out << "property float opacity\n"; + out << "property float scale_0\n"; + out << "property float scale_1\n"; + out << "property float scale_2\n"; + out << "property float rot_0\n"; + out << "property float rot_1\n"; + out << "property float rot_2\n"; + out << "property float rot_3\n"; + +#ifdef SPZ_BUILD_EXTENSIONS + writeExtensionsToPlyHeader(data.extensions, out); +#endif + + out << "end_header\n"; + + // Write in chunks of 1M points to bound peak memory regardless of cloud size. + const int64_t kChunkSize = 1'000'000; + std::vector values(std::min(kChunkSize, N) * D); + std::array bufQuat = {}; + for (int64_t start = 0; start < N; start += kChunkSize) { + const int64_t end = std::min(start + kChunkSize, N); + const int64_t rowCount = end - start; + if (rowCount * D != static_cast(values.size())) { + values.resize(rowCount * D); + } + int64_t outIdx = 0; + for (int64_t i = start; i < end; i++) { + const int64_t i3 = i * 3; + const int64_t i4 = i * 4; + // Position (x, y, z) + for (size_t j = 0; j < 3; j++) { values[outIdx + j] = data.positions[i3 + j]; } + if (c.rotFlipPFunc) { + c.rotFlipPFunc(values.data() + outIdx); + } else { + for (size_t j = 0; j < 3; j++) { values[outIdx + j] *= c.flipP[j]; } + } + outIdx += 3; + // Normals (nx, ny, nz): always zero, but some viewers expect them present + values[outIdx] = 0.0f; values[outIdx + 1] = 0.0f; values[outIdx + 2] = 0.0f; + outIdx += 3; + // Color (r, g, b): DC component for spherical harmonics + values[outIdx++] = data.colors[i3 + 0]; + values[outIdx++] = data.colors[i3 + 1]; + values[outIdx++] = data.colors[i3 + 2]; + // Spherical harmonics: Interleave so the coefficients are the fastest-changing axis and + // the channel (r, g, b) is slower-changing axis. + for (int64_t k = 0; k < 3; k++) { + for (int64_t j = 0; j < shDim; j++) { + values[outIdx + j] = data.sh[(i * shDim + j) * 3 + k]; + } + for (int32_t band = 0; band < data.shDegree && band < SH_MAX_DEGREE; ++band) { + if (c.rotFlipShFuncs[static_cast(band)]) { + c.rotFlipShFuncs[static_cast(band)](values.data() + outIdx + static_cast(band * (band + 2))); + } + } + for (int64_t j = 0; j < shDim; j++) { values[outIdx + j] *= c.flipSh[j]; } + outIdx += shDim; + } + // Alpha + values[outIdx++] = data.alphas[i]; + // Scale (sx, sy, sz) + values[outIdx++] = data.scales[i3 + 0]; + values[outIdx++] = data.scales[i3 + 1]; + values[outIdx++] = data.scales[i3 + 2]; + // Rotation (qw, qx, qy, qz) + for (int32_t j = 0; j < 4; j++) { bufQuat[j] = data.rotations[i4 + j]; } + if (c.rotFlipQFunc) { + c.rotFlipQFunc(bufQuat.data()); + } else { + for (int32_t j = 0; j < 3; j++) { bufQuat[j] *= c.flipQ[j]; } + } + // data.rotations are x,y,z,w per point; PLY expects w then x,y,z (see property order below). + values[outIdx++] = bufQuat[3]; + values[outIdx++] = bufQuat[0]; + values[outIdx++] = bufQuat[1]; + values[outIdx++] = bufQuat[2]; + } + CHECK_EQ(outIdx, rowCount * D); + out.write(reinterpret_cast(values.data()), rowCount * D * sizeof(float)); + } + +#ifdef SPZ_BUILD_EXTENSIONS + writeExtensionsToPlyData(data.extensions, out); +#endif + + out.close(); + if (!out.good()) { + SpzLog("[SPZ ERROR] Failed to write to: %s", filename.c_str()); + return false; + } + return true; +} + +bool hasExtensionSupport() { +#ifdef SPZ_BUILD_EXTENSIONS + return true; +#else + return false; +#endif +} + +} // namespace spz diff --git a/vendor/spz/load-spz.h b/vendor/spz/load-spz.h new file mode 100644 index 0000000..a20c941 --- /dev/null +++ b/vendor/spz/load-spz.h @@ -0,0 +1,195 @@ +/* +MIT License + +Copyright (c) 2025 Niantic Labs +Copyright (c) 2025 Adobe Inc. + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +*/ + +#pragma once +#include +#include +#include +#include +#include + +#include "splat-types.h" +#ifdef SPZ_BUILD_EXTENSIONS +#include "splat-extensions.h" +#endif + +#ifdef ANDROID +#include +#endif + +namespace spz { +#ifdef ANDROID +static constexpr char LOG_TAG[] = "SPZ"; +template +static void SpzLog(const char *fmt, Args &&...args) { + __android_log_print(ANDROID_LOG_INFO, LOG_TAG, fmt, std::forward(args)...); +} +#else +template +static void SpzLog(const char *fmt, Args &&...args) { + printf(fmt, std::forward(args)...); + printf("\n"); + fflush(stdout); +} +#endif // ANDROID + +template +static void SpzLog(const char *fmt) { + SpzLog("%s", fmt); +} + +constexpr int SH_MAX_COEFFS = 24; // Maximum number of SH coefficients (degree 0..4) +constexpr int DEFAULT_SH1_BITS = 5; +constexpr int DEFAULT_SH_REST_BITS = 4; + +// Latest version of the packed format, update this when changing the format. +constexpr int LATEST_SPZ_HEADER_VERSION = 4; + +// Minimum version of the ZSTD-compressed SPZ format. +constexpr int MIN_ZSTD_SPZ_HEADER_VERSION = 4; + +// Minimum version of that uses SmallestThree quaternions. +constexpr int MIN_SMALLEST_THREE_QUATERNIONS_VERSION = 3; + +// NGSP Magic Number used for spz file identification. +constexpr uint32_t NGSP_MAGIC = 0x5053474e; + +// Represents a single inflated gaussian. Each gaussian has 344 bytes (position, rotation, scale, +// color, alpha, and 24 SH coeffs x 3 channels). Although the data is easier to interpret in this +// format, it is not more precise than the packed format, since it was inflated. +struct UnpackedGaussian { + std::array position; // x, y, z + std::array rotation; // x, y, z, w + std::array scale; // std::log(scale) + std::array color; // rgb sh0 encoding + float alpha; // inverse logistic + std::array shR; + std::array shG; + std::array shB; +}; + +// Represents a single low precision gaussian. Each gaussian has exactly 92 bytes (for degree 4 +// spherical harmonics); the struct layout is fixed so that at() can index into non-interleaved buffers. +struct PackedGaussian { + std::array position{}; + std::array rotation{}; + std::array scale{}; + std::array color{}; + uint8_t alpha = 0; + std::array shR{}; + std::array shG{}; + std::array shB{}; + + UnpackedGaussian unpack( + bool usesFloat16, bool usesQuaternionSmallestThree, int32_t fractionalBits, const CoordinateConverter &c) const; +}; + +// Represents a full splat with lower precision. Each splat has at most 92 bytes (for degree 4 +// spherical harmonics), although splats with fewer spherical harmonics degrees will have less. +// The data is stored non-interleaved. +struct PackedGaussians { + uint32_t version = LATEST_SPZ_HEADER_VERSION; // Version of the packed format + int32_t numPoints = 0; // Total number of points (gaussians) + int32_t shDegree = 0; // Degree of spherical harmonics + int32_t fractionalBits = 0; // Number of bits used for fractional part of fixed-point coords + bool antialiased = false; // Whether gaussians should be rendered with mip-splat antialiasing + bool usesQuaternionSmallestThree = true; // Whether gaussians use the smallest three method to store quaternions + bool hadSkippedExtensions = false; // True when extensions were present in the file but ignored at load time + + std::vector positions; + std::vector scales; + std::vector rotations; + std::vector alphas; + std::vector colors; + std::vector sh; + +#ifdef SPZ_BUILD_EXTENSIONS + std::vector extensions; // List of extensions, if any +#endif + + bool usesFloat16() const; + PackedGaussian at(int32_t i) const; + UnpackedGaussian unpack(int32_t i, const CoordinateConverter &c) const; +}; + +struct PackOptions { + uint32_t version = LATEST_SPZ_HEADER_VERSION; // Version of the packed format + + CoordinateSystem from = CoordinateSystem::UNSPECIFIED; + + // Quantization bits are only used during packing to reduce information entropy for g-zipping. + // Unpacking doesn't need these values since g-unzipping already fills zero bits for quantized data. + uint8_t sh1Bits = DEFAULT_SH1_BITS; // Bits for SH degree 1 coefficients + uint8_t shRestBits = DEFAULT_SH_REST_BITS; // Bits for SH degree 2+ coefficients +}; + +struct UnpackOptions { + CoordinateSystem to = CoordinateSystem::UNSPECIFIED; +}; + +// Structure for PLY extra elements (non-vertex elements) +struct PlyExtraElement { + std::string name; + int32_t count; + size_t bytesPerElement; + bool isKnown; // true for elements we explicitly handle (like safe_orbit) +}; + +// Saves Gaussian splat in packed format, returning a vector of bytes. +bool saveSpz( + const GaussianCloud &gaussians, const PackOptions &options, std::vector *output); + +// Loads Gaussian splat from a vector of bytes in packed format. +GaussianCloud loadSpz(const std::vector &data, const UnpackOptions &options); + +// Loads Gaussian splat from a file / byte pointer / vector in packed format. +PackedGaussians loadSpzPacked(const std::string &filename); +PackedGaussians loadSpzPacked(const uint8_t *data, size_t size); +PackedGaussians loadSpzPacked(const std::vector &data); + +// Saves Gaussian splat in packed format to a file +bool saveSpz( + const GaussianCloud &gaussians, const PackOptions &options, const std::string &filename); + +// Loads Gaussian splat from a file in packed format +GaussianCloud loadSpz(const std::string &filename, const UnpackOptions &o); + +// Loads Gaussian splat from a byte pointer in packed format. +GaussianCloud loadSpz(const uint8_t *data, size_t size, const UnpackOptions &options); + +// Saves Gaussian splat data in .ply format +bool saveSplatToPly( + const spz::GaussianCloud &gaussians, const PackOptions &options, const std::string &filename); + +// Loads Gaussian splat data in .ply format +GaussianCloud loadSplatFromPly(const std::string &filename, const UnpackOptions &options); + +void serializePackedGaussians(const PackedGaussians &packed, std::ostream *out); + +bool compressGzipped(const uint8_t *data, size_t size, std::vector *out); + +// Returns true if the build has extension support enabled, false otherwise +bool hasExtensionSupport(); +} // namespace spz diff --git a/vendor/spz/splat-c-types.h b/vendor/spz/splat-c-types.h new file mode 100644 index 0000000..97ab5e9 --- /dev/null +++ b/vendor/spz/splat-c-types.h @@ -0,0 +1,61 @@ +/* +MIT License + +Copyright (c) 2025 Niantic Labs +Copyright (c) 2025 Adobe Inc. + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +*/ + +#ifndef SPZ_SPLAT_C_TYPES_H_ +#define SPZ_SPLAT_C_TYPES_H_ + +#define _USE_MATH_DEFINES +#include +#include +#include + +// These types are used to bridge between the C++ API and C (to interop with Swift and C#). + +typedef struct { + size_t count; + float *data; +} SpzFloatBuffer; + +// Forward declaration - full definition is in splat-extensions.h +#ifdef SPZ_BUILD_EXTENSIONS +#include "splat-extensions.h" +#else +typedef struct SpzExtensionNode SpzExtensionNode; +#endif + +typedef struct { + int32_t numPoints; + int32_t shDegree; + bool antialiased; + SpzFloatBuffer positions; + SpzFloatBuffer scales; + SpzFloatBuffer rotations; + SpzFloatBuffer alphas; + SpzFloatBuffer colors; + SpzFloatBuffer sh; + SpzExtensionNode* extensions; +} GaussianCloudData; + +#endif // SPZ_SPLAT_C_TYPES_H_ diff --git a/vendor/spz/splat-types.cc b/vendor/spz/splat-types.cc new file mode 100644 index 0000000..7369bc4 --- /dev/null +++ b/vendor/spz/splat-types.cc @@ -0,0 +1,117 @@ +/* +MIT License + +Copyright (c) 2024 Niantic Labs + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +*/ + +#include "splat-types.h" + +#include +#include + +namespace spz { + +float halfToFloat(Half h) { + auto sgn = ((h >> 15) & 0x1); + auto exponent = ((h >> 10) & 0x1f); + auto mantissa = h & 0x3ff; + + float signMul = sgn == 1 ? -1.0 : 1.0; + if (exponent == 0) { + // Subnormal numbers (no exponent, 0 in the mantissa decimal). + return signMul * std::pow(2.0f, -14.0f) * static_cast(mantissa) / 1024.0f; + } + + if (exponent == 31) { + // Infinity or NaN. + return mantissa != 0 ? std::numeric_limits::quiet_NaN() : signMul * std::numeric_limits::infinity(); + } + + // non-zero exponent implies 1 in the mantissa decimal. + return signMul * std::pow(2.0f, static_cast(exponent) - 15.0f) + * (1.0f + static_cast(mantissa) / 1024.0f); +} + +Half floatToHalf(float f) { + uint32_t f32 = *reinterpret_cast(&f); + int32_t sign = (f32 >> 31) & 0x01; // 1 bit -> 1 bit + int32_t exponent = ((f32 >> 23) & 0xff); // 8 bits -> 5 bits + int32_t mantissa = f32 & 0x7fffff; // 23 bits -> 10 bits + + // Handle inf and nan from float. + if (exponent == 0xFF) { + if (mantissa == 0) { + return (sign << 15) | 0x7C00; // Inf + } + + return (sign << 15) | 0x7C01; // Nan + } + + // If the exponent is greater than the range of half, return +/- Inf. + int32_t centeredExp = exponent - 127; + if (centeredExp > 15) { + return (sign << 15) | 0x7C00; + } + + // Normal numbers. centeredExp = [-15, 15] + if (centeredExp > -15) { + return (sign << 15) | ((centeredExp + 15) << 10) | (mantissa >> 13); + } + + // Subnormal numbers. + int32_t fullMantissa = 0x800000 | mantissa; + int32_t shift = -(centeredExp + 14); // Shift is in [-1 to -113] + int32_t newMantissa = fullMantissa >> shift; + return (sign << 15) | (newMantissa >> 13); +} + +float norm(const Vec3f &a) { + return std::sqrt(squaredNorm(a)); +} + +Vec3f normalized(const Vec3f &v) { + float n = norm(v); + return {v[0] / n, v[1] / n, v[2] / n}; +} + +float norm(const Quat4f &q) { + return std::sqrt(q[0] * q[0] + q[1] * q[1] + q[2] * q[2] + q[3] * q[3]); +} + +Quat4f axisAngleQuat(const Vec3f &scaledAxis) { + const float &a0 = scaledAxis[0]; + const float &a1 = scaledAxis[1]; + const float &a2 = scaledAxis[2]; + const float thetaSquared = a0 * a0 + a1 * a1 + a2 * a2; + // For points not at the origin, the full conversion is numerically stable. + if (thetaSquared > 0.0f) { + const float theta = std::sqrt(thetaSquared); + const float halfTheta = theta * 0.5f; + const float k = std::sin(halfTheta) / theta; + return normalized(std::array{std::cos(halfTheta), a0 * k, a1 * k, a2 * k}); + } + // If thetaSquared is 0, then we will get NaNs when dividing by theta. By approximating with a + // Taylor series, and truncating at one term, the value will be computed correctly. + const float k = 0.5f; + return normalized(Quat4f{1.0f, a0 * k, a1 * k, a2 * k}); +} + +} // namespace spz diff --git a/vendor/spz/splat-types.h b/vendor/spz/splat-types.h new file mode 100644 index 0000000..4ac5606 --- /dev/null +++ b/vendor/spz/splat-types.h @@ -0,0 +1,574 @@ +/* +MIT License + +Copyright (c) 2025 Niantic Labs +Copyright (c) 2025 Adobe Inc. + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +*/ + +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +#include "splat-c-types.h" +#ifdef SPZ_BUILD_EXTENSIONS +#include "splat-extensions.h" +#endif + +namespace spz { +constexpr int SH_MAX_DEGREE = 4; + +inline SpzFloatBuffer copyFloatBuffer(const std::vector &vector) { + SpzFloatBuffer buffer = {0, nullptr}; + if (!vector.empty()) { + buffer.count = vector.size(); + buffer.data = new float[buffer.count]; + std::memcpy(buffer.data, vector.data(), buffer.count * sizeof(float)); + } + return buffer; +} + +enum class CoordinateSystem : uint32_t { + UNSPECIFIED = 0, + LDB = 1, // Left Down Back + RDB = 2, // Right Down Back + LUB = 3, // Left Up Back + RUB = 4, // Right Up Back, Three.js coordinate system + LDF = 5, // Left Down Front + RDF = 6, // Right Down Front, PLY coordinate system + LUF = 7, // Left Up Front, GLB coordinate system + RUF = 8, // Right Up Front, Unity coordinate system + LFD = 9, // Left Front Down + RFD = 10, // Right Front Down + LFU = 11, // Left Front Up + RFU = 12, // Right Front Up + LBD = 13, // Left Back Down + RBD = 14, // Right Back Down + LBU = 15, // Left Back Up + RBU = 16, // Right Back Up +}; + +using AnalyticRotateShFn = void (*)(float*); + +enum class SplatAttribute : int32_t { + Positions, + Alphas, + Colors, + Scales, + Rotations, + Sh, +}; + +// Iteration order for the six per-point attributes. The emscripten bindings +// (`slots[]` in spz-bindings.cc) index a parallel array in lockstep — keep the +// two ordered identically. +inline constexpr std::array kAllSplatAttributes = { + SplatAttribute::Positions, SplatAttribute::Alphas, SplatAttribute::Colors, + SplatAttribute::Scales, SplatAttribute::Rotations, SplatAttribute::Sh, +}; + +struct CoordinateConverter { + std::array flipP = {1.0f, 1.0f, 1.0f}; // x, y, z flips. + std::array flipQ = {1.0f, 1.0f, 1.0f}; // x, y, z flips, w is never flipped. + std::array flipSh = // Flips for the 24 spherical harmonics coefficients. + {1.0f, 1.0f, 1.0f, 1.0f, 1.0f, + 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, + 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, + 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, + 1.0f, 1.0f, 1.0f, 1.0f}; + // Per-band SH transform; for cross-family, rotation and flip are baked in together. + // For within-family, these are empty and flipSh above is used instead. + std::array, SH_MAX_DEGREE> rotFlipShFuncs = {}; + std::function rotFlipPFunc = nullptr; + std::function rotFlipQFunc = nullptr; +}; + +constexpr std::array axesMatch(CoordinateSystem a, CoordinateSystem b) { + auto aNum = static_cast(a) - 1; + auto bNum = static_cast(b) - 1; + if (aNum < 0 || bNum < 0) { + return {true, true, true}; + } + return { + ((aNum >> 0) & 1) == ((bNum >> 0) & 1), + ((aNum >> 1) & 1) == ((bNum >> 1) & 1), + ((aNum >> 2) & 1) == ((bNum >> 2) & 1)}; +} + +constexpr bool needRotation(CoordinateSystem a, CoordinateSystem b) { + auto aNum = static_cast(a) - 1; + auto bNum = static_cast(b) - 1; + if (aNum < 0 || bNum < 0) { + return false; + } + return ((aNum >> 3) & 1) != ((bNum >> 3) & 1); +} + +// SH rotation matrices for R_x(+π/2) and R_x(-π/2). Each entry operates on the 2l+1 +// coefficients of band l (l = 1..4). The minus table is the transpose (= inverse) of the plus +// table since both are orthogonal matrices. +inline constexpr std::array kAnalyticRotatePlusPiHalfAboutXTable = { + [](float* p) { + const float t0 = p[0], t1 = p[1], t2 = p[2]; + p[0] = t1; + p[1] = -t0; + p[2] = t2; + }, + [](float* p) { + std::array s{}; + for (int i = 0; i < 5; ++i) { + s[static_cast(i)] = p[i]; + } + const float s3 = std::sqrt(3.f); + p[0] = s[3]; + p[1] = -s[1]; + p[2] = -0.5f * s[2] - (s3 / 2.f) * s[4]; + p[3] = -s[0]; + p[4] = -(s3 / 2.f) * s[2] + 0.5f * s[4]; + }, + [](float* p) { + std::array s{}; + for (int i = 0; i < 7; ++i) { + s[static_cast(i)] = p[i]; + } + const float s15 = std::sqrt(15.f); + p[0] = -std::sqrt(5.f / 8.f) * s[3] + std::sqrt(3.f / 8.f) * s[5]; + p[1] = -s[1]; + p[2] = -std::sqrt(3.f / 8.f) * s[3] - std::sqrt(5.f / 8.f) * s[5]; + p[3] = std::sqrt(5.f / 8.f) * s[0] + std::sqrt(3.f / 8.f) * s[2]; + p[4] = -0.25f * s[4] - (s15 / 4.f) * s[6]; + p[5] = -std::sqrt(3.f / 8.f) * s[0] + std::sqrt(5.f / 8.f) * s[2]; + p[6] = -(s15 / 4.f) * s[4] + 0.25f * s[6]; + }, + [](float* p) { + std::array s{}; + for (int i = 0; i < 9; ++i) { + s[static_cast(i)] = p[i]; + } + const float s2 = std::sqrt(2.f); + const float s5 = std::sqrt(5.f); + const float s7 = std::sqrt(7.f); + const float s14 = std::sqrt(14.f); + const float s35 = std::sqrt(35.f); + p[0] = -(s14 / 4.f) * s[5] + (s2 / 4.f) * s[7]; + p[1] = -0.75f * s[1] + (s7 / 4.f) * s[3]; + p[2] = -(s2 / 4.f) * s[5] - (s14 / 4.f) * s[7]; + p[3] = (s7 / 4.f) * s[1] + 0.75f * s[3]; + p[4] = (3.f / 8.f) * s[4] + (s5 / 4.f) * s[6] + (s35 / 8.f) * s[8]; + p[5] = (s14 / 4.f) * s[0] + (s2 / 4.f) * s[2]; + p[6] = (s5 / 4.f) * s[4] + 0.5f * s[6] - (s7 / 4.f) * s[8]; + p[7] = -(s2 / 4.f) * s[0] + (s14 / 4.f) * s[2]; + p[8] = (s35 / 8.f) * s[4] - (s7 / 4.f) * s[6] + 0.125f * s[8]; + }, +}; + +inline constexpr std::array kAnalyticRotateMinusPiHalfAboutXTable = { + [](float* p) { + const float t0 = p[0], t1 = p[1], t2 = p[2]; + p[0] = -t1; + p[1] = t0; + p[2] = t2; + }, + [](float* p) { + std::array s{}; + for (int i = 0; i < 5; ++i) s[static_cast(i)] = p[i]; + const float s3 = std::sqrt(3.f); + p[0] = -s[3]; + p[1] = -s[1]; + p[2] = -0.5f * s[2] - (s3 / 2.f) * s[4]; + p[3] = s[0]; + p[4] = -(s3 / 2.f) * s[2] + 0.5f * s[4]; + }, + [](float* p) { + std::array s{}; + for (int i = 0; i < 7; ++i) s[static_cast(i)] = p[i]; + const float s15 = std::sqrt(15.f); + p[0] = std::sqrt(5.f / 8.f) * s[3] - std::sqrt(3.f / 8.f) * s[5]; + p[1] = -s[1]; + p[2] = std::sqrt(3.f / 8.f) * s[3] + std::sqrt(5.f / 8.f) * s[5]; + p[3] = -std::sqrt(5.f / 8.f) * s[0] - std::sqrt(3.f / 8.f) * s[2]; + p[4] = -0.25f * s[4] - (s15 / 4.f) * s[6]; + p[5] = std::sqrt(3.f / 8.f) * s[0] - std::sqrt(5.f / 8.f) * s[2]; + p[6] = -(s15 / 4.f) * s[4] + 0.25f * s[6]; + }, + [](float* p) { + std::array s{}; + for (int i = 0; i < 9; ++i) s[static_cast(i)] = p[i]; + const float s2 = std::sqrt(2.f); + const float s5 = std::sqrt(5.f); + const float s7 = std::sqrt(7.f); + const float s14 = std::sqrt(14.f); + const float s35 = std::sqrt(35.f); + p[0] = (s14 / 4.f) * s[5] - (s2 / 4.f) * s[7]; + p[1] = -0.75f * s[1] + (s7 / 4.f) * s[3]; + p[2] = (s2 / 4.f) * s[5] + (s14 / 4.f) * s[7]; + p[3] = (s7 / 4.f) * s[1] + 0.75f * s[3]; + p[4] = (3.f / 8.f) * s[4] + (s5 / 4.f) * s[6] + (s35 / 8.f) * s[8]; + p[5] = -(s14 / 4.f) * s[0] - (s2 / 4.f) * s[2]; + p[6] = (s5 / 4.f) * s[4] + 0.5f * s[6] - (s7 / 4.f) * s[8]; + p[7] = (s2 / 4.f) * s[0] - (s14 / 4.f) * s[2]; + p[8] = (s35 / 8.f) * s[4] - (s7 / 4.f) * s[6] + 0.125f * s[8]; + }, +}; + +inline CoordinateConverter coordinateConverter(CoordinateSystem from, CoordinateSystem to, int shDegree) { + CoordinateConverter result; + + if (needRotation(from, to)) { + // A cross-family conversion decomposes into R_x(±π/2) plus a within-family flip. + // Shift `from` into the same family as `to` so the recursive call resolves to a flip only. + const bool backward = ((static_cast(from) - 1) >> 3) & 1; + const CoordinateSystem innerFrom = backward + ? static_cast(static_cast(from) - 8) // rotated→standard: shift from into standard family + : static_cast(static_cast(from) + 8); // standard→rotated: shift from into rotated family + result = coordinateConverter(innerFrom, to, shDegree); + + // Bake the inner flip into each transform function so callers never need to know the order. + // Forward (standard→rotated): rotate then flip. + // Backward (rotated→standard): flip then rotate. + const auto fp = result.flipP; + const auto fq = result.flipQ; + if (backward) { + result.rotFlipPFunc = [fp](float* p) { + p[0] *= fp[0]; p[1] *= fp[1]; p[2] *= fp[2]; + const float y = p[1], z = p[2]; p[1] = z; p[2] = -y; // R_x(-π/2) + }; + result.rotFlipQFunc = [fq](float* p) { + p[0] *= fq[0]; p[1] *= fq[1]; p[2] *= fq[2]; + const float s = std::sqrt(2.f) / 2.f; + const float x = p[0], y = p[1], z = p[2], w = p[3]; + p[0] = s * (-w + x); p[1] = s * (y + z); p[2] = s * (-y + z); p[3] = s * (w + x); + }; + } else { + result.rotFlipPFunc = [fp](float* p) { + const float y = p[1], z = p[2]; p[1] = -z; p[2] = y; // R_x(+π/2) + p[0] *= fp[0]; p[1] *= fp[1]; p[2] *= fp[2]; + }; + result.rotFlipQFunc = [fq](float* p) { + const float s = std::sqrt(2.f) / 2.f; + const float x = p[0], y = p[1], z = p[2], w = p[3]; + p[0] = s * (w + x); p[1] = s * (y - z); p[2] = s * (y + z); p[3] = s * (w - x); + p[0] *= fq[0]; p[1] *= fq[1]; p[2] *= fq[2]; + }; + } + result.flipP = {1.0f, 1.0f, 1.0f}; + result.flipQ = {1.0f, 1.0f, 1.0f}; + + const AnalyticRotateShFn* shTable = backward + ? kAnalyticRotateMinusPiHalfAboutXTable.data() + : kAnalyticRotatePlusPiHalfAboutXTable.data(); + for (int b = 0; b < shDegree && b < SH_MAX_DEGREE; ++b) { + const AnalyticRotateShFn rotFn = shTable[b]; + const size_t bandStart = static_cast(b * (b + 2)); + const size_t bandSize = static_cast(2 * b + 3); + std::array bandFlip{}; + for (size_t k = 0; k < bandSize; ++k) bandFlip[k] = result.flipSh[bandStart + k]; + if (backward) { + result.rotFlipShFuncs[static_cast(b)] = [rotFn, bandFlip, bandSize](float* p) { + for (size_t k = 0; k < bandSize; ++k) p[k] *= bandFlip[k]; + rotFn(p); + }; + } else { + result.rotFlipShFuncs[static_cast(b)] = [rotFn, bandFlip, bandSize](float* p) { + rotFn(p); + for (size_t k = 0; k < bandSize; ++k) p[k] *= bandFlip[k]; + }; + } + } + result.flipSh.fill(1.0f); + return result; + } + + auto [xMatch, yMatch, zMatch] = axesMatch(from, to); + float x = xMatch ? 1.0f : -1.0f; + float y = yMatch ? 1.0f : -1.0f; + float z = zMatch ? 1.0f : -1.0f; + + result.flipP = {x, y, z}; + result.flipQ = {y * z, x * z, x * y}; + result.flipSh = { + y, // 0 + z, // 1 + x, // 2 + x * y, // 3 + y * z, // 4 + 1.0f, // 5 + x * z, // 6 + 1.0f, // 7 + y, // 8 + x * y * z, // 9 + y, // 10 + z, // 11 + x, // 12 + z, // 13 + x, // 14 + // Used https://github.com/nerfstudio-project/gsplat/blob/main/gsplat/cuda/csrc/SphericalHarmonicsCUDA.cu + // to compute these values. + x * y, // 15 + y * z, // 16 + x * y, // 17 + y * z, // 18 + 1.0f, // 19 + x * z, // 20 + 1.0f, // 21 + x * z, // 22 + y, // 23 + }; + return result; +} + +// A point cloud composed of Gaussians. Each gaussian is represented by: +// - xyz position +// - xyz scales (on log scale, compute exp(x) to get scale factor) +// - xyzw quaternion +// - alpha (before sigmoid activation, compute sigmoid(a) to get alpha value between 0 and 1) +// - rgb color (as SH DC component, compute 0.5 + 0.282095 * x to get color value between 0 and 1) +// - 0 to 71 spherical harmonics coefficients (see comment below) +struct GaussianCloud { + // Total number of points (gaussians) in this splat. + int32_t numPoints = 0; + + // Degree of spherical harmonics for this splat. + int32_t shDegree = 0; + + // Whether the gaussians should be rendered in antialiased mode (mip splatting) + bool antialiased = false; + + // See block comment above for details + std::vector positions; + std::vector scales; + std::vector rotations; + std::vector alphas; + std::vector colors; + + // Spherical harmonics coefficients. The number of coefficients per point depends on shDegree: + // 0 -> 0 + // 1 -> 9 (3 coeffs x 3 channels) + // 2 -> 24 (8 coeffs x 3 channels) + // 3 -> 45 (15 coeffs x 3 channels) + // 4 -> 72 (24 coeffs x 3 channels) + // The color channel is the inner (fastest varying) axis, and the coefficient is the outer + // (slower varying) axis, i.e. for degree 1, the order of the 9 values is: + // sh1n1_r, sh1n1_g, sh1n1_b, sh10_r, sh10_g, sh10_b, sh1p1_r, sh1p1_g, sh1p1_b + std::vector sh; + +#ifdef SPZ_BUILD_EXTENSIONS + std::vector extensions; // List of extensions, if any +#endif + + // The caller is responsible for freeing the pointers in the returned GaussianCloudData + GaussianCloudData data() const { + GaussianCloudData data; + data.numPoints = numPoints; + data.shDegree = shDegree; + data.antialiased = antialiased; + data.positions = copyFloatBuffer(positions); + data.scales = copyFloatBuffer(scales); + data.rotations = copyFloatBuffer(rotations); + data.alphas = copyFloatBuffer(alphas); + data.colors = copyFloatBuffer(colors); + data.sh = copyFloatBuffer(sh); +#ifdef SPZ_BUILD_EXTENSIONS + data.extensions = copyExtensions(extensions); +#else + data.extensions = nullptr; +#endif + return data; + } + + // Convert between two coordinate systems, for example from RDF (ply format) to RUB (used by spz). + // This is performed in-place. + void convertCoordinates(CoordinateSystem from, CoordinateSystem to) { + if (numPoints == 0) { + // There is nothing to convert. + return; + } + CoordinateConverter c = coordinateConverter(from, to, shDegree); + // Positions: call transform (which includes flip for cross-family), else apply flipP. + if (c.rotFlipPFunc) { + for (size_t i = 0; i < positions.size(); i += 3) { + c.rotFlipPFunc(positions.data() + i); + } + } else { + for (size_t i = 0; i < positions.size(); i += 3) { + positions[i + 0] *= c.flipP[0]; + positions[i + 1] *= c.flipP[1]; + positions[i + 2] *= c.flipP[2]; + } + } + // Rotations: same pattern. Within-family only flips x,y,z; w is untouched. Cross-family rotates all four. + if (c.rotFlipQFunc) { + for (size_t i = 0; i < rotations.size(); i += 4) { + c.rotFlipQFunc(rotations.data() + i); + } + } else { + for (size_t i = 0; i < rotations.size(); i += 4) { + rotations[i + 0] *= c.flipQ[0]; + rotations[i + 1] *= c.flipQ[1]; + rotations[i + 2] *= c.flipQ[2]; + // Don't modify rotations[i + 3] (w component) + } + } + const size_t numCoeffs = sh.size() / 3; + const size_t numCoeffsPerPoint = numCoeffs / numPoints; + if (c.rotFlipShFuncs[0]) { + // Cross-family: rotation+flip baked into rotFlipShFuncs per band. + for (size_t coeffBase = 0; coeffBase < numCoeffs; coeffBase += numCoeffsPerPoint) { + for (int band = 0; band < shDegree && band < SH_MAX_DEGREE; ++band) { + const size_t bandStart = static_cast(band * (band + 2)); + const size_t bandSize = static_cast(2 * band + 3); + if (bandStart + bandSize > numCoeffsPerPoint) { break; } + std::array tmp{}; + for (int channel = 0; channel < 3; ++channel) { + for (size_t k = 0; k < bandSize; ++k) { + tmp[k] = sh[(coeffBase + bandStart + k) * 3 + static_cast(channel)]; + } + c.rotFlipShFuncs[static_cast(band)](tmp.data()); + for (size_t k = 0; k < bandSize; ++k) { + sh[(coeffBase + bandStart + k) * 3 + static_cast(channel)] = tmp[k]; + } + } + } + } + } else { + // Within-family: rotate spherical harmonics by inverting coefficients that reference the + // y and z axes, for each RGB channel. See spherical_harmonics_kernel_impl.h for spherical + // harmonics formulas. + for (size_t coeffBase = 0; coeffBase < numCoeffs; coeffBase += numCoeffsPerPoint) { + for (size_t j = 0; j < numCoeffsPerPoint; ++j) { + const size_t base = (coeffBase + j) * 3; + sh[base + 0] *= c.flipSh[j]; + sh[base + 1] *= c.flipSh[j]; + sh[base + 2] *= c.flipSh[j]; + } + } + } + } + + // Rotates the GaussianCloud by 180 degrees about the x axis (converts from RUB to RDF coordinates + // and vice versa. This is performed in-place. + void rotate180DegAboutX() { convertCoordinates(CoordinateSystem::RUB, CoordinateSystem::RDF); } + + float medianVolume() const { + if (numPoints == 0) { + return 0.01f; + } + // The volume of an ellipsoid is 4/3 * pi * x * y * z, where x, y, and z are the radii on each + // axis. Scales are stored on a log scale, and exp(x) * exp(y) * exp(z) = exp(x + y + z). So we + // can sort by value = (x + y + z) and compute volume = 4/3 * pi * exp(value) later. + std::vector scaleSums; + scaleSums.reserve(scales.size() / 3); + for (size_t i = 0; i < scales.size(); i += 3) { + scaleSums.push_back(scales[i] + scales[i + 1] + scales[i + 2]); + } + const auto mid = scaleSums.begin() + scaleSums.size() / 2; + std::nth_element(scaleSums.begin(), mid, scaleSums.end()); + return (M_PI * 4 / 3) * exp(*mid); + } +}; + +// SPZ Splat math helpers, lightweight implementations of vector and quaternion math. +using Vec3f = std::array; // x, y, z +using Quat4f = std::array; // w, x, y, z +using Half = uint16_t; + +// Half-precision helpers. +float halfToFloat(Half h); +Half floatToHalf(float f); + +// Vector helpers. +Vec3f normalized(const Vec3f &v); +float norm(const Vec3f &a); + +// Quaternion helpers. +float norm(const Quat4f &q); +inline Quat4f normalized(const Quat4f &v) { + float norm = std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2] + v[3] * v[3]); + return {v[0] / norm, v[1] / norm, v[2] / norm, v[3] / norm}; +} +Quat4f axisAngleQuat(const Vec3f &scaledAxis); + +// Constexpr helpers. +constexpr Vec3f vec3f(const float *data) { return {data[0], data[1], data[2]}; } + +constexpr float dot(const Vec3f &a, const Vec3f &b) { + return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; +} + +constexpr float squaredNorm(const Vec3f &v) { return dot(v, v); } + +constexpr Quat4f quat4f(const float *data) { return {data[0], data[1], data[2], data[3]}; } + +constexpr Vec3f times(const Quat4f &q, const Vec3f &p) { + auto [w, x, y, z] = q; + auto [vx, vy, vz] = p; + auto x2 = x + x; + auto y2 = y + y; + auto z2 = z + z; + auto wx2 = w * x2; + auto wy2 = w * y2; + auto wz2 = w * z2; + auto xx2 = x * x2; + auto xy2 = x * y2; + auto xz2 = x * z2; + auto yy2 = y * y2; + auto yz2 = y * z2; + auto zz2 = z * z2; + return { + vx * (1.0f - (yy2 + zz2)) + vy * (xy2 - wz2) + vz * (xz2 + wy2), + vx * (xy2 + wz2) + vy * (1.0f - (xx2 + zz2)) + vz * (yz2 - wx2), + vx * (xz2 - wy2) + vy * (yz2 + wx2) + vz * (1.0f - (xx2 + yy2))}; +} + +inline Quat4f times(const Quat4f &a, const Quat4f &b) { + auto [w, x, y, z] = a; + auto [qw, qx, qy, qz] = b; + return normalized(std::array{ + w * qw - x * qx - y * qy - z * qz, + w * qx + x * qw + y * qz - z * qy, + w * qy - x * qz + y * qw + z * qx, + w * qz + x * qy - y * qx + z * qw}); +} + +constexpr Quat4f times(const Quat4f &a, float s) { + return {a[0] * s, a[1] * s, a[2] * s, a[3] * s}; +} + +constexpr Quat4f plus(const Quat4f &a, const Quat4f &b) { + return {a[0] + b[0], a[1] + b[1], a[2] + b[2], a[3] + b[3]}; +} + +constexpr Vec3f times(const Vec3f &v, float s) { return {v[0] * s, v[1] * s, v[2] * s}; } + +constexpr Vec3f plus(const Vec3f &a, const Vec3f &b) { + return {a[0] + b[0], a[1] + b[1], a[2] + b[2]}; +} + +constexpr Vec3f times(const Vec3f &a, const Vec3f &b) { + return {a[0] * b[0], a[1] * b[1], a[2] * b[2]}; +} + +} // namespace spz diff --git a/vendor/spz/splat-utils.h b/vendor/spz/splat-utils.h new file mode 100644 index 0000000..443d5c3 --- /dev/null +++ b/vendor/spz/splat-utils.h @@ -0,0 +1,143 @@ +#pragma once + +#include +#include +#include + +#include "load-spz.h" +#include "splat-types.h" + +namespace spz { + +// Scale factor for DC color components. To convert to RGB, we should multiply by 0.282, but it can +// be useful to represent base colors that are out of range if the higher spherical harmonics bands +// bring them back into range so we multiply by a smaller value. +constexpr float colorScale = 0.15f; +constexpr float sqrt1_2 = 0.707106781186547524401f; + +constexpr int32_t dimForDegree(int32_t degree) { + switch (degree) { + case 0: return 0; + case 1: return 3; + case 2: return 8; + case 3: return 15; + case 4: return 24; + default: return 0; + } +} + +constexpr int32_t degreeForDim(int32_t dim) { + if (dim < 3) + return 0; + if (dim < 8) + return 1; + if (dim < 15) + return 2; + if (dim < 24) + return 3; + return 4; +} + +constexpr size_t floatsPerPoint(SplatAttribute attr, int32_t shDegree) { + switch (attr) { + case SplatAttribute::Positions: return 3; + case SplatAttribute::Alphas: return 1; + case SplatAttribute::Colors: return 3; + case SplatAttribute::Scales: return 3; + case SplatAttribute::Rotations: return 4; + case SplatAttribute::Sh: return static_cast(dimForDegree(shDegree)) * 3; + } + return 0; +} + +inline float unquantizeSH(uint8_t x) { + return (static_cast(x) - 128.0f) / 128.0f; +} + +inline float invSigmoid(float x) { + return std::log(x / (1.0f - x)); +} + +// Decode a quaternion stored as its first three components in 8-bit fixed point; +// w is reconstructed from the constraint that the quaternion is unit-length and +// non-negative. Within-family: `c.flipQ` scales xyz (fused into the decode). +// Cross-family: coordinateConverter sets `c.flipQ` to identity (so the fused +// multiply is a no-op) and `c.rotFlipQFunc` rotates the full xyzw afterwards. +inline void unpackQuaternionFirstThree( + float rotation[4], const uint8_t r[3], + const CoordinateConverter &c = CoordinateConverter()) { + Vec3f xyz = times( + plus( + times( + Vec3f{static_cast(r[0]), static_cast(r[1]), static_cast(r[2])}, + 1.0f / 127.5f), + Vec3f{-1, -1, -1}), + c.flipQ); + std::copy(xyz.data(), xyz.data() + 3, &rotation[0]); + rotation[3] = std::sqrt(std::max(0.0f, 1.0f - squaredNorm(xyz))); + if (c.rotFlipQFunc) { + c.rotFlipQFunc(rotation); + } +} + +// Decode a quaternion stored in the "smallest three" encoding: 2 bits identify +// the largest component (reconstructed from the unit-length constraint) and the +// other three are 10-bit signed magnitudes scaled by 1/sqrt(2). +inline void unpackQuaternionSmallestThree( + float rotation[4], const uint8_t r[4], + const CoordinateConverter &c = CoordinateConverter()) { + uint32_t comp = + r[0] + (r[1] << 8) + (r[2] << 16) + (r[3] << 24); + + constexpr uint32_t c_mask = (1u << 9u) - 1u; + + const int i_largest = comp >> 30; + float sum_squares = 0; + for (int i = 3; i >= 0; --i) { + if (i != i_largest) { + uint32_t mag = comp & c_mask; + uint32_t negbit = (comp >> 9u) & 0x1u; + comp = comp >> 10u; + rotation[i] = sqrt1_2 * static_cast(mag) / static_cast(c_mask); + if (negbit == 1) { + rotation[i] = -rotation[i]; + } + sum_squares += rotation[i] * rotation[i]; + } + } + rotation[i_largest] = std::sqrt(1.0f - sum_squares); + + // Within-family flip via c.flipQ (identity in the cross-family case, where + // c.rotFlipQFunc carries the full flip+rotation on the xyzw quaternion). + for (int i = 0; i < 3; i++) { + rotation[i] *= c.flipQ[i]; + } + if (c.rotFlipQFunc) { + c.rotFlipQFunc(rotation); + } +} + +// Access a PackedGaussians attribute buffer by enum tag. Lets callers iterate +// attributes (see kAllSplatAttributes) instead of hard-coding the six member +// names in every loop. +inline std::vector &packedBuffer(PackedGaussians &p, SplatAttribute attr) { + switch (attr) { + case SplatAttribute::Positions: return p.positions; + case SplatAttribute::Alphas: return p.alphas; + case SplatAttribute::Colors: return p.colors; + case SplatAttribute::Scales: return p.scales; + case SplatAttribute::Rotations: return p.rotations; + case SplatAttribute::Sh: return p.sh; + } + return p.positions; // unreachable +} + +inline constexpr size_t floatCount(const PackedGaussians &p, SplatAttribute attr) { + return floatsPerPoint(attr, p.shDegree) * static_cast(p.numPoints); +} + +inline void releasePacked(PackedGaussians &p, SplatAttribute attr) { + std::vector().swap(packedBuffer(p, attr)); +} + +} // namespace spz