chore: import upstream snapshot with attribution

2026-07-13 13:18:05 +08:00
commit acd8e21031
297 changed files with 56514 additions and 0 deletions
@@ -0,0 +1,23 @@
+Language: Cpp
+BasedOnStyle: LLVM
+UseTab: Never
+IndentWidth: 4
+TabWidth: 4
+BreakBeforeBraces: Attach
+AllowShortIfStatementsOnASingleLine: false
+IndentCaseLabels: false
+ColumnLimit: 0
+AccessModifierOffset: -4
+PointerAlignment: Left
+SpaceAfterTemplateKeyword: false
+Cpp11BracedListStyle: false
+SpaceBeforeCpp11BracedList: true
+AlignAfterOpenBracket: DontAlign
+ContinuationIndentWidth: 4
+InsertNewlineAtEOF: true
+# BreakAfterOpenBracketBracedList: true
+# BreakAfterOpenBracketFunction: true
+# BreakAfterOpenBracketIf: true
+# BreakBeforeCloseBracketBracedList: true
+# BreakBeforeCloseBracketFunction: true
+# BreakBeforeCloseBracketIf: true
@@ -0,0 +1,74 @@
+cmake_minimum_required(VERSION 3.31)
+cmake_policy(SET CMP0091 NEW)
+cmake_policy(SET CMP0042 NEW)
+
+project(binding)
+
+find_package(Threads REQUIRED)
+find_package(Eigen3 CONFIG REQUIRED)
+find_package(nanoflann CONFIG REQUIRED)
+find_package(cmake-js CONFIG REQUIRED)
+find_package(WebP CONFIG REQUIRED)
+find_package(libavif CONFIG REQUIRED)
+
+set(BINDING_BINARY_DIR "${CMAKE_CURRENT_BINARY_DIR}" CACHE PATH "node binding output dir")
+
+add_library(${PROJECT_NAME} SHARED)
+list(APPEND BINDING_SOURCE_FILES
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/thread_pool.cpp"
+    # GAUSSIAN SPLAT
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/splat/splat.cpp"
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/splat/splat_lod.cpp"
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/splat/splat_block.cpp"
+    # NODE API FILES
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/node/api_splat.cpp"
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/node/api_spatial.cpp"
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/node/api_thread_pool.cpp"
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/node/api_webp.cpp"
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/node/api_avif.cpp"
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/node/binding.cpp"
+)
+list(APPEND BINDING_HEADER_FILES
+    "${CMAKE_CURRENT_SOURCE_DIR}/include/future_helpers.h"
+    "${CMAKE_CURRENT_SOURCE_DIR}/include/container_helpers.h"
+    "${CMAKE_CURRENT_SOURCE_DIR}/include/thread_pool.h"
+    # GAUSSIAN SPLAT
+    "${CMAKE_CURRENT_SOURCE_DIR}/include/splat/splat.h"
+    "${CMAKE_CURRENT_SOURCE_DIR}/include/splat/splat_lod.h"
+    "${CMAKE_CURRENT_SOURCE_DIR}/include/splat/splat_block.h"
+    # NODE API FILES
+    "${CMAKE_CURRENT_SOURCE_DIR}/include/node/api_buffer.h"
+    "${CMAKE_CURRENT_SOURCE_DIR}/include/node/api_splat.h"
+    "${CMAKE_CURRENT_SOURCE_DIR}/include/node/api_spatial.h"
+    "${CMAKE_CURRENT_SOURCE_DIR}/include/node/api_thread_pool.h"
+    "${CMAKE_CURRENT_SOURCE_DIR}/include/node/api_webp.h"
+    "${CMAKE_CURRENT_SOURCE_DIR}/include/node/api_avif.h"
+)
+target_sources(${PROJECT_NAME}
+PRIVATE
+    ${BINDING_SOURCE_FILES}
+
+PRIVATE FILE_SET binding_private_headers
+    TYPE HEADERS
+    BASE_DIRS "${CMAKE_CURRENT_SOURCE_DIR}/include"
+    FILES ${BINDING_HEADER_FILES}
+)
+target_link_libraries(${PROJECT_NAME}
+    PRIVATE cmake_js::node Threads::Threads Eigen3::Eigen nanoflann::nanoflann WebP::webp avif
+)
+set_target_properties(
+    ${PROJECT_NAME} PROPERTIES
+    PREFIX ""
+    OUTPUT_NAME "splat-transform"
+    SUFFIX ".node"
+    CXX_EXTENSIONS OFF
+    COMPILE_FEATURES cxx_std_20
+    COMPILE_OPTIONS $<$<CXX_COMPILER_ID:MSVC>:/Zc:__cplusplus>
+    LINK_OPTIONS $<$<AND:$<CXX_COMPILER_ID:GNU,Clang>,$<NOT:$<PLATFORM_ID:Apple>>>:-Wl,--exclude-libs,ALL>
+    COMPILE_DEFINITIONS "UNICODE;_UNICODE"
+    INTERPROCEDURAL_OPTIMIZATION ON
+    LIBRARY_OUTPUT_DIRECTORY ${BINDING_BINARY_DIR}
+    RUNTIME_OUTPUT_DIRECTORY ${BINDING_BINARY_DIR}
+    RUNTIME_OUTPUT_DIRECTORY_RELEASE ${BINDING_BINARY_DIR}
+    LIBRARY_OUTPUT_DIRECTORY_RELEASE ${BINDING_BINARY_DIR}
+)
@@ -0,0 +1,155 @@
+{
+  "version": 10,
+  "cmakeMinimumRequired": {
+    "major": 3,
+    "minor": 31,
+    "patch": 0
+  },
+  "configurePresets": [
+    {
+      "$comment": [
+        "default build config for any env",
+        "needs: vcpkg",
+        "vcpkg install root should be in environments as VCPKG_ROOT"
+      ],
+      "name": "default",
+      "toolchainFile": "$env{VCPKG_ROOT}/scripts/buildsystems/vcpkg.cmake",
+      "cacheVariables": {
+        "CMAKE_PREFIX_PATH": "${sourceDir}/cmake"
+      }
+    },
+    {
+      "name": "x64-windows-static",
+      "$comment": ["build config for x64 windows"],
+      "inherits": "default",
+      "cacheVariables": {
+        "VCPKG_TARGET_TRIPLET": "x64-windows-static"
+      },
+      "condition": {
+        "type": "equals",
+        "lhs": "${hostSystemName}",
+        "rhs": "Windows"
+      }
+    },
+    {
+      "name": "arm64-windows-static",
+      "$comment": ["build config for arm64 windows"],
+      "inherits": "default",
+      "cacheVariables": {
+        "VCPKG_TARGET_TRIPLET": "arm64-windows-static"
+      },
+      "condition": {
+        "type": "equals",
+        "lhs": "${hostSystemName}",
+        "rhs": "Windows"
+      }
+    },
+    {
+      "name": "x64-linux",
+      "$comment": ["build config for x64 linux"],
+      "inherits": "default",
+      "cacheVariables": {
+        "VCPKG_TARGET_TRIPLET": "x64-linux",
+        "CMAKE_MAKE_PROGRAM": "/usr/bin/ninja",
+        "CMAKE_C_COMPILER": "/usr/bin/cc",
+        "CMAKE_CXX_COMPILER": "/usr/bin/c++"
+      },
+      "condition": {
+        "type": "equals",
+        "lhs": "${hostSystemName}",
+        "rhs": "Linux"
+      }
+    },
+    {
+      "name": "arm64-linux",
+      "$comment": ["build config for arm64 linux"],
+      "inherits": "default",
+      "cacheVariables": {
+        "VCPKG_TARGET_TRIPLET": "arm64-linux",
+        "CMAKE_ASM_NASM_COMPILER": "/usr/bin/nasm",
+        "CMAKE_MAKE_PROGRAM": "/usr/bin/ninja",
+        "CMAKE_C_COMPILER": "/usr/bin/cc",
+        "CMAKE_CXX_COMPILER": "/usr/bin/c++"
+      },
+      "condition": {
+        "type": "equals",
+        "lhs": "${hostSystemName}",
+        "rhs": "Linux"
+      }
+    },
+    {
+      "name": "arm64-osx",
+      "$comment": ["build config for arm64 osx"],
+      "inherits": "default",
+      "cacheVariables": {
+        "VCPKG_TARGET_TRIPLET": "arm64-osx",
+        "CMAKE_MAKE_PROGRAM": "/opt/homebrew/bin/ninja",
+        "CMAKE_ASM_NASM_COMPILER": "/opt/homebrew/bin/nasm",
+        "CMAKE_C_COMPILER": "/usr/bin/cc",
+        "CMAKE_CXX_COMPILER": "/usr/bin/c++"
+      },
+      "condition": {
+        "type": "equals",
+        "lhs": "${hostSystemName}",
+        "rhs": "Darwin"
+      }
+    }
+  ],
+  "buildPresets": [
+    {
+      "name": "default",
+      "configurePreset": "default",
+      "jobs": 8
+    },
+    {
+      "name": "x64-windows-static",
+      "inherits": "default",
+      "configurePreset": "x64-windows-static",
+      "condition": {
+        "type": "equals",
+        "lhs": "${hostSystemName}",
+        "rhs": "Windows"
+      }
+    },
+    {
+      "name": "arm64-windows-static",
+      "inherits": "default",
+      "configurePreset": "arm64-windows-static",
+      "condition": {
+        "type": "equals",
+        "lhs": "${hostSystemName}",
+        "rhs": "Windows"
+      }
+    },
+    {
+      "name": "x64-linux",
+      "inherits": "default",
+      "configurePreset": "x64-linux",
+      "condition": {
+        "type": "equals",
+        "lhs": "${hostSystemName}",
+        "rhs": "Linux"
+      }
+    },
+    {
+      "name": "arm64-linux",
+      "inherits": "default",
+      "configurePreset": "arm64-linux",
+      "condition": {
+        "type": "equals",
+        "lhs": "${hostSystemName}",
+        "rhs": "Linux"
+      }
+    },
+    {
+      "name": "arm64-osx",
+      "inherits": "default",
+      "configurePreset": "arm64-osx",
+      "condition": {
+        "type": "equals",
+        "lhs": "${hostSystemName}",
+        "rhs": "Darwin"
+      }
+    }
+  ]
+}
@@ -0,0 +1,47 @@
+import path from 'node:path';
+import { createRequire } from 'node:module';
+import yargs from 'yargs/yargs';
+import { hideBin } from 'yargs/helpers';
+import { spawnProcess } from '@internal/utils/process.js';
+
+const __dirname = import.meta.dirname;
+const require = createRequire(import.meta.url);
+
+await yargs(hideBin(process.argv))
+    .scriptName('build')
+    .showHelpOnFail(false)
+    .strict()
+    .command({
+        command: ['build', '$0'],
+        builder(argv) {
+            argv.option('preset', {
+                describe: 'CMake preset',
+                default: 'default',
+            }).option('target', {
+                describe: 'Platform target',
+                demandOption: true,
+            });
+        },
+        async handler(argv) {
+            await spawnProcess(
+                'node',
+                [
+                    require.resolve('cmake-js/bin/cmake-js'),
+                    'build',
+                    `--CDBINDING_BINARY_DIR=${path.resolve(`../splat-transform-${argv.target}`)}`,
+                    '--',
+                    '--preset',
+                    argv.preset,
+                ],
+                {
+                    cwd: process.cwd(),
+                    env: {
+                        ...process.env,
+                        FORCE_COLOR: 1,
+                        PATH: path.join(__dirname, '.`/node_modules/.bin') + path.delimiter + process.env.PATH,
+                    },
+                },
+            ).promise;
+        },
+    })
+    .parseAsync();
@@ -0,0 +1,25 @@
+# cmake-js config for find_package
+add_library(cmake_js::node INTERFACE IMPORTED)
+
+if(MSVC)
+    if(CMAKE_JS_NODELIB_DEF AND CMAKE_JS_NODELIB_TARGET)
+        # Generate node.lib
+        add_custom_command(OUTPUT ${CMAKE_JS_NODELIB_TARGET}
+            COMMAND ${CMAKE_AR} /def:${CMAKE_JS_NODELIB_DEF} /out:${CMAKE_JS_NODELIB_TARGET} ${CMAKE_STATIC_LINKER_FLAGS}
+        )
+        add_custom_target(node_lib DEPENDS ${CMAKE_JS_NODELIB_TARGET})
+        target_link_libraries(cmake_js::node INTERFACE ${CMAKE_JS_NODELIB_TARGET})
+        add_dependencies(cmake_js::node node_lib)
+    endif()
+endif()
+
+if(CMAKE_SYSTEM_NAME STREQUAL "Windows")
+    target_link_libraries(cmake_js::node INTERFACE "delayimp.lib")
+endif()
+
+set_target_properties(cmake_js::node PROPERTIES
+    INTERFACE_COMPILE_FEATURES cxx_std_14
+    INTERFACE_COMPILE_DEFINITIONS "NAPI_VERSION=8"
+    INTERFACE_INCLUDE_DIRECTORIES "${CMAKE_JS_INC}"
+    INTERFACE_SOURCES "${CMAKE_JS_SRC}"
+)
@@ -0,0 +1,39 @@
+#pragma once
+#include <initializer_list>
+#include <ranges>
+#include <vector>
+
+namespace helpers::container {
+
+template<typename R, typename V>
+concept range_of = std::ranges::range<R> && std::same_as<std::ranges::range_value_t<R>, V>;
+
+template<std::ranges::range T, std::ranges::range U>
+    requires std::same_as<std::ranges::range_value_t<T>, std::ranges::range_value_t<U>> &&
+             requires(T&& dst, U&& src) {
+#if __cpp_lib_containers_ranges
+                 dst.append_range(std::forward<U>(src));
+#else
+                 dst.insert(dst.end(), src.begin(), src.end());
+#endif
+             }
+inline void append_range(T&& dst, U&& src) {
+#if __cpp_lib_containers_ranges
+    dst.append_range(std::forward<U>(src));
+#else
+    dst.insert(dst.end(), src.begin(), src.end());
+#endif
+}
+
+template<std::ranges::range T>
+    requires requires(T&& dst, std::initializer_list<std::ranges::range_value_t<T>>&& src) {
+#if __cpp_lib_containers_ranges
+        dst.append_range(std::move(src));
+#else
+        dst.insert(dst.end(), src.begin(), src.end());
+#endif
+    }
+inline void append_range(T&& dst, std::initializer_list<std::ranges::range_value_t<T>>&& src) {
+    append_range<T, std::initializer_list<std::ranges::range_value_t<T>>&&>(std::forward<T>(dst), std::move(src));
+}
+} // namespace helpers::container
@@ -0,0 +1,135 @@
+#pragma once
+#include <concepts>
+#include <cstddef>
+#include <exception>
+#include <functional>
+#include <future>
+#include <ranges>
+#include <type_traits>
+
+namespace helpers::future {
+namespace detail {
+template<typename T>
+struct future_value;
+
+template<typename T>
+struct future_value<std::future<T>> {
+    using type = T;
+};
+
+template<typename T>
+concept future_range = std::ranges::range<T> && requires {
+    typename future_value<std::ranges::range_value_t<std::remove_cvref_t<T>>>::type;
+};
+
+template<future_range T>
+using future_range_value_t = future_value<std::ranges::range_value_t<std::remove_cvref_t<T>>>::type;
+} // namespace detail
+
+struct FailurePolicyDrainAll;
+// Stops draining after the first observed failure. This does not cancel remaining
+// work; use only when unfinished tasks do not depend on this stack frame.
+struct FailurePolicyStopOnFirstFailure;
+
+template<typename T>
+concept future_failure_policy = std::same_as<T, FailurePolicyDrainAll> || std::same_as<T, FailurePolicyStopOnFirstFailure>;
+
+template<future_failure_policy P = FailurePolicyDrainAll, detail::future_range T>
+    requires std::same_as<detail::future_range_value_t<T>, void>
+void drain_futures(T&& futures) {
+    std::exception_ptr first_exception;
+    for (auto& future : futures) {
+        try {
+            future.get();
+        } catch (...) {
+            if (!first_exception) {
+                first_exception = std::current_exception();
+            }
+
+            if constexpr (std::is_same_v<P, FailurePolicyStopOnFirstFailure>) {
+                break;
+            }
+        }
+    }
+
+    if (first_exception) {
+        std::rethrow_exception(first_exception);
+    }
+}
+
+template<future_failure_policy P = FailurePolicyDrainAll, detail::future_range T, typename F>
+    requires std::same_as<detail::future_range_value_t<T>, void> &&
+             std::invocable<F&, size_t>
+void drain_futures(T&& futures, F&& f) {
+    std::exception_ptr first_exception;
+    size_t index = 0;
+    for (auto& future : futures) {
+        try {
+            future.get();
+            std::invoke(f, index);
+        } catch (...) {
+            if (!first_exception) {
+                first_exception = std::current_exception();
+            }
+
+            if constexpr (std::is_same_v<P, FailurePolicyStopOnFirstFailure>) {
+                break;
+            }
+        }
+        index++;
+    }
+
+    if (first_exception) {
+        std::rethrow_exception(first_exception);
+    }
+}
+
+template<future_failure_policy P = FailurePolicyDrainAll, detail::future_range T, typename F>
+    requires(!std::same_as<detail::future_range_value_t<T>, void>) &&
+            std::invocable<F&, detail::future_range_value_t<T>&&>
+void drain_futures(T&& futures, F&& f) {
+    std::exception_ptr first_exception;
+    for (auto& future : futures) {
+        try {
+            std::invoke(f, future.get());
+        } catch (...) {
+            if (!first_exception) {
+                first_exception = std::current_exception();
+            }
+            if constexpr (std::is_same_v<P, FailurePolicyStopOnFirstFailure>) {
+                break;
+            }
+        }
+    }
+
+    if (first_exception) {
+        std::rethrow_exception(first_exception);
+    }
+}
+
+template<future_failure_policy P = FailurePolicyDrainAll, detail::future_range T, typename F>
+    requires(!std::same_as<detail::future_range_value_t<T>, void>) &&
+            std::invocable<F&, detail::future_range_value_t<T>&&, size_t>
+void drain_futures(T&& futures, F&& f) {
+    std::exception_ptr first_exception;
+    size_t index = 0;
+    for (auto& future : futures) {
+        try {
+            std::invoke(f, future.get(), index);
+        } catch (...) {
+            if (!first_exception) {
+                first_exception = std::current_exception();
+            }
+
+            if constexpr (std::is_same_v<P, FailurePolicyStopOnFirstFailure>) {
+                break;
+            }
+        }
+        index++;
+    }
+
+    if (first_exception) {
+        std::rethrow_exception(first_exception);
+    }
+}
+} // namespace helpers::future
@@ -0,0 +1,9 @@
+#pragma once
+#include <napi.h>
+
+namespace node_api::imaging {
+Napi::Value avif_encode_rgba(const Napi::CallbackInfo& info);
+Napi::Value avif_encode_rgba_batched(const Napi::CallbackInfo& info);
+Napi::Value avif_decode_rgba(const Napi::CallbackInfo& info);
+Napi::Value avif_decode_rgba_batched(const Napi::CallbackInfo& info);
+} // namespace node_api::imaging
@@ -0,0 +1,92 @@
+#pragma once
+#include <memory>
+#include <napi.h>
+#include <span>
+#include <vector>
+
+namespace node_api::buffer {
+template<typename T, typename D = std::default_delete<T>>
+class UniqueBufferFinalizer {
+public:
+    UniqueBufferFinalizer() noexcept = default;
+    UniqueBufferFinalizer(const UniqueBufferFinalizer<T, D>& other) = delete;
+    UniqueBufferFinalizer(UniqueBufferFinalizer<T, D>&& other) noexcept = default;
+    UniqueBufferFinalizer(std::unique_ptr<T, D>&& ptr) noexcept : ptr(std::move(ptr)) {}
+    UniqueBufferFinalizer(T* ptr) noexcept : ptr(ptr) {}
+    void operator()(Napi::Env env, T* data) noexcept {}
+
+    inline static Napi::Buffer<T> make_buffer(Napi::Env& env, std::unique_ptr<T, D>&& ptr, size_t size) noexcept {
+        auto finalizer = UniqueBufferFinalizer(std::move(ptr));
+        auto data_ptr = finalizer.ptr.get();
+        return Napi::Buffer<T>::NewOrCopy(env, data_ptr, size, std::move(finalizer));
+    }
+
+protected:
+    std::unique_ptr<T, D> ptr;
+};
+
+template<typename T, typename D = std::default_delete<std::vector<T>>>
+class UniqueVecBufferFinalizer : public UniqueBufferFinalizer<std::vector<T>, D> {
+public:
+    void operator()(Napi::Env env, T* data) noexcept {}
+
+    template<typename U>
+    inline static Napi::Buffer<T> make_buffer(Napi::Env& env, U&& data) noexcept {
+        auto finalizer = UniqueVecBufferFinalizer(std::forward<U>(data));
+        auto data_ptr = finalizer.ptr->data();
+        auto data_size = finalizer.ptr->size();
+        return Napi::Buffer<T>::NewOrCopy(env, data_ptr, data_size, std::move(finalizer));
+    }
+};
+
+template<typename T>
+class SharedBufferFinalizer {
+public:
+    SharedBufferFinalizer() noexcept = default;
+    SharedBufferFinalizer(const SharedBufferFinalizer<T>& other) noexcept = default;
+    SharedBufferFinalizer(SharedBufferFinalizer<T>&& other) noexcept = default;
+    SharedBufferFinalizer(const std::shared_ptr<T>& ptr) noexcept : ptr(ptr) {}
+    SharedBufferFinalizer(std::shared_ptr<T>&& ptr) noexcept : ptr(std::move(ptr)) {}
+    void operator()(Napi::Env env, T* data) noexcept {}
+
+    template<typename U>
+        requires std::is_same_v<std::remove_reference_t<U>, std::shared_ptr<T>>
+    inline static Napi::Buffer<T> make_buffer(Napi::Env& env, U&& ptr, size_t size) noexcept {
+        auto finalizer = SharedBufferFinalizer(std::forward<U>(ptr));
+        auto data_ptr = finalizer.ptr.get();
+        return Napi::Buffer<T>::NewOrCopy(env, data_ptr, size, std::move(finalizer));
+    }
+
+    template<typename U>
+        requires std::is_same_v<std::remove_reference_t<U>, std::shared_ptr<T>>
+    inline static Napi::Buffer<T> make_buffer(Napi::Env& env, U&& ptr, size_t offset, size_t size) noexcept {
+        auto finalizer = SharedBufferFinalizer(std::forward<U>(ptr));
+        auto data_ptr = finalizer.ptr.get();
+        return Napi::Buffer<T>::NewOrCopy(env, data_ptr + offset, size, std::move(finalizer));
+    }
+
+protected:
+    std::shared_ptr<T> ptr;
+};
+
+template<typename T>
+class SharedVecBufferFinalizer : public SharedBufferFinalizer<std::vector<T>> {
+public:
+    void operator()(Napi::Env env, T* data) noexcept {}
+
+    template<typename U>
+    inline static Napi::Buffer<T> make_buffer(Napi::Env& env, U&& data) noexcept {
+        auto finalizer = SharedVecBufferFinalizer(std::forward<U>(data));
+        auto data_ptr = finalizer.ptr->data();
+        auto data_size = finalizer.ptr->size();
+        return Napi::Buffer<T>::NewOrCopy(env, data_ptr, data_size, std::move(finalizer));
+    }
+
+    template<typename U>
+    inline static Napi::Buffer<T> make_buffer(Napi::Env& env, U&& data, size_t offset, size_t size) noexcept {
+        auto finalizer = SharedVecBufferFinalizer(std::forward<U>(data));
+        auto data_ptr = finalizer.ptr->data();
+        return Napi::Buffer<T>::NewOrCopy(env, data_ptr + offset, size, std::move(finalizer));
+    }
+};
+} // namespace node_api::buffer
@@ -0,0 +1,6 @@
+#pragma once
+#include <napi.h>
+
+namespace node_api::spatial {
+Napi::Value cluster_average(const Napi::CallbackInfo& info);
+}
@@ -0,0 +1,7 @@
+#pragma once
+#include <napi.h>
+
+namespace node_api::splat {
+Napi::Value generate_lod(const Napi::CallbackInfo& info);
+Napi::Value split(const Napi::CallbackInfo& info);
+} // namespace node_api::splat
@@ -0,0 +1,26 @@
+#pragma once
+#include <algorithm>
+#include <bit>
+#include <cstddef>
+#include <napi.h>
+#include <thread>
+#include <thread_pool.h>
+
+namespace node_api::threading {
+class ThreadPool : public Napi::ObjectWrap<ThreadPool> {
+public:
+    static Napi::FunctionReference Init(Napi::Env env, Napi::Object exports);
+
+    ThreadPool(const Napi::CallbackInfo& info);
+    ::threading::ThreadPool& impl() noexcept;
+
+    ThreadPool() = delete;
+
+private:
+    static size_t calc_thread_count(const Napi::CallbackInfo& info);
+    Napi::Value thread_count(const Napi::CallbackInfo& info);
+    Napi::Value task_count(const Napi::CallbackInfo& info);
+
+    ::threading::ThreadPool impl_;
+};
+} // namespace node_api::threading
@@ -0,0 +1,8 @@
+#pragma once
+#include <napi.h>
+
+namespace node_api::imaging {
+Napi::Value webp_encode_rgba_lossless(const Napi::CallbackInfo& info);
+Napi::Value webp_encode_rgba(const Napi::CallbackInfo& info);
+Napi::Value webp_decode_rgba(const Napi::CallbackInfo& info);
+} // namespace node_api::imaging
@@ -0,0 +1,60 @@
+#pragma once
+#include <array>
+#include <eigen3/Eigen/Dense>
+#include <memory>
+#include <vector>
+
+namespace splat {
+class SH {
+public:
+    SH() noexcept;
+    SH(size_t size);
+    SH(const SH& other);
+    SH(SH&& other) noexcept;
+
+    SH& operator=(const SH& other);
+    SH& operator=(SH&& other) noexcept;
+
+    void swap(SH& other) noexcept;
+
+    SH& set_zero() noexcept;
+    SH& add_multiplied(const SH& other, float scalar) noexcept;
+
+    const float& operator[](size_t index) const;
+    float& operator[](size_t index);
+
+    float* data() noexcept;
+    size_t size() const noexcept;
+
+    ~SH() noexcept;
+
+private:
+    void reset() noexcept;
+
+    size_t size_;
+    std::unique_ptr<float[]> ptr;
+};
+
+struct Gaussian {
+    Eigen::Vector3f mean;
+    Eigen::Vector3f scale;
+    Eigen::Vector4f rotation;
+    Eigen::Matrix3f covariance;
+    SH sh;
+    float opacity;
+    Eigen::AlignedBox3f bounding_box;
+
+    void compute_covariance();
+    void decompose_covariance();
+    void compute_bounding_box(float k = 3.0f);
+    float area() const;
+};
+
+struct Splat {
+    std::vector<Gaussian> gaussians;
+    Eigen::AlignedBox3f bounding_box;
+
+    void compute_bounding_box();
+    void compute_compact_bounding_box();
+};
+} // namespace splat
@@ -0,0 +1,17 @@
+#pragma once
+#include <splat/splat.h>
+#include <type_traits>
+#include <vector>
+
+namespace splat::block {
+namespace detail {
+std::vector<Splat> split(const Splat& splat, size_t max_block_size);
+} // namespace detail
+
+template<typename T>
+    requires std::is_same_v<std::remove_reference_t<T>, Splat>
+std::vector<Splat> split(T&& input, double precision) {
+    auto max_block_size = static_cast<size_t>(static_cast<double>(input.gaussians.size()) * precision);
+    return detail::split(input, max_block_size);
+}
+} // namespace splat::block
@@ -0,0 +1,43 @@
+#pragma once
+#include <container_helpers.h>
+#include <cstdint>
+#include <splat/splat.h>
+#include <vector>
+
+namespace splat::lod {
+namespace detail {
+Splat reduce_gaussians(size_t id, const Splat& splat, size_t target_count, float scale_boost, size_t max_step);
+}
+
+struct SplatLod {
+    std::vector<Splat> splats;
+    std::vector<uint32_t> levels;
+};
+
+struct SplatLevelParameters {
+    float precision;
+    float scale_boost;
+};
+
+template<typename T>
+    requires requires(T) {
+        requires helpers::container::range_of<T, SplatLevelParameters>;
+    }
+inline void generate_lod(size_t id, SplatLod& lod, T&& parameters, size_t min_size, size_t max_step) {
+    auto base_count = lod.splats.back().gaussians.size();
+    for (auto& parameter : parameters) {
+        auto target_count = std::max<size_t>(static_cast<size_t>(std::floor(static_cast<double>(base_count) * parameter.precision)), 1);
+        auto& current = lod.splats.back();
+        if (current.gaussians.size() > min_size && target_count < current.gaussians.size()) {
+            auto generated = detail::reduce_gaussians(
+                id,
+                current,
+                target_count,
+                parameter.scale_boost,
+                max_step);
+            lod.splats.push_back(std::move(generated));
+        }
+        lod.levels.push_back(lod.splats.size() - 1);
+    }
+}
+} // namespace splat::lod
@@ -0,0 +1,122 @@
+#pragma once
+#include <atomic>
+#include <concepts>
+#include <condition_variable>
+#include <functional>
+#include <future>
+#include <memory>
+#include <mutex>
+#include <queue>
+#include <shared_mutex>
+#include <stdexcept>
+#include <thread>
+#include <type_traits>
+#include <utility>
+#include <vector>
+
+namespace threading {
+class ThreadPool {
+public:
+    ThreadPool(size_t thread_count = std::thread::hardware_concurrency());
+    ThreadPool(ThreadPool&& other) noexcept;
+
+    template<typename F, typename... Args,
+        typename ReturnType = std::invoke_result_t<std::decay_t<F>, std::decay_t<Args>...>>
+        requires std::invocable<std::decay_t<F>, std::decay_t<Args>...>
+    std::future<ReturnType> submit_task(F&& f, Args&&... args);
+    size_t thread_count() const noexcept;
+    size_t task_count() const noexcept;
+    ThreadPool& operator=(ThreadPool&& other) noexcept;
+    void stop() noexcept;
+    ~ThreadPool() noexcept;
+
+    ThreadPool(const ThreadPool& other) = delete;
+    ThreadPool& operator=(const ThreadPool& other) = delete;
+
+private:
+#ifdef __cpp_lib_move_only_function
+    using Task = std::move_only_function<void()>;
+#else
+    using Task = std::function<void()>;
+#endif
+
+    struct Worker {
+        std::atomic<bool> stopped { false };
+        std::thread thread;
+
+        void stop() noexcept;
+        void clean_up() noexcept;
+        ~Worker() noexcept;
+    };
+
+    struct State {
+        std::vector<std::unique_ptr<Worker>> workers;
+        std::queue<Task> tasks;
+        mutable std::shared_mutex queue_mutex;
+        mutable std::shared_mutex worker_mutex;
+        std::condition_variable_any cv;
+        std::atomic<bool> stopped { false };
+
+        void stop() noexcept;
+        ~State() noexcept;
+    };
+
+    static void worker_loop(State* state, Worker* worker);
+
+    std::unique_ptr<State> state;
+};
+
+template<typename F, typename... Args, typename ReturnType>
+    requires std::invocable<std::decay_t<F>, std::decay_t<Args>...>
+std::future<ReturnType> ThreadPool::submit_task(F&& f, Args&&... args) {
+    auto* state = this->state.get();
+    if (state == nullptr || state->stopped.load(std::memory_order_relaxed)) {
+        throw std::runtime_error("Cannot submit tasks to stopped thread pool");
+    }
+#ifdef __cpp_lib_move_only_function
+    auto promise = std::promise<ReturnType>();
+    auto future = promise.get_future();
+    {
+        auto lk = std::lock_guard(state->queue_mutex);
+        if (state->stopped.load(std::memory_order_relaxed)) {
+            throw std::runtime_error("Cannot submit tasks to stopped thread pool");
+        }
+        state->tasks.emplace([f = std::forward<F>(f), ... args = std::forward<Args>(args), promise = std::move(promise)]() mutable noexcept -> void {
+            try {
+                if constexpr (std::is_void_v<ReturnType>) {
+                    std::invoke(std::move(f), std::move(args)...);
+                    promise.set_value();
+                } else {
+                    promise.set_value(std::invoke(std::move(f), std::move(args)...));
+                }
+            } catch (...) {
+                promise.set_exception(std::current_exception());
+            }
+        });
+    }
+#else
+    auto promise = std::make_shared<std::promise<ReturnType>>();
+    auto future = promise->get_future();
+    {
+        auto lk = std::lock_guard(state->queue_mutex);
+        if (state->stopped.load(std::memory_order_relaxed)) {
+            throw std::runtime_error("Cannot submit tasks to stopped thread pool");
+        }
+        state->tasks.emplace([f = std::forward<F>(f), ... args = std::forward<Args>(args), promise = promise]() noexcept -> void {
+            try {
+                if constexpr (std::is_void_v<ReturnType>) {
+                    std::invoke(std::move(f), std::move(args)...);
+                    promise->set_value();
+                } else {
+                    promise->set_value(std::invoke(std::move(f), std::move(args)...));
+                }
+            } catch (...) {
+                promise->set_exception(std::current_exception());
+            }
+        });
+    }
+#endif
+    state->cv.notify_one();
+    return future;
+}
+} // namespace threading
@@ -0,0 +1,23 @@
+{
+  "name": "@internal/splat-transform-native-source",
+  "version": "0.0.1",
+  "private": true,
+  "type": "module",
+  "dependencies": {
+    "cmake-js": "^8.0.0",
+    "node-addon-api": "^8.9.0",
+    "yargs": "^18.0.0"
+  },
+  "scripts": {
+    "build:native:win32-arm64": "node ./build.js build --preset arm64-windows-static --target win32-arm64-msvc",
+    "build:native:win32-x64": "node ./build.js build --preset x64-windows-static --target win32-x64-msvc",
+    "build:native:linux-arm64-gnu": "node ./build.js build --preset arm64-linux --target linux-arm64-gnu",
+    "build:native:linux-x64-gnu": "node ./build.js build --preset x64-linux --target linux-x64-gnu",
+    "build:native:darwin-arm64": "node ./build.js build --preset arm64-osx --target darwin-arm64"
+  },
+  "binary": {
+    "napi_versions": [
+      8
+    ]
+  }
+}
@@ -0,0 +1,301 @@
+#include <avif/avif_cxx.h>
+#include <bit>
+#include <future_helpers.h>
+#include <memory>
+#include <node/api_avif.h>
+#include <node/api_buffer.h>
+#include <node/api_thread_pool.h>
+#include <span>
+#include <thread>
+#include <thread_pool.h>
+#include <vector>
+
+namespace {
+class RawData {
+public:
+    RawData() noexcept : data_(AVIF_DATA_EMPTY) {}
+    RawData(RawData&& other) noexcept : data_(other.data_) {
+        other.data_ = AVIF_DATA_EMPTY;
+    }
+    RawData(avifRWData&& data) noexcept : data_(data) {
+        data = AVIF_DATA_EMPTY;
+    }
+    RawData(const RawData& other) = delete;
+    void operator()(Napi::Env env, uint8_t* data) noexcept {}
+
+    RawData& operator=(RawData&& other) noexcept {
+        if (this != &other) {
+            if (this->data_.data != nullptr) {
+                avifRWDataFree(&this->data_);
+            }
+            this->data_ = other.data_;
+            other.data_ = AVIF_DATA_EMPTY;
+        }
+        return *this;
+    }
+
+    RawData& operator=(const RawData& other) = delete;
+
+    avifRWData& data() {
+        return this->data_;
+    }
+
+    const avifRWData& data() const {
+        return this->data_;
+    }
+
+    Napi::Buffer<uint8_t> make_buffer(Napi::Env& env) {
+        auto ptr = this->data_.data;
+        auto size = this->data_.size;
+        return Napi::Buffer<uint8_t>::NewOrCopy(env, ptr, size, std::move(*this));
+    }
+
+    ~RawData() {
+        if (this->data_.data != nullptr) {
+            avifRWDataFree(&this->data_);
+        }
+    }
+
+private:
+    avifRWData data_;
+};
+
+class RGBImageData {
+public:
+    RGBImageData() noexcept : data_({ 0 }) {}
+    RGBImageData(RGBImageData&& other) noexcept : data_(other.data_) {
+        other.data_ = { 0 };
+    }
+    RGBImageData(avifRGBImage&& data) noexcept : data_(data) {
+        data = { 0 };
+    }
+    RGBImageData(const RGBImageData& other) = delete;
+    void operator()(Napi::Env env, uint8_t* data) noexcept {}
+
+    RGBImageData& operator=(RGBImageData&& other) noexcept {
+        if (this != &other) {
+            if (this->data_.pixels != nullptr) {
+                avifRGBImageFreePixels(&this->data_);
+            }
+            this->data_ = other.data_;
+            other.data_ = { 0 };
+        }
+        return *this;
+    }
+
+    RGBImageData& operator=(const RGBImageData& other) = delete;
+
+    avifRGBImage& data() {
+        return this->data_;
+    }
+
+    const avifRGBImage& data() const {
+        return this->data_;
+    }
+
+    Napi::Buffer<uint8_t> make_buffer(Napi::Env& env) {
+        auto ptr = this->data_.pixels;
+        auto size = this->data_.rowBytes * this->data_.height;
+        return Napi::Buffer<uint8_t>::NewOrCopy(env, ptr, size, std::move(*this));
+    }
+
+    ~RGBImageData() {
+        if (this->data_.pixels != nullptr) {
+            avifRGBImageFreePixels(&this->data_);
+        }
+    }
+
+private:
+    avifRGBImage data_;
+};
+
+RawData avif_encode_rga_impl(std::span<uint8_t> pixels, int32_t width, int32_t height, int32_t quality, int32_t max_threads) {
+    avif::EncoderPtr encoder = nullptr;
+    avifRGBImage rgb;
+    std::memset(&rgb, 0, sizeof(rgb));
+    avif::ImagePtr image(avifImageCreate(width, height, 8, AVIF_PIXEL_FORMAT_YUV444));
+
+    avifRGBImageSetDefaults(&rgb, image.get());
+    rgb.pixels = pixels.data();
+    rgb.rowBytes = 4 * rgb.width;
+    avifResult convertResult = avifImageRGBToYUV(image.get(), &rgb);
+
+    if (convertResult != AVIF_RESULT_OK) {
+        throw std::runtime_error("Failed to convert to YUV(A)");
+    }
+
+    encoder.reset(avifEncoderCreate());
+    if (encoder == nullptr) {
+        throw std::runtime_error("Out of memory");
+    }
+
+    encoder->quality = quality;
+    encoder->qualityAlpha = AVIF_QUALITY_LOSSLESS;
+    encoder->maxThreads = max_threads;
+
+    avifResult addImageResult = avifEncoderAddImage(encoder.get(), image.get(), 1, AVIF_ADD_IMAGE_FLAG_SINGLE);
+
+    if (addImageResult != AVIF_RESULT_OK) {
+        throw std::runtime_error("Failed to add image to encoder");
+    }
+
+    RawData avifOutput;
+    avifResult finishResult = avifEncoderFinish(encoder.get(), &avifOutput.data());
+
+    if (finishResult != AVIF_RESULT_OK) {
+        throw std::runtime_error("Failed to finish encode");
+    }
+
+    return avifOutput;
+}
+
+RGBImageData avif_decode_rgba_impl(std::span<uint8_t> data, int32_t max_threads) {
+    avif::DecoderPtr decoder(avifDecoderCreate());
+
+    if (decoder == nullptr) {
+        throw std::runtime_error("Out of memory");
+    }
+
+    decoder->maxThreads = max_threads;
+
+    avifResult result = avifDecoderSetIOMemory(decoder.get(), data.data(), data.size());
+
+    if (result != AVIF_RESULT_OK) {
+        throw std::runtime_error("Cannot set IO on avifDecoder");
+    }
+
+    result = avifDecoderParse(decoder.get());
+
+    if (result != AVIF_RESULT_OK) {
+        throw std::runtime_error("Failed to decode image");
+    }
+
+    result = avifDecoderNextImage(decoder.get());
+
+    if (result != AVIF_RESULT_OK) {
+        throw std::runtime_error("Failed to get image");
+    }
+
+    RGBImageData rgb;
+    avifRGBImageSetDefaults(&rgb.data(), decoder->image);
+    result = avifRGBImageAllocatePixels(&rgb.data());
+
+    if (result != AVIF_RESULT_OK) {
+        throw std::runtime_error("Allocation of RGB samples failed");
+    }
+
+    result = avifImageYUVToRGB(decoder->image, &rgb.data());
+
+    if (result != AVIF_RESULT_OK) {
+        throw std::runtime_error("Conversion from YUV failed");
+    }
+
+    return rgb;
+}
+} // namespace
+
+namespace node_api::imaging {
+Napi::Value avif_encode_rgba(const Napi::CallbackInfo& info) {
+    auto env = info.Env();
+    if (info.Length() < 4 || !info[0].IsBuffer() || !info[1].IsNumber() || !info[2].IsNumber() || !info[3].IsNumber()) {
+        Napi::TypeError::New(env, "Wrong Arguments").ThrowAsJavaScriptException();
+        return env.Null();
+    }
+    auto buffer = info[0].As<Napi::Buffer<uint8_t>>();
+    auto width = info[1].As<Napi::Number>().Int32Value();
+    auto height = info[2].As<Napi::Number>().Int32Value();
+    auto quality = info[3].As<Napi::Number>().Int32Value();
+
+    return avif_encode_rga_impl(
+        std::span(buffer.Data(), buffer.Length()),
+        width, height,
+        quality, static_cast<int32_t>(std::thread::hardware_concurrency()))
+        .make_buffer(env);
+}
+Napi::Value avif_encode_rgba_batched(const Napi::CallbackInfo& info) {
+    auto env = info.Env();
+    if (info.Length() < 2 || !info[0].IsArray() || !info[1].IsObject()) {
+        Napi::TypeError::New(env, "Wrong Arguments").ThrowAsJavaScriptException();
+        return env.Null();
+    }
+
+    constexpr auto max_thread = 16;
+    auto& pool = node_api::threading::ThreadPool::Unwrap(info[1].As<Napi::Object>())->impl();
+    auto thread_count = pool.thread_count();
+    auto inputs = info[0].As<Napi::Array>();
+    auto input_count = inputs.Length();
+    auto outputs = Napi::Array::New(env, input_count);
+
+    {
+        auto futures = std::vector<std::future<RawData>>();
+        futures.reserve(inputs.Length());
+        for (auto i = 0; i < input_count; i++) {
+            auto input = inputs[i].AsValue().As<Napi::Object>();
+            auto buffer = input.Get("data").As<Napi::Buffer<uint8_t>>();
+            auto width = input.Get("width").As<Napi::Number>().Int32Value();
+            auto height = input.Get("height").As<Napi::Number>().Int32Value();
+            auto quality = input.Get("quality").As<Napi::Number>().Int32Value();
+            futures.push_back(pool.submit_task(avif_encode_rga_impl, std::span(buffer.Data(), buffer.Length()), width, height, quality, max_thread));
+        }
+
+        helpers::future::drain_futures(futures, [&](RawData&& data, size_t i) {
+            outputs[i] = data.make_buffer(env);
+        });
+    }
+
+    return outputs;
+}
+Napi::Value avif_decode_rgba(const Napi::CallbackInfo& info) {
+    auto env = info.Env();
+    if (info.Length() < 1 || !info[0].IsBuffer()) {
+        Napi::TypeError::New(env, "Wrong Arguments").ThrowAsJavaScriptException();
+        return env.Null();
+    }
+    auto buffer = info[0].As<Napi::Buffer<uint8_t>>();
+
+    auto&& rgb = avif_decode_rgba_impl(std::span(buffer.Data(), buffer.Length()),
+        static_cast<int32_t>(std::thread::hardware_concurrency()));
+
+    auto object = Napi::Object::New(env);
+    {
+        object.Set("width", Napi::Number::New(env, rgb.data().width));
+        object.Set("height", Napi::Number::New(env, rgb.data().height));
+        object.Set("data", rgb.make_buffer(env));
+    }
+    return object;
+}
+
+Napi::Value avif_decode_rgba_batched(const Napi::CallbackInfo& info) {
+    auto env = info.Env();
+    if (info.Length() < 2 || !info[0].IsArray() || !info[1].IsObject()) {
+        Napi::TypeError::New(env, "Wrong Arguments").ThrowAsJavaScriptException();
+        return env.Null();
+    }
+
+    constexpr auto max_thread = 16;
+    auto& pool = node_api::threading::ThreadPool::Unwrap(info[1].As<Napi::Object>())->impl();
+    auto thread_count = pool.thread_count();
+    auto inputs = info[0].As<Napi::Array>();
+    auto input_count = inputs.Length();
+    auto outputs = Napi::Array::New(env, input_count);
+
+    {
+        auto futures = std::vector<std::future<RGBImageData>>();
+        futures.reserve(inputs.Length());
+        for (auto i = 0; i < input_count; i++) {
+            auto input = inputs[i].AsValue().As<Napi::Buffer<uint8_t>>();
+            futures.push_back(pool.submit_task(avif_decode_rgba_impl, std::span(input.Data(), input.Length()), max_thread));
+        }
+
+        helpers::future::drain_futures(futures, [&](RGBImageData&& rgb, size_t i) {
+            auto object = Napi::Object::New(env);
+            object.Set("width", Napi::Number::New(env, rgb.data().width));
+            object.Set("height", Napi::Number::New(env, rgb.data().height));
+            object.Set("data", rgb.make_buffer(env));
+            outputs[i] = object;
+        });
+    }
+
+    return outputs;
+}
+} // namespace node_api::imaging
@@ -0,0 +1,150 @@
+#include <algorithm>
+#include <cstdint>
+#include <future>
+#include <future_helpers.h>
+#include <memory>
+#include <node/api_spatial.h>
+#include <node/api_thread_pool.h>
+#include <random>
+#include <span>
+#include <thread>
+#include <thread_pool.h>
+#include <vector>
+
+namespace {
+struct PartialAverage {
+    std::unique_ptr<float[]> sums;
+    std::unique_ptr<size_t[]> counts;
+};
+} // namespace
+
+namespace node_api::spatial {
+Napi::Value cluster_average(const Napi::CallbackInfo& info) {
+    auto env = info.Env();
+    if (info.Length() < 5 || !info[0].IsArray() || !info[1].IsBuffer() || !info[2].IsNumber() || !info[3].IsArray() || !info[4].IsObject()) {
+        Napi::TypeError::New(env, "Wrong Arguments").ThrowAsJavaScriptException();
+        return env.Null();
+    }
+
+    auto input_table = info[0].As<Napi::Array>();
+    auto input_labels = info[1].As<Napi::Buffer<uint32_t>>();
+    auto cluster_count = static_cast<size_t>(info[2].As<Napi::Number>().Uint32Value());
+    auto output_table = info[3].As<Napi::Array>();
+    auto thread_pool_wrapper = node_api::threading::ThreadPool::Unwrap(info[4].As<Napi::Object>());
+
+    auto col_size = input_table.Length();
+
+    auto data_table = std::vector<std::span<float>>();
+    auto output = std::vector<std::span<float>>();
+    data_table.reserve(col_size);
+    output.reserve(col_size);
+
+    for (auto i = 0; i < col_size; i++) {
+        auto data = input_table.Get(i).As<Napi::Buffer<float>>();
+        auto out = output_table.Get(i).As<Napi::Buffer<float>>();
+        data_table.emplace_back(data.Data(), data.Length());
+        output.emplace_back(out.Data(), out.Length());
+    }
+
+    auto labels = std::span<uint32_t>(input_labels.Data(), input_labels.Length());
+    auto row_size = data_table[0].size();
+
+    auto& pool = thread_pool_wrapper->impl();
+    auto partials = std::vector<PartialAverage>();
+    auto fallback_rows = std::vector<size_t>(cluster_count, 0);
+    if (row_size > 0) {
+        std::random_device random_device;
+        auto seed = (static_cast<uint64_t>(random_device()) << 32) ^ static_cast<uint64_t>(random_device());
+        auto rng = std::mt19937_64(seed);
+        auto row_distribution = std::uniform_int_distribution<size_t>(0, row_size - 1);
+        for (auto& row_index : fallback_rows) {
+            row_index = row_distribution(rng);
+        }
+    }
+
+    auto process_rows = [&](size_t begin, size_t end) -> PartialAverage {
+        auto partial = PartialAverage {
+            .sums = std::make_unique<float[]>(cluster_count * col_size),
+            .counts = std::make_unique<size_t[]>(cluster_count),
+        };
+        for (auto row_index = begin; row_index < end; row_index++) {
+            auto cluster_index = labels[row_index];
+            if (cluster_index >= cluster_count) {
+                continue;
+            }
+
+            partial.counts[cluster_index]++;
+            auto sum_offset = static_cast<size_t>(cluster_index) * col_size;
+            for (auto column_index = 0; column_index < col_size; column_index++) {
+                partial.sums[sum_offset + column_index] += data_table[column_index][row_index];
+            }
+        }
+        return partial;
+    };
+
+    auto reduce_clusters = [&](size_t begin, size_t end) {
+        for (auto cluster_index = begin; cluster_index < end; cluster_index++) {
+            size_t count = 0;
+            for (auto column_index = 0; column_index < col_size; column_index++) {
+                output[column_index][cluster_index] = 0.0f;
+            }
+
+            auto sum_offset = cluster_index * col_size;
+            for (auto& partial : partials) {
+                count += partial.counts[cluster_index];
+                for (auto column_index = 0; column_index < col_size; column_index++) {
+                    output[column_index][cluster_index] += partial.sums[sum_offset + column_index];
+                }
+            }
+
+            if (count == 0) {
+                auto row_index = fallback_rows[cluster_index];
+                for (auto column_index = 0; column_index < col_size; column_index++) {
+                    output[column_index][cluster_index] = data_table[column_index][row_index];
+                }
+            } else {
+                auto inv_count = 1.0f / static_cast<float>(count);
+                for (auto column_index = 0; column_index < col_size; column_index++) {
+                    output[column_index][cluster_index] *= inv_count;
+                }
+            }
+        }
+    };
+
+    {
+        size_t used_thread_count = std::min(pool.thread_count(), row_size);
+
+        auto row_tasks = std::vector<std::future<PartialAverage>>();
+        row_tasks.reserve(used_thread_count);
+
+        auto rows_per_thread = (row_size + used_thread_count - 1) / used_thread_count;
+        // process rows, calculate sum per thread.
+        for (size_t thread_index = 0; thread_index < used_thread_count; thread_index++) {
+            auto begin = thread_index * rows_per_thread;
+            auto end = std::min(begin + rows_per_thread, row_size);
+            row_tasks.push_back(pool.submit_task(process_rows, begin, end));
+        }
+
+        partials.reserve(row_tasks.size());
+        helpers::future::drain_futures(row_tasks, [&](PartialAverage&& data) {
+            partials.push_back(std::move(data));
+        });
+    }
+
+    {
+        size_t used_thread_count = std::min(pool.thread_count(), cluster_count);
+        auto tasks = std::vector<std::future<void>>();
+        tasks.reserve(used_thread_count);
+        auto clusters_per_thread = (cluster_count + used_thread_count - 1) / used_thread_count;
+        // reduce threaded cluster result
+        for (size_t thread_index = 0; thread_index < used_thread_count; thread_index++) {
+            auto begin = thread_index * clusters_per_thread;
+            auto end = std::min(begin + clusters_per_thread, cluster_count);
+            tasks.push_back(pool.submit_task(reduce_clusters, begin, end));
+        }
+        helpers::future::drain_futures(tasks);
+    }
+
+    return env.Null();
+}
+} // namespace node_api::spatial
@@ -0,0 +1,357 @@
+#include <algorithm>
+#include <atomic>
+#include <cassert>
+#include <container_helpers.h>
+#include <future_helpers.h>
+#include <node/api_buffer.h>
+#include <node/api_splat.h>
+#include <node/api_thread_pool.h>
+#include <numeric>
+#include <ranges>
+#include <span>
+#include <splat/splat.h>
+#include <splat/splat_block.h>
+#include <splat/splat_lod.h>
+#include <thread_pool.h>
+
+namespace {
+enum SplatTable {
+    SPLAT_TABLE_MEAN_X = 0,
+    SPLAT_TABLE_MEAN_Y = 1,
+    SPLAT_TABLE_MEAN_Z = 2,
+    SPLAT_TABLE_SCALE_X = 3,
+    SPLAT_TABLE_SCALE_Y = 4,
+    SPLAT_TABLE_SCALE_Z = 5,
+    SPLAT_TABLE_QUAT_X = 6,
+    SPLAT_TABLE_QUAT_Y = 7,
+    SPLAT_TABLE_QUAT_Z = 8,
+    SPLAT_TABLE_QUAT_W = 9,
+    SPLAT_TABLE_COLOR_R = 10,
+    SPLAT_TABLE_COLOR_G = 11,
+    SPLAT_TABLE_COLOR_B = 12,
+    SPLAT_TABLE_OPACITY = 13,
+    SPLAT_TABLE_SH_OFFSET = 14
+};
+
+inline splat::Splat read_splat(const std::vector<std::span<float>>& buffers, size_t sh_size, threading::ThreadPool& pool) {
+    auto thread_count = pool.thread_count();
+    auto read_count = std::atomic<size_t>(0);
+    auto gaussian_count = buffers[0].size();
+    auto gaussians = std::vector<splat::Gaussian>(gaussian_count);
+
+    auto gaussian_per_thread = std::max(gaussian_count / thread_count, static_cast<size_t>(1));
+    auto futures = std::vector<std::future<void>>();
+    auto used_threads = std::min(gaussian_count, thread_count);
+
+    futures.reserve(used_threads);
+    read_count.store(0, std::memory_order_release);
+
+    auto process_gaussian = [&, buffers](size_t start, size_t count) {
+        for (auto i = 0; i < count; i++) {
+            auto index = start + i;
+
+            auto& gaussian = gaussians[index];
+            gaussian.mean = Eigen::Vector3f(
+                buffers[SPLAT_TABLE_MEAN_X][index],
+                buffers[SPLAT_TABLE_MEAN_Y][index],
+                buffers[SPLAT_TABLE_MEAN_Z][index]);
+            gaussian.scale = Eigen::Vector3f(
+                buffers[SPLAT_TABLE_SCALE_X][index],
+                buffers[SPLAT_TABLE_SCALE_Y][index],
+                buffers[SPLAT_TABLE_SCALE_Z][index]);
+            gaussian.rotation = Eigen::Vector4f(
+                buffers[SPLAT_TABLE_QUAT_X][index],
+                buffers[SPLAT_TABLE_QUAT_Y][index],
+                buffers[SPLAT_TABLE_QUAT_Z][index],
+                buffers[SPLAT_TABLE_QUAT_W][index]);
+            gaussian.opacity = buffers[SPLAT_TABLE_OPACITY][index];
+
+            // SH
+            gaussian.sh = ::splat::SH(sh_size + 3);
+            gaussian.sh[0] = buffers[SPLAT_TABLE_COLOR_R][index];
+            gaussian.sh[1] = buffers[SPLAT_TABLE_COLOR_G][index];
+            gaussian.sh[2] = buffers[SPLAT_TABLE_COLOR_B][index];
+            for (auto j = 3; j < gaussian.sh.size(); j++) {
+                gaussian.sh[j] = buffers[SPLAT_TABLE_SH_OFFSET + j - 3][index];
+            }
+
+            gaussian.rotation.normalize();
+            gaussian.compute_covariance();
+            gaussian.compute_bounding_box();
+        }
+        read_count.fetch_add(count, std::memory_order_acq_rel);
+    };
+
+    for (auto i = 0; i < used_threads; i++) {
+        futures.push_back(pool.submit_task(process_gaussian, gaussian_per_thread * i, gaussian_per_thread));
+    }
+
+    process_gaussian(gaussian_per_thread * used_threads, gaussian_count - gaussian_per_thread * used_threads);
+
+    helpers::future::drain_futures(futures);
+
+    assert(read_count.load(std::memory_order_acquire) == gaussian_count);
+
+    auto splat = splat::Splat {
+        .gaussians = std::move(gaussians),
+        .bounding_box = Eigen::AlignedBox3f {},
+    };
+    splat.compute_bounding_box();
+    return splat;
+}
+
+inline std::vector<Napi::Float32Array> write_splat(
+    Napi::Env& env,
+    size_t buffer_count,
+    const splat::Splat& splat,
+    threading::ThreadPool& pool) {
+    auto buffers = std::vector<Napi::Float32Array>();
+    auto buffer_spans = std::vector<std::span<float>>();
+    auto thread_count = pool.thread_count();
+    auto written_count = std::atomic<size_t>(0);
+    buffers.reserve(buffer_count);
+    buffer_spans.reserve(buffer_count);
+
+    for (auto i = 0; i < buffer_count; i++) {
+        buffers.push_back(Napi::Float32Array::New(env, splat.gaussians.size()));
+        buffer_spans.emplace_back(buffers.back().Data(), splat.gaussians.size());
+    }
+
+    auto gaussian_per_thread = std::max(splat.gaussians.size() / thread_count, static_cast<size_t>(1));
+    auto futures = std::vector<std::future<void>>();
+    auto used_threads = std::min(splat.gaussians.size(), thread_count);
+
+    futures.reserve(used_threads);
+    written_count.store(0, std::memory_order_release);
+
+    // transform gaussian to js readable struct
+    auto process_gaussian = [&](size_t start, size_t count) {
+        for (auto i = 0; i < count; i++) {
+            auto index = start + i;
+            auto write_index = index;
+            auto& gaussian = splat.gaussians[index];
+
+            buffer_spans[SPLAT_TABLE_MEAN_X][write_index] = gaussian.mean.x();
+            buffer_spans[SPLAT_TABLE_MEAN_Y][write_index] = gaussian.mean.y();
+            buffer_spans[SPLAT_TABLE_MEAN_Z][write_index] = gaussian.mean.z();
+            buffer_spans[SPLAT_TABLE_SCALE_X][write_index] = gaussian.scale.x();
+            buffer_spans[SPLAT_TABLE_SCALE_Y][write_index] = gaussian.scale.y();
+            buffer_spans[SPLAT_TABLE_SCALE_Z][write_index] = gaussian.scale.z();
+            buffer_spans[SPLAT_TABLE_QUAT_X][write_index] = gaussian.rotation.x();
+            buffer_spans[SPLAT_TABLE_QUAT_Y][write_index] = gaussian.rotation.y();
+            buffer_spans[SPLAT_TABLE_QUAT_Z][write_index] = gaussian.rotation.z();
+            buffer_spans[SPLAT_TABLE_QUAT_W][write_index] = gaussian.rotation.w();
+            buffer_spans[SPLAT_TABLE_COLOR_R][write_index] = gaussian.sh[0];
+            buffer_spans[SPLAT_TABLE_COLOR_G][write_index] = gaussian.sh[1];
+            buffer_spans[SPLAT_TABLE_COLOR_B][write_index] = gaussian.sh[2];
+            buffer_spans[SPLAT_TABLE_OPACITY][write_index] = gaussian.opacity;
+
+            for (auto j = 3; j < gaussian.sh.size(); j++) {
+                buffer_spans[SPLAT_TABLE_SH_OFFSET + j - 3][write_index] = gaussian.sh[j];
+            }
+        }
+
+        written_count.fetch_add(count, std::memory_order_acq_rel);
+    };
+
+    for (auto i = 0; i < used_threads; i++) {
+        futures.push_back(pool.submit_task(process_gaussian, gaussian_per_thread * i, gaussian_per_thread));
+    }
+
+    process_gaussian(gaussian_per_thread * used_threads, splat.gaussians.size() - gaussian_per_thread * used_threads);
+
+    helpers::future::drain_futures(futures);
+
+    assert(written_count.load(std::memory_order_acquire) == splat.gaussians.size());
+
+    return buffers;
+}
+} // namespace
+namespace node_api::splat {
+Napi::Value generate_lod(const Napi::CallbackInfo& info) {
+    auto env = info.Env();
+    if (info.Length() < 7 || !info[0].IsArray() || !info[1].IsNumber() || !info[2].IsBuffer() || !info[3].IsNumber() || !info[4].IsNumber() || !info[5].IsNumber() || !info[6].IsObject()) {
+        Napi::TypeError::New(env, "Wrong Arguments").ThrowAsJavaScriptException();
+        return env.Null();
+    }
+
+    auto& pool = node_api::threading::ThreadPool::Unwrap(info[6].As<Napi::Object>())->impl();
+    size_t thread_count = pool.thread_count();
+
+    auto sh_size = info[1].As<Napi::Number>().Uint32Value();
+    std::vector<::splat::Splat> blocks;
+    // read & create blocks.
+    {
+        auto buffers = std::vector<std::span<float>>();
+        auto array = info[0].As<Napi::Array>();
+        buffers.reserve(sh_size + SPLAT_TABLE_SH_OFFSET);
+        for (auto i = 0; i < sh_size + SPLAT_TABLE_SH_OFFSET; i++) {
+            auto buffer = array[i].AsValue().As<Napi::Buffer<float>>();
+            buffers.emplace_back(buffer.Data(), buffer.Length());
+        }
+        blocks = ::splat::block::split(
+            read_splat(buffers, sh_size, pool),
+            info[3].As<Napi::Number>().DoubleValue());
+    }
+
+    auto level_parameters = info[2].As<Napi::Buffer<::splat::lod::SplatLevelParameters>>();
+    auto levels = level_parameters.Length();
+    auto level_parameters_span = std::span(level_parameters.Data() + 1, level_parameters.Length() - 1);
+
+    auto results = std::vector<::splat::Splat>();
+    auto gaussian_count = std::make_unique<std::vector<uint32_t>>();
+    auto block_boxes = std::make_unique<std::vector<float>>();
+    auto block_refs = std::make_unique<std::vector<uint32_t>>();
+
+    // reverse data, small.
+    results.reserve(blocks.size() * levels);
+    gaussian_count->reserve(blocks.size() * levels);
+    block_boxes->reserve(blocks.size() * 6);
+    block_refs->reserve(blocks.size() * levels);
+
+    {
+        auto futures = std::vector<std::future<::splat::lod::SplatLod>>();
+        auto used_threads = std::min(blocks.size(), thread_count);
+        auto min_size = info[4].As<Napi::Number>().Uint32Value();
+        auto max_step = info[5].As<Napi::Number>().Uint32Value();
+
+        futures.reserve(blocks.size());
+
+        auto process_block = [&, levels, max_step](size_t index) -> ::splat::lod::SplatLod {
+            auto& block = blocks[index];
+            auto lod = ::splat::lod::SplatLod();
+            lod.levels.reserve(levels);
+            lod.splats.reserve(levels);
+            lod.levels.push_back(0);
+            lod.splats.push_back(std::move(block));
+            if (level_parameters_span.size() > 0) {
+                ::splat::lod::generate_lod(index, lod, level_parameters_span, min_size, max_step);
+            }
+            return lod;
+        };
+
+        for (auto i = 0; i < blocks.size(); i++) {
+            futures.push_back(pool.submit_task(process_block, i));
+        }
+
+        helpers::future::drain_futures(futures, [&](::splat::lod::SplatLod block) {
+            auto base_offset = static_cast<uint32_t>(results.size());
+            {
+                auto& bbx = block.splats.front().bounding_box;
+                helpers::container::append_range(*block_boxes, bbx.min());
+                helpers::container::append_range(*block_boxes, bbx.max());
+            }
+            helpers::container::append_range(*gaussian_count, block.splats | std::views::transform([](::splat::Splat& splat) -> uint32_t {
+                return static_cast<uint32_t>(splat.gaussians.size());
+            }));
+            helpers::container::append_range(*block_refs, block.levels | std::views::transform([base_offset](uint32_t ref) -> uint32_t {
+                return ref + base_offset;
+            }));
+            for (auto& splat : block.splats) {
+                results.push_back(std::move(splat));
+            }
+        });
+
+        // release blocks
+        {
+            auto _ = std::move(blocks);
+        }
+    }
+
+    // auto total = std::reduce(gaussian_count->begin(), gaussian_count->end());
+    auto buffers_per_splat = sh_size + SPLAT_TABLE_SH_OFFSET;
+    auto buffers = std::vector<Napi::Float32Array>();
+
+    buffers.reserve(buffers_per_splat * results.size());
+
+    for (auto& splat : results) {
+        helpers::container::append_range(buffers, write_splat(env, buffers_per_splat, splat, pool));
+
+        // free data already transformed.
+        {
+            auto _ = std::move(splat);
+        }
+    }
+
+    auto object = Napi::Object::New(env);
+    {
+        auto data = Napi::Array::New(env, buffers.size());
+        auto offset = 0;
+        object.Set("data", data);
+        for (auto i = 0; i < buffers.size(); i++) {
+            data[i] = buffers[i];
+        }
+        auto _ = std::move(buffers);
+    }
+    object.Set("blockBoxes", node_api::buffer::UniqueVecBufferFinalizer<float>::make_buffer(env, std::move(block_boxes)));
+    object.Set("blockRefs", node_api::buffer::UniqueVecBufferFinalizer<uint32_t>::make_buffer(env, std::move(block_refs)));
+    object.Set("gaussianCount", node_api::buffer::UniqueVecBufferFinalizer<uint32_t>::make_buffer(env, std::move(gaussian_count)));
+    return object;
+}
+
+Napi::Value split(const Napi::CallbackInfo& info) {
+    auto env = info.Env();
+    if (info.Length() < 4 || !info[0].IsArray() || !info[1].IsNumber() || !info[2].IsNumber() || !info[3].IsObject()) {
+        Napi::TypeError::New(env, "Wrong Arguments").ThrowAsJavaScriptException();
+        return env.Null();
+    }
+
+    auto& pool = node_api::threading::ThreadPool::Unwrap(info[3].As<Napi::Object>())->impl();
+    size_t thread_count = pool.thread_count();
+
+    auto sh_size = info[1].As<Napi::Number>().Uint32Value();
+    std::vector<::splat::Splat> blocks;
+    // read & create blocks.
+    {
+        auto buffers = std::vector<std::span<float>>();
+        auto array = info[0].As<Napi::Array>();
+        buffers.reserve(sh_size + SPLAT_TABLE_SH_OFFSET);
+        for (auto i = 0; i < sh_size + SPLAT_TABLE_SH_OFFSET; i++) {
+            auto buffer = array[i].AsValue().As<Napi::Buffer<float>>();
+            buffers.emplace_back(buffer.Data(), buffer.Length());
+        }
+        blocks = ::splat::block::split(
+            read_splat(buffers, sh_size, pool),
+            info[2].As<Napi::Number>().DoubleValue());
+    }
+
+    auto buffers_per_splat = sh_size + SPLAT_TABLE_SH_OFFSET;
+    auto block_boxes = std::make_unique<std::vector<float>>();
+    auto gaussian_count = std::make_unique<std::vector<uint32_t>>();
+    auto buffers = std::vector<Napi::Float32Array>();
+
+    block_boxes->reserve(blocks.size() * 6);
+    gaussian_count->reserve(blocks.size());
+    buffers.reserve(buffers_per_splat * blocks.size());
+
+    for (auto& splat : blocks) {
+
+        helpers::container::append_range(buffers, write_splat(env, buffers_per_splat, splat, pool));
+        gaussian_count->push_back(static_cast<uint32_t>(splat.gaussians.size()));
+        {
+            auto& bbx = splat.bounding_box;
+            helpers::container::append_range(*block_boxes, bbx.min());
+            helpers::container::append_range(*block_boxes, bbx.max());
+        }
+
+        // free data already transformed.
+        {
+            auto _ = std::move(splat);
+        }
+    }
+
+    auto object = Napi::Object::New(env);
+    {
+        auto data = Napi::Array::New(env, buffers.size());
+        auto offset = 0;
+        object.Set("data", data);
+        for (auto i = 0; i < buffers.size(); i++) {
+            data[i] = buffers[i];
+        }
+        auto _ = std::move(buffers);
+    }
+    object.Set("blockBoxes", node_api::buffer::UniqueVecBufferFinalizer<float>::make_buffer(env, std::move(block_boxes)));
+    object.Set("gaussianCount", node_api::buffer::UniqueVecBufferFinalizer<uint32_t>::make_buffer(env, std::move(gaussian_count)));
+    return object;
+}
+} // namespace node_api::splat
@@ -0,0 +1,42 @@
+#include <napi.h>
+#include <node/api_thread_pool.h>
+#include <thread_pool.h>
+
+namespace node_api::threading {
+Napi::FunctionReference ThreadPool::Init(Napi::Env env, Napi::Object exports) {
+    auto ctor = DefineClass(env, "ThreadPool", {
+                                                   InstanceAccessor("threadCount", &ThreadPool::thread_count, nullptr, static_cast<napi_property_attributes>(napi_enumerable | napi_configurable)),
+                                                   InstanceAccessor("taskCount", &ThreadPool::task_count, nullptr, static_cast<napi_property_attributes>(napi_enumerable | napi_configurable)),
+                                               });
+    exports.Set("ThreadPool", ctor);
+    return Napi::Persistent(ctor);
+}
+
+ThreadPool::ThreadPool(const Napi::CallbackInfo& info)
+    : Napi::ObjectWrap<ThreadPool>(info), impl_(calc_thread_count(info)) {
+}
+
+Napi::Value ThreadPool::thread_count(const Napi::CallbackInfo& info) {
+    return Napi::Number::New(info.Env(), static_cast<double>(this->impl_.thread_count()));
+}
+
+Napi::Value ThreadPool::task_count(const Napi::CallbackInfo& info) {
+    return Napi::Number::New(info.Env(), static_cast<double>(this->impl_.task_count()));
+}
+
+::threading::ThreadPool& ThreadPool::impl() noexcept {
+    return this->impl_;
+}
+
+size_t ThreadPool::calc_thread_count(const Napi::CallbackInfo& info) {
+    size_t thread_count = 0u;
+    if (info.Length() > 0 && info[0].IsNumber()) {
+        thread_count = info[0].As<Napi::Number>().Uint32Value();
+    }
+    if (thread_count == 0) {
+        thread_count = std::min(std::max(std::bit_ceil(std::thread::hardware_concurrency()) / 2u, 1u), 16u);
+    }
+
+    return thread_count;
+}
+} // namespace node_api::threading
@@ -0,0 +1,108 @@
+#include <node/api_buffer.h>
+#include <node/api_webp.h>
+#include <span>
+#include <webp/decode.h>
+#include <webp/encode.h>
+
+namespace {
+size_t webp_encode_rgba(const std::span<uint8_t>& data, int32_t width, int32_t height, bool lossless, float quality, uint8_t** output) {
+    WebPPicture pic;
+    WebPConfig config;
+    WebPMemoryWriter wrt;
+    int ok;
+
+    if (output == nullptr) {
+        return 0;
+    }
+
+    if (!WebPConfigPreset(&config, WEBP_PRESET_DEFAULT, quality) ||
+        !WebPPictureInit(&pic)) {
+        return 0; // shouldn't happen, except if system installation is broken
+    }
+
+    config.lossless = static_cast<int32_t>(lossless);
+    config.exact = 1;
+    pic.use_argb = static_cast<int32_t>(lossless);
+    pic.width = width;
+    pic.height = height;
+    pic.writer = WebPMemoryWrite;
+    pic.custom_ptr = &wrt;
+    WebPMemoryWriterInit(&wrt);
+
+    ok = WebPPictureImportRGBA(&pic, data.data(), 4 * width) && WebPEncode(&config, &pic);
+    WebPPictureFree(&pic);
+    if (!ok) {
+        WebPMemoryWriterClear(&wrt);
+        return 0;
+    }
+    *output = wrt.mem;
+    return wrt.size;
+}
+} // namespace
+
+namespace node_api::imaging {
+Napi::Value webp_encode_rgba(const Napi::CallbackInfo& info) {
+    auto env = info.Env();
+    if (info.Length() < 4 || !info[0].IsBuffer() || !info[1].IsNumber() || !info[2].IsNumber() || !info[3].IsNumber()) {
+        Napi::TypeError::New(env, "Wrong Arguments").ThrowAsJavaScriptException();
+        return env.Null();
+    }
+    auto buffer = info[0].As<Napi::Buffer<uint8_t>>();
+    auto width = info[1].As<Napi::Number>().Int32Value();
+    auto height = info[2].As<Napi::Number>().Int32Value();
+    auto quality = info[3].As<Napi::Number>().FloatValue();
+    auto span = std::span<uint8_t>(buffer.Data(), buffer.Length());
+
+    uint8_t* output = nullptr;
+    auto compressed = ::webp_encode_rgba(span, width, height, false, quality, &output);
+    if (compressed == 0) {
+        Napi::TypeError::New(env, "WebP Encode Failed").ThrowAsJavaScriptException();
+        return env.Null();
+    }
+    return node_api::buffer::UniqueBufferFinalizer<uint8_t, decltype(&WebPFree)>::make_buffer(env, std::unique_ptr<uint8_t, decltype(&WebPFree)>(output, WebPFree), compressed);
+}
+
+Napi::Value webp_encode_rgba_lossless(const Napi::CallbackInfo& info) {
+    auto env = info.Env();
+    if (info.Length() < 3 || !info[0].IsBuffer() || !info[1].IsNumber() || !info[2].IsNumber()) {
+        Napi::TypeError::New(env, "Wrong Arguments").ThrowAsJavaScriptException();
+        return env.Null();
+    }
+
+    auto buffer = info[0].As<Napi::Buffer<uint8_t>>();
+    auto width = info[1].As<Napi::Number>().Int32Value();
+    auto height = info[2].As<Napi::Number>().Int32Value();
+    auto span = std::span<uint8_t>(buffer.Data(), buffer.Length());
+
+    uint8_t* output = nullptr;
+    auto compressed = ::webp_encode_rgba(span, width, height, true, 70, &output);
+    if (compressed == 0) {
+        Napi::TypeError::New(env, "WebP Encode Failed").ThrowAsJavaScriptException();
+        return env.Null();
+    }
+    return node_api::buffer::UniqueBufferFinalizer<uint8_t, decltype(&WebPFree)>::make_buffer(env, std::unique_ptr<uint8_t, decltype(&WebPFree)>(output, WebPFree), compressed);
+}
+
+Napi::Value webp_decode_rgba(const Napi::CallbackInfo& info) {
+    auto env = info.Env();
+    if (info.Length() < 1 || !info[0].IsBuffer()) {
+        Napi::TypeError::New(env, "Wrong Arguments").ThrowAsJavaScriptException();
+        return env.Null();
+    }
+    auto buffer = info[0].As<Napi::Buffer<uint8_t>>();
+    int32_t width = 0;
+    int32_t height = 0;
+    auto rgba = WebPDecodeRGBA(buffer.Data(), buffer.Length(), &width, &height);
+    if (rgba == nullptr) {
+        Napi::TypeError::New(env, "WebP Decode Failed").ThrowAsJavaScriptException();
+        return env.Null();
+    }
+    auto object = Napi::Object::New(env);
+    {
+        object.Set("data", node_api::buffer::UniqueBufferFinalizer<uint8_t, decltype(&WebPFree)>::make_buffer(env, std::unique_ptr<uint8_t, decltype(&WebPFree)>(rgba, WebPFree), width * height * 4));
+        object.Set("width", Napi::Number::New(env, width));
+        object.Set("height", Napi::Number::New(env, height));
+    }
+    return object;
+}
+} // namespace node_api::imaging
@@ -0,0 +1,57 @@
+#include <concepts>
+#include <napi.h>
+#include <node/api_avif.h>
+#include <node/api_spatial.h>
+#include <node/api_splat.h>
+#include <node/api_thread_pool.h>
+#include <node/api_webp.h>
+#include <type_traits>
+#include <utility>
+
+namespace {
+struct InstanceData {
+    Napi::FunctionReference thread_pool;
+};
+
+template<typename F>
+    requires requires(F f, const Napi::CallbackInfo& info) {
+        { f(info) } -> std::same_as<Napi::Value>;
+    }
+auto wrap_call(F&& f) {
+    return [f = std::forward<F>(f)](const Napi::CallbackInfo& info) {
+        try {
+            return f(info);
+        } catch (const Napi::Error&) {
+            throw;
+        } catch (const std::exception& ex) {
+            auto env = info.Env();
+            Napi::Error::New(env, ex.what()).ThrowAsJavaScriptException();
+            return env.Null();
+        } catch (...) {
+            auto env = info.Env();
+            Napi::Error::New(env, "Unknown runtime exception").ThrowAsJavaScriptException();
+            return env.Null();
+        }
+    };
+}
+} // namespace
+
+static Napi::Object Init(Napi::Env env, Napi::Object exports) {
+    exports.Set("generate_splat_lod", Napi::Function::New(env, wrap_call(node_api::splat::generate_lod)));
+    exports.Set("split_splat", Napi::Function::New(env, wrap_call(node_api::splat::split)));
+    exports.Set("webp_encode_rgba", Napi::Function::New(env, wrap_call(node_api::imaging::webp_encode_rgba)));
+    exports.Set("webp_encode_rgba_lossless", Napi::Function::New(env, wrap_call(node_api::imaging::webp_encode_rgba_lossless)));
+    exports.Set("webp_decode_rgba", Napi::Function::New(env, wrap_call(node_api::imaging::webp_decode_rgba)));
+    exports.Set("avif_encode_rgba", Napi::Function::New(env, wrap_call(node_api::imaging::avif_encode_rgba)));
+    exports.Set("avif_encode_rgba_batched", Napi::Function::New(env, wrap_call(node_api::imaging::avif_encode_rgba_batched)));
+    exports.Set("avif_decode_rgba", Napi::Function::New(env, wrap_call(node_api::imaging::avif_decode_rgba)));
+    exports.Set("avif_decode_rgba_batched", Napi::Function::New(env, wrap_call(node_api::imaging::avif_decode_rgba_batched)));
+    exports.Set("cluster_average", Napi::Function::New(env, wrap_call(node_api::spatial::cluster_average)));
+    auto instance_data = new InstanceData {
+        .thread_pool = node_api::threading::ThreadPool::Init(env, exports)
+    };
+    env.SetInstanceData<InstanceData>(instance_data);
+    return exports;
+}
+
+NODE_API_MODULE(addon_napi, Init)
@@ -0,0 +1,188 @@
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <cstring>
+#include <eigen3/Eigen/Dense>
+#include <numbers>
+#include <splat/splat.h>
+#include <stdexcept>
+#include <tuple>
+#include <utility>
+
+namespace splat {
+SH::SH() noexcept : size_(0), ptr(nullptr) {
+}
+
+SH::SH(size_t size) : size_(size), ptr(size > 0 ? std::make_unique<float[]>(size) : nullptr) {
+}
+
+SH::SH(const SH& other) : size_(other.size_), ptr(other.size_ > 0 ? std::make_unique<float[]>(other.size_) : nullptr) {
+    if (this->ptr != nullptr) {
+        std::copy(other.ptr.get(), other.ptr.get() + this->size_, this->ptr.get());
+    }
+}
+
+SH::SH(SH&& other) noexcept : size_(std::exchange(other.size_, 0)),
+                              ptr(std::exchange(other.ptr, nullptr)) {
+}
+
+SH& SH::operator=(const SH& other) {
+    if (this != &other) {
+        this->reset();
+        if (other.size_ > 0) {
+            this->ptr = std::make_unique<float[]>(other.size_);
+            std::copy(other.ptr.get(), other.ptr.get() + other.size_, this->ptr.get());
+        }
+        this->size_ = other.size_;
+    }
+    return *this;
+}
+
+SH& SH::operator=(SH&& other) noexcept {
+    if (this != &other) {
+        this->size_ = std::exchange(other.size_, 0);
+        this->ptr = std::exchange(other.ptr, nullptr);
+    }
+    return *this;
+}
+
+void SH::swap(SH& other) noexcept {
+    std::swap(this->ptr, other.ptr);
+    std::swap(this->size_, other.size_);
+}
+
+SH& SH::set_zero() noexcept {
+    if (this->ptr) {
+        std::memset(this->ptr.get(), 0, sizeof(float) * this->size_);
+    }
+
+    return *this;
+}
+
+SH& SH::add_multiplied(const SH& other, float scalar) noexcept {
+    assert(other.size_ == this->size_);
+
+    for (auto i = 0; i < this->size_; i++) {
+        this->ptr[i] += other.ptr[i] * scalar;
+    }
+
+    return *this;
+}
+
+const float& SH::operator[](size_t index) const {
+    return this->ptr[index];
+}
+
+float& SH::operator[](size_t index) {
+    return this->ptr[index];
+}
+
+void SH::reset() noexcept {
+    this->ptr.reset();
+    this->size_ = 0;
+}
+
+float* SH::data() noexcept {
+    return this->ptr.get();
+}
+
+size_t SH::size() const noexcept {
+    return this->size_;
+}
+
+SH::~SH() noexcept {
+}
+
+void Gaussian::compute_bounding_box(float k) {
+    auto sigma = this->scale.cwiseAbs();
+    this->bounding_box = Eigen::AlignedBox3f(this->mean - k * sigma, this->mean + k * sigma);
+}
+
+void Gaussian::compute_covariance() {
+    Eigen::Matrix3f s = Eigen::Matrix3f::Identity();
+
+    s(0, 0) = this->scale.x();
+    s(1, 1) = this->scale.y();
+    s(2, 2) = this->scale.z();
+
+    auto r = Eigen::Matrix3f(Eigen::Quaternionf(this->rotation));
+
+    Eigen::Matrix3f t = r * s;
+
+    this->covariance = t * t.transpose();
+}
+
+void Gaussian::decompose_covariance() {
+    auto matrix = this->covariance;
+    Eigen::SelfAdjointEigenSolver<Eigen::Matrix3f> eigen_solver(matrix);
+    auto eigen_values = eigen_solver.eigenvalues();
+    auto eigen_vectors = eigen_solver.eigenvectors();
+
+    if (eigen_values.hasNaN()) {
+        throw std::runtime_error("Found Nans in covariance decomposing!");
+    }
+
+    int loops = 0;
+    while (eigen_values[0] == 0 || eigen_values[1] == 0 || eigen_values[2] == 0) {
+        matrix(0, 0) += std::max(matrix(0, 0) * 0.0001f, std::numeric_limits<float>::epsilon());
+        matrix(1, 1) += std::max(matrix(1, 1) * 0.0001f, std::numeric_limits<float>::epsilon());
+        matrix(2, 2) += std::max(matrix(2, 2) * 0.0001f, std::numeric_limits<float>::epsilon());
+
+        Eigen::SelfAdjointEigenSolver<Eigen::Matrix3f> eigen_solver2(matrix);
+        eigen_values = eigen_solver2.eigenvalues();
+        eigen_vectors = eigen_solver2.eigenvectors();
+
+        loops++;
+    }
+
+    auto v1 = eigen_vectors.col(0);
+    auto v2 = eigen_vectors.col(1);
+    auto v3 = eigen_vectors.col(2);
+
+    auto test = v1.cross(v2);
+
+    if (test.dot(v3) < 0) {
+        eigen_vectors.col(2) *= -1;
+    }
+
+    float a = std::sqrt(std::abs(eigen_values.x()));
+    float b = std::sqrt(std::abs(eigen_values.y()));
+    float c = std::sqrt(std::abs(eigen_values.z()));
+
+    this->scale = { a, b, c };
+    auto q = Eigen::Quaternionf(eigen_vectors);
+    this->rotation = { q.x(), q.y(), q.z(), q.w() };
+    this->rotation.normalize();
+}
+
+float Gaussian::area() const {
+    constexpr float P = 1.6075f;
+    auto numerator = std::pow(this->scale.x() * this->scale.y(), P) + std::pow(this->scale.x() * this->scale.z(), P) +
+                     std::pow(this->scale.y() * this->scale.z(), P);
+    return 4.0f * std::numbers::pi_v<float> * std::pow(numerator / 3.0f, 1.0f / P);
+}
+
+void Splat::compute_bounding_box() {
+    this->bounding_box = Eigen::AlignedBox3f();
+    for (auto& gaussian : this->gaussians) {
+        this->bounding_box.extend(gaussian.bounding_box);
+    }
+}
+
+void Splat::compute_compact_bounding_box() {
+    this->bounding_box = Eigen::AlignedBox3f();
+    auto ref = std::vector<std::tuple<size_t, float>>();
+    ref.reserve(this->gaussians.size());
+    for (auto i = 0; i < this->gaussians.size(); i++) {
+        ref.emplace_back(i, this->gaussians[i].area());
+    }
+    std::sort(ref.begin(), ref.end(), [](auto& a, auto& b) -> bool { return std::get<1>(a) < std::get<1>(b); });
+    auto needed = static_cast<size_t>(static_cast<double>(this->gaussians.size()) * 0.95);
+    if (needed == 0) {
+        needed = this->gaussians.size();
+    }
+    for (auto i = 0; i < needed; i++) {
+        this->bounding_box.extend(this->gaussians[std::get<0>(ref[i])].bounding_box);
+    }
+}
+} // namespace splat
@@ -0,0 +1,137 @@
+#include <array>
+#include <container_helpers.h>
+#include <eigen3/Eigen/Dense>
+#include <queue>
+#include <ranges>
+#include <splat/splat_block.h>
+#include <utility>
+#include <vector>
+
+namespace {
+struct RefSplat {
+    const splat::Splat& source;
+    std::vector<size_t> gaussians;
+    Eigen::AlignedBox3f box;
+    bool need_split;
+
+    RefSplat(const splat::Splat& splat) : source(splat), box(splat.bounding_box), need_split(true) {
+        this->gaussians.reserve(splat.gaussians.size());
+        for (auto i = 0; i < splat.gaussians.size(); i++) {
+            this->gaussians.push_back(i);
+        }
+    }
+    RefSplat(RefSplat& splat) : source(splat.source), need_split(true) {}
+    RefSplat(RefSplat&& splat) = default;
+
+    void compute_bounding_box() {
+        this->box = Eigen::AlignedBox3f();
+        for (auto index : this->gaussians) {
+            this->box.extend(this->source.gaussians[index].bounding_box);
+        }
+    }
+
+    void compute_need_split(const Eigen::AlignedBox3f& parent_box, size_t max_block_size) {
+        this->need_split = this->gaussians.size() > max_block_size &&
+                           ((this->box.min() - parent_box.min()).norm() > 0.001 ||
+                               (this->box.max() - parent_box.max()).norm() > 0.001);
+    }
+
+    splat::Splat to_owned() {
+        std::vector<splat::Gaussian> gaussians;
+        std::vector<float> sh;
+        gaussians.reserve(this->gaussians.size());
+        for (auto index : this->gaussians) {
+            gaussians.push_back(this->source.gaussians[index]);
+        }
+        {
+            auto _ = std::move(this->gaussians);
+        }
+        auto result = splat::Splat {
+            .gaussians = std::move(gaussians),
+            .bounding_box = this->box,
+        };
+        result.compute_compact_bounding_box();
+        return result;
+    }
+};
+
+std::array<Eigen::AlignedBox3f, 8> split_box(const Eigen::AlignedBox3f& box) {
+    auto center = box.center();
+    return {
+        Eigen::AlignedBox3f().extend(center).extend(box.corner(Eigen::AlignedBox3f::CornerType::BottomLeftFloor)),
+        Eigen::AlignedBox3f().extend(center).extend(box.corner(Eigen::AlignedBox3f::CornerType::BottomRightFloor)),
+        Eigen::AlignedBox3f().extend(center).extend(box.corner(Eigen::AlignedBox3f::CornerType::TopLeftFloor)),
+        Eigen::AlignedBox3f().extend(center).extend(box.corner(Eigen::AlignedBox3f::CornerType::TopRightFloor)),
+        Eigen::AlignedBox3f().extend(center).extend(box.corner(Eigen::AlignedBox3f::CornerType::BottomLeftCeil)),
+        Eigen::AlignedBox3f().extend(center).extend(box.corner(Eigen::AlignedBox3f::CornerType::BottomRightCeil)),
+        Eigen::AlignedBox3f().extend(center).extend(box.corner(Eigen::AlignedBox3f::CornerType::TopLeftCeil)),
+        Eigen::AlignedBox3f().extend(center).extend(box.corner(Eigen::AlignedBox3f::CornerType::TopRightCeil)),
+    };
+}
+
+std::array<RefSplat, 8> split_block(RefSplat& splat, size_t max_block_size) {
+    auto boxes = split_box(splat.box);
+    std::array<RefSplat, 8> result { RefSplat(splat),
+        RefSplat(splat),
+        RefSplat(splat),
+        RefSplat(splat),
+        RefSplat(splat),
+        RefSplat(splat),
+        RefSplat(splat),
+        RefSplat(splat) };
+
+    for (auto& r : result) {
+        r.gaussians.reserve(splat.gaussians.size() / 4);
+    }
+
+    for (auto index : splat.gaussians) {
+        auto& box = splat.source.gaussians[index].bounding_box;
+        auto max_interaction = 0.0;
+        auto max_interaction_index = 0;
+        for (auto i = 0; i < 8; i++) {
+            auto intersection = boxes[i].intersection(box);
+            if (!intersection.isEmpty()) {
+                auto current = intersection.volume();
+                if (current > max_interaction) {
+                    max_interaction = current;
+                    max_interaction_index = i;
+                }
+            }
+        }
+        result[max_interaction_index].gaussians.push_back(index);
+    }
+
+    for (auto i = 0; i < result.size(); i++) {
+        result[i].box = boxes[i];
+        result[i].compute_need_split(splat.box, max_block_size);
+    }
+
+    return result;
+}
+} // namespace
+
+namespace splat::block::detail {
+std::vector<Splat> split(const Splat& splat, size_t max_block_size) {
+    if (splat.gaussians.size() <= max_block_size) {
+        return std::vector({ splat });
+    }
+
+    std::queue<RefSplat> queue;
+    std::vector<Splat> results;
+    queue.emplace(splat);
+
+    while (!queue.empty()) {
+        auto r = std::move(queue.front());
+        queue.pop();
+        if (r.need_split) {
+            for (auto& splitted : split_block(r, max_block_size)) {
+                queue.emplace(std::move(splitted));
+            }
+        } else if (r.gaussians.size() > 0) {
+            results.emplace_back(r.to_owned());
+        }
+    }
+
+    return results;
+}
+} // namespace splat::block::detail
@@ -0,0 +1,396 @@
+#include <algorithm>
+#include <bit>
+#include <cmath>
+#include <container_helpers.h>
+#include <eigen3/Eigen/Dense>
+#include <initializer_list>
+#include <nanoflann.hpp>
+#include <numbers>
+#include <ranges>
+#include <splat/splat.h>
+#include <splat/splat_lod.h>
+#include <thread>
+#include <tuple>
+#include <vector>
+
+namespace {
+constexpr double TWO_PI_POW_1_5 = 0x1.f7fccep+3;
+constexpr double EPS_COV = 1e-8;
+constexpr double LOG_2PI = 1.8378771;
+constexpr float OPACITY_PRUNE_THRESHOLD = 0.005f;
+
+class GaussianCloud {
+public:
+    GaussianCloud(const splat::Splat& splat) : splat(splat) {}
+
+    inline size_t kdtree_get_point_count() const { return splat.gaussians.size(); }
+    inline float kdtree_get_pt(size_t i, int d) const { return splat.gaussians[i].mean[d]; }
+    template<class BBOX>
+    bool kdtree_get_bbox(BBOX&) const { return false; }
+
+private:
+    const splat::Splat& splat;
+};
+
+struct GaussianCache {
+    Eigen::Matrix3d rotation;
+    Eigen::Matrix3d rotation_transpose;
+    Eigen::Matrix3d sigma;
+    Eigen::Vector3d variances;
+    Eigen::Vector3d inv_diag;
+    double log_det;
+    double mass;
+};
+
+using KDTree = nanoflann::KDTreeSingleIndexAdaptor<
+    nanoflann::L2_Simple_Adaptor<float, GaussianCloud>,
+    GaussianCloud, 3>;
+
+// https://github.com/graphdeco-inria/gaussian-hierarchy/
+static splat::Gaussian merge_gaussians(const splat::Splat& input, const std::vector<uint32_t>& idx, float scale_boost) {
+    splat::Gaussian out;
+    out.mean.setZero();
+    out.covariance.setZero();
+    out.scale.setZero();
+    out.rotation.setZero();
+    out.sh = splat::SH(input.gaussians[0].sh.size());
+    out.sh.set_zero();
+    out.opacity = 0.0;
+
+    std::vector<float> weights(idx.size());
+    float total_weight = 0.0f;
+
+    for (auto i = 0; i < idx.size(); i++) {
+        const auto& g = input.gaussians[idx[i]];
+        weights[i] = g.area() * g.opacity;
+        total_weight += weights[i];
+    }
+
+    for (auto& weight : weights) {
+        weight /= total_weight;
+    }
+
+    for (auto i = 0; i < idx.size(); i++) {
+        const auto& g = input.gaussians[idx[i]];
+        out.mean += g.mean * weights[i];
+        out.sh.add_multiplied(g.sh, weights[i]);
+    }
+
+    for (auto i = 0; i < idx.size(); i++) {
+        const auto& g = input.gaussians[idx[i]];
+        auto d = g.mean - out.mean;
+        out.covariance += weights[i] * (g.covariance + d * d.transpose());
+    }
+
+    out.decompose_covariance();
+    out.scale *= scale_boost;
+    out.opacity = total_weight / out.area();
+
+    if (out.opacity > 1.0f) {
+        out.scale *= std::pow(out.opacity, 1.0f / 3.0f);
+        out.opacity = 1.0;
+    }
+    out.compute_covariance();
+
+    return out;
+}
+
+static Eigen::Matrix3d sigma_from_rotation_variances(const Eigen::Matrix3d& rotation, const Eigen::Vector3d& variances) {
+    Eigen::Matrix3d s = Eigen::Matrix3d::Identity();
+
+    s(0, 0) = variances.x();
+    s(1, 1) = variances.y();
+    s(2, 2) = variances.z();
+
+    return rotation * s * rotation.transpose();
+}
+
+static void build_cache(const splat::Splat& splat, std::vector<GaussianCache>& caches) {
+    static Eigen::Vector3d VECTOR_ONE(1.0, 1.0, 1.0);
+    static Eigen::Vector3d VECTOR_EPS(EPS_COV, EPS_COV, EPS_COV);
+
+    for (auto& gaussian : splat.gaussians) {
+        auto rotation = Eigen::Matrix3d(Eigen::Quaterniond(gaussian.rotation.w(), gaussian.rotation.x(), gaussian.rotation.y(), gaussian.rotation.z()));
+        auto t = rotation.transpose();
+        Eigen::Vector3f vf = gaussian.scale.cwiseMax(1e-12f).cwiseSquare();
+        Eigen::Vector3d v = Eigen::Vector3d(vf.x(), vf.y(), vf.z()) + VECTOR_EPS;
+        Eigen::Vector3d safe_v = v.cwiseMax(1e-30);
+
+        caches.push_back(GaussianCache {
+            .rotation = rotation,
+            .rotation_transpose = t,
+            .sigma = sigma_from_rotation_variances(rotation, v),
+            .variances = v,
+            .inv_diag = VECTOR_ONE.array() / safe_v.array(),
+            .log_det = safe_v.array().log().sum(),
+            .mass = TWO_PI_POW_1_5 * gaussian.opacity * gaussian.scale.prod() + 1e-12 });
+    }
+}
+
+std::vector<Eigen::Vector3d> make_gaussian_samples(size_t n, uint32_t seed) {
+    std::vector<Eigen::Vector3d> result;
+    result.reserve(n);
+    auto rand = [&seed]() -> double {
+        auto t = (seed += 0x6d2b79f5);
+        t = (t ^ (t >> 15)) * (t | 1);
+        t ^= t + (t ^ (t >> 7)) * (t | 61);
+        return static_cast<double>((t ^ (t >> 14)) >> 0) / 4294967296.0;
+    };
+
+    for (auto i = 0; i < n; i++) {
+        auto u1 = std::max(rand(), 1e-12);
+        auto u2 = rand();
+        auto u3 = std::max(rand(), 1e-12);
+        auto u4 = rand();
+        auto r1 = std::sqrt(-2 * std::log(u1));
+        auto t1 = 2 * std::numbers::pi_v<double> * u2;
+        auto r2 = std::sqrt(-2 * std::log(u3));
+        auto t2 = 2 * std::numbers::pi_v<double> * u4;
+        result.emplace_back(r1 * std::cos(t1), r1 * std::sin(t1), r2 * std::cos(t2));
+    }
+
+    return result;
+}
+
+static double gauss_log_pdf_diagrot(const Eigen::Vector3d& v, const Eigen::Vector3d& m, const Eigen::Matrix3d& r, const Eigen::Vector3d& inv, double log_det) {
+    Eigen::Vector3d d = v - m;
+    Eigen::Vector3d y = r * d;
+    double quad = y.cwiseSquare().transpose() * inv;
+    return -0.5 * (3.0 * LOG_2PI + log_det + quad);
+}
+
+static double log_add_exp(double a, double b) {
+    if (-std::numeric_limits<double>::infinity() == a) {
+        return b;
+    }
+    if (-std::numeric_limits<double>::infinity() == b) {
+        return a;
+    }
+    auto m = std::max(a, b);
+    return m + std::log(std::exp(a - m) + std::exp(b - m));
+};
+
+static double compute_edge_cost(const splat::Splat& splat, const std::vector<GaussianCache>& caches, const std::vector<Eigen::Vector3d>& samples, std::tuple<uint32_t, uint32_t>& edge) {
+    auto [i, j] = edge;
+    auto& u = splat.gaussians[i];
+    auto& v = splat.gaussians[j];
+    auto& cache_i = caches[i];
+    auto& cache_j = caches[j];
+    auto wi = cache_i.mass;
+    auto wj = cache_j.mass;
+    auto w = wi + wj;
+    auto w_safe = w > 0 ? w : 1;
+    auto pi = wi / w_safe;
+    pi = std::max(1e-12, std::min(1.0 - 1e-12, pi));
+    auto pj = 1 - pi;
+    auto log_pi = std::log(pi);
+    auto log_pj = std::log(pj);
+    Eigen::Vector3d u_mean(u.mean.x(), u.mean.y(), u.mean.z());
+    Eigen::Vector3d v_mean(v.mean.x(), v.mean.y(), v.mean.z());
+
+    auto mm = pi * u_mean + pj * v_mean;
+    auto di = u_mean - mm;
+    auto dj = v_mean - mm;
+    Eigen::Matrix3d sigma = pi * cache_i.sigma + pj * cache_j.sigma;
+
+    sigma += pi * (di * di.transpose()) + pj * (dj * dj.transpose());
+
+    // Force symmetry + regularize
+    sigma(0, 1) = sigma(1, 0) = 0.5 * (sigma(0, 1) + sigma(1, 0));
+    sigma(0, 2) = sigma(2, 0) = 0.5 * (sigma(0, 2) + sigma(2, 0));
+    sigma(1, 2) = sigma(2, 1) = 0.5 * (sigma(1, 2) + sigma(2, 1));
+
+    sigma(0, 0) += EPS_COV;
+    sigma(1, 1) += EPS_COV;
+    sigma(2, 2) += EPS_COV;
+
+    // det should never be neg.
+    auto det_m = std::max(sigma.determinant(), 1e-30);
+    auto log_det_m = std::log(det_m);
+    // E_p[-log q_m] computed analytically as entropy of merged Gaussian
+    auto ep_neg_log_q = 0.5 * (3.0 * LOG_2PI + log_det_m + 3.0);
+
+    // Sample from each component separately with same z-vectors
+    Eigen::Vector3d std_i = cache_i.variances.cwiseMax(0).cwiseSqrt();
+    Eigen::Vector3d std_j = cache_j.variances.cwiseMax(0).cwiseSqrt();
+
+    auto sum_log_p_on_i = 0.0;
+    auto sum_log_p_on_j = 0.0;
+
+    for (auto& sample : samples) {
+        Eigen::Vector3d ti = std_i.array() * sample.array();
+        Eigen::Vector3d tj = std_j.array() * sample.array();
+        Eigen::Vector3d xi = u_mean + cache_i.rotation * ti;
+        Eigen::Vector3d xj = v_mean + cache_j.rotation * tj;
+
+        // Evaluate log p_ij at samples from component i
+        auto log_ni_on_i = gauss_log_pdf_diagrot(xi, u_mean, cache_i.rotation_transpose, cache_i.inv_diag, cache_i.log_det);
+        auto log_nj_on_i = gauss_log_pdf_diagrot(xi, v_mean, cache_j.rotation_transpose, cache_j.inv_diag, cache_j.log_det);
+        sum_log_p_on_i += log_add_exp(log_pi + log_ni_on_i, log_pj + log_nj_on_i);
+
+        // Evaluate log p_ij at samples from component j
+        auto log_ni_on_j = gauss_log_pdf_diagrot(xj, u_mean, cache_i.rotation_transpose, cache_i.inv_diag, cache_i.log_det);
+        auto log_nj_on_j = gauss_log_pdf_diagrot(xj, v_mean, cache_j.rotation_transpose, cache_j.inv_diag, cache_j.log_det);
+        sum_log_p_on_j += log_add_exp(log_pi + log_ni_on_j, log_pj + log_nj_on_j);
+    }
+
+    auto ei = sum_log_p_on_i / samples.size();
+    auto ej = sum_log_p_on_j / samples.size();
+    auto ep_log_p = pi * ei + pj * ej;
+    auto geo = ep_log_p + ep_neg_log_q;
+
+    auto c_sh = 0.0;
+
+    for (auto i = 0; i < u.sh.size(); i++) {
+        auto d = u.sh[i] - v.sh[i];
+        c_sh += (d * d);
+    }
+
+    return geo + c_sh;
+}
+
+inline static bool validate_gaussian(const splat::Gaussian& gaussian) {
+    return gaussian.scale.cwiseNotEqual(0.0f).all() &&
+           gaussian.opacity >= OPACITY_PRUNE_THRESHOLD;
+}
+} // namespace
+
+namespace splat::lod::detail {
+Splat reduce_gaussians(size_t id, const Splat& input, size_t target_count, float scale_boost, size_t max_step) {
+    assert(input.gaussians.size() > target_count);
+
+    Splat current;
+    uint32_t knn_k = 16;
+    auto samples = make_gaussian_samples(1, 0);
+    std::vector<uint32_t> out_indices;
+    std::vector<float> out_distances;
+    std::vector<size_t> order;
+    std::vector<std::tuple<uint32_t, uint32_t>> edges;
+    std::vector<std::tuple<uint32_t, uint32_t>> pairs;
+    std::vector<bool> used;
+    std::vector<double> costs;
+    std::vector<GaussianCache> caches;
+    auto step = 0;
+
+    // prune gaussians
+    current.gaussians.reserve(input.gaussians.size());
+    helpers::container::append_range(current.gaussians, input.gaussians | std::views::filter(validate_gaussian));
+
+    while (current.gaussians.size() > target_count && current.gaussians.size() > 0 && step < max_step) {
+        GaussianCloud cloud(current);
+        KDTree tree(3, cloud,
+            nanoflann::KDTreeSingleIndexAdaptorParams(
+                10, nanoflann::KDTreeSingleIndexAdaptorFlags::None,
+                std::max(std::bit_floor(std::thread::hardware_concurrency()) / 2u, 1u)));
+
+        auto k_eff = std::min<uint32_t>(std::max<uint32_t>(1, knn_k), std::max<uint32_t>(1, current.gaussians.size() - 1));
+        size_t query_count = static_cast<size_t>(std::min<size_t>(current.gaussians.size(), k_eff + 1));
+        auto max_pairs = current.gaussians.size() - target_count;
+
+        assert(max_pairs > 0);
+
+        Splat generated;
+
+        edges.clear();
+        edges.reserve(current.gaussians.size() * k_eff);
+        caches.clear();
+        caches.reserve(current.gaussians.size());
+        generated.gaussians.reserve(current.gaussians.size() / 2);
+
+        build_cache(current, caches);
+
+        // knn search
+        for (auto i = 0; i < current.gaussians.size(); i++) {
+            out_indices.assign(query_count, 0);
+            out_distances.assign(query_count, 0.0);
+            auto result_count = tree.knnSearch(current.gaussians[i].mean.data(), query_count, out_indices.data(), out_distances.data());
+            const size_t take = std::min<size_t>(k_eff, result_count > 0 ? result_count - 1 : size_t { 0 });
+            for (auto j = 0; j < take; j++) {
+                auto neighbor = static_cast<int>(out_indices[j + 1]);
+                if (neighbor <= i) {
+                    continue;
+                }
+                edges.emplace_back(i, neighbor);
+            }
+        }
+
+        assert(edges.size() > 0);
+
+        // compute edge costs.
+        {
+
+            costs.clear();
+            costs.reserve(edges.size());
+            for (auto& edge : edges) {
+                costs.push_back(compute_edge_cost(current, caches, samples, edge));
+            }
+        }
+
+        // take pairs
+        {
+            order.clear();
+            order.reserve(edges.size());
+            for (auto i = 0; i < costs.size(); ++i) {
+                if (std::isfinite(costs[i])) {
+                    order.push_back(i);
+                }
+            }
+
+            std::stable_sort(order.begin(), order.end(), [&](const size_t lhs, const size_t rhs) {
+                return costs[lhs] < costs[rhs];
+            });
+
+            used.assign(current.gaussians.size(), false);
+            pairs.clear();
+            pairs.reserve(max_pairs);
+
+            for (auto edge_idx : order) {
+                const auto [u, v] = edges[edge_idx];
+                if (used[u] || used[v]) {
+                    continue;
+                }
+                used[u] = true;
+                used[v] = true;
+                pairs.emplace_back(u, v);
+                if (pairs.size() >= max_pairs) {
+                    break;
+                }
+            }
+        }
+
+        // merge pairs
+        for (auto& pair : pairs) {
+            auto [u, v] = pair;
+            auto g = merge_gaussians(current, { u, v }, scale_boost);
+            if (validate_gaussian(g)) {
+                generated.gaussians.push_back(std::move(g));
+            }
+        }
+
+        // keep unused gaussian
+        for (auto i = 0; i < current.gaussians.size(); i++) {
+            if (!used[i]) {
+                generated.gaussians.push_back(current.gaussians[i]);
+            }
+        }
+
+        current = std::move(generated);
+
+        // trace hard working on current block.
+        if (step >= max_step / 2) {
+            printf("block: %zu, step: %d, required: %zu, current: %zu(%.2f%%)\n",
+                id, step, target_count, current.gaussians.size(),
+                static_cast<double>(current.gaussians.size()) / static_cast<double>(target_count) * 100.0);
+        }
+
+        // cannot reduce any more...
+        if (pairs.size() == 0) {
+            break;
+        }
+
+        step++;
+    }
+    current.gaussians.shrink_to_fit();
+    return current;
+}
+} // namespace splat::lod::detail
@@ -0,0 +1,134 @@
+#include <functional>
+#include <memory>
+#include <mutex>
+#include <shared_mutex>
+#include <thread>
+#include <thread_pool.h>
+#include <utility>
+
+namespace threading {
+ThreadPool::ThreadPool(size_t thread_count)
+    : state(std::make_unique<State>()) {
+    if (thread_count == 0) {
+        thread_count = 1;
+    }
+
+    this->state->workers.reserve(thread_count);
+
+    for (auto i = 0; i < thread_count; i++) {
+        auto worker = std::make_unique<Worker>();
+        worker->thread = std::thread(ThreadPool::worker_loop, this->state.get(), worker.get());
+        this->state->workers.push_back(std::move(worker));
+    }
+}
+
+ThreadPool::ThreadPool(ThreadPool&& other) noexcept
+    : state(std::move(other.state)) {
+}
+
+ThreadPool& ThreadPool::operator=(ThreadPool&& other) noexcept {
+    if (this != &other) {
+        this->stop();
+        this->state = std::move(other.state);
+    }
+    return *this;
+}
+
+size_t ThreadPool::thread_count() const noexcept {
+    if (this->state == nullptr) {
+        return 0;
+    }
+
+    auto lk = std::shared_lock(this->state->worker_mutex);
+    return this->state->workers.size();
+}
+
+size_t ThreadPool::task_count() const noexcept {
+    if (this->state == nullptr) {
+        return 0;
+    }
+
+    auto lk = std::shared_lock(this->state->queue_mutex);
+    return this->state->tasks.size();
+}
+
+void ThreadPool::stop() noexcept {
+    if (this->state == nullptr) {
+        return;
+    }
+
+    this->state->stop();
+}
+
+ThreadPool::~ThreadPool() noexcept {
+    this->stop();
+}
+
+void ThreadPool::State::stop() noexcept {
+    bool expected = false;
+    if (!this->stopped.compare_exchange_strong(expected, true, std::memory_order_relaxed)) {
+        return;
+    }
+
+    {
+        auto lk = std::unique_lock(this->queue_mutex);
+        auto _ = std::move(this->tasks);
+        lk.unlock();
+    }
+
+    {
+        auto lk = std::lock_guard(this->worker_mutex);
+        for (auto& worker : this->workers) {
+            worker->stop();
+        }
+    }
+
+    this->cv.notify_all();
+
+    for (auto& worker : this->workers) {
+        worker->clean_up();
+    }
+}
+
+ThreadPool::State::~State() noexcept {
+    this->stop();
+}
+
+void ThreadPool::Worker::stop() noexcept {
+    this->stopped.store(true, std::memory_order_relaxed);
+}
+
+void ThreadPool::Worker::clean_up() noexcept {
+    this->stop();
+
+    if (this->thread.joinable()) {
+        try {
+            this->thread.join();
+        } catch (...) {
+        }
+    }
+}
+
+ThreadPool::Worker::~Worker() noexcept {
+    this->clean_up();
+}
+
+void ThreadPool::worker_loop(ThreadPool::State* state, ThreadPool::Worker* worker) {
+    while (true) {
+        ThreadPool::Task task;
+        {
+            auto lk = std::unique_lock(state->queue_mutex);
+            state->cv.wait(lk, [&]() {
+                return state->stopped.load(std::memory_order_relaxed) || worker->stopped.load(std::memory_order_relaxed) || !state->tasks.empty();
+            });
+
+            if (state->stopped.load(std::memory_order_relaxed) || worker->stopped.load(std::memory_order_relaxed) || state->tasks.empty()) {
+                return;
+            }
+            task = std::move(state->tasks.front());
+            state->tasks.pop();
+        }
+        task();
+    }
+}
+} // namespace threading
@@ -0,0 +1,14 @@
+{
+  "default-registry": {
+    "kind": "git",
+    "baseline": "ac7af7424cbaf9057cb246b620f455303dccd6ed",
+    "repository": "https://github.com/microsoft/vcpkg"
+  },
+  "registries": [
+    {
+      "kind": "artifact",
+      "location": "https://github.com/microsoft/vcpkg-ce-catalog/archive/refs/heads/main.zip",
+      "name": "microsoft"
+    }
+  ]
+}
@@ -0,0 +1,11 @@
+{
+  "dependencies": [
+    "eigen3",
+    "nanoflann",
+    "libwebp",
+    {
+      "name": "libavif",
+      "features": ["aom"]
+    }
+  ]
+}