diff --git a/.clang-tidy b/.clang-tidy
index 1ad12ab..d243816 100644
--- a/.clang-tidy
+++ b/.clang-tidy
@@ -18,6 +18,7 @@ Checks: >
   -readability-implicit-bool-conversion,
   -readability-avoid-const-params-in-decls,
   -readability-braces-around-statements,
+  -readability-convert-member-functions-to-static,
   -readability-isolate-declaration,
   -readability-math-missing-parentheses,
   -cppcoreguidelines-avoid-magic-numbers,
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 07d6ab7..b5e61d3 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -20,6 +20,8 @@ add_library(genmetaballs_core
   genmetaballs/src/cuda/core/camera.cuh
   genmetaballs/src/cuda/core/utils.cu
   genmetaballs/src/cuda/core/utils.cuh
+  genmetaballs/src/cuda/core/fmb.cuh
+  genmetaballs/src/cuda/core/fmb.cu
   genmetaballs/src/cuda/core/geometry.cuh
   genmetaballs/src/cuda/core/geometry.cu
   genmetaballs/src/cuda/core/confidence.cuh
diff --git a/genmetaballs/src/cuda/bindings.cu b/genmetaballs/src/cuda/bindings.cu
index 3ae2652..ff3850b 100644
--- a/genmetaballs/src/cuda/bindings.cu
+++ b/genmetaballs/src/cuda/bindings.cu
@@ -2,11 +2,13 @@
 #include <nanobind/nanobind.h>
 #include <nanobind/ndarray.h>
 #include <nanobind/operators.h>
+#include <nanobind/stl/tuple.h>
 #include <nanobind/stl/vector.h>
 
 #include "core/blender.cuh"
 #include "core/camera.cuh"
 #include "core/confidence.cuh"
+#include "core/fmb.cuh"
 #include "core/geometry.cuh"
 #include "core/image.cuh"
 #include "core/utils.cuh"
@@ -22,6 +24,23 @@ void bind_image_view(nb::module_& m, const char* name);
 
 NB_MODULE(_genmetaballs_bindings, m) {
 
+    /*
+     * FMB module bindings
+     */
+
+    nb::module_ fmb = m.def_submodule("fmb", "Fuzzy meta ball data types");
+
+    nb::class_<FMB>(fmb, "FMB")
+        .def(nb::init<Pose, float, float, float>())
+        .def_prop_ro("pose", &FMB::get_pose)
+        .def_prop_ro("extent",
+                     [](const FMB& self) {
+                         auto extent = self.get_extent();
+                         return std::tuple{extent.x, extent.y, extent.z};
+                     })
+        .def("quadratic_form", &FMB::quadratic_form,
+             "Evaluate the associated quadratic form at the given vector", nb::arg("vec"));
+
     /*
      * Geometry module bindings
      */
diff --git a/genmetaballs/src/cuda/core/fmb.cu b/genmetaballs/src/cuda/core/fmb.cu
new file mode 100644
index 0000000..83d9cdc
--- /dev/null
+++ b/genmetaballs/src/cuda/core/fmb.cu
@@ -0,0 +1,33 @@
+#include "fmb.cuh"
+#include "geometry.cuh"
+#include "utils.cuh"
+
+CUDA_CALLABLE float FMB::quadratic_form(const Vec3D vec) const {
+    const auto shftd_vec = vec - pose_.get_tran();
+    const auto rot_shftd_vec = pose_.get_rot().apply(shftd_vec);
+    const auto scaled_rot_shftd_vec = Vec3D(
+        rot_shftd_vec.x / extent_.x, rot_shftd_vec.y / extent_.y, rot_shftd_vec.z / extent_.z);
+    return dot(rot_shftd_vec, scaled_rot_shftd_vec);
+}
+
+template <>
+__host__ FMBScene<MemoryLocation::HOST>::FMBScene(size_t size)
+    : fmbs_{new FMB[size]}, log_weights_{new float[size]}, size_{size} {}
+
+template <>
+__host__ FMBScene<MemoryLocation::DEVICE>::FMBScene(size_t size) : size_{size} {
+    CUDA_CHECK(cudaMalloc(&fmbs_, size * sizeof(FMB)));
+    CUDA_CHECK(cudaMalloc(&log_weights_, size * sizeof(float)));
+}
+
+template <>
+__host__ FMBScene<MemoryLocation::HOST>::~FMBScene() {
+    delete[] fmbs_;
+    delete[] log_weights_;
+}
+
+template <>
+__host__ FMBScene<MemoryLocation::DEVICE>::~FMBScene() {
+    CUDA_CHECK(cudaFree(fmbs_));
+    CUDA_CHECK(cudaFree(log_weights_));
+}
diff --git a/genmetaballs/src/cuda/core/fmb.cuh b/genmetaballs/src/cuda/core/fmb.cuh
index bf2f162..94242be 100644
--- a/genmetaballs/src/cuda/core/fmb.cuh
+++ b/genmetaballs/src/cuda/core/fmb.cuh
@@ -1,26 +1,116 @@
 #pragma once
 
+#include <cuda/std/span>
+#include <cuda/std/tuple>
+#include <stdexcept>
+
 #include "geometry.cuh"
+#include "utils.cuh"
+
+class FMB {
+private:
+    // In Gaussian terms:
+    // - mean: pose.tran
+    // - cov: pose.rot.mat().inv() * diag(extent) * pose.rot.mat()
+    Pose pose_;
+    float3 extent_;
+
+public:
+    FMB() : pose_{}, extent_{1.0f, 1.0f, 1.0f} {};
+
+    FMB(const Pose& pose, float x_extent, float y_extent, float z_extent) noexcept(false)
+        : pose_{pose} {
+        if (x_extent <= 0 || y_extent <= 0 || z_extent <= 0)
+            throw std::domain_error("a metaball cannot have negative extent");
+        extent_ = {x_extent, y_extent, z_extent};
+    }
+
+    CUDA_CALLABLE Pose get_pose() const {
+        return pose_;
+    }
+    CUDA_CALLABLE float3 get_extent() const {
+        return extent_;
+    }
 
-struct FMB {
-    Pose pose; // mean + orientation
-    float3 extent;
+    CUDA_CALLABLE float quadratic_form(const Vec3D) const;
 };
 
-template <template <typename> class containter_template>
-class FMBs {
+template <MemoryLocation location>
+class FMBScene {
 private:
-    containter_template<FMB> fmbs_;
-    containter_template<float> log_weights_;
+    FMB* fmbs_;
+    float* log_weights_;
+    size_t size_;
 
 public:
-    FMBs(uint32_t size) : fmbs_(size), log_weights_(size) {
-        // TODO: set all log_weights_ to 0
+    __host__ FMBScene(size_t size);
+
+    __host__ ~FMBScene();
+
+    CUDA_CALLABLE cuda::std::tuple<FMB&, float&> operator[](const uint32_t i) {
+        return cuda::std::tie(fmbs_[i], log_weights_[i]);
     }
-    CUDA_CALLABLE const containter_template<FMB>& get_all_fmbs() const {
-        return fmbs_;
+
+    CUDA_CALLABLE cuda::std::tuple<const FMB&, const float&> operator[](const uint32_t i) const {
+        return cuda::std::tie(fmbs_[i], log_weights_[i]);
+    }
+
+    class Iterator {
+    private:
+        FMB* fmb_ptr_;
+        float* log_weight_ptr_;
+
+    public:
+        CUDA_CALLABLE Iterator(FMB* const fmb_ptr, float* const log_weight_ptr)
+            : fmb_ptr_{fmb_ptr}, log_weight_ptr_{log_weight_ptr} {}
+        CUDA_CALLABLE cuda::std::tuple<FMB&, float&> operator*() {
+            return cuda::std::tie(*fmb_ptr_, *log_weight_ptr_);
+        }
+        CUDA_CALLABLE bool operator!=(const Iterator& other) const {
+            return fmb_ptr_ != other.fmb_ptr_ || log_weight_ptr_ != other.log_weight_ptr_;
+        }
+        CUDA_CALLABLE Iterator& operator++() {
+            fmb_ptr_++, log_weight_ptr_++;
+            return *this;
+        }
+    };
+
+    class ConstIterator {
+    private:
+        const FMB* fmb_ptr_;
+        const float* log_weight_ptr_;
+
+    public:
+        CUDA_CALLABLE ConstIterator(const FMB* const fmb_ptr, const float* const log_weight_ptr)
+            : fmb_ptr_{fmb_ptr}, log_weight_ptr_{log_weight_ptr} {}
+        CUDA_CALLABLE cuda::std::tuple<const FMB&, const float&> operator*() const {
+            return cuda::std::tie(*fmb_ptr_, *log_weight_ptr_);
+        }
+        CUDA_CALLABLE bool operator!=(const ConstIterator& other) const {
+            return fmb_ptr_ != other.fmb_ptr_ || log_weight_ptr_ != other.log_weight_ptr_;
+        }
+        CUDA_CALLABLE ConstIterator& operator++() {
+            fmb_ptr_++, log_weight_ptr_++;
+            return *this;
+        }
+    };
+
+    CUDA_CALLABLE Iterator begin() {
+        return Iterator(fmbs_, log_weights_);
+    }
+    CUDA_CALLABLE Iterator end() {
+        return Iterator(fmbs_ + size_, log_weights_ + size_);
+    }
+    CUDA_CALLABLE ConstIterator begin() const {
+        return ConstIterator(fmbs_, log_weights_);
+    }
+    CUDA_CALLABLE ConstIterator end() const {
+        return ConstIterator(fmbs_ + size_, log_weights_ + size_);
     }
     CUDA_CALLABLE const FMB& get_fmb(uint32_t idx) const {
         return fmbs_[idx];
     }
+    CUDA_CALLABLE size_t size() const {
+        return size_;
+    }
 };
diff --git a/genmetaballs/src/cuda/core/forward.cu b/genmetaballs/src/cuda/core/forward.cu
index 309729c..d75ce21 100644
--- a/genmetaballs/src/cuda/core/forward.cu
+++ b/genmetaballs/src/cuda/core/forward.cu
@@ -34,8 +34,8 @@ __global__ render_kernel(const Getter fmb_getter, const Blender blender,
         float w0 = 0.0f, tf = 0.0f, sumexpd = 0.0f;
         for (const auto& fmb : fmb_getter->get_metaballs(ray)) {
             const auto& [t, d] = Intersector::intersect(fmb, ray);
-            w = blender->blend(t, d);
-            sumexpd += exp(d);
+            w = blender->blend(t, d, fmb, ray);
+            sumexpd += exp(d); // numerically unstable. use logsumexp
             tf += t;
             w0 += w;
         }
diff --git a/genmetaballs/src/cuda/core/getter.cuh b/genmetaballs/src/cuda/core/getter.cuh
index 86590f1..4f17090 100644
--- a/genmetaballs/src/cuda/core/getter.cuh
+++ b/genmetaballs/src/cuda/core/getter.cuh
@@ -9,16 +9,15 @@
 #include "utils.cuh"
 
 // This is the dummy version of getter, where all FMBs are relevant to any ray
-template <template <typename> class containter_template>
+template <MemoryLocation location>
 struct AllGetter {
-    const FMBs<containter_template>& fmbs;
-    Pose extr; // Current assumption: rays are in camera frame
+    FMBScene<location>& scene;
+    Pose& extr; // Current assumption: rays are in camera frame
 
-    CUDA_CALLABLE AllGetter(const FMBs<containter_template>& fmbs, const Pose& extr)
-        : fmbs(fmbs), extr(extr) {}
+    CUDA_CALLABLE AllGetter(FMBScene<location>& scene, Pose& extr) : scene(scene), extr(extr) {}
 
     // It does not bother using the ray, because it simply returns all FMBs
-    CUDA_CALLABLE const containter_template<FMB>& get_metaballs(const Ray& ray) const {
-        return fmbs.get_all_fmbs();
+    CUDA_CALLABLE FMBScene<location>& get_metaballs(const Ray& ray) const {
+        return scene;
     }
 };
diff --git a/genmetaballs/src/genmetaballs/__init__.py b/genmetaballs/src/genmetaballs/__init__.py
index 9c40d8c..e69de29 100644
--- a/genmetaballs/src/genmetaballs/__init__.py
+++ b/genmetaballs/src/genmetaballs/__init__.py
@@ -1,3 +0,0 @@
-from .gpu_add import gpu_add
-
-__all__ = ["gpu_add"]
diff --git a/genmetaballs/src/genmetaballs/core/__init__.py b/genmetaballs/src/genmetaballs/core/__init__.py
index ac2a6fd..eae5c9c 100644
--- a/genmetaballs/src/genmetaballs/core/__init__.py
+++ b/genmetaballs/src/genmetaballs/core/__init__.py
@@ -1,6 +1,6 @@
 from typing import Literal
 
-from genmetaballs._genmetaballs_bindings import geometry
+from genmetaballs._genmetaballs_bindings import fmb, geometry
 from genmetaballs._genmetaballs_bindings.blender import (
     FourParameterBlender,
     ThreeParameterBlender,
@@ -51,6 +51,7 @@ def make_image(height: int, width: int, device: DeviceType) -> CPUImage | GPUIma
     "array2d_float",
     "ZeroParameterConfidence",
     "TwoParameterConfidence",
+    "fmb",
     "geometry",
     "Camera",
     "Intrinsics",
diff --git a/genmetaballs/src/genmetaballs/gpu_add.py b/genmetaballs/src/genmetaballs/gpu_add.py
deleted file mode 100644
index 214b75f..0000000
--- a/genmetaballs/src/genmetaballs/gpu_add.py
+++ /dev/null
@@ -1,5 +0,0 @@
-from . import _genmetaballs_bindings as _gmbb
-
-
-def gpu_add(a: list[float], b: list[float]) -> list[float]:
-    return _gmbb.gpu_add(a, b)
diff --git a/tests/cpp_tests/test_fmb.cu b/tests/cpp_tests/test_fmb.cu
new file mode 100644
index 0000000..73d6602
--- /dev/null
+++ b/tests/cpp_tests/test_fmb.cu
@@ -0,0 +1,31 @@
+#include <cstdint>
+#include <cuda_runtime.h>
+#include <gtest/gtest.h>
+#include <thrust/device_vector.h>
+#include <thrust/host_vector.h>
+
+#include "core/fmb.cuh"
+#include "core/utils.cuh"
+
+__global__ void dummy_kernel(FMBScene<MemoryLocation::DEVICE>& scene, int* num_fmbs) {
+
+    int _num_fmbs = 0;
+
+    for (auto [fmb, w] : scene) {
+        _num_fmbs += 1;
+    }
+
+    *num_fmbs = _num_fmbs;
+}
+
+TEST(FMBTests, KernelRangeBasedForLoopSmokeTest) {
+
+    FMBScene<MemoryLocation::DEVICE> dummy_scene(10);
+    thrust::device_vector<int> device_res(1);
+
+    dummy_kernel<<<1, 1>>>(dummy_scene, thrust::raw_pointer_cast(device_res.data()));
+
+    thrust::host_vector<int> host_res = device_res;
+
+    EXPECT_EQ(host_res[0], 10);
+}
diff --git a/tests/cpp_tests/test_getter.cu b/tests/cpp_tests/test_getter.cu
index 096a499..f245cdd 100644
--- a/tests/cpp_tests/test_getter.cu
+++ b/tests/cpp_tests/test_getter.cu
@@ -11,55 +11,12 @@
 #include "core/getter.cuh"
 #include "core/utils.cuh"
 
-// a whacky helper to get the container type from the template parameter
-template <typename Container>
-struct template_container_of;
-
-template <template <typename...> class Template, typename... Args>
-struct template_container_of<Template<Args...>> {
-    template <typename T>
-    using type = Template<T>;
-};
-
-template <typename InstantiatedContainer>
-struct GetterTestTypes;
-
-template <template <typename...> class Template, typename... Args>
-struct GetterTestTypes<Template<Args...>> {
-    template <typename T>
-    using ContainerTemplate = typename template_container_of<Template<Args...>>::template type<T>;
-
-    using FMBsType = FMBs<ContainerTemplate>;
-    using AllGetterType = AllGetter<ContainerTemplate>;
-};
-
-// Template fixture for all container types
-template <typename Container>
-class AllGetterTestFixture : public ::testing::Test {};
-
-using ContainerTypes = ::testing::Types<std::vector<FMB>, thrust::device_vector<FMB>>;
-TYPED_TEST_SUITE(AllGetterTestFixture, ContainerTypes);
-
-// CUDA kernel that constructs AllGetter and calls get_metaballs
-template <typename AllGetterType, typename FMBsType>
-__global__ void test_get_metaballs_kernel_device(const FMBsType* fmbs, const Pose* extr,
-                                                 const Ray* rays, int num_rays, int* out_sizes) {
-    int idx = threadIdx.x + blockIdx.x * blockDim.x;
-    AllGetterType getter(*fmbs, *extr);
-    const auto& fmbs_returned = getter.get_metaballs(rays[idx]);
-    out_sizes[idx] = static_cast<int>(fmbs_returned.size());
-}
-
-TYPED_TEST(AllGetterTestFixture, ReturnsAllFMBsForAnyRay) {
-    // Extract types from TypeParam (std::vector<FMB> or thrust::device_vector<FMB>)
-    using FMBsType = typename GetterTestTypes<TypeParam>::FMBsType;
-    using AllGetterType = typename GetterTestTypes<TypeParam>::AllGetterType;
-
+TEST(AllGetterTest, AllGetterHostTest) {
     constexpr uint32_t num_fmbs = 40;
-    FMBsType fmbs(num_fmbs);
+    FMBScene<MemoryLocation::HOST> scene(num_fmbs);
 
-    Pose extr = Pose();
-    AllGetterType getter(fmbs, extr);
+    Pose extr;
+    AllGetter<MemoryLocation::HOST> getter(scene, extr);
 
     // Create test rays
     std::vector<Ray> rays = {
@@ -76,64 +33,65 @@ TYPED_TEST(AllGetterTestFixture, ReturnsAllFMBsForAnyRay) {
     };
 
     // Get reference to all FMBs from the original FMBs object
-    const auto& all_fmbs_ref = fmbs.get_all_fmbs();
 
     // Test on host - AllGetter should return the same container for all rays
     for (const auto& ray : rays) {
         const auto& fmbs_returned = getter.get_metaballs(ray);
 
         // Verify that we get the same container reference (both are host objects)
-        EXPECT_EQ(&fmbs_returned, &all_fmbs_ref)
+        EXPECT_EQ(&fmbs_returned, &scene)
             << "AllGetter should return the same FMBs container for all rays";
 
-        // For thrust::device_vector, also verify device pointers match
-        if constexpr (std::is_same_v<TypeParam, thrust::device_vector<FMB>>) {
-            EXPECT_EQ(thrust::raw_pointer_cast(fmbs_returned.data()),
-                      thrust::raw_pointer_cast(all_fmbs_ref.data()))
-                << "Device pointers should match for thrust::device_vector";
-        }
-
         // Verify container sizes match
-        EXPECT_EQ(fmbs_returned.size(), all_fmbs_ref.size())
+        EXPECT_EQ(fmbs_returned.size(), scene.size())
             << "Returned FMBs container size must match all_fmbs size";
     }
+}
+
+__global__ void test_get_metaballs_kernel_device(const AllGetter<MemoryLocation::DEVICE> fmb_getter,
+                                                 const Ray* rays, int num_rays, int* out_sizes) {
+    int idx = threadIdx.x + blockIdx.x * blockDim.x;
+    const auto& fmbs_returned = fmb_getter.get_metaballs(rays[idx]);
+    out_sizes[idx] = static_cast<int>(fmbs_returned.size());
+}
+
+TEST(AllGetterTest, AllGetterDeviceTest) {
+    constexpr uint32_t num_fmbs = 40;
+    FMBScene<MemoryLocation::DEVICE> device_scene(num_fmbs);
+    Pose extr;
+
+    AllGetter<MemoryLocation::DEVICE> getter(device_scene, extr);
+
+    // Create test rays
+    std::vector<Ray> rays = {
+        Ray{Vec3D{0.0f, 0.0f, 0.0f}, Vec3D{1.0f, 0.0f, 0.0f}},
+        Ray{Vec3D{1.0f, 1.0f, 1.0f}, Vec3D{0.0f, 1.0f, 0.0f}},
+        Ray{Vec3D{-1.0f, -1.0f, -1.0f}, Vec3D{0.0f, 0.0f, 1.0f}},
+        Ray{Vec3D{2.5f, -3.1f, 0.2f}, Vec3D{-0.5f, 0.6f, 0.0f}},
+        Ray{Vec3D{4.4f, 0.0f, -0.9f}, Vec3D{0.3f, -0.2f, 1.0f}},
+        Ray{Vec3D{5.0f, 2.2f, 1.1f}, Vec3D{-1.0f, 2.0f, 0.2f}},
+        Ray{Vec3D{0.0f, 7.0f, 6.0f}, Vec3D{0.0f, -1.0f, -1.0f}},
+        Ray{Vec3D{-2.0f, 0.0f, 0.0f}, Vec3D{0.2f, 1.1f, 0.7f}},
+        Ray{Vec3D{9.1f, -0.3f, 2.7f}, Vec3D{-0.3f, 0.1f, 0.0f}},
+        Ray{Vec3D{1.2f, 8.8f, -4.5f}, Vec3D{1.0f, 0.0f, 1.0f}},
+    };
 
     // Test on GPU for device containers
-    // Construct AllGetter on device and call get_metaballs
-    if constexpr (std::is_same_v<TypeParam, thrust::device_vector<FMB>>) {
-        Ray* d_rays = nullptr;
-        FMBsType* d_fmbs = nullptr;
-        Pose* d_extr = nullptr;
-        int* d_sizes = nullptr;
-        int num_rays = static_cast<int>(rays.size());
-
-        CUDA_CHECK(cudaMalloc(&d_rays, num_rays * sizeof(Ray)));
-        CUDA_CHECK(cudaMalloc(&d_fmbs, sizeof(FMBsType)));
-        CUDA_CHECK(cudaMalloc(&d_extr, sizeof(Pose)));
-        CUDA_CHECK(cudaMalloc(&d_sizes, num_rays * sizeof(int)));
-        CUDA_CHECK(cudaMemcpy(d_rays, rays.data(), num_rays * sizeof(Ray), cudaMemcpyHostToDevice));
-        CUDA_CHECK(cudaMemcpy(d_fmbs, &fmbs, sizeof(FMBsType), cudaMemcpyHostToDevice));
-        CUDA_CHECK(cudaMemcpy(d_extr, &extr, sizeof(Pose), cudaMemcpyHostToDevice));
-
-        // Launch kernel that constructs getter and calls get_metaballs on the device
-        test_get_metaballs_kernel_device<AllGetterType, FMBsType>
-            <<<1, num_rays>>>(d_fmbs, d_extr, d_rays, num_rays, d_sizes);
-        CUDA_CHECK(cudaGetLastError());
-        CUDA_CHECK(cudaDeviceSynchronize());
-
-        std::vector<int> host_sizes(num_rays);
-        CUDA_CHECK(
-            cudaMemcpy(host_sizes.data(), d_sizes, num_rays * sizeof(int), cudaMemcpyDeviceToHost));
-
-        // Verify that get_metaballs returns the correct size for all rays
-        for (int i = 0; i < num_rays; ++i) {
-            EXPECT_EQ(static_cast<size_t>(host_sizes[i]), all_fmbs_ref.size())
-                << "Device get_metaballs returned correct size for ray " << i;
-        }
-
-        CUDA_CHECK(cudaFree(d_rays));
-        CUDA_CHECK(cudaFree(d_fmbs));
-        CUDA_CHECK(cudaFree(d_extr));
-        CUDA_CHECK(cudaFree(d_sizes));
+    int num_rays = static_cast<int>(rays.size());
+    thrust::device_vector<int> device_sizes(num_rays);
+
+    // Launch kernel that constructs getter and calls get_metaballs on the device
+    test_get_metaballs_kernel_device<<<1, num_rays>>>(
+        getter, thrust::raw_pointer_cast(rays.data()), num_rays,
+        thrust::raw_pointer_cast(device_sizes.data()));
+    CUDA_CHECK(cudaGetLastError());
+    CUDA_CHECK(cudaDeviceSynchronize());
+
+    thrust::host_vector<int> host_sizes = device_sizes;
+
+    // Verify that get_metaballs returns the correct size for all rays
+    for (int i = 0; i < num_rays; ++i) {
+        EXPECT_EQ(static_cast<size_t>(host_sizes[i]), device_scene.size())
+            << "Device get_metaballs returned correct size for ray " << i;
     }
 }
diff --git a/tests/python_tests/test_fmb.py b/tests/python_tests/test_fmb.py
new file mode 100644
index 0000000..211f99a
--- /dev/null
+++ b/tests/python_tests/test_fmb.py
@@ -0,0 +1,27 @@
+import numpy as np
+import pytest
+from scipy.spatial.distance import mahalanobis
+from scipy.spatial.transform import Rotation as Rot
+
+from genmetaballs.core import fmb, geometry
+
+FMB = fmb.FMB
+Pose, Vec3D, Rotation = geometry.Pose, geometry.Vec3D, geometry.Rotation
+
+
+@pytest.fixture
+def rng() -> np.random.Generator:
+    return np.random.default_rng(0)
+
+
+def test_fmb_quadratic_form(rng):
+    for _ in range(100):
+        quat = rng.uniform(size=4)
+        tran, extent, vec = rng.uniform(size=(3, 3))
+        pose = Pose.from_components(Rotation.from_quat(*quat), Vec3D(*tran))
+        scipy_rot_mat = Rot.from_quat(quat).as_matrix()
+        cov = scipy_rot_mat.T @ np.diag(extent) @ scipy_rot_mat
+        assert np.isclose(
+            FMB(pose, *extent).quadratic_form(Vec3D(*vec)),
+            mahalanobis(vec, tran, np.linalg.inv(cov)) ** 2,
+        )