added files

added aabb line trace

.claude/settings.local.json

@@ -1,7 +0,0 @@
-{
-  "permissions": {
-    "allow": [
-      "Bash(ls:*)"
-    ]
-  }
-}

.idea/editor.xml

@@ -17,7 +17,7 @@
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppBoostFormatTooManyArgs/@EntryIndexedValue" value="WARNING" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppCStyleCast/@EntryIndexedValue" value="SUGGESTION" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppCVQualifierCanNotBeAppliedToReference/@EntryIndexedValue" value="WARNING" type="string" />
-    <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppClassCanBeFinal/@EntryIndexedValue" value="WARNING" type="string" />
+    <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppClassCanBeFinal/@EntryIndexedValue" value="DO_NOT_SHOW" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppClassIsIncomplete/@EntryIndexedValue" value="WARNING" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppClassNeedsConstructorBecauseOfUninitializedMember/@EntryIndexedValue" value="WARNING" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppClassNeverUsed/@EntryIndexedValue" value="WARNING" type="string" />
@@ -103,14 +103,14 @@
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppImplicitDefaultConstructorNotAvailable/@EntryIndexedValue" value="WARNING" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppIncompatiblePointerConversion/@EntryIndexedValue" value="WARNING" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppIncompleteSwitchStatement/@EntryIndexedValue" value="WARNING" type="string" />
-    <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppInconsistentNaming/@EntryIndexedValue" value="HINT" type="string" />
+    <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppInconsistentNaming/@EntryIndexedValue" value="DO_NOT_SHOW" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppIntegralToPointerConversion/@EntryIndexedValue" value="WARNING" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppInvalidLineContinuation/@EntryIndexedValue" value="WARNING" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppJoinDeclarationAndAssignment/@EntryIndexedValue" value="SUGGESTION" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppLambdaCaptureNeverUsed/@EntryIndexedValue" value="WARNING" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppLocalVariableMayBeConst/@EntryIndexedValue" value="HINT" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppLocalVariableMightNotBeInitialized/@EntryIndexedValue" value="WARNING" type="string" />
-    <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppLocalVariableWithNonTrivialDtorIsNeverUsed/@EntryIndexedValue" value="WARNING" type="string" />
+    <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppLocalVariableWithNonTrivialDtorIsNeverUsed/@EntryIndexedValue" value="DO_NOT_SHOW" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppLongFloat/@EntryIndexedValue" value="WARNING" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppMemberFunctionMayBeConst/@EntryIndexedValue" value="SUGGESTION" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppMemberFunctionMayBeStatic/@EntryIndexedValue" value="SUGGESTION" type="string" />
@@ -202,7 +202,7 @@
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppStaticDataMemberInUnnamedStruct/@EntryIndexedValue" value="WARNING" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppStaticSpecifierOnAnonymousNamespaceMember/@EntryIndexedValue" value="SUGGESTION" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppStringLiteralToCharPointerConversion/@EntryIndexedValue" value="WARNING" type="string" />
-    <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppTabsAreDisallowed/@EntryIndexedValue" value="WARNING" type="string" />
+    <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppTabsAreDisallowed/@EntryIndexedValue" value="DO_NOT_SHOW" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppTemplateArgumentsCanBeDeduced/@EntryIndexedValue" value="HINT" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppTemplateParameterNeverUsed/@EntryIndexedValue" value="HINT" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppTemplateParameterShadowing/@EntryIndexedValue" value="WARNING" type="string" />
@@ -216,7 +216,7 @@
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppUnmatchedPragmaEndRegionDirective/@EntryIndexedValue" value="WARNING" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppUnmatchedPragmaRegionDirective/@EntryIndexedValue" value="WARNING" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppUnnamedNamespaceInHeaderFile/@EntryIndexedValue" value="WARNING" type="string" />
-    <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppUnnecessaryWhitespace/@EntryIndexedValue" value="WARNING" type="string" />
+    <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppUnnecessaryWhitespace/@EntryIndexedValue" value="DO_NOT_SHOW" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppUnsignedZeroComparison/@EntryIndexedValue" value="WARNING" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppUnusedIncludeDirective/@EntryIndexedValue" value="WARNING" type="string" />
     <option name="/Default/CodeInspection/Highlighting/InspectionSeverities/=CppUseAlgorithmWithCount/@EntryIndexedValue" value="SUGGESTION" type="string" />

CONTRIBUTING.md

@@ -1,32 +1,36 @@
-## 🤝 Contributing to OMath or other Orange's Projects
+# Contributing
 
-### ❕ Prerequisites
+## Prerequisites
 
-- A working up-to-date OMath installation
-- C++ knowledge
-- Git knowledge
-- Ability to ask for help (Feel free to create empty pull-request or PM a maintainer
-  in [Telegram](https://t.me/orange_cpp))
+- C++ compiler with C++23 support (Clang 18+, GCC 14+, MSVC 19.38+)
+- CMake 3.25+
+- Git
+- Familiarity with the codebase (see `INSTALL.md` for setup)
 
-### ⏬ Setting up OMath
+For questions, create a draft PR or reach out via [Telegram](https://t.me/orange_cpp).
 
-Please read INSTALL.md file in repository
+## Workflow
 
-### 🔀 Pull requests and Branches
+1. [Fork](https://docs.github.com/en/pull-requests/collaborating-with-pull-requests/working-with-forks/fork-a-repo) the repository.
+2. Create a feature branch from `main`.
+3. Make your changes, ensuring tests pass.
+4. Open a [pull request](https://docs.github.com/en/pull-requests/collaborating-with-pull-requests/proposing-changes-to-your-work-with-pull-requests/creating-a-pull-request-from-a-fork) against `main`.
 
-In order to send code back to the official OMath repository, you must first create a copy of OMath on your github
-account ([fork](https://help.github.com/articles/creating-a-pull-request-from-a-fork/)) and
-then [create a pull request](https://help.github.com/articles/creating-a-pull-request-from-a-fork/) back to OMath.
+## Code Style
 
-OMath development is performed on multiple branches. Changes are then pull requested into master. By default, changes
-merged into master will not roll out to stable build users unless the `stable` tag is updated.
+Follow the project `.clang-format`. Run `clang-format` before committing.
 
-### 📜 Code-Style
+## Building
 
-The orange code-style can be found in `.clang-format`.
+Use one of the CMake presets defined in `CMakePresets.json`:
 
-### 📦 Building
+```bash
+cmake --preset <preset-name> -DOMATH_BUILD_TESTS=ON
+cmake --build --preset <preset-name>
+```
 
-OMath has already created the  `cmake-build` and `out` directories where cmake/bin files are located. By default, you
-can build OMath by running `cmake --build cmake-build/build/windows-release --target omath -j 6` in the source
-directory.
+Run `cmake --list-presets` to see available configurations.
+
+## Tests
+
+All new functionality must include unit tests. Run the test binary after building to verify nothing is broken.

examples/example_barycentric/example_barycentric.cpp

@@ -71,18 +71,18 @@ void drawChar(char c, float x, float y, float scale, const Color& color, std::ve
         lines.push_back(x + x1 * w);
         lines.push_back(y + y1 * h);
         lines.push_back(0.0f);
-        lines.push_back(color.x);
-        lines.push_back(color.y);
-        lines.push_back(color.z);
+        lines.push_back(color.value().x);
+        lines.push_back(color.value().y);
+        lines.push_back(color.value().z);
         lines.push_back(1.0f); // size
         lines.push_back(1.0f); // isLine
 
         lines.push_back(x + x2 * w);
         lines.push_back(y + y2 * h);
         lines.push_back(0.0f);
-        lines.push_back(color.x);
-        lines.push_back(color.y);
-        lines.push_back(color.z);
+        lines.push_back(color.value().x);
+        lines.push_back(color.value().y);
+        lines.push_back(color.value().z);
         lines.push_back(1.0f); // size
         lines.push_back(1.0f); // isLine
     };

include/omath/3d_primitives/aabb.hpp

@@ -0,0 +1,28 @@
+//
+// Created by Vladislav on 24.03.2026.
+//
+
+#pragma once
+#include "omath/linear_algebra/vector3.hpp"
+
+namespace omath::primitives
+{
+    template<class Type>
+    struct Aabb final
+    {
+        Vector3<Type> min;
+        Vector3<Type> max;
+
+        [[nodiscard]]
+        constexpr Vector3<Type> center() const noexcept
+        {
+            return (min + max) / static_cast<Type>(2);
+        }
+
+        [[nodiscard]]
+        constexpr Vector3<Type> extents() const noexcept
+        {
+            return (max - min) / static_cast<Type>(2);
+        }
+    };
+} // namespace omath::primitives

include/omath/collision/bvh_tree.hpp

@@ -0,0 +1,412 @@
+//
+// Created by Orange on 04/08/2026.
+//
+
+#pragma once
+#include "omath/3d_primitives/aabb.hpp"
+#include "omath/collision/line_tracer.hpp"
+#include <algorithm>
+#include <cstdint>
+#include <numeric>
+#include <span>
+#include <vector>
+
+namespace omath::collision
+{
+    template<class Type = float>
+    class BvhTree final
+    {
+    public:
+        using AabbType = primitives::Aabb<Type>;
+
+        struct HitResult
+        {
+            std::size_t object_index;
+            Type distance_sqr;
+        };
+
+        BvhTree() = default;
+
+        explicit BvhTree(std::span<const AabbType> aabbs)
+            : m_aabbs(aabbs.begin(), aabbs.end())
+        {
+            if (aabbs.empty())
+                return;
+
+            m_indices.resize(aabbs.size());
+            std::iota(m_indices.begin(), m_indices.end(), std::size_t{0});
+
+            m_nodes.reserve(aabbs.size() * 2);
+
+            build(m_aabbs, 0, aabbs.size());
+        }
+
+        [[nodiscard]]
+        std::vector<std::size_t> query_overlaps(const AabbType& query_aabb) const
+        {
+            std::vector<std::size_t> results;
+
+            if (m_nodes.empty())
+                return results;
+
+            query_overlaps_impl(0, query_aabb, results);
+            return results;
+        }
+
+        template<class RayType = Ray<>>
+        [[nodiscard]]
+        std::vector<HitResult> query_ray(const RayType& ray) const
+        {
+            std::vector<HitResult> results;
+
+            if (m_nodes.empty())
+                return results;
+
+            query_ray_impl(0, ray, results);
+
+            std::ranges::sort(results, [](const HitResult& a, const HitResult& b)
+                              { return a.distance_sqr < b.distance_sqr; });
+            return results;
+        }
+
+        [[nodiscard]]
+        std::size_t node_count() const noexcept
+        {
+            return m_nodes.size();
+        }
+
+        [[nodiscard]]
+        bool empty() const noexcept
+        {
+            return m_nodes.empty();
+        }
+
+    private:
+        static constexpr std::size_t k_sah_bucket_count = 12;
+        static constexpr std::size_t k_leaf_threshold = 1;
+        static constexpr std::size_t k_null_index = std::numeric_limits<std::size_t>::max();
+
+        struct Node
+        {
+            AabbType bounds;
+            std::size_t left = k_null_index;
+            std::size_t right = k_null_index;
+
+            // For leaf nodes: index range into m_indices
+            std::size_t first_index = 0;
+            std::size_t index_count = 0;
+
+            [[nodiscard]]
+            bool is_leaf() const noexcept
+            {
+                return left == k_null_index;
+            }
+        };
+
+        struct SahBucket
+        {
+            AabbType bounds = {
+                {std::numeric_limits<Type>::max(), std::numeric_limits<Type>::max(),
+                 std::numeric_limits<Type>::max()},
+                {std::numeric_limits<Type>::lowest(), std::numeric_limits<Type>::lowest(),
+                 std::numeric_limits<Type>::lowest()}
+            };
+            std::size_t count = 0;
+        };
+
+        [[nodiscard]]
+        static constexpr Type surface_area(const AabbType& aabb) noexcept
+        {
+            const auto d = aabb.max - aabb.min;
+            return static_cast<Type>(2) * (d.x * d.y + d.y * d.z + d.z * d.x);
+        }
+
+        [[nodiscard]]
+        static constexpr AabbType merge(const AabbType& a, const AabbType& b) noexcept
+        {
+            return {
+                {std::min(a.min.x, b.min.x), std::min(a.min.y, b.min.y), std::min(a.min.z, b.min.z)},
+                {std::max(a.max.x, b.max.x), std::max(a.max.y, b.max.y), std::max(a.max.z, b.max.z)}
+            };
+        }
+
+        [[nodiscard]]
+        static constexpr bool overlaps(const AabbType& a, const AabbType& b) noexcept
+        {
+            return a.min.x <= b.max.x && a.max.x >= b.min.x
+                && a.min.y <= b.max.y && a.max.y >= b.min.y
+                && a.min.z <= b.max.z && a.max.z >= b.min.z;
+        }
+
+        std::size_t build(std::span<const AabbType> aabbs, std::size_t begin, std::size_t end)
+        {
+            const auto node_idx = m_nodes.size();
+            m_nodes.emplace_back();
+
+            auto& node = m_nodes[node_idx];
+            node.bounds = compute_bounds(aabbs, begin, end);
+
+            const auto count = end - begin;
+
+            if (count <= k_leaf_threshold)
+            {
+                node.first_index = begin;
+                node.index_count = count;
+                return node_idx;
+            }
+
+            // Find best SAH split
+            const auto centroid_bounds = compute_centroid_bounds(aabbs, begin, end);
+            const auto split = find_best_split(aabbs, begin, end, node.bounds, centroid_bounds);
+
+            // If SAH says don't split, make a leaf
+            if (!split.has_value())
+            {
+                node.first_index = begin;
+                node.index_count = count;
+                return node_idx;
+            }
+
+            const auto [axis, split_pos] = split.value();
+
+            // Partition indices around the split
+            const auto mid = partition_indices(aabbs, begin, end, axis, split_pos);
+
+            // Degenerate partition fallback: split in the middle
+            const auto actual_mid = (mid == begin || mid == end) ? begin + count / 2 : mid;
+
+            // Build children — careful: m_nodes may reallocate, so don't hold references across build calls
+            const auto left_idx = build(aabbs, begin, actual_mid);
+            const auto right_idx = build(aabbs, actual_mid, end);
+
+            m_nodes[node_idx].left = left_idx;
+            m_nodes[node_idx].right = right_idx;
+            m_nodes[node_idx].index_count = 0;
+
+            return node_idx;
+        }
+
+        [[nodiscard]]
+        AabbType compute_bounds(std::span<const AabbType> aabbs, std::size_t begin, std::size_t end) const
+        {
+            AabbType bounds = {
+                {std::numeric_limits<Type>::max(), std::numeric_limits<Type>::max(),
+                 std::numeric_limits<Type>::max()},
+                {std::numeric_limits<Type>::lowest(), std::numeric_limits<Type>::lowest(),
+                 std::numeric_limits<Type>::lowest()}
+            };
+
+            for (auto i = begin; i < end; ++i)
+                bounds = merge(bounds, aabbs[m_indices[i]]);
+
+            return bounds;
+        }
+
+        [[nodiscard]]
+        AabbType compute_centroid_bounds(std::span<const AabbType> aabbs, std::size_t begin, std::size_t end) const
+        {
+            AabbType bounds = {
+                {std::numeric_limits<Type>::max(), std::numeric_limits<Type>::max(),
+                 std::numeric_limits<Type>::max()},
+                {std::numeric_limits<Type>::lowest(), std::numeric_limits<Type>::lowest(),
+                 std::numeric_limits<Type>::lowest()}
+            };
+
+            for (auto i = begin; i < end; ++i)
+            {
+                const auto c = aabbs[m_indices[i]].center();
+                bounds.min.x = std::min(bounds.min.x, c.x);
+                bounds.min.y = std::min(bounds.min.y, c.y);
+                bounds.min.z = std::min(bounds.min.z, c.z);
+                bounds.max.x = std::max(bounds.max.x, c.x);
+                bounds.max.y = std::max(bounds.max.y, c.y);
+                bounds.max.z = std::max(bounds.max.z, c.z);
+            }
+
+            return bounds;
+        }
+
+        struct SplitResult
+        {
+            int axis;
+            Type position;
+        };
+
+        [[nodiscard]]
+        std::optional<SplitResult> find_best_split(std::span<const AabbType> aabbs, std::size_t begin,
+                                                   std::size_t end, const AabbType& node_bounds,
+                                                   const AabbType& centroid_bounds) const
+        {
+            const auto count = end - begin;
+            const auto leaf_cost = static_cast<Type>(count);
+            auto best_cost = leaf_cost;
+            std::optional<SplitResult> best_split;
+
+            for (int axis = 0; axis < 3; ++axis)
+            {
+                const auto axis_min = get_component(centroid_bounds.min, axis);
+                const auto axis_max = get_component(centroid_bounds.max, axis);
+
+                if (axis_max - axis_min < std::numeric_limits<Type>::epsilon())
+                    continue;
+
+                SahBucket buckets[k_sah_bucket_count] = {};
+
+                const auto inv_extent = static_cast<Type>(k_sah_bucket_count) / (axis_max - axis_min);
+
+                // Fill buckets
+                for (auto i = begin; i < end; ++i)
+                {
+                    const auto centroid = get_component(aabbs[m_indices[i]].center(), axis);
+                    auto bucket_idx = static_cast<std::size_t>((centroid - axis_min) * inv_extent);
+                    bucket_idx = std::min(bucket_idx, k_sah_bucket_count - 1);
+
+                    buckets[bucket_idx].count++;
+                    if (buckets[bucket_idx].count == 1)
+                        buckets[bucket_idx].bounds = aabbs[m_indices[i]];
+                    else
+                        buckets[bucket_idx].bounds = merge(buckets[bucket_idx].bounds, aabbs[m_indices[i]]);
+                }
+
+                // Evaluate split costs using prefix/suffix sweeps
+                AabbType prefix_bounds[k_sah_bucket_count - 1];
+                std::size_t prefix_count[k_sah_bucket_count - 1];
+
+                prefix_bounds[0] = buckets[0].bounds;
+                prefix_count[0] = buckets[0].count;
+                for (std::size_t i = 1; i < k_sah_bucket_count - 1; ++i)
+                {
+                    prefix_bounds[i] = (buckets[i].count > 0)
+                                           ? merge(prefix_bounds[i - 1], buckets[i].bounds)
+                                           : prefix_bounds[i - 1];
+                    prefix_count[i] = prefix_count[i - 1] + buckets[i].count;
+                }
+
+                AabbType suffix_bounds = buckets[k_sah_bucket_count - 1].bounds;
+                std::size_t suffix_count = buckets[k_sah_bucket_count - 1].count;
+
+                const auto parent_area = surface_area(node_bounds);
+                const auto inv_parent_area = static_cast<Type>(1) / parent_area;
+
+                for (auto i = static_cast<int>(k_sah_bucket_count) - 2; i >= 0; --i)
+                {
+                    const auto left_count = prefix_count[i];
+                    const auto right_count = suffix_count;
+
+                    if (left_count == 0 || right_count == 0)
+                    {
+                        if (i > 0)
+                        {
+                            suffix_bounds = (buckets[i].count > 0)
+                                                ? merge(suffix_bounds, buckets[i].bounds)
+                                                : suffix_bounds;
+                            suffix_count += buckets[i].count;
+                        }
+                        continue;
+                    }
+
+                    const auto cost = static_cast<Type>(1)
+                                      + (static_cast<Type>(left_count) * surface_area(prefix_bounds[i])
+                                         + static_cast<Type>(right_count) * surface_area(suffix_bounds))
+                                            * inv_parent_area;
+
+                    if (cost < best_cost)
+                    {
+                        best_cost = cost;
+                        best_split = SplitResult{
+                            axis,
+                            axis_min + static_cast<Type>(i + 1) * (axis_max - axis_min)
+                                           / static_cast<Type>(k_sah_bucket_count)
+                        };
+                    }
+
+                    suffix_bounds = (buckets[i].count > 0)
+                                        ? merge(suffix_bounds, buckets[i].bounds)
+                                        : suffix_bounds;
+                    suffix_count += buckets[i].count;
+                }
+            }
+
+            return best_split;
+        }
+
+        std::size_t partition_indices(std::span<const AabbType> aabbs, std::size_t begin, std::size_t end,
+                                      int axis, Type split_pos)
+        {
+            auto it = std::partition(m_indices.begin() + static_cast<std::ptrdiff_t>(begin),
+                                     m_indices.begin() + static_cast<std::ptrdiff_t>(end),
+                                     [&](std::size_t idx)
+                                     { return get_component(aabbs[idx].center(), axis) < split_pos; });
+
+            return static_cast<std::size_t>(std::distance(m_indices.begin(), it));
+        }
+
+        [[nodiscard]]
+        static constexpr Type get_component(const Vector3<Type>& v, int axis) noexcept
+        {
+            switch (axis)
+            {
+            case 0:
+                return v.x;
+            case 1:
+                return v.y;
+            default:
+                return v.z;
+            }
+        }
+
+        void query_overlaps_impl(std::size_t node_idx, const AabbType& query_aabb,
+                                 std::vector<std::size_t>& results) const
+        {
+            const auto& node = m_nodes[node_idx];
+
+            if (!overlaps(node.bounds, query_aabb))
+                return;
+
+            if (node.is_leaf())
+            {
+                for (auto i = node.first_index; i < node.first_index + node.index_count; ++i)
+                    if (overlaps(query_aabb, m_aabbs[m_indices[i]]))
+                        results.push_back(m_indices[i]);
+                return;
+            }
+
+            query_overlaps_impl(node.left, query_aabb, results);
+            query_overlaps_impl(node.right, query_aabb, results);
+        }
+
+        template<class RayType>
+        void query_ray_impl(std::size_t node_idx, const RayType& ray,
+                            std::vector<HitResult>& results) const
+        {
+            const auto& node = m_nodes[node_idx];
+
+            // Quick AABB-ray rejection using the slab method
+            const auto hit = LineTracer<RayType>::get_ray_hit_point(ray, node.bounds);
+            if (hit == ray.end)
+                return;
+
+            if (node.is_leaf())
+            {
+                for (auto i = node.first_index; i < node.first_index + node.index_count; ++i)
+                {
+                    const auto leaf_hit = LineTracer<RayType>::get_ray_hit_point(
+                        ray, m_aabbs[m_indices[i]]);
+                    if (leaf_hit != ray.end)
+                    {
+                        const auto diff = leaf_hit - ray.start;
+                        results.push_back({m_indices[i], diff.dot(diff)});
+                    }
+                }
+                return;
+            }
+
+            query_ray_impl(node.left, ray, results);
+            query_ray_impl(node.right, ray, results);
+        }
+
+        std::vector<Node> m_nodes;
+        std::vector<std::size_t> m_indices;
+        std::vector<AabbType> m_aabbs;
+    };
+} // namespace omath::collision

include/omath/collision/line_tracer.hpp

@@ -3,6 +3,7 @@
 //
 #pragma once
 
+#include "omath/3d_primitives/aabb.hpp"
 #include "omath/linear_algebra/triangle.hpp"
 #include "omath/linear_algebra/vector3.hpp"
 
@@ -34,6 +35,7 @@ namespace omath::collision
     class LineTracer final
     {
         using TriangleType = Triangle<typename RayType::VectorType>;
+        using AABBType = primitives::Aabb<typename RayType::VectorType::ContainedType>;
 
     public:
         LineTracer() = delete;
@@ -87,6 +89,54 @@ namespace omath::collision
             return ray.start + ray_dir * t_hit;
         }
 
+        // Slab method ray-AABB intersection
+        // Returns the hit point on the AABB surface, or ray.end if no intersection
+        [[nodiscard]]
+        constexpr static auto get_ray_hit_point(const RayType& ray, const AABBType& aabb) noexcept
+        {
+            using T = typename RayType::VectorType::ContainedType;
+            const auto dir = ray.direction_vector();
+
+            auto t_min = -std::numeric_limits<T>::infinity();
+            auto t_max = std::numeric_limits<T>::infinity();
+
+            const auto process_axis = [&](const T& d, const T& origin, const T& box_min,
+                                          const T& box_max) -> bool
+            {
+                constexpr T k_epsilon = std::numeric_limits<T>::epsilon();
+                if (std::abs(d) < k_epsilon)
+                    return origin >= box_min && origin <= box_max;
+
+                const T inv = T(1) / d;
+                T t0 = (box_min - origin) * inv;
+                T t1 = (box_max - origin) * inv;
+                if (t0 > t1)
+                    std::swap(t0, t1);
+
+                t_min = std::max(t_min, t0);
+                t_max = std::min(t_max, t1);
+                return t_min <= t_max;
+            };
+
+            if (!process_axis(dir.x, ray.start.x, aabb.min.x, aabb.max.x))
+                return ray.end;
+            if (!process_axis(dir.y, ray.start.y, aabb.min.y, aabb.max.y))
+                return ray.end;
+            if (!process_axis(dir.z, ray.start.z, aabb.min.z, aabb.max.z))
+                return ray.end;
+
+            // t_hit: use entry point if in front of origin, otherwise 0 (started inside)
+            const T t_hit = std::max(T(0), t_min);
+
+            if (t_max < T(0))
+                return ray.end; // box entirely behind origin
+
+            if (!ray.infinite_length && t_hit > T(1))
+                return ray.end; // box beyond ray endpoint
+
+            return ray.start + dir * t_hit;
+        }
+
         template<class MeshType>
         [[nodiscard]]
         constexpr static auto get_ray_hit_point(const RayType& ray, const MeshType& mesh) noexcept

include/omath/engines/cry_engine/camera.hpp

@@ -9,5 +9,5 @@
 
 namespace omath::cry_engine
 {
-    using Camera = projection::Camera<Mat4X4, ViewAngles, CameraTrait>;
+    using Camera = projection::Camera<Mat4X4, ViewAngles, CameraTrait, false, NDCDepthRange::ZERO_TO_ONE>;
 } // namespace omath::cry_engine

include/omath/engines/cry_engine/formulas.hpp

@@ -22,7 +22,8 @@ namespace omath::cry_engine
     Mat4X4 rotation_matrix(const ViewAngles& angles) noexcept;
 
     [[nodiscard]]
-    Mat4X4 calc_perspective_projection_matrix(float field_of_view, float aspect_ratio, float near, float far) noexcept;
+    Mat4X4 calc_perspective_projection_matrix(float field_of_view, float aspect_ratio, float near, float far,
+                                             NDCDepthRange ndc_depth_range = NDCDepthRange::NEGATIVE_ONE_TO_ONE) noexcept;
 
     template<class FloatingType>
     requires std::is_floating_point_v<FloatingType>

include/omath/engines/cry_engine/traits/camera_trait.hpp

@@ -18,7 +18,7 @@ namespace omath::cry_engine
         static Mat4X4 calc_view_matrix(const ViewAngles& angles, const Vector3<float>& cam_origin) noexcept;
         [[nodiscard]]
         static Mat4X4 calc_projection_matrix(const projection::FieldOfView& fov, const projection::ViewPort& view_port,
-                                             float near, float far) noexcept;
+                                             float near, float far, NDCDepthRange ndc_depth_range) noexcept;
     };
 
 } // namespace omath::cry_engine

include/omath/engines/frostbite_engine/camera.hpp

@@ -9,5 +9,5 @@
 
 namespace omath::frostbite_engine
 {
-    using Camera =  projection::Camera<Mat4X4, ViewAngles, CameraTrait>;
+    using Camera =  projection::Camera<Mat4X4, ViewAngles, CameraTrait, false, NDCDepthRange::ZERO_TO_ONE>;
 } // namespace omath::unity_engine

include/omath/engines/frostbite_engine/formulas.hpp

@@ -22,7 +22,8 @@ namespace omath::frostbite_engine
     Mat4X4 rotation_matrix(const ViewAngles& angles) noexcept;
 
     [[nodiscard]]
-    Mat4X4 calc_perspective_projection_matrix(float field_of_view, float aspect_ratio, float near, float far) noexcept;
+    Mat4X4 calc_perspective_projection_matrix(float field_of_view, float aspect_ratio, float near, float far,
+                                             NDCDepthRange ndc_depth_range = NDCDepthRange::NEGATIVE_ONE_TO_ONE) noexcept;
 
     template<class FloatingType>
     requires std::is_floating_point_v<FloatingType>

include/omath/engines/frostbite_engine/traits/camera_trait.hpp

@@ -18,7 +18,7 @@ namespace omath::frostbite_engine
         static Mat4X4 calc_view_matrix(const ViewAngles& angles, const Vector3<float>& cam_origin) noexcept;
         [[nodiscard]]
         static Mat4X4 calc_projection_matrix(const projection::FieldOfView& fov, const projection::ViewPort& view_port,
-                                             float near, float far) noexcept;
+                                             float near, float far, NDCDepthRange ndc_depth_range) noexcept;
     };
 
 } // namespace omath::unreal_engine

include/omath/engines/iw_engine/camera.hpp

@@ -9,5 +9,5 @@
 
 namespace omath::iw_engine
 {
-    using Camera =  projection::Camera<Mat4X4, ViewAngles, CameraTrait>;
+    using Camera =  projection::Camera<Mat4X4, ViewAngles, CameraTrait, false, NDCDepthRange::ZERO_TO_ONE>;
 } // namespace omath::iw_engine

include/omath/engines/iw_engine/formulas.hpp

@@ -22,7 +22,8 @@ namespace omath::iw_engine
     [[nodiscard]] Mat4X4 calc_view_matrix(const ViewAngles& angles, const Vector3<float>& cam_origin) noexcept;
 
     [[nodiscard]]
-    Mat4X4 calc_perspective_projection_matrix(float field_of_view, float aspect_ratio, float near, float far) noexcept;
+    Mat4X4 calc_perspective_projection_matrix(float field_of_view, float aspect_ratio, float near, float far,
+                                             NDCDepthRange ndc_depth_range = NDCDepthRange::NEGATIVE_ONE_TO_ONE) noexcept;
 
     template<class FloatingType>
     requires std::is_floating_point_v<FloatingType>

include/omath/engines/iw_engine/traits/camera_trait.hpp

@@ -18,7 +18,7 @@ namespace omath::iw_engine
         static Mat4X4 calc_view_matrix(const ViewAngles& angles, const Vector3<float>& cam_origin) noexcept;
         [[nodiscard]]
         static Mat4X4 calc_projection_matrix(const projection::FieldOfView& fov, const projection::ViewPort& view_port,
-                                             float near, float far) noexcept;
+                                             float near, float far, NDCDepthRange ndc_depth_range) noexcept;
     };
 
 } // namespace omath::iw_engine

include/omath/engines/opengl_engine/camera.hpp

@@ -8,5 +8,5 @@
 
 namespace omath::opengl_engine
 {
-    using Camera =  projection::Camera<Mat4X4, ViewAngles, CameraTrait, true>;
+    using Camera =  projection::Camera<Mat4X4, ViewAngles, CameraTrait, true, NDCDepthRange::NEGATIVE_ONE_TO_ONE>;
 } // namespace omath::opengl_engine

include/omath/engines/opengl_engine/formulas.hpp

@@ -21,7 +21,8 @@ namespace omath::opengl_engine
     Mat4X4 rotation_matrix(const ViewAngles& angles) noexcept;
 
     [[nodiscard]]
-    Mat4X4 calc_perspective_projection_matrix(float field_of_view, float aspect_ratio, float near, float far) noexcept;
+    Mat4X4 calc_perspective_projection_matrix(float field_of_view, float aspect_ratio, float near, float far,
+                                             NDCDepthRange ndc_depth_range = NDCDepthRange::NEGATIVE_ONE_TO_ONE) noexcept;
 
     template<class FloatingType>
     requires std::is_floating_point_v<FloatingType>

include/omath/engines/opengl_engine/traits/camera_trait.hpp

@@ -18,7 +18,7 @@ namespace omath::opengl_engine
         static Mat4X4 calc_view_matrix(const ViewAngles& angles, const Vector3<float>& cam_origin) noexcept;
         [[nodiscard]]
         static Mat4X4 calc_projection_matrix(const projection::FieldOfView& fov, const projection::ViewPort& view_port,
-                                             float near, float far) noexcept;
+                                             float near, float far, NDCDepthRange ndc_depth_range) noexcept;
     };
 
 } // namespace omath::opengl_engine

include/omath/engines/source_engine/camera.hpp

@@ -7,5 +7,5 @@
 #include "traits/camera_trait.hpp"
 namespace omath::source_engine
 {
-    using Camera =  projection::Camera<Mat4X4, ViewAngles, CameraTrait>;
+    using Camera =  projection::Camera<Mat4X4, ViewAngles, CameraTrait, false, NDCDepthRange::ZERO_TO_ONE>;
 } // namespace omath::source_engine

include/omath/engines/source_engine/formulas.hpp

@@ -21,7 +21,8 @@ namespace omath::source_engine
     [[nodiscard]] Mat4X4 calc_view_matrix(const ViewAngles& angles, const Vector3<float>& cam_origin) noexcept;
 
     [[nodiscard]]
-    Mat4X4 calc_perspective_projection_matrix(float field_of_view, float aspect_ratio, float near, float far) noexcept;
+    Mat4X4 calc_perspective_projection_matrix(float field_of_view, float aspect_ratio, float near, float far,
+                                             NDCDepthRange ndc_depth_range = NDCDepthRange::NEGATIVE_ONE_TO_ONE) noexcept;
 
     template<class FloatingType>
     requires std::is_floating_point_v<FloatingType>

include/omath/engines/source_engine/traits/camera_trait.hpp

@@ -18,7 +18,7 @@ namespace omath::source_engine
         static Mat4X4 calc_view_matrix(const ViewAngles& angles, const Vector3<float>& cam_origin) noexcept;
         [[nodiscard]]
         static Mat4X4 calc_projection_matrix(const projection::FieldOfView& fov, const projection::ViewPort& view_port,
-                                             float near, float far) noexcept;
+                                             float near, float far, NDCDepthRange ndc_depth_range) noexcept;
     };
 
 } // namespace omath::source_engine

include/omath/engines/unity_engine/camera.hpp

@@ -9,5 +9,5 @@
 
 namespace omath::unity_engine
 {
-    using Camera = projection::Camera<Mat4X4, ViewAngles, CameraTrait>;
+    using Camera = projection::Camera<Mat4X4, ViewAngles, CameraTrait, false, NDCDepthRange::ZERO_TO_ONE>;
 } // namespace omath::unity_engine

include/omath/engines/unity_engine/formulas.hpp

@@ -22,7 +22,8 @@ namespace omath::unity_engine
     Mat4X4 rotation_matrix(const ViewAngles& angles) noexcept;
 
     [[nodiscard]]
-    Mat4X4 calc_perspective_projection_matrix(float field_of_view, float aspect_ratio, float near, float far) noexcept;
+    Mat4X4 calc_perspective_projection_matrix(float field_of_view, float aspect_ratio, float near, float far,
+                                             NDCDepthRange ndc_depth_range = NDCDepthRange::NEGATIVE_ONE_TO_ONE) noexcept;
 
     template<class FloatingType>
     requires std::is_floating_point_v<FloatingType>

include/omath/engines/unity_engine/traits/camera_trait.hpp

@@ -18,7 +18,7 @@ namespace omath::unity_engine
         static Mat4X4 calc_view_matrix(const ViewAngles& angles, const Vector3<float>& cam_origin) noexcept;
         [[nodiscard]]
         static Mat4X4 calc_projection_matrix(const projection::FieldOfView& fov, const projection::ViewPort& view_port,
-                                             float near, float far) noexcept;
+                                             float near, float far, NDCDepthRange ndc_depth_range) noexcept;
     };
 
 } // namespace omath::unity_engine

include/omath/engines/unreal_engine/camera.hpp

@@ -9,5 +9,5 @@
 
 namespace omath::unreal_engine
 {
-    using Camera = projection::Camera<Mat4X4, ViewAngles, CameraTrait>;
+    using Camera = projection::Camera<Mat4X4, ViewAngles, CameraTrait, false, NDCDepthRange::ZERO_TO_ONE>;
 } // namespace omath::unreal_engine

include/omath/engines/unreal_engine/formulas.hpp

@@ -22,7 +22,8 @@ namespace omath::unreal_engine
     Mat4X4 rotation_matrix(const ViewAngles& angles) noexcept;
 
     [[nodiscard]]
-    Mat4X4 calc_perspective_projection_matrix(float field_of_view, float aspect_ratio, float near, float far) noexcept;
+    Mat4X4 calc_perspective_projection_matrix(float field_of_view, float aspect_ratio, float near, float far,
+                                             NDCDepthRange ndc_depth_range = NDCDepthRange::NEGATIVE_ONE_TO_ONE) noexcept;
 
     template<class FloatingType>
     requires std::is_floating_point_v<FloatingType>

include/omath/engines/unreal_engine/traits/camera_trait.hpp

@@ -18,7 +18,7 @@ namespace omath::unreal_engine
         static Mat4X4 calc_view_matrix(const ViewAngles& angles, const Vector3<float>& cam_origin) noexcept;
         [[nodiscard]]
         static Mat4X4 calc_projection_matrix(const projection::FieldOfView& fov, const projection::ViewPort& view_port,
-                                             float near, float far) noexcept;
+                                             float near, float far, NDCDepthRange ndc_depth_range) noexcept;
     };
 
 } // namespace omath::unreal_engine

include/omath/linear_algebra/mat.hpp

@@ -37,6 +37,12 @@ namespace omath
         COLUMN_MAJOR
     };
 
+    enum class NDCDepthRange : uint8_t
+    {
+        NEGATIVE_ONE_TO_ONE = 0, // OpenGL: [-1.0, 1.0]
+        ZERO_TO_ONE              // DirectX / Vulkan: [0.0, 1.0]
+    };
+
     template<typename M1, typename M2> concept MatTemplateEqual
             = (M1::rows == M2::rows) && (M1::columns == M2::columns)
               && std::is_same_v<typename M1::value_type, typename M2::value_type> && (M1::store_type == M2::store_type);
@@ -658,56 +664,98 @@ namespace omath
         } * mat_translation<Type, St>(-camera_origin);
     }
 
-    template<class Type = float, MatStoreType St = MatStoreType::ROW_MAJOR>
+    template<class Type = float, MatStoreType St = MatStoreType::ROW_MAJOR,
+             NDCDepthRange DepthRange = NDCDepthRange::NEGATIVE_ONE_TO_ONE>
     [[nodiscard]]
     Mat<4, 4, Type, St> mat_perspective_left_handed(const float field_of_view, const float aspect_ratio,
                                                     const float near, const float far) noexcept
     {
         const float fov_half_tan = std::tan(angles::degrees_to_radians(field_of_view) / 2.f);
 
-        return {{1.f / (aspect_ratio * fov_half_tan), 0.f, 0.f, 0.f},
-                {0.f, 1.f / fov_half_tan, 0.f, 0.f},
-                {0.f, 0.f, (far + near) / (far - near), -(2.f * near * far) / (far - near)},
-                {0.f, 0.f, 1.f, 0.f}};
+        if constexpr (DepthRange == NDCDepthRange::ZERO_TO_ONE)
+            return {{1.f / (aspect_ratio * fov_half_tan), 0.f, 0.f, 0.f},
+                    {0.f, 1.f / fov_half_tan, 0.f, 0.f},
+                    {0.f, 0.f, far / (far - near), -(near * far) / (far - near)},
+                    {0.f, 0.f, 1.f, 0.f}};
+        else if constexpr (DepthRange == NDCDepthRange::NEGATIVE_ONE_TO_ONE)
+            return {{1.f / (aspect_ratio * fov_half_tan), 0.f, 0.f, 0.f},
+                    {0.f, 1.f / fov_half_tan, 0.f, 0.f},
+                    {0.f, 0.f, (far + near) / (far - near), -(2.f * near * far) / (far - near)},
+                    {0.f, 0.f, 1.f, 0.f}};
+        else
+            std::unreachable();
     }
 
-    template<class Type = float, MatStoreType St = MatStoreType::ROW_MAJOR>
+    template<class Type = float, MatStoreType St = MatStoreType::ROW_MAJOR,
+             NDCDepthRange DepthRange = NDCDepthRange::NEGATIVE_ONE_TO_ONE>
     [[nodiscard]]
     Mat<4, 4, Type, St> mat_perspective_right_handed(const float field_of_view, const float aspect_ratio,
                                                      const float near, const float far) noexcept
     {
         const float fov_half_tan = std::tan(angles::degrees_to_radians(field_of_view) / 2.f);
 
-        return {{1.f / (aspect_ratio * fov_half_tan), 0.f, 0.f, 0.f},
-                {0.f, 1.f / fov_half_tan, 0.f, 0.f},
-                {0.f, 0.f, -(far + near) / (far - near), -(2.f * near * far) / (far - near)},
-                {0.f, 0.f, -1.f, 0.f}};
+        if constexpr (DepthRange == NDCDepthRange::ZERO_TO_ONE)
+            return {{1.f / (aspect_ratio * fov_half_tan), 0.f, 0.f, 0.f},
+                    {0.f, 1.f / fov_half_tan, 0.f, 0.f},
+                    {0.f, 0.f, -far / (far - near), -(near * far) / (far - near)},
+                    {0.f, 0.f, -1.f, 0.f}};
+        else if constexpr (DepthRange == NDCDepthRange::NEGATIVE_ONE_TO_ONE)
+            return {{1.f / (aspect_ratio * fov_half_tan), 0.f, 0.f, 0.f},
+                    {0.f, 1.f / fov_half_tan, 0.f, 0.f},
+                    {0.f, 0.f, -(far + near) / (far - near), -(2.f * near * far) / (far - near)},
+                    {0.f, 0.f, -1.f, 0.f}};
+        else
+            std::unreachable();
     }
-    template<class Type = float, MatStoreType St = MatStoreType::ROW_MAJOR>
+    template<class Type = float, MatStoreType St = MatStoreType::ROW_MAJOR,
+             NDCDepthRange DepthRange = NDCDepthRange::NEGATIVE_ONE_TO_ONE>
     [[nodiscard]]
     Mat<4, 4, Type, St> mat_ortho_left_handed(const Type left, const Type right, const Type bottom, const Type top,
                                               const Type near, const Type far) noexcept
     {
-        return
-        {
-            { static_cast<Type>(2) / (right - left), 0.f,       0.f,    -(right + left) / (right - left)},
-            { 0.f,      static_cast<Type>(2) / (top - bottom),  0.f,    -(top + bottom) / (top - bottom)},
-            { 0.f,      0.f,       static_cast<Type>(2) / (far - near), -(far + near) / (far - near)    },
-            { 0.f,      0.f,       0.f,                                 1.f                             }
-        };
+        if constexpr (DepthRange == NDCDepthRange::ZERO_TO_ONE)
+            return
+            {
+                { static_cast<Type>(2) / (right - left), 0.f,       0.f,    -(right + left) / (right - left)},
+                { 0.f,      static_cast<Type>(2) / (top - bottom),  0.f,    -(top + bottom) / (top - bottom)},
+                { 0.f,      0.f,       static_cast<Type>(1) / (far - near), -near / (far - near)            },
+                { 0.f,      0.f,       0.f,                                 1.f                             }
+            };
+        else if constexpr (DepthRange == NDCDepthRange::NEGATIVE_ONE_TO_ONE)
+            return
+            {
+                { static_cast<Type>(2) / (right - left), 0.f,       0.f,    -(right + left) / (right - left)},
+                { 0.f,      static_cast<Type>(2) / (top - bottom),  0.f,    -(top + bottom) / (top - bottom)},
+                { 0.f,      0.f,       static_cast<Type>(2) / (far - near), -(far + near) / (far - near)    },
+                { 0.f,      0.f,       0.f,                                 1.f                             }
+            };
+        else
+            std::unreachable();
     }
-    template<class Type = float, MatStoreType St = MatStoreType::ROW_MAJOR>
+    template<class Type = float, MatStoreType St = MatStoreType::ROW_MAJOR,
+             NDCDepthRange DepthRange = NDCDepthRange::NEGATIVE_ONE_TO_ONE>
     [[nodiscard]]
     Mat<4, 4, Type, St> mat_ortho_right_handed(const Type left, const Type right, const Type bottom, const Type top,
                                                const Type near, const Type far) noexcept
     {
-        return
-        {
-                { static_cast<Type>(2) / (right - left), 0.f,       0.f,     -(right + left) / (right - left)},
-                { 0.f,      static_cast<Type>(2) / (top - bottom),  0.f,     -(top + bottom) / (top - bottom)},
-                { 0.f,      0.f,       -static_cast<Type>(2) / (far - near), -(far + near) / (far - near)    },
-                { 0.f,      0.f,       0.f,                                  1.f                             }
-        };
+        if constexpr (DepthRange == NDCDepthRange::ZERO_TO_ONE)
+            return
+            {
+                    { static_cast<Type>(2) / (right - left), 0.f,       0.f,     -(right + left) / (right - left)},
+                    { 0.f,      static_cast<Type>(2) / (top - bottom),  0.f,     -(top + bottom) / (top - bottom)},
+                    { 0.f,      0.f,       -static_cast<Type>(1) / (far - near), -near / (far - near)            },
+                    { 0.f,      0.f,       0.f,                                  1.f                             }
+            };
+        else if constexpr (DepthRange == NDCDepthRange::NEGATIVE_ONE_TO_ONE)
+            return
+            {
+                    { static_cast<Type>(2) / (right - left), 0.f,       0.f,     -(right + left) / (right - left)},
+                    { 0.f,      static_cast<Type>(2) / (top - bottom),  0.f,     -(top + bottom) / (top - bottom)},
+                    { 0.f,      0.f,       -static_cast<Type>(2) / (far - near), -(far + near) / (far - near)    },
+                    { 0.f,      0.f,       0.f,                                  1.f                             }
+            };
+        else
+            std::unreachable();
     }
     template<class T = float, MatStoreType St = MatStoreType::COLUMN_MAJOR>
    Mat<4, 4, T, St> mat_look_at_left_handed(const Vector3<T>& eye, const Vector3<T>& center, const Vector3<T>& up)

include/omath/omath.hpp

@@ -35,6 +35,7 @@
 #include "omath/collision/line_tracer.hpp"
 #include "omath/collision/gjk_algorithm.hpp"
 #include "omath/collision/epa_algorithm.hpp"
+#include "omath/collision/bvh_tree.hpp"
 // Pathfinding algorithms
 #include "omath/pathfinding/a_star.hpp"
 #include "omath/pathfinding/navigation_mesh.hpp"

include/omath/projection/camera.hpp

@@ -4,6 +4,7 @@
 
 #pragma once
 
+#include "omath/3d_primitives/aabb.hpp"
 #include "omath/linear_algebra/mat.hpp"
 #include "omath/linear_algebra/triangle.hpp"
 #include "omath/linear_algebra/vector3.hpp"
@@ -44,19 +45,21 @@ namespace omath::projection
     template<class T, class MatType, class ViewAnglesType>
     concept CameraEngineConcept =
             requires(const Vector3<float>& cam_origin, const Vector3<float>& look_at, const ViewAnglesType& angles,
-                     const FieldOfView& fov, const ViewPort& viewport, float znear, float zfar) {
+                     const FieldOfView& fov, const ViewPort& viewport, float znear, float zfar,
+                     NDCDepthRange ndc_depth_range) {
                 // Presence + return types
                 { T::calc_look_at_angle(cam_origin, look_at) } -> std::same_as<ViewAnglesType>;
                 { T::calc_view_matrix(angles, cam_origin) } -> std::same_as<MatType>;
-                { T::calc_projection_matrix(fov, viewport, znear, zfar) } -> std::same_as<MatType>;
+                { T::calc_projection_matrix(fov, viewport, znear, zfar, ndc_depth_range) } -> std::same_as<MatType>;
 
                 // Enforce noexcept as in the trait declaration
                 requires noexcept(T::calc_look_at_angle(cam_origin, look_at));
                 requires noexcept(T::calc_view_matrix(angles, cam_origin));
-                requires noexcept(T::calc_projection_matrix(fov, viewport, znear, zfar));
+                requires noexcept(T::calc_projection_matrix(fov, viewport, znear, zfar, ndc_depth_range));
             };
 
-    template<class Mat4X4Type, class ViewAnglesType, class TraitClass, bool inverted_z = false>
+    template<class Mat4X4Type, class ViewAnglesType, class TraitClass, bool inverted_z = false,
+             NDCDepthRange depth_range = NDCDepthRange::NEGATIVE_ONE_TO_ONE>
     requires CameraEngineConcept<TraitClass, Mat4X4Type, ViewAnglesType>
     class Camera final
     {
@@ -135,7 +138,8 @@ namespace omath::projection
         {
             if (!m_projection_matrix.has_value())
                 m_projection_matrix = TraitClass::calc_projection_matrix(m_field_of_view, m_view_port,
-                                                                         m_near_plane_distance, m_far_plane_distance);
+                                                                         m_near_plane_distance, m_far_plane_distance,
+                                                                         depth_range);
 
             return m_projection_matrix.value();
         }
@@ -271,18 +275,94 @@ namespace omath::projection
                 return a[axis] < -a[3] && b[axis] < -b[3] && c[axis] < -c[3];
             };
 
-            // Clip volume in clip space (OpenGL-style):
+            // Clip volume in clip space:
             // -w <= x <= w
             // -w <= y <= w
-            // -w <= z <= w
+            // z_min <= z <= w  (z_min = -w for [-1,1], 0 for [0,1])
 
-            for (int i = 0; i < 3; i++)
+            // x and y planes
+            for (int i = 0; i < 2; i++)
             {
                 if (all_outside_plane(i, c0, c1, c2, false))
-                    return true; // x < -w (left)
+                    return true;
                 if (all_outside_plane(i, c0, c1, c2, true))
-                    return true; // x >  w (right)
+                    return true;
             }
+
+            // z far plane: z > w
+            if (all_outside_plane(2, c0, c1, c2, true))
+                return true;
+
+            // z near plane
+            if constexpr (depth_range == NDCDepthRange::ZERO_TO_ONE)
+            {
+                // 0 <= z, so reject if z < 0 for all vertices
+                if (c0[2] < 0.f && c1[2] < 0.f && c2[2] < 0.f)
+                    return true;
+            }
+            else
+            {
+                // -w <= z
+                if (all_outside_plane(2, c0, c1, c2, false))
+                    return true;
+            }
+            return false;
+        }
+
+        [[nodiscard]] bool is_aabb_culled_by_frustum(const primitives::Aabb<float>& aabb) const noexcept
+        {
+            const auto& m = get_view_projection_matrix();
+
+            // Gribb-Hartmann: extract 6 frustum planes from the view-projection matrix.
+            // Each plane is (a, b, c, d) such that ax + by + cz + d >= 0 means inside.
+            // For a 4x4 matrix with rows r0..r3:
+            //   Left   = r3 + r0
+            //   Right  = r3 - r0
+            //   Bottom = r3 + r1
+            //   Top    = r3 - r1
+            //   Near   = r3 + r2  ([-1,1]) or r2 ([0,1])
+            //   Far    = r3 - r2
+            struct Plane final
+            {
+                float a, b, c, d;
+            };
+
+            const auto extract_plane = [&m](const int sign, const int row) -> Plane
+            {
+                return {
+                        m.at(3, 0) + static_cast<float>(sign) * m.at(row, 0),
+                        m.at(3, 1) + static_cast<float>(sign) * m.at(row, 1),
+                        m.at(3, 2) + static_cast<float>(sign) * m.at(row, 2),
+                        m.at(3, 3) + static_cast<float>(sign) * m.at(row, 3),
+                };
+            };
+
+            std::array<Plane, 6> planes = {
+                    extract_plane(1, 0), // left
+                    extract_plane(-1, 0), // right
+                    extract_plane(1, 1), // bottom
+                    extract_plane(-1, 1), // top
+                    extract_plane(-1, 2), // far
+            };
+
+            // Near plane depends on NDC depth range
+            if constexpr (depth_range == NDCDepthRange::ZERO_TO_ONE)
+                planes[5] = {m.at(2, 0), m.at(2, 1), m.at(2, 2), m.at(2, 3)};
+            else
+                planes[5] = extract_plane(1, 2);
+
+            // For each plane, find the AABB corner most in the direction of the plane normal
+            // (the "positive vertex"). If it's outside, the entire AABB is outside.
+            for (const auto& [a, b, c, d] : planes)
+            {
+                const float px = a >= 0.f ? aabb.max.x : aabb.min.x;
+                const float py = b >= 0.f ? aabb.max.y : aabb.min.y;
+                const float pz = c >= 0.f ? aabb.max.z : aabb.min.z;
+
+                if (a * px + b * py + c * pz + d < 0.f)
+                    return true;
+            }
+
             return false;
         }
 
@@ -306,7 +386,8 @@ namespace omath::projection
                 return std::unexpected(Error::WORLD_POSITION_IS_OUT_OF_SCREEN_BOUNDS);
 
             // ReSharper disable once CppTooWideScope
-            const auto clipped_manually = clipping == ViewPortClipping::MANUAL && (projected.at(2, 0) < 0.0f - eps
+            constexpr auto z_min = depth_range == NDCDepthRange::ZERO_TO_ONE ? 0.0f : -1.0f;
+            const auto clipped_manually = clipping == ViewPortClipping::MANUAL && (projected.at(2, 0) < z_min - eps
                                           || projected.at(2, 0) > 1.0f + eps);
             if (clipped_manually)
                 return std::unexpected(Error::WORLD_POSITION_IS_OUT_OF_SCREEN_BOUNDS);
@@ -368,8 +449,26 @@ namespace omath::projection
         [[nodiscard]] constexpr static bool is_ndc_out_of_bounds(const Type& ndc) noexcept
         {
             constexpr auto eps = std::numeric_limits<float>::epsilon();
-            return std::ranges::any_of(ndc.raw_array(),
-                                       [](const auto& val) { return val < -1.0f - eps || val > 1.0f + eps; });
+
+            const auto& data = ndc.raw_array();
+            // x and y are always in [-1, 1]
+            if (data[0] < -1.0f - eps || data[0] > 1.0f + eps)
+                return true;
+            if (data[1] < -1.0f - eps || data[1] > 1.0f + eps)
+                return true;
+            return is_ndc_z_value_out_of_bounds(data[2]);
+        }
+        template<class ZType>
+         [[nodiscard]]
+        constexpr static bool is_ndc_z_value_out_of_bounds(const ZType& z_ndc) noexcept
+        {
+            constexpr auto eps = std::numeric_limits<float>::epsilon();
+            if constexpr (depth_range == NDCDepthRange::NEGATIVE_ONE_TO_ONE)
+                return z_ndc < -1.0f - eps || z_ndc > 1.0f + eps;
+            if constexpr (depth_range == NDCDepthRange::ZERO_TO_ONE)
+                return z_ndc < 0.0f - eps || z_ndc > 1.0f + eps;
+
+            std::unreachable();
         }
 
         // NDC REPRESENTATION:

source/engines/cry_engine/formulas.cpp

@@ -35,8 +35,15 @@ namespace omath::cry_engine
                * mat_rotation_axis_x<float, MatStoreType::ROW_MAJOR>(angles.pitch);
     }
     Mat4X4 calc_perspective_projection_matrix(const float field_of_view, const float aspect_ratio, const float near,
-                                              const float far) noexcept
+                                              const float far, const NDCDepthRange ndc_depth_range) noexcept
     {
-        return mat_perspective_left_handed(field_of_view, aspect_ratio, near, far);
+        if (ndc_depth_range == NDCDepthRange::ZERO_TO_ONE)
+            return mat_perspective_left_handed<float, MatStoreType::ROW_MAJOR, NDCDepthRange::ZERO_TO_ONE>(
+                    field_of_view, aspect_ratio, near, far);
+
+        if (ndc_depth_range == NDCDepthRange::NEGATIVE_ONE_TO_ONE)
+            return mat_perspective_left_handed<float, MatStoreType::ROW_MAJOR, NDCDepthRange::NEGATIVE_ONE_TO_ONE>(
+                    field_of_view, aspect_ratio, near, far);
+        std::unreachable();
     }
 } // namespace omath::unity_engine

source/engines/cry_engine/traits/camera_trait.cpp

@@ -19,8 +19,9 @@ namespace omath::cry_engine
     }
     Mat4X4 CameraTrait::calc_projection_matrix(const projection::FieldOfView& fov,
                                                const projection::ViewPort& view_port, const float near,
-                                               const float far) noexcept
+                                               const float far, const NDCDepthRange ndc_depth_range) noexcept
     {
-        return calc_perspective_projection_matrix(fov.as_degrees(), view_port.aspect_ratio(), near, far);
+        return calc_perspective_projection_matrix(fov.as_degrees(), view_port.aspect_ratio(), near, far,
+                                                  ndc_depth_range);
     }
 } // namespace omath::unity_engine

source/engines/frostbite_engine/formulas.cpp

@@ -35,8 +35,16 @@ namespace omath::frostbite_engine
                * mat_rotation_axis_x<float, MatStoreType::ROW_MAJOR>(angles.pitch);
     }
     Mat4X4 calc_perspective_projection_matrix(const float field_of_view, const float aspect_ratio, const float near,
-                                              const float far) noexcept
+                                              const float far, const NDCDepthRange ndc_depth_range) noexcept
     {
-        return mat_perspective_left_handed(field_of_view, aspect_ratio, near, far);
+        if (ndc_depth_range == NDCDepthRange::ZERO_TO_ONE)
+            return mat_perspective_left_handed<float, MatStoreType::ROW_MAJOR, NDCDepthRange::ZERO_TO_ONE>(
+                    field_of_view, aspect_ratio, near, far);
+
+        if (ndc_depth_range == NDCDepthRange::NEGATIVE_ONE_TO_ONE)
+            return mat_perspective_left_handed<float, MatStoreType::ROW_MAJOR, NDCDepthRange::NEGATIVE_ONE_TO_ONE>(
+                    field_of_view, aspect_ratio, near, far);
+
+        std::unreachable();
     }
 } // namespace omath::unity_engine

source/engines/frostbite_engine/traits/camera_trait.cpp

@@ -19,8 +19,9 @@ namespace omath::frostbite_engine
     }
     Mat4X4 CameraTrait::calc_projection_matrix(const projection::FieldOfView& fov,
                                                const projection::ViewPort& view_port, const float near,
-                                               const float far) noexcept
+                                               const float far, const NDCDepthRange ndc_depth_range) noexcept
     {
-        return calc_perspective_projection_matrix(fov.as_degrees(), view_port.aspect_ratio(), near, far);
+        return calc_perspective_projection_matrix(fov.as_degrees(), view_port.aspect_ratio(), near, far,
+                                                  ndc_depth_range);
     }
 } // namespace omath::unity_engine

source/engines/iw_engine/formulas.cpp

@@ -36,18 +36,27 @@ namespace omath::iw_engine
     }
 
     Mat4X4 calc_perspective_projection_matrix(const float field_of_view, const float aspect_ratio, const float near,
-                                              const float far) noexcept
+                                              const float far, const NDCDepthRange ndc_depth_range) noexcept
     {
         // NOTE: Need magic number to fix fov calculation, since IW engine inherit Quake proj matrix calculation
         constexpr auto k_multiply_factor = 0.75f;
 
         const float fov_half_tan = std::tan(angles::degrees_to_radians(field_of_view) / 2.f) * k_multiply_factor;
 
-        return {
-                {1.f / (aspect_ratio * fov_half_tan), 0, 0, 0},
-                {0, 1.f / (fov_half_tan), 0, 0},
-                {0, 0, (far + near) / (far - near), -(2.f * far * near) / (far - near)},
-                {0, 0, 1, 0},
-        };
+        if (ndc_depth_range == NDCDepthRange::ZERO_TO_ONE)
+            return {
+                    {1.f / (aspect_ratio * fov_half_tan), 0, 0, 0},
+                    {0, 1.f / (fov_half_tan), 0, 0},
+                    {0, 0, far / (far - near), -(near * far) / (far - near)},
+                    {0, 0, 1, 0},
+            };
+        if (ndc_depth_range == NDCDepthRange::NEGATIVE_ONE_TO_ONE)
+            return {
+                    {1.f / (aspect_ratio * fov_half_tan), 0, 0, 0},
+                    {0, 1.f / (fov_half_tan), 0, 0},
+                    {0, 0, (far + near) / (far - near), -(2.f * far * near) / (far - near)},
+                    {0, 0, 1, 0},
+            };
+        std::unreachable();
     };
 } // namespace omath::iw_engine

source/engines/iw_engine/traits/camera_trait.cpp

@@ -19,8 +19,9 @@ namespace omath::iw_engine
     }
     Mat4X4 CameraTrait::calc_projection_matrix(const projection::FieldOfView& fov,
                                                const projection::ViewPort& view_port, const float near,
-                                               const float far) noexcept
+                                               const float far, const NDCDepthRange ndc_depth_range) noexcept
     {
-        return calc_perspective_projection_matrix(fov.as_degrees(), view_port.aspect_ratio(), near, far);
+        return calc_perspective_projection_matrix(fov.as_degrees(), view_port.aspect_ratio(), near, far,
+                                                  ndc_depth_range);
     }
 } // namespace omath::iw_engine

source/engines/opengl_engine/formulas.cpp

@@ -8,15 +8,15 @@ namespace omath::opengl_engine
 
     Vector3<float> forward_vector(const ViewAngles& angles) noexcept
     {
-        const auto vec
-                = rotation_matrix(angles) * mat_column_from_vector<float, MatStoreType::COLUMN_MAJOR>(k_abs_forward);
+        const auto vec =
+                rotation_matrix(angles) * mat_column_from_vector<float, MatStoreType::COLUMN_MAJOR>(k_abs_forward);
 
         return {vec.at(0, 0), vec.at(1, 0), vec.at(2, 0)};
     }
     Vector3<float> right_vector(const ViewAngles& angles) noexcept
     {
-        const auto vec
-                = rotation_matrix(angles) * mat_column_from_vector<float, MatStoreType::COLUMN_MAJOR>(k_abs_right);
+        const auto vec =
+                rotation_matrix(angles) * mat_column_from_vector<float, MatStoreType::COLUMN_MAJOR>(k_abs_right);
 
         return {vec.at(0, 0), vec.at(1, 0), vec.at(2, 0)};
     }
@@ -28,7 +28,7 @@ namespace omath::opengl_engine
     }
     Mat4X4 calc_view_matrix(const ViewAngles& angles, const Vector3<float>& cam_origin) noexcept
     {
-        return mat_look_at_right_handed(cam_origin, cam_origin+forward_vector(angles), up_vector(angles));
+        return mat_look_at_right_handed(cam_origin, cam_origin + forward_vector(angles), up_vector(angles));
     }
     Mat4X4 rotation_matrix(const ViewAngles& angles) noexcept
     {
@@ -37,15 +37,16 @@ namespace omath::opengl_engine
                * mat_rotation_axis_x<float, MatStoreType::COLUMN_MAJOR>(angles.pitch);
     }
     Mat4X4 calc_perspective_projection_matrix(const float field_of_view, const float aspect_ratio, const float near,
-                                              const float far) noexcept
+                                              const float far, const NDCDepthRange ndc_depth_range) noexcept
     {
-        const float fov_half_tan = std::tan(angles::degrees_to_radians(field_of_view) / 2.f);
+        if (ndc_depth_range == NDCDepthRange::NEGATIVE_ONE_TO_ONE)
+            return mat_perspective_right_handed<float, MatStoreType::COLUMN_MAJOR, NDCDepthRange::NEGATIVE_ONE_TO_ONE>(
+                    field_of_view, aspect_ratio, near, far);
 
-        return {
-                {1.f / (aspect_ratio * fov_half_tan), 0, 0, 0},
-                {0, 1.f / (fov_half_tan), 0, 0},
-                {0, 0, -(far + near) / (far - near), -(2.f * far * near) / (far - near)},
-                {0, 0, -1, 0},
-        };
+        if (ndc_depth_range == NDCDepthRange::ZERO_TO_ONE)
+            return mat_perspective_right_handed<float, MatStoreType::COLUMN_MAJOR, NDCDepthRange::ZERO_TO_ONE>(
+                        field_of_view, aspect_ratio, near, far);
+
+        std::unreachable();
     }
 } // namespace omath::opengl_engine

source/engines/opengl_engine/traits/camera_trait.cpp

@@ -20,8 +20,9 @@ namespace omath::opengl_engine
     }
     Mat4X4 CameraTrait::calc_projection_matrix(const projection::FieldOfView& fov,
                                                const projection::ViewPort& view_port, const float near,
-                                               const float far) noexcept
+                                               const float far, const NDCDepthRange ndc_depth_range) noexcept
     {
-        return calc_perspective_projection_matrix(fov.as_degrees(), view_port.aspect_ratio(), near, far);
+        return calc_perspective_projection_matrix(fov.as_degrees(), view_port.aspect_ratio(), near, far,
+                                                  ndc_depth_range);
     }
 } // namespace omath::opengl_engine

source/engines/source_engine/formulas.cpp

@@ -36,18 +36,27 @@ namespace omath::source_engine
     }
 
     Mat4X4 calc_perspective_projection_matrix(const float field_of_view, const float aspect_ratio, const float near,
-                                              const float far) noexcept
+                                              const float far, const NDCDepthRange ndc_depth_range) noexcept
     {
         // NOTE: Need magic number to fix fov calculation, since source inherit Quake proj matrix calculation
         constexpr auto k_multiply_factor = 0.75f;
 
         const float fov_half_tan = std::tan(angles::degrees_to_radians(field_of_view) / 2.f) * k_multiply_factor;
 
-        return {
-                {1.f / (aspect_ratio * fov_half_tan), 0, 0, 0},
-                {0, 1.f / (fov_half_tan), 0, 0},
-                {0, 0, (far + near) / (far - near), -(2.f * far * near) / (far - near)},
-                {0, 0, 1, 0},
-        };
+        if (ndc_depth_range == NDCDepthRange::ZERO_TO_ONE)
+            return {
+                    {1.f / (aspect_ratio * fov_half_tan), 0, 0, 0},
+                    {0, 1.f / (fov_half_tan), 0, 0},
+                    {0, 0, far / (far - near), -(near * far) / (far - near)},
+                    {0, 0, 1, 0},
+            };
+        if (ndc_depth_range == NDCDepthRange::NEGATIVE_ONE_TO_ONE)
+            return {
+                    {1.f / (aspect_ratio * fov_half_tan), 0, 0, 0},
+                    {0, 1.f / (fov_half_tan), 0, 0},
+                    {0, 0, (far + near) / (far - near), -(2.f * far * near) / (far - near)},
+                    {0, 0, 1, 0},
+            };
+        std::unreachable();
     }
 } // namespace omath::source_engine

source/engines/source_engine/traits/camera_trait.cpp

@@ -20,8 +20,9 @@ namespace omath::source_engine
     }
     Mat4X4 CameraTrait::calc_projection_matrix(const projection::FieldOfView& fov,
                                                const projection::ViewPort& view_port, const float near,
-                                               const float far) noexcept
+                                               const float far, const NDCDepthRange ndc_depth_range) noexcept
     {
-        return calc_perspective_projection_matrix(fov.as_degrees(), view_port.aspect_ratio(), near, far);
+        return calc_perspective_projection_matrix(fov.as_degrees(), view_port.aspect_ratio(), near, far,
+                                                  ndc_depth_range);
     }
 } // namespace omath::source_engine

source/engines/unity_engine/formulas.cpp

@@ -35,8 +35,15 @@ namespace omath::unity_engine
                * mat_rotation_axis_x<float, MatStoreType::ROW_MAJOR>(angles.pitch);
     }
     Mat4X4 calc_perspective_projection_matrix(const float field_of_view, const float aspect_ratio, const float near,
-                                              const float far) noexcept
+                                              const float far, const NDCDepthRange ndc_depth_range) noexcept
     {
-        return omath::mat_perspective_right_handed(field_of_view, aspect_ratio, near, far);
+        if (ndc_depth_range == NDCDepthRange::ZERO_TO_ONE)
+            return omath::mat_perspective_right_handed<float, MatStoreType::ROW_MAJOR, NDCDepthRange::ZERO_TO_ONE>(
+                    field_of_view, aspect_ratio, near, far);
+        if (ndc_depth_range == NDCDepthRange::NEGATIVE_ONE_TO_ONE)
+            return omath::mat_perspective_right_handed<float, MatStoreType::ROW_MAJOR,
+                                                       NDCDepthRange::NEGATIVE_ONE_TO_ONE>(field_of_view, aspect_ratio,
+                                                                                           near, far);
+        std::unreachable();
     }
 } // namespace omath::unity_engine

source/engines/unity_engine/traits/camera_trait.cpp

@@ -19,8 +19,9 @@ namespace omath::unity_engine
     }
     Mat4X4 CameraTrait::calc_projection_matrix(const projection::FieldOfView& fov,
                                                const projection::ViewPort& view_port, const float near,
-                                               const float far) noexcept
+                                               const float far, const NDCDepthRange ndc_depth_range) noexcept
     {
-        return calc_perspective_projection_matrix(fov.as_degrees(), view_port.aspect_ratio(), near, far);
+        return calc_perspective_projection_matrix(fov.as_degrees(), view_port.aspect_ratio(), near, far,
+                                                  ndc_depth_range);
     }
 } // namespace omath::unity_engine

source/engines/unreal_engine/formulas.cpp

@@ -35,8 +35,12 @@ namespace omath::unreal_engine
                * mat_rotation_axis_y<float, MatStoreType::ROW_MAJOR>(angles.pitch);
     }
     Mat4X4 calc_perspective_projection_matrix(const float field_of_view, const float aspect_ratio, const float near,
-                                              const float far) noexcept
+                                              const float far, const NDCDepthRange ndc_depth_range) noexcept
     {
+        if (ndc_depth_range == NDCDepthRange::ZERO_TO_ONE)
+            return mat_perspective_left_handed<float, MatStoreType::ROW_MAJOR, NDCDepthRange::ZERO_TO_ONE>(
+                    field_of_view, aspect_ratio, near, far);
+
         return mat_perspective_left_handed(field_of_view, aspect_ratio, near, far);
     }
 } // namespace omath::unreal_engine

source/engines/unreal_engine/traits/camera_trait.cpp

@@ -19,8 +19,9 @@ namespace omath::unreal_engine
     }
     Mat4X4 CameraTrait::calc_projection_matrix(const projection::FieldOfView& fov,
                                                const projection::ViewPort& view_port, const float near,
-                                               const float far) noexcept
+                                               const float far, const NDCDepthRange ndc_depth_range) noexcept
     {
-        return calc_perspective_projection_matrix(fov.as_degrees(), view_port.aspect_ratio(), near, far);
+        return calc_perspective_projection_matrix(fov.as_degrees(), view_port.aspect_ratio(), near, far,
+                                                  ndc_depth_range);
     }
 } // namespace omath::unreal_engine

tests/engines/unit_test_traits_engines.cpp

@@ -21,6 +21,9 @@
 #include <omath/engines/unreal_engine/traits/camera_trait.hpp>
 
 #include <omath/engines/source_engine/traits/pred_engine_trait.hpp>
+#include <omath/engines/source_engine/traits/camera_trait.hpp>
+
+#include <omath/engines/cry_engine/traits/camera_trait.hpp>
 
 #include <omath/projectile_prediction/projectile.hpp>
 #include <omath/projectile_prediction/target.hpp>
@@ -218,9 +221,14 @@ TEST(TraitTests, Frostbite_Pred_And_Mesh_And_Camera)
     // CameraTrait look at should be callable
     const auto angles = e::CameraTrait::calc_look_at_angle({0, 0, 0}, {0, 1, 1});
     (void)angles;
-    const auto proj = e::CameraTrait::calc_projection_matrix(projection::FieldOfView::from_degrees(60.f), {1280.f, 720.f}, 0.1f, 1000.f);
+    const auto proj = e::CameraTrait::calc_projection_matrix(projection::FieldOfView::from_degrees(60.f), {1280.f, 720.f}, 0.1f, 1000.f, NDCDepthRange::NEGATIVE_ONE_TO_ONE);
     const auto expected = e::calc_perspective_projection_matrix(60.f, 1280.f / 720.f, 0.1f, 1000.f);
     expect_matrix_near(proj, expected);
+
+    const auto proj_zo = e::CameraTrait::calc_projection_matrix(projection::FieldOfView::from_degrees(60.f), {1280.f, 720.f}, 0.1f, 1000.f, NDCDepthRange::ZERO_TO_ONE);
+    const auto expected_zo = e::calc_perspective_projection_matrix(60.f, 1280.f / 720.f, 0.1f, 1000.f, NDCDepthRange::ZERO_TO_ONE);
+    expect_matrix_near(proj_zo, expected_zo);
+    EXPECT_NE(proj, proj_zo);
 }
 
 TEST(TraitTests, IW_Pred_And_Mesh_And_Camera)
@@ -264,10 +272,15 @@ TEST(TraitTests, IW_Pred_And_Mesh_And_Camera)
     e::ViewAngles va;
     expect_matrix_near(e::MeshTrait::rotation_matrix(va), e::rotation_matrix(va));
 
-    const auto proj = e::CameraTrait::calc_projection_matrix(projection::FieldOfView::from_degrees(45.f), {1920.f, 1080.f}, 0.1f, 1000.f);
+    const auto proj = e::CameraTrait::calc_projection_matrix(projection::FieldOfView::from_degrees(45.f), {1920.f, 1080.f}, 0.1f, 1000.f, NDCDepthRange::NEGATIVE_ONE_TO_ONE);
     const auto expected = e::calc_perspective_projection_matrix(45.f, 1920.f / 1080.f, 0.1f, 1000.f);
     expect_matrix_near(proj, expected);
 
+    const auto proj_zo = e::CameraTrait::calc_projection_matrix(projection::FieldOfView::from_degrees(45.f), {1920.f, 1080.f}, 0.1f, 1000.f, NDCDepthRange::ZERO_TO_ONE);
+    const auto expected_zo = e::calc_perspective_projection_matrix(45.f, 1920.f / 1080.f, 0.1f, 1000.f, NDCDepthRange::ZERO_TO_ONE);
+    expect_matrix_near(proj_zo, expected_zo);
+    EXPECT_NE(proj, proj_zo);
+
     // non-airborne
     t.m_is_airborne = false;
     const auto pred_ground_iw = e::PredEngineTrait::predict_target_position(t, 2.f, 9.81f);
@@ -314,10 +327,15 @@ TEST(TraitTests, OpenGL_Pred_And_Mesh_And_Camera)
     e::ViewAngles va;
     expect_matrix_near(e::MeshTrait::rotation_matrix(va), e::rotation_matrix(va));
 
-    const auto proj = e::CameraTrait::calc_projection_matrix(projection::FieldOfView::from_degrees(60.f), {1280.f, 720.f}, 0.1f, 1000.f);
+    const auto proj = e::CameraTrait::calc_projection_matrix(projection::FieldOfView::from_degrees(60.f), {1280.f, 720.f}, 0.1f, 1000.f, NDCDepthRange::NEGATIVE_ONE_TO_ONE);
     const auto expected = e::calc_perspective_projection_matrix(60.f, 1280.f / 720.f, 0.1f, 1000.f);
     expect_matrix_near(proj, expected);
 
+    const auto proj_zo = e::CameraTrait::calc_projection_matrix(projection::FieldOfView::from_degrees(60.f), {1280.f, 720.f}, 0.1f, 1000.f, NDCDepthRange::ZERO_TO_ONE);
+    const auto expected_zo = e::calc_perspective_projection_matrix(60.f, 1280.f / 720.f, 0.1f, 1000.f, NDCDepthRange::ZERO_TO_ONE);
+    expect_matrix_near(proj_zo, expected_zo);
+    EXPECT_NE(proj, proj_zo);
+
     // non-airborne
     t.m_is_airborne = false;
     const auto pred_ground_gl = e::PredEngineTrait::predict_target_position(t, 2.f, 9.81f);
@@ -364,10 +382,15 @@ TEST(TraitTests, Unity_Pred_And_Mesh_And_Camera)
     e::ViewAngles va;
     expect_matrix_near(e::MeshTrait::rotation_matrix(va), e::rotation_matrix(va));
 
-    const auto proj = e::CameraTrait::calc_projection_matrix(projection::FieldOfView::from_degrees(60.f), {1280.f, 720.f}, 0.1f, 1000.f);
+    const auto proj = e::CameraTrait::calc_projection_matrix(projection::FieldOfView::from_degrees(60.f), {1280.f, 720.f}, 0.1f, 1000.f, NDCDepthRange::NEGATIVE_ONE_TO_ONE);
     const auto expected = e::calc_perspective_projection_matrix(60.f, 1280.f / 720.f, 0.1f, 1000.f);
     expect_matrix_near(proj, expected);
 
+    const auto proj_zo = e::CameraTrait::calc_projection_matrix(projection::FieldOfView::from_degrees(60.f), {1280.f, 720.f}, 0.1f, 1000.f, NDCDepthRange::ZERO_TO_ONE);
+    const auto expected_zo = e::calc_perspective_projection_matrix(60.f, 1280.f / 720.f, 0.1f, 1000.f, NDCDepthRange::ZERO_TO_ONE);
+    expect_matrix_near(proj_zo, expected_zo);
+    EXPECT_NE(proj, proj_zo);
+
     // non-airborne
     t.m_is_airborne = false;
     const auto pred_ground_unity = e::PredEngineTrait::predict_target_position(t, 2.f, 9.81f);
@@ -414,12 +437,237 @@ TEST(TraitTests, Unreal_Pred_And_Mesh_And_Camera)
     e::ViewAngles va;
     expect_matrix_near(e::MeshTrait::rotation_matrix(va), e::rotation_matrix(va));
 
-    const auto proj = e::CameraTrait::calc_projection_matrix(projection::FieldOfView::from_degrees(60.f), {1280.f, 720.f}, 0.1f, 1000.f);
+    const auto proj = e::CameraTrait::calc_projection_matrix(projection::FieldOfView::from_degrees(60.f), {1280.f, 720.f}, 0.1f, 1000.f, NDCDepthRange::NEGATIVE_ONE_TO_ONE);
     const auto expected = e::calc_perspective_projection_matrix(60.f, 1280.f / 720.f, 0.1f, 1000.f);
     expect_matrix_near(proj, expected);
 
+    const auto proj_zo = e::CameraTrait::calc_projection_matrix(projection::FieldOfView::from_degrees(60.f), {1280.f, 720.f}, 0.1f, 1000.f, NDCDepthRange::ZERO_TO_ONE);
+    const auto expected_zo = e::calc_perspective_projection_matrix(60.f, 1280.f / 720.f, 0.1f, 1000.f, NDCDepthRange::ZERO_TO_ONE);
+    expect_matrix_near(proj_zo, expected_zo);
+    EXPECT_NE(proj, proj_zo);
+
     // non-airborne
     t.m_is_airborne = false;
     const auto pred_ground_unreal = e::PredEngineTrait::predict_target_position(t, 2.f, 9.81f);
     EXPECT_NEAR(pred_ground_unreal.x, 4.f, 1e-6f);
 }
+
+// ── NDC Depth Range tests for Source and CryEngine camera traits ────────────
+
+TEST(NDCDepthRangeTests, Source_BothDepthRanges)
+{
+    namespace e = omath::source_engine;
+
+    const auto proj_no = e::CameraTrait::calc_projection_matrix(
+            projection::FieldOfView::from_degrees(90.f), {1920.f, 1080.f}, 0.1f, 1000.f,
+            NDCDepthRange::NEGATIVE_ONE_TO_ONE);
+    const auto expected_no = e::calc_perspective_projection_matrix(
+            90.f, 1920.f / 1080.f, 0.1f, 1000.f, NDCDepthRange::NEGATIVE_ONE_TO_ONE);
+    expect_matrix_near(proj_no, expected_no);
+
+    const auto proj_zo = e::CameraTrait::calc_projection_matrix(
+            projection::FieldOfView::from_degrees(90.f), {1920.f, 1080.f}, 0.1f, 1000.f,
+            NDCDepthRange::ZERO_TO_ONE);
+    const auto expected_zo = e::calc_perspective_projection_matrix(
+            90.f, 1920.f / 1080.f, 0.1f, 1000.f, NDCDepthRange::ZERO_TO_ONE);
+    expect_matrix_near(proj_zo, expected_zo);
+
+    EXPECT_NE(proj_no, proj_zo);
+}
+
+TEST(NDCDepthRangeTests, CryEngine_BothDepthRanges)
+{
+    namespace e = omath::cry_engine;
+
+    const auto proj_no = e::CameraTrait::calc_projection_matrix(
+            projection::FieldOfView::from_degrees(90.f), {1920.f, 1080.f}, 0.1f, 1000.f,
+            NDCDepthRange::NEGATIVE_ONE_TO_ONE);
+    const auto expected_no = e::calc_perspective_projection_matrix(
+            90.f, 1920.f / 1080.f, 0.1f, 1000.f, NDCDepthRange::NEGATIVE_ONE_TO_ONE);
+    expect_matrix_near(proj_no, expected_no);
+
+    const auto proj_zo = e::CameraTrait::calc_projection_matrix(
+            projection::FieldOfView::from_degrees(90.f), {1920.f, 1080.f}, 0.1f, 1000.f,
+            NDCDepthRange::ZERO_TO_ONE);
+    const auto expected_zo = e::calc_perspective_projection_matrix(
+            90.f, 1920.f / 1080.f, 0.1f, 1000.f, NDCDepthRange::ZERO_TO_ONE);
+    expect_matrix_near(proj_zo, expected_zo);
+
+    EXPECT_NE(proj_no, proj_zo);
+}
+
+// ── Verify Z mapping for ZERO_TO_ONE across all engines ─────────────────────
+
+// Helper: projects a point at given z through a left-handed projection matrix and returns NDC z
+static float project_z_lh(const Mat<4, 4>& proj, float z)
+{
+    auto clip = proj * mat_column_from_vector<float>({0, 0, z});
+    return clip.at(2, 0) / clip.at(3, 0);
+}
+
+TEST(NDCDepthRangeTests, Source_ZeroToOne_ZRange)
+{
+    namespace e = omath::source_engine;
+    // Source is left-handed
+    const auto proj = e::calc_perspective_projection_matrix(90.f, 16.f / 9.f, 0.1f, 1000.f, NDCDepthRange::ZERO_TO_ONE);
+
+    EXPECT_NEAR(project_z_lh(proj, 0.1f), 0.0f, 1e-4f);
+    EXPECT_NEAR(project_z_lh(proj, 1000.f), 1.0f, 1e-4f);
+    EXPECT_GT(project_z_lh(proj, 500.f), 0.0f);
+    EXPECT_LT(project_z_lh(proj, 500.f), 1.0f);
+}
+
+TEST(NDCDepthRangeTests, IW_ZeroToOne_ZRange)
+{
+    namespace e = omath::iw_engine;
+    // IW is left-handed
+    const auto proj = e::calc_perspective_projection_matrix(90.f, 16.f / 9.f, 0.1f, 1000.f, NDCDepthRange::ZERO_TO_ONE);
+
+    EXPECT_NEAR(project_z_lh(proj, 0.1f), 0.0f, 1e-4f);
+    EXPECT_NEAR(project_z_lh(proj, 1000.f), 1.0f, 1e-4f);
+    EXPECT_GT(project_z_lh(proj, 500.f), 0.0f);
+    EXPECT_LT(project_z_lh(proj, 500.f), 1.0f);
+}
+
+TEST(NDCDepthRangeTests, OpenGL_ZeroToOne_ZRange)
+{
+    namespace e = omath::opengl_engine;
+    // OpenGL is right-handed (negative z forward), column-major
+    const auto proj = e::calc_perspective_projection_matrix(90.f, 16.f / 9.f, 0.1f, 1000.f, NDCDepthRange::ZERO_TO_ONE);
+
+    // OpenGL engine uses column-major matrices, project manually
+    auto proj_z = [&](float z) {
+        auto clip = proj * mat_column_from_vector<float, MatStoreType::COLUMN_MAJOR>({0, 0, z});
+        return clip.at(2, 0) / clip.at(3, 0);
+    };
+
+    EXPECT_NEAR(proj_z(-0.1f), 0.0f, 1e-4f);
+    EXPECT_NEAR(proj_z(-1000.f), 1.0f, 1e-4f);
+    EXPECT_GT(proj_z(-500.f), 0.0f);
+    EXPECT_LT(proj_z(-500.f), 1.0f);
+}
+
+TEST(NDCDepthRangeTests, Frostbite_ZeroToOne_ZRange)
+{
+    namespace e = omath::frostbite_engine;
+    const auto proj = e::calc_perspective_projection_matrix(90.f, 16.f / 9.f, 0.1f, 1000.f, NDCDepthRange::ZERO_TO_ONE);
+
+    EXPECT_NEAR(project_z_lh(proj, 0.1f), 0.0f, 1e-4f);
+    EXPECT_NEAR(project_z_lh(proj, 1000.f), 1.0f, 1e-4f);
+    EXPECT_GT(project_z_lh(proj, 500.f), 0.0f);
+    EXPECT_LT(project_z_lh(proj, 500.f), 1.0f);
+}
+
+TEST(NDCDepthRangeTests, Unity_ZeroToOne_ZRange)
+{
+    namespace e = omath::unity_engine;
+    // Unity is right-handed, row-major
+    const auto proj = e::calc_perspective_projection_matrix(90.f, 16.f / 9.f, 0.1f, 1000.f, NDCDepthRange::ZERO_TO_ONE);
+
+    auto proj_z = [&](float z) {
+        auto clip = proj * mat_column_from_vector<float>({0, 0, z});
+        return clip.at(2, 0) / clip.at(3, 0);
+    };
+
+    EXPECT_NEAR(proj_z(-0.1f), 0.0f, 1e-4f);
+    EXPECT_NEAR(proj_z(-1000.f), 1.0f, 1e-4f);
+    EXPECT_GT(proj_z(-500.f), 0.0f);
+    EXPECT_LT(proj_z(-500.f), 1.0f);
+}
+
+TEST(NDCDepthRangeTests, Unreal_ZeroToOne_ZRange)
+{
+    namespace e = omath::unreal_engine;
+    const auto proj = e::calc_perspective_projection_matrix(90.f, 16.f / 9.f, 0.1f, 1000.f, NDCDepthRange::ZERO_TO_ONE);
+
+    EXPECT_NEAR(project_z_lh(proj, 0.1f), 0.0f, 1e-4f);
+    EXPECT_NEAR(project_z_lh(proj, 1000.f), 1.0f, 1e-4f);
+    EXPECT_GT(project_z_lh(proj, 500.f), 0.0f);
+    EXPECT_LT(project_z_lh(proj, 500.f), 1.0f);
+}
+
+TEST(NDCDepthRangeTests, CryEngine_ZeroToOne_ZRange)
+{
+    namespace e = omath::cry_engine;
+    const auto proj = e::calc_perspective_projection_matrix(90.f, 16.f / 9.f, 0.1f, 1000.f, NDCDepthRange::ZERO_TO_ONE);
+
+    EXPECT_NEAR(project_z_lh(proj, 0.1f), 0.0f, 1e-4f);
+    EXPECT_NEAR(project_z_lh(proj, 1000.f), 1.0f, 1e-4f);
+    EXPECT_GT(project_z_lh(proj, 500.f), 0.0f);
+    EXPECT_LT(project_z_lh(proj, 500.f), 1.0f);
+}
+
+// ── Verify Z mapping for NEGATIVE_ONE_TO_ONE across all engines ─────────────
+
+TEST(NDCDepthRangeTests, Source_NegativeOneToOne_ZRange)
+{
+    namespace e = omath::source_engine;
+    const auto proj = e::calc_perspective_projection_matrix(90.f, 16.f / 9.f, 0.1f, 1000.f, NDCDepthRange::NEGATIVE_ONE_TO_ONE);
+
+    EXPECT_NEAR(project_z_lh(proj, 0.1f), -1.0f, 1e-3f);
+    EXPECT_NEAR(project_z_lh(proj, 1000.f), 1.0f, 1e-3f);
+}
+
+TEST(NDCDepthRangeTests, IW_NegativeOneToOne_ZRange)
+{
+    namespace e = omath::iw_engine;
+    const auto proj = e::calc_perspective_projection_matrix(90.f, 16.f / 9.f, 0.1f, 1000.f, NDCDepthRange::NEGATIVE_ONE_TO_ONE);
+
+    EXPECT_NEAR(project_z_lh(proj, 0.1f), -1.0f, 1e-3f);
+    EXPECT_NEAR(project_z_lh(proj, 1000.f), 1.0f, 1e-3f);
+}
+
+TEST(NDCDepthRangeTests, Frostbite_NegativeOneToOne_ZRange)
+{
+    namespace e = omath::frostbite_engine;
+    const auto proj = e::calc_perspective_projection_matrix(90.f, 16.f / 9.f, 0.1f, 1000.f, NDCDepthRange::NEGATIVE_ONE_TO_ONE);
+
+    EXPECT_NEAR(project_z_lh(proj, 0.1f), -1.0f, 1e-3f);
+    EXPECT_NEAR(project_z_lh(proj, 1000.f), 1.0f, 1e-3f);
+}
+
+TEST(NDCDepthRangeTests, Unreal_NegativeOneToOne_ZRange)
+{
+    namespace e = omath::unreal_engine;
+    const auto proj = e::calc_perspective_projection_matrix(90.f, 16.f / 9.f, 0.1f, 1000.f, NDCDepthRange::NEGATIVE_ONE_TO_ONE);
+
+    EXPECT_NEAR(project_z_lh(proj, 0.1f), -1.0f, 1e-3f);
+    EXPECT_NEAR(project_z_lh(proj, 1000.f), 1.0f, 1e-3f);
+}
+
+TEST(NDCDepthRangeTests, CryEngine_NegativeOneToOne_ZRange)
+{
+    namespace e = omath::cry_engine;
+    const auto proj = e::calc_perspective_projection_matrix(90.f, 16.f / 9.f, 0.1f, 1000.f, NDCDepthRange::NEGATIVE_ONE_TO_ONE);
+
+    EXPECT_NEAR(project_z_lh(proj, 0.1f), -1.0f, 1e-3f);
+    EXPECT_NEAR(project_z_lh(proj, 1000.f), 1.0f, 1e-3f);
+}
+
+TEST(NDCDepthRangeTests, OpenGL_NegativeOneToOne_ZRange)
+{
+    namespace e = omath::opengl_engine;
+    const auto proj = e::calc_perspective_projection_matrix(90.f, 16.f / 9.f, 0.1f, 1000.f, NDCDepthRange::NEGATIVE_ONE_TO_ONE);
+
+    auto proj_z = [&](float z) {
+        auto clip = proj * mat_column_from_vector<float, MatStoreType::COLUMN_MAJOR>({0, 0, z});
+        return clip.at(2, 0) / clip.at(3, 0);
+    };
+
+    EXPECT_NEAR(proj_z(-0.1f), -1.0f, 1e-3f);
+    EXPECT_NEAR(proj_z(-1000.f), 1.0f, 1e-3f);
+}
+
+TEST(NDCDepthRangeTests, Unity_NegativeOneToOne_ZRange)
+{
+    namespace e = omath::unity_engine;
+    const auto proj = e::calc_perspective_projection_matrix(90.f, 16.f / 9.f, 0.1f, 1000.f, NDCDepthRange::NEGATIVE_ONE_TO_ONE);
+
+    auto proj_z = [&](float z) {
+        auto clip = proj * mat_column_from_vector<float>({0, 0, z});
+        return clip.at(2, 0) / clip.at(3, 0);
+    };
+
+    EXPECT_NEAR(proj_z(-0.1f), -1.0f, 1e-3f);
+    EXPECT_NEAR(proj_z(-1000.f), 1.0f, 1e-3f);
+}

tests/general/unit_test_bvh_tree.cpp

@@ -0,0 +1,876 @@
+//
+// Created by Orange on 04/08/2026.
+//
+#include <gtest/gtest.h>
+#include <omath/collision/bvh_tree.hpp>
+#include <algorithm>
+#include <random>
+#include <set>
+
+using Aabb = omath::primitives::Aabb<float>;
+using BvhTree = omath::collision::BvhTree<float>;
+using Ray = omath::collision::Ray<>;
+using HitResult = BvhTree::HitResult;
+
+using AabbD = omath::primitives::Aabb<double>;
+using BvhTreeD = omath::collision::BvhTree<double>;
+
+// ============================================================================
+// Helper: brute-force overlap query for verification
+// ============================================================================
+static std::set<std::size_t> brute_force_overlaps(const std::vector<Aabb>& aabbs, const Aabb& query)
+{
+    std::set<std::size_t> result;
+    for (std::size_t i = 0; i < aabbs.size(); ++i)
+    {
+        if (query.min.x <= aabbs[i].max.x && query.max.x >= aabbs[i].min.x
+            && query.min.y <= aabbs[i].max.y && query.max.y >= aabbs[i].min.y
+            && query.min.z <= aabbs[i].max.z && query.max.z >= aabbs[i].min.z)
+            result.insert(i);
+    }
+    return result;
+}
+
+// ============================================================================
+// Construction tests
+// ============================================================================
+
+TEST(UnitTestBvhTree, EmptyTree)
+{
+    const BvhTree tree;
+    EXPECT_TRUE(tree.empty());
+    EXPECT_EQ(tree.node_count(), 0);
+    EXPECT_TRUE(tree.query_overlaps({}).empty());
+}
+
+TEST(UnitTestBvhTree, EmptySpan)
+{
+    const std::vector<Aabb> empty;
+    const BvhTree tree(empty);
+    EXPECT_TRUE(tree.empty());
+    EXPECT_EQ(tree.node_count(), 0);
+}
+
+TEST(UnitTestBvhTree, SingleElement)
+{
+    const std::vector<Aabb> aabbs = {
+        {{0.f, 0.f, 0.f}, {1.f, 1.f, 1.f}}
+    };
+
+    const BvhTree tree(aabbs);
+    EXPECT_FALSE(tree.empty());
+    EXPECT_EQ(tree.node_count(), 1);
+
+    const auto results = tree.query_overlaps({{0.5f, 0.5f, 0.5f}, {1.5f, 1.5f, 1.5f}});
+    ASSERT_EQ(results.size(), 1);
+    EXPECT_EQ(results[0], 0);
+
+    const auto miss = tree.query_overlaps({{5.f, 5.f, 5.f}, {6.f, 6.f, 6.f}});
+    EXPECT_TRUE(miss.empty());
+}
+
+TEST(UnitTestBvhTree, TwoElements)
+{
+    const std::vector<Aabb> aabbs = {
+        {{0.f, 0.f, 0.f}, {1.f, 1.f, 1.f}},
+        {{5.f, 5.f, 5.f}, {6.f, 6.f, 6.f}},
+    };
+
+    const BvhTree tree(aabbs);
+    EXPECT_FALSE(tree.empty());
+
+    // Hit first only
+    auto r = tree.query_overlaps({{-0.5f, -0.5f, -0.5f}, {0.5f, 0.5f, 0.5f}});
+    ASSERT_EQ(r.size(), 1);
+    EXPECT_EQ(r[0], 0);
+
+    // Hit second only
+    r = tree.query_overlaps({{5.5f, 5.5f, 5.5f}, {7.f, 7.f, 7.f}});
+    ASSERT_EQ(r.size(), 1);
+    EXPECT_EQ(r[0], 1);
+
+    // Hit both
+    r = tree.query_overlaps({{-1.f, -1.f, -1.f}, {10.f, 10.f, 10.f}});
+    EXPECT_EQ(r.size(), 2);
+}
+
+TEST(UnitTestBvhTree, ThreeElements)
+{
+    const std::vector<Aabb> aabbs = {
+        {{0.f, 0.f, 0.f}, {1.f, 1.f, 1.f}},
+        {{2.f, 2.f, 2.f}, {3.f, 3.f, 3.f}},
+        {{10.f, 10.f, 10.f}, {11.f, 11.f, 11.f}},
+    };
+
+    const BvhTree tree(aabbs);
+
+    const auto results = tree.query_overlaps({{0.5f, 0.5f, 0.5f}, {2.5f, 2.5f, 2.5f}});
+    EXPECT_EQ(results.size(), 2);
+
+    const auto far = tree.query_overlaps({{9.5f, 9.5f, 9.5f}, {10.5f, 10.5f, 10.5f}});
+    ASSERT_EQ(far.size(), 1);
+    EXPECT_EQ(far[0], 2);
+}
+
+TEST(UnitTestBvhTree, NodeCountGrowsSublinearly)
+{
+    // For N objects, node count should be at most 2N-1
+    std::vector<Aabb> aabbs;
+    for (int i = 0; i < 100; ++i)
+    {
+        const auto f = static_cast<float>(i) * 3.f;
+        aabbs.push_back({{f, 0.f, 0.f}, {f + 1.f, 1.f, 1.f}});
+    }
+    const BvhTree tree(aabbs);
+    EXPECT_LE(tree.node_count(), 2 * aabbs.size());
+}
+
+// ============================================================================
+// Overlap query tests
+// ============================================================================
+
+TEST(UnitTestBvhTree, OverlapExactTouch)
+{
+    // Two boxes share exactly one face
+    const std::vector<Aabb> aabbs = {
+        {{0.f, 0.f, 0.f}, {1.f, 1.f, 1.f}},
+        {{1.f, 0.f, 0.f}, {2.f, 1.f, 1.f}},
+    };
+
+    const BvhTree tree(aabbs);
+
+    // Query exactly at the shared face — should overlap both
+    const auto r = tree.query_overlaps({{0.5f, 0.f, 0.f}, {1.5f, 1.f, 1.f}});
+    EXPECT_EQ(r.size(), 2);
+}
+
+TEST(UnitTestBvhTree, OverlapEdgeTouch)
+{
+    // Query AABB edge-touches an object AABB
+    const std::vector<Aabb> aabbs = {
+        {{0.f, 0.f, 0.f}, {1.f, 1.f, 1.f}},
+    };
+
+    const BvhTree tree(aabbs);
+
+    // Touching at corner point (1,1,1)
+    const auto r = tree.query_overlaps({{1.f, 1.f, 1.f}, {2.f, 2.f, 2.f}});
+    EXPECT_EQ(r.size(), 1);
+}
+
+TEST(UnitTestBvhTree, OverlapQueryInsideObject)
+{
+    // Query is fully inside an object
+    const std::vector<Aabb> aabbs = {
+        {{-10.f, -10.f, -10.f}, {10.f, 10.f, 10.f}},
+    };
+
+    const BvhTree tree(aabbs);
+    const auto r = tree.query_overlaps({{-1.f, -1.f, -1.f}, {1.f, 1.f, 1.f}});
+    ASSERT_EQ(r.size(), 1);
+    EXPECT_EQ(r[0], 0);
+}
+
+TEST(UnitTestBvhTree, OverlapObjectInsideQuery)
+{
+    // Object is fully inside the query
+    const std::vector<Aabb> aabbs = {
+        {{4.f, 4.f, 4.f}, {5.f, 5.f, 5.f}},
+    };
+
+    const BvhTree tree(aabbs);
+    const auto r = tree.query_overlaps({{0.f, 0.f, 0.f}, {10.f, 10.f, 10.f}});
+    ASSERT_EQ(r.size(), 1);
+    EXPECT_EQ(r[0], 0);
+}
+
+TEST(UnitTestBvhTree, OverlapMissOnSingleAxis)
+{
+    // Overlap on X and Y but not Z
+    const std::vector<Aabb> aabbs = {
+        {{0.f, 0.f, 0.f}, {1.f, 1.f, 1.f}},
+    };
+
+    const BvhTree tree(aabbs);
+
+    const auto r = tree.query_overlaps({{0.f, 0.f, 5.f}, {1.f, 1.f, 6.f}});
+    EXPECT_TRUE(r.empty());
+}
+
+TEST(UnitTestBvhTree, OverlapNegativeCoordinates)
+{
+    const std::vector<Aabb> aabbs = {
+        {{-5.f, -5.f, -5.f}, {-3.f, -3.f, -3.f}},
+        {{-2.f, -2.f, -2.f}, {0.f, 0.f, 0.f}},
+        {{1.f, 1.f, 1.f}, {3.f, 3.f, 3.f}},
+    };
+
+    const BvhTree tree(aabbs);
+
+    const auto r = tree.query_overlaps({{-6.f, -6.f, -6.f}, {-4.f, -4.f, -4.f}});
+    ASSERT_EQ(r.size(), 1);
+    EXPECT_EQ(r[0], 0);
+}
+
+TEST(UnitTestBvhTree, OverlapMixedNegativePositive)
+{
+    const std::vector<Aabb> aabbs = {
+        {{-1.f, -1.f, -1.f}, {1.f, 1.f, 1.f}},
+    };
+
+    const BvhTree tree(aabbs);
+
+    // Query spans negative and positive
+    const auto r = tree.query_overlaps({{-0.5f, -0.5f, -0.5f}, {0.5f, 0.5f, 0.5f}});
+    ASSERT_EQ(r.size(), 1);
+}
+
+TEST(UnitTestBvhTree, OverlapNoHitsAmongMany)
+{
+    std::vector<Aabb> aabbs;
+    for (int i = 0; i < 50; ++i)
+    {
+        const auto f = static_cast<float>(i) * 5.f;
+        aabbs.push_back({{f, 0.f, 0.f}, {f + 1.f, 1.f, 1.f}});
+    }
+
+    const BvhTree tree(aabbs);
+
+    // Query far from all objects
+    const auto r = tree.query_overlaps({{-100.f, -100.f, -100.f}, {-90.f, -90.f, -90.f}});
+    EXPECT_TRUE(r.empty());
+}
+
+TEST(UnitTestBvhTree, OverlapAllObjects)
+{
+    std::vector<Aabb> aabbs;
+    for (int i = 0; i < 64; ++i)
+    {
+        const auto f = static_cast<float>(i);
+        aabbs.push_back({{f, f, f}, {f + 0.5f, f + 0.5f, f + 0.5f}});
+    }
+
+    const BvhTree tree(aabbs);
+
+    const auto r = tree.query_overlaps({{-1.f, -1.f, -1.f}, {100.f, 100.f, 100.f}});
+    EXPECT_EQ(r.size(), 64);
+}
+
+TEST(UnitTestBvhTree, OverlapReturnsCorrectIndices)
+{
+    // Specific index verification
+    const std::vector<Aabb> aabbs = {
+        {{0.f, 0.f, 0.f}, {1.f, 1.f, 1.f}},     // 0
+        {{10.f, 0.f, 0.f}, {11.f, 1.f, 1.f}},    // 1
+        {{20.f, 0.f, 0.f}, {21.f, 1.f, 1.f}},    // 2
+        {{30.f, 0.f, 0.f}, {31.f, 1.f, 1.f}},    // 3
+        {{40.f, 0.f, 0.f}, {41.f, 1.f, 1.f}},    // 4
+    };
+
+    const BvhTree tree(aabbs);
+
+    // Hit only index 2
+    auto r = tree.query_overlaps({{19.5f, -1.f, -1.f}, {20.5f, 2.f, 2.f}});
+    ASSERT_EQ(r.size(), 1);
+    EXPECT_EQ(r[0], 2);
+
+    // Hit only index 4
+    r = tree.query_overlaps({{39.5f, -1.f, -1.f}, {40.5f, 2.f, 2.f}});
+    ASSERT_EQ(r.size(), 1);
+    EXPECT_EQ(r[0], 4);
+}
+
+// ============================================================================
+// Spatial distribution tests
+// ============================================================================
+
+TEST(UnitTestBvhTree, ObjectsAlongXAxis)
+{
+    // All objects on a line along X
+    std::vector<Aabb> aabbs;
+    for (int i = 0; i < 20; ++i)
+    {
+        const auto f = static_cast<float>(i) * 4.f;
+        aabbs.push_back({{f, 0.f, 0.f}, {f + 1.f, 1.f, 1.f}});
+    }
+
+    const BvhTree tree(aabbs);
+
+    const auto r = tree.query_overlaps({{0.f, 0.f, 0.f}, {1.f, 1.f, 1.f}});
+    EXPECT_EQ(r.size(), 1);
+
+    const auto mid = tree.query_overlaps({{7.5f, -1.f, -1.f}, {8.5f, 2.f, 2.f}});
+    EXPECT_EQ(mid.size(), 1);
+}
+
+TEST(UnitTestBvhTree, ObjectsAlongYAxis)
+{
+    std::vector<Aabb> aabbs;
+    for (int i = 0; i < 20; ++i)
+    {
+        const auto f = static_cast<float>(i) * 4.f;
+        aabbs.push_back({{0.f, f, 0.f}, {1.f, f + 1.f, 1.f}});
+    }
+
+    const BvhTree tree(aabbs);
+
+    const auto r = tree.query_overlaps({{-1.f, 38.f, -1.f}, {2.f, 40.f, 2.f}});
+    EXPECT_EQ(r.size(), 1);
+}
+
+TEST(UnitTestBvhTree, ObjectsAlongZAxis)
+{
+    std::vector<Aabb> aabbs;
+    for (int i = 0; i < 20; ++i)
+    {
+        const auto f = static_cast<float>(i) * 4.f;
+        aabbs.push_back({{0.f, 0.f, f}, {1.f, 1.f, f + 1.f}});
+    }
+
+    const BvhTree tree(aabbs);
+
+    const auto r = tree.query_overlaps({{-1.f, -1.f, 38.f}, {2.f, 2.f, 40.f}});
+    EXPECT_EQ(r.size(), 1);
+}
+
+TEST(UnitTestBvhTree, ObjectsInPlaneXY)
+{
+    // Grid in the XY plane
+    std::vector<Aabb> aabbs;
+    for (int x = 0; x < 10; ++x)
+        for (int y = 0; y < 10; ++y)
+        {
+            const auto fx = static_cast<float>(x) * 3.f;
+            const auto fy = static_cast<float>(y) * 3.f;
+            aabbs.push_back({{fx, fy, 0.f}, {fx + 1.f, fy + 1.f, 1.f}});
+        }
+
+    const BvhTree tree(aabbs);
+    EXPECT_EQ(tree.query_overlaps({{-1.f, -1.f, -1.f}, {100.f, 100.f, 2.f}}).size(), 100);
+
+    // Single cell query
+    const auto r = tree.query_overlaps({{0.f, 0.f, 0.f}, {0.5f, 0.5f, 0.5f}});
+    EXPECT_EQ(r.size(), 1);
+}
+
+TEST(UnitTestBvhTree, ClusteredObjects)
+{
+    // Two clusters far apart
+    std::vector<Aabb> aabbs;
+    for (int i = 0; i < 25; ++i)
+    {
+        const auto f = static_cast<float>(i) * 0.5f;
+        aabbs.push_back({{f, f, f}, {f + 0.4f, f + 0.4f, f + 0.4f}});
+    }
+    for (int i = 0; i < 25; ++i)
+    {
+        const auto f = 100.f + static_cast<float>(i) * 0.5f;
+        aabbs.push_back({{f, f, f}, {f + 0.4f, f + 0.4f, f + 0.4f}});
+    }
+
+    const BvhTree tree(aabbs);
+
+    // Query near first cluster
+    const auto r1 = tree.query_overlaps({{-1.f, -1.f, -1.f}, {15.f, 15.f, 15.f}});
+    EXPECT_EQ(r1.size(), 25);
+
+    // Query near second cluster
+    const auto r2 = tree.query_overlaps({{99.f, 99.f, 99.f}, {115.f, 115.f, 115.f}});
+    EXPECT_EQ(r2.size(), 25);
+
+    // Query between clusters — should find nothing
+    const auto gap = tree.query_overlaps({{50.f, 50.f, 50.f}, {60.f, 60.f, 60.f}});
+    EXPECT_TRUE(gap.empty());
+}
+
+TEST(UnitTestBvhTree, OverlappingObjects)
+{
+    // Objects that overlap each other
+    const std::vector<Aabb> aabbs = {
+        {{0.f, 0.f, 0.f}, {2.f, 2.f, 2.f}},
+        {{1.f, 1.f, 1.f}, {3.f, 3.f, 3.f}},
+        {{1.5f, 1.5f, 1.5f}, {4.f, 4.f, 4.f}},
+    };
+
+    const BvhTree tree(aabbs);
+
+    // Query at the overlap region of all three
+    const auto r = tree.query_overlaps({{1.5f, 1.5f, 1.5f}, {2.f, 2.f, 2.f}});
+    EXPECT_EQ(r.size(), 3);
+}
+
+TEST(UnitTestBvhTree, IdenticalObjects)
+{
+    // All objects at the same position
+    std::vector<Aabb> aabbs;
+    for (int i = 0; i < 10; ++i)
+        aabbs.push_back({{0.f, 0.f, 0.f}, {1.f, 1.f, 1.f}});
+
+    const BvhTree tree(aabbs);
+
+    const auto r = tree.query_overlaps({{0.f, 0.f, 0.f}, {1.f, 1.f, 1.f}});
+    EXPECT_EQ(r.size(), 10);
+}
+
+TEST(UnitTestBvhTree, DegenerateThickPlanes)
+{
+    // Very flat AABBs (thickness ~0 in one axis)
+    const std::vector<Aabb> aabbs = {
+        {{0.f, 0.f, 0.f}, {10.f, 10.f, 0.001f}},
+        {{0.f, 0.f, 5.f}, {10.f, 10.f, 5.001f}},
+    };
+
+    const BvhTree tree(aabbs);
+
+    const auto r = tree.query_overlaps({{0.f, 0.f, -0.01f}, {10.f, 10.f, 0.01f}});
+    ASSERT_EQ(r.size(), 1);
+}
+
+TEST(UnitTestBvhTree, VaryingSizes)
+{
+    // Objects of wildly different sizes
+    const std::vector<Aabb> aabbs = {
+        {{0.f, 0.f, 0.f}, {0.01f, 0.01f, 0.01f}},         // tiny
+        {{-500.f, -500.f, -500.f}, {500.f, 500.f, 500.f}}, // huge
+        {{10.f, 10.f, 10.f}, {11.f, 11.f, 11.f}},          // normal
+    };
+
+    const BvhTree tree(aabbs);
+
+    // The huge box should overlap almost any query
+    auto r = tree.query_overlaps({{200.f, 200.f, 200.f}, {201.f, 201.f, 201.f}});
+    ASSERT_EQ(r.size(), 1);
+    EXPECT_EQ(r[0], 1);
+
+    // Query at origin hits the tiny and the huge
+    r = tree.query_overlaps({{-0.1f, -0.1f, -0.1f}, {0.1f, 0.1f, 0.1f}});
+    EXPECT_EQ(r.size(), 2);
+}
+
+// ============================================================================
+// Ray query tests
+// ============================================================================
+
+TEST(UnitTestBvhTree, RayQueryBasic)
+{
+    const std::vector<Aabb> aabbs = {
+        {{0.f, 0.f, 0.f}, {1.f, 1.f, 1.f}},
+        {{5.f, 0.f, 0.f}, {6.f, 1.f, 1.f}},
+        {{0.f, 5.f, 0.f}, {1.f, 6.f, 1.f}},
+    };
+
+    const BvhTree tree(aabbs);
+
+    Ray ray;
+    ray.start = {-1.f, 0.5f, 0.5f};
+    ray.end = {10.f, 0.5f, 0.5f};
+    ray.infinite_length = true;
+
+    const auto hits = tree.query_ray(ray);
+    EXPECT_GE(hits.size(), 2);
+
+    if (hits.size() >= 2)
+        EXPECT_LE(hits[0].distance_sqr, hits[1].distance_sqr);
+}
+
+TEST(UnitTestBvhTree, RayQueryMissesAll)
+{
+    const std::vector<Aabb> aabbs = {
+        {{0.f, 0.f, 0.f}, {1.f, 1.f, 1.f}},
+        {{5.f, 0.f, 0.f}, {6.f, 1.f, 1.f}},
+    };
+
+    const BvhTree tree(aabbs);
+
+    // Ray above everything
+    Ray ray;
+    ray.start = {-1.f, 100.f, 0.5f};
+    ray.end = {10.f, 100.f, 0.5f};
+    ray.infinite_length = true;
+
+    const auto hits = tree.query_ray(ray);
+    EXPECT_TRUE(hits.empty());
+}
+
+TEST(UnitTestBvhTree, RayQueryAlongY)
+{
+    const std::vector<Aabb> aabbs = {
+        {{0.f, 0.f, 0.f}, {1.f, 1.f, 1.f}},
+        {{0.f, 5.f, 0.f}, {1.f, 6.f, 1.f}},
+        {{0.f, 10.f, 0.f}, {1.f, 11.f, 1.f}},
+    };
+
+    const BvhTree tree(aabbs);
+
+    Ray ray;
+    ray.start = {0.5f, -1.f, 0.5f};
+    ray.end = {0.5f, 20.f, 0.5f};
+    ray.infinite_length = true;
+
+    const auto hits = tree.query_ray(ray);
+    EXPECT_EQ(hits.size(), 3);
+
+    // Verify sorted by distance
+    for (std::size_t i = 1; i < hits.size(); ++i)
+        EXPECT_LE(hits[i - 1].distance_sqr, hits[i].distance_sqr);
+}
+
+TEST(UnitTestBvhTree, RayQueryAlongZ)
+{
+    const std::vector<Aabb> aabbs = {
+        {{0.f, 0.f, 0.f}, {1.f, 1.f, 1.f}},
+        {{0.f, 0.f, 10.f}, {1.f, 1.f, 11.f}},
+    };
+
+    const BvhTree tree(aabbs);
+
+    Ray ray;
+    ray.start = {0.5f, 0.5f, -5.f};
+    ray.end = {0.5f, 0.5f, 20.f};
+    ray.infinite_length = true;
+
+    const auto hits = tree.query_ray(ray);
+    EXPECT_EQ(hits.size(), 2);
+}
+
+TEST(UnitTestBvhTree, RayQueryDiagonal)
+{
+    const std::vector<Aabb> aabbs = {
+        {{0.f, 0.f, 0.f}, {1.f, 1.f, 1.f}},
+        {{5.f, 5.f, 5.f}, {6.f, 6.f, 6.f}},
+        {{10.f, 10.f, 10.f}, {11.f, 11.f, 11.f}},
+    };
+
+    const BvhTree tree(aabbs);
+
+    // Diagonal ray through all three
+    Ray ray;
+    ray.start = {-1.f, -1.f, -1.f};
+    ray.end = {15.f, 15.f, 15.f};
+    ray.infinite_length = true;
+
+    const auto hits = tree.query_ray(ray);
+    EXPECT_EQ(hits.size(), 3);
+}
+
+TEST(UnitTestBvhTree, RayQueryOnEmptyTree)
+{
+    const BvhTree tree;
+
+    Ray ray;
+    ray.start = {0.f, 0.f, 0.f};
+    ray.end = {10.f, 0.f, 0.f};
+    ray.infinite_length = true;
+
+    const auto hits = tree.query_ray(ray);
+    EXPECT_TRUE(hits.empty());
+}
+
+TEST(UnitTestBvhTree, RayQuerySortedByDistance)
+{
+    // Many boxes along a line
+    std::vector<Aabb> aabbs;
+    for (int i = 0; i < 20; ++i)
+    {
+        const auto f = static_cast<float>(i) * 3.f;
+        aabbs.push_back({{f, 0.f, 0.f}, {f + 1.f, 1.f, 1.f}});
+    }
+
+    const BvhTree tree(aabbs);
+
+    Ray ray;
+    ray.start = {-1.f, 0.5f, 0.5f};
+    ray.end = {100.f, 0.5f, 0.5f};
+    ray.infinite_length = true;
+
+    const auto hits = tree.query_ray(ray);
+    EXPECT_EQ(hits.size(), 20);
+
+    for (std::size_t i = 1; i < hits.size(); ++i)
+        EXPECT_LE(hits[i - 1].distance_sqr, hits[i].distance_sqr);
+}
+
+// ============================================================================
+// Brute-force verification tests
+// ============================================================================
+
+TEST(UnitTestBvhTree, BruteForceVerificationGrid)
+{
+    std::vector<Aabb> aabbs;
+    for (int x = 0; x < 10; ++x)
+        for (int y = 0; y < 10; ++y)
+            for (int z = 0; z < 10; ++z)
+            {
+                const auto fx = static_cast<float>(x) * 3.f;
+                const auto fy = static_cast<float>(y) * 3.f;
+                const auto fz = static_cast<float>(z) * 3.f;
+                aabbs.push_back({{fx, fy, fz}, {fx + 1.f, fy + 1.f, fz + 1.f}});
+            }
+
+    const BvhTree tree(aabbs);
+
+    // Test several queries and compare to brute force
+    const std::vector<Aabb> queries = {
+        {{0.f, 0.f, 0.f}, {1.5f, 1.5f, 1.5f}},
+        {{-1.f, -1.f, -1.f}, {100.f, 100.f, 100.f}},
+        {{13.f, 13.f, 13.f}, {14.f, 14.f, 14.f}},
+        {{-50.f, -50.f, -50.f}, {-40.f, -40.f, -40.f}},
+        {{5.5f, 5.5f, 5.5f}, {7.5f, 7.5f, 7.5f}},
+    };
+
+    for (const auto& q : queries)
+    {
+        const auto bvh_results = tree.query_overlaps(q);
+        const auto brute_results = brute_force_overlaps(aabbs, q);
+
+        const std::set<std::size_t> bvh_set(bvh_results.begin(), bvh_results.end());
+        EXPECT_EQ(bvh_set, brute_results)
+            << "Mismatch for query [(" << q.min.x << "," << q.min.y << "," << q.min.z
+            << ") -> (" << q.max.x << "," << q.max.y << "," << q.max.z << ")]";
+    }
+}
+
+TEST(UnitTestBvhTree, BruteForceVerificationRandom)
+{
+    std::mt19937 rng(42);
+    std::uniform_real_distribution<float> pos_dist(-50.f, 50.f);
+    std::uniform_real_distribution<float> size_dist(0.5f, 3.f);
+
+    std::vector<Aabb> aabbs;
+    for (int i = 0; i < 200; ++i)
+    {
+        const auto x = pos_dist(rng);
+        const auto y = pos_dist(rng);
+        const auto z = pos_dist(rng);
+        const auto sx = size_dist(rng);
+        const auto sy = size_dist(rng);
+        const auto sz = size_dist(rng);
+        aabbs.push_back({{x, y, z}, {x + sx, y + sy, z + sz}});
+    }
+
+    const BvhTree tree(aabbs);
+
+    // Run 50 random queries
+    for (int i = 0; i < 50; ++i)
+    {
+        const auto qx = pos_dist(rng);
+        const auto qy = pos_dist(rng);
+        const auto qz = pos_dist(rng);
+        const auto qsx = size_dist(rng);
+        const auto qsy = size_dist(rng);
+        const auto qsz = size_dist(rng);
+        const Aabb query = {{qx, qy, qz}, {qx + qsx, qy + qsy, qz + qsz}};
+
+        const auto bvh_results = tree.query_overlaps(query);
+        const auto brute_results = brute_force_overlaps(aabbs, query);
+
+        const std::set<std::size_t> bvh_set(bvh_results.begin(), bvh_results.end());
+        EXPECT_EQ(bvh_set, brute_results) << "Mismatch on random query iteration " << i;
+    }
+}
+
+// ============================================================================
+// Large dataset tests
+// ============================================================================
+
+TEST(UnitTestBvhTree, LargeGridDataset)
+{
+    std::vector<Aabb> aabbs;
+    for (int x = 0; x < 10; ++x)
+        for (int y = 0; y < 10; ++y)
+            for (int z = 0; z < 10; ++z)
+            {
+                const auto fx = static_cast<float>(x) * 3.f;
+                const auto fy = static_cast<float>(y) * 3.f;
+                const auto fz = static_cast<float>(z) * 3.f;
+                aabbs.push_back({{fx, fy, fz}, {fx + 1.f, fy + 1.f, fz + 1.f}});
+            }
+
+    const BvhTree tree(aabbs);
+    EXPECT_FALSE(tree.empty());
+
+    const auto results = tree.query_overlaps({{0.f, 0.f, 0.f}, {1.5f, 1.5f, 1.5f}});
+    EXPECT_EQ(results.size(), 1);
+
+    const auto all_results = tree.query_overlaps({{-1.f, -1.f, -1.f}, {100.f, 100.f, 100.f}});
+    EXPECT_EQ(all_results.size(), 1000);
+}
+
+TEST(UnitTestBvhTree, FiveThousandObjects)
+{
+    std::vector<Aabb> aabbs;
+    for (int i = 0; i < 5000; ++i)
+    {
+        const auto f = static_cast<float>(i) * 2.f;
+        aabbs.push_back({{f, 0.f, 0.f}, {f + 1.f, 1.f, 1.f}});
+    }
+
+    const BvhTree tree(aabbs);
+    EXPECT_FALSE(tree.empty());
+
+    // Query that should hit exactly 1
+    const auto r = tree.query_overlaps({{0.f, 0.f, 0.f}, {0.5f, 0.5f, 0.5f}});
+    EXPECT_EQ(r.size(), 1);
+
+    // Query that misses
+    const auto miss = tree.query_overlaps({{-100.f, -100.f, -100.f}, {-90.f, -90.f, -90.f}});
+    EXPECT_TRUE(miss.empty());
+}
+
+// ============================================================================
+// Double precision tests
+// ============================================================================
+
+TEST(UnitTestBvhTree, DoublePrecision)
+{
+    const std::vector<AabbD> aabbs = {
+        {{0.0, 0.0, 0.0}, {1.0, 1.0, 1.0}},
+        {{5.0, 5.0, 5.0}, {6.0, 6.0, 6.0}},
+        {{10.0, 10.0, 10.0}, {11.0, 11.0, 11.0}},
+    };
+
+    const BvhTreeD tree(aabbs);
+    EXPECT_FALSE(tree.empty());
+
+    const auto r = tree.query_overlaps({{0.5, 0.5, 0.5}, {1.5, 1.5, 1.5}});
+    ASSERT_EQ(r.size(), 1);
+    EXPECT_EQ(r[0], 0);
+
+    const auto r2 = tree.query_overlaps({{4.5, 4.5, 4.5}, {5.5, 5.5, 5.5}});
+    ASSERT_EQ(r2.size(), 1);
+    EXPECT_EQ(r2[0], 1);
+}
+
+TEST(UnitTestBvhTree, DoublePrecisionLargeCoordinates)
+{
+    const std::vector<AabbD> aabbs = {
+        {{1e10, 1e10, 1e10}, {1e10 + 1.0, 1e10 + 1.0, 1e10 + 1.0}},
+        {{-1e10, -1e10, -1e10}, {-1e10 + 1.0, -1e10 + 1.0, -1e10 + 1.0}},
+    };
+
+    const BvhTreeD tree(aabbs);
+
+    const auto r = tree.query_overlaps({{1e10 - 0.5, 1e10 - 0.5, 1e10 - 0.5},
+                                         {1e10 + 0.5, 1e10 + 0.5, 1e10 + 0.5}});
+    ASSERT_EQ(r.size(), 1);
+    EXPECT_EQ(r[0], 0);
+}
+
+// ============================================================================
+// Edge case tests
+// ============================================================================
+
+TEST(UnitTestBvhTree, ZeroSizeQuery)
+{
+    const std::vector<Aabb> aabbs = {
+        {{0.f, 0.f, 0.f}, {1.f, 1.f, 1.f}},
+    };
+
+    const BvhTree tree(aabbs);
+
+    // Point query inside the box
+    const auto r = tree.query_overlaps({{0.5f, 0.5f, 0.5f}, {0.5f, 0.5f, 0.5f}});
+    EXPECT_EQ(r.size(), 1);
+
+    // Point query outside the box
+    const auto miss = tree.query_overlaps({{5.f, 5.f, 5.f}, {5.f, 5.f, 5.f}});
+    EXPECT_TRUE(miss.empty());
+}
+
+TEST(UnitTestBvhTree, ZeroSizeObjects)
+{
+    // Point-like AABBs
+    const std::vector<Aabb> aabbs = {
+        {{1.f, 1.f, 1.f}, {1.f, 1.f, 1.f}},
+        {{5.f, 5.f, 5.f}, {5.f, 5.f, 5.f}},
+    };
+
+    const BvhTree tree(aabbs);
+
+    const auto r = tree.query_overlaps({{0.f, 0.f, 0.f}, {2.f, 2.f, 2.f}});
+    ASSERT_EQ(r.size(), 1);
+    EXPECT_EQ(r[0], 0);
+}
+
+TEST(UnitTestBvhTree, NoDuplicateResults)
+{
+    std::vector<Aabb> aabbs;
+    for (int i = 0; i < 50; ++i)
+    {
+        const auto f = static_cast<float>(i) * 2.f;
+        aabbs.push_back({{f, 0.f, 0.f}, {f + 1.f, 1.f, 1.f}});
+    }
+
+    const BvhTree tree(aabbs);
+
+    const auto r = tree.query_overlaps({{-1.f, -1.f, -1.f}, {200.f, 2.f, 2.f}});
+
+    // Check for duplicates
+    const std::set<std::size_t> unique_results(r.begin(), r.end());
+    EXPECT_EQ(unique_results.size(), r.size());
+    EXPECT_EQ(r.size(), 50);
+}
+
+TEST(UnitTestBvhTree, LargeSpread)
+{
+    // Objects with huge gaps between them
+    const std::vector<Aabb> aabbs = {
+        {{0.f, 0.f, 0.f}, {1.f, 1.f, 1.f}},
+        {{1000.f, 0.f, 0.f}, {1001.f, 1.f, 1.f}},
+        {{-1000.f, 0.f, 0.f}, {-999.f, 1.f, 1.f}},
+        {{0.f, 1000.f, 0.f}, {1.f, 1001.f, 1.f}},
+        {{0.f, -1000.f, 0.f}, {1.f, -999.f, 1.f}},
+    };
+
+    const BvhTree tree(aabbs);
+
+    auto r = tree.query_overlaps({{999.f, -1.f, -1.f}, {1002.f, 2.f, 2.f}});
+    ASSERT_EQ(r.size(), 1);
+    EXPECT_EQ(r[0], 1);
+
+    r = tree.query_overlaps({{-1001.f, -1.f, -1.f}, {-998.f, 2.f, 2.f}});
+    ASSERT_EQ(r.size(), 1);
+    EXPECT_EQ(r[0], 2);
+}
+
+TEST(UnitTestBvhTree, AllObjectsSameCenter)
+{
+    // All AABBs centered at origin but different sizes
+    std::vector<Aabb> aabbs;
+    for (int i = 1; i <= 10; ++i)
+    {
+        const auto s = static_cast<float>(i);
+        aabbs.push_back({{-s, -s, -s}, {s, s, s}});
+    }
+
+    const BvhTree tree(aabbs);
+
+    // Small query at origin should hit all
+    const auto r = tree.query_overlaps({{-0.1f, -0.1f, -0.1f}, {0.1f, 0.1f, 0.1f}});
+    EXPECT_EQ(r.size(), 10);
+
+    // Query touching only the largest box
+    const auto r2 = tree.query_overlaps({{9.5f, 9.5f, 9.5f}, {10.5f, 10.5f, 10.5f}});
+    ASSERT_EQ(r2.size(), 1);
+    EXPECT_EQ(r2[0], 9);
+}
+
+TEST(UnitTestBvhTree, MultipleQueriesSameTree)
+{
+    std::vector<Aabb> aabbs;
+    for (int i = 0; i < 100; ++i)
+    {
+        const auto f = static_cast<float>(i) * 2.f;
+        aabbs.push_back({{f, 0.f, 0.f}, {f + 1.f, 1.f, 1.f}});
+    }
+
+    const BvhTree tree(aabbs);
+
+    // Run many queries on the same tree
+    for (int i = 0; i < 100; ++i)
+    {
+        const auto f = static_cast<float>(i) * 2.f;
+        const auto r = tree.query_overlaps({{f, 0.f, 0.f}, {f + 1.f, 1.f, 1.f}});
+        ASSERT_GE(r.size(), 1) << "Query for object " << i << " should find at least itself";
+    }
+}

tests/general/unit_test_line_tracer_aabb.cpp

@@ -0,0 +1,125 @@
+//
+// Created by Vlad on 3/25/2025.
+//
+#include "omath/collision/line_tracer.hpp"
+#include "omath/3d_primitives/aabb.hpp"
+#include <gtest/gtest.h>
+
+using Vec3 = omath::Vector3<float>;
+using Ray = omath::collision::Ray<>;
+using LineTracer = omath::collision::LineTracer<>;
+using AABB = omath::primitives::Aabb<float>;
+
+static Ray make_ray(Vec3 start, Vec3 end, bool infinite = false)
+{
+    Ray r;
+    r.start = start;
+    r.end = end;
+    r.infinite_length = infinite;
+    return r;
+}
+
+// Ray passing straight through the center along Z axis
+TEST(LineTracerAABBTests, HitCenterAlongZ)
+{
+    const AABB box{{-1.f, -1.f, -1.f}, {1.f, 1.f, 1.f}};
+    const auto ray = make_ray({0.f, 0.f, -5.f}, {0.f, 0.f, 5.f});
+
+    const auto hit = LineTracer::get_ray_hit_point(ray, box);
+    EXPECT_NE(hit, ray.end);
+    EXPECT_NEAR(hit.z, -1.f, 1e-4f);
+    EXPECT_NEAR(hit.x, 0.f, 1e-4f);
+    EXPECT_NEAR(hit.y, 0.f, 1e-4f);
+}
+
+// Ray passing straight through the center along X axis
+TEST(LineTracerAABBTests, HitCenterAlongX)
+{
+    const AABB box{{-1.f, -1.f, -1.f}, {1.f, 1.f, 1.f}};
+    const auto ray = make_ray({-5.f, 0.f, 0.f}, {5.f, 0.f, 0.f});
+
+    const auto hit = LineTracer::get_ray_hit_point(ray, box);
+    EXPECT_NE(hit, ray.end);
+    EXPECT_NEAR(hit.x, -1.f, 1e-4f);
+}
+
+// Ray that misses entirely (too far in Y)
+TEST(LineTracerAABBTests, MissReturnsEnd)
+{
+    const AABB box{{-1.f, -1.f, -1.f}, {1.f, 1.f, 1.f}};
+    const auto ray = make_ray({0.f, 5.f, -5.f}, {0.f, 5.f, 5.f});
+
+    const auto hit = LineTracer::get_ray_hit_point(ray, box);
+    EXPECT_EQ(hit, ray.end);
+}
+
+// Ray that stops short before reaching the box
+TEST(LineTracerAABBTests, RayTooShortReturnsEnd)
+{
+    const AABB box{{3.f, -1.f, -1.f}, {5.f, 1.f, 1.f}};
+    const auto ray = make_ray({0.f, 0.f, 0.f}, {2.f, 0.f, 0.f});
+
+    const auto hit = LineTracer::get_ray_hit_point(ray, box);
+    EXPECT_EQ(hit, ray.end);
+}
+
+// Infinite ray that starts before the box should hit
+TEST(LineTracerAABBTests, InfiniteRayHits)
+{
+    const AABB box{{3.f, -1.f, -1.f}, {5.f, 1.f, 1.f}};
+    const auto ray = make_ray({0.f, 0.f, 0.f}, {2.f, 0.f, 0.f}, true);
+
+    const auto hit = LineTracer::get_ray_hit_point(ray, box);
+    EXPECT_NE(hit, ray.end);
+    EXPECT_NEAR(hit.x, 3.f, 1e-4f);
+}
+
+// Ray starting inside the box — t_min=0, so hit point equals ray.start
+TEST(LineTracerAABBTests, RayStartsInsideBox)
+{
+    const AABB box{{-1.f, -1.f, -1.f}, {1.f, 1.f, 1.f}};
+    const auto ray = make_ray({0.f, 0.f, 0.f}, {0.f, 0.f, 5.f});
+
+    const auto hit = LineTracer::get_ray_hit_point(ray, box);
+    EXPECT_NE(hit, ray.end);
+    // t_min is clamped to 0, so hit == start
+    EXPECT_NEAR(hit.x, 0.f, 1e-4f);
+    EXPECT_NEAR(hit.y, 0.f, 1e-4f);
+    EXPECT_NEAR(hit.z, 0.f, 1e-4f);
+}
+
+// Ray parallel to XY plane, pointing along X, at Z outside the box
+TEST(LineTracerAABBTests, ParallelRayOutsideSlabMisses)
+{
+    const AABB box{{-1.f, -1.f, -1.f}, {1.f, 1.f, 1.f}};
+    // Z component of ray is 3.0 — outside box's Z slab
+    const auto ray = make_ray({-5.f, 0.f, 3.f}, {5.f, 0.f, 3.f});
+
+    const auto hit = LineTracer::get_ray_hit_point(ray, box);
+    EXPECT_EQ(hit, ray.end);
+}
+
+// Ray parallel to XY plane, pointing along X, at Z inside the box
+TEST(LineTracerAABBTests, ParallelRayInsideSlabHits)
+{
+    const AABB box{{-1.f, -1.f, -1.f}, {1.f, 1.f, 1.f}};
+    const auto ray = make_ray({-5.f, 0.f, 0.f}, {5.f, 0.f, 0.f});
+
+    const auto hit = LineTracer::get_ray_hit_point(ray, box);
+    EXPECT_NE(hit, ray.end);
+    EXPECT_NEAR(hit.x, -1.f, 1e-4f);
+}
+
+// Diagonal ray hitting a corner region
+TEST(LineTracerAABBTests, DiagonalRayHits)
+{
+    const AABB box{{0.f, 0.f, 0.f}, {2.f, 2.f, 2.f}};
+    const auto ray = make_ray({-1.f, -1.f, -1.f}, {3.f, 3.f, 3.f});
+
+    const auto hit = LineTracer::get_ray_hit_point(ray, box);
+    EXPECT_NE(hit, ray.end);
+    // Entry point should be at (0,0,0)
+    EXPECT_NEAR(hit.x, 0.f, 1e-4f);
+    EXPECT_NEAR(hit.y, 0.f, 1e-4f);
+    EXPECT_NEAR(hit.z, 0.f, 1e-4f);
+}

tests/general/unit_test_mat.cpp

@@ -240,4 +240,126 @@ TEST(UnitTestMatStandalone, MatPerspectiveLeftHanded)
     projected /= projected.at(3, 0);
 
     EXPECT_TRUE(projected.at(2, 0) > -1.0f && projected.at(2, 0) < 0.f);
+}
+
+TEST(UnitTestMatStandalone, MatPerspectiveLeftHandedZeroToOne)
+{
+    const auto proj = mat_perspective_left_handed<float, MatStoreType::ROW_MAJOR, NDCDepthRange::ZERO_TO_ONE>(
+            90.f, 16.f / 9.f, 0.1f, 1000.f);
+
+    // Near plane point should map to z ~ 0
+    auto near_pt = proj * mat_column_from_vector<float>({0, 0, 0.1f});
+    near_pt /= near_pt.at(3, 0);
+    EXPECT_NEAR(near_pt.at(2, 0), 0.0f, 1e-4f);
+
+    // Far plane point should map to z ~ 1
+    auto far_pt = proj * mat_column_from_vector<float>({0, 0, 1000.f});
+    far_pt /= far_pt.at(3, 0);
+    EXPECT_NEAR(far_pt.at(2, 0), 1.0f, 1e-4f);
+
+    // Mid-range point should be in [0, 1]
+    auto mid_pt = proj * mat_column_from_vector<float>({0, 0, 500.f});
+    mid_pt /= mid_pt.at(3, 0);
+    EXPECT_GT(mid_pt.at(2, 0), 0.0f);
+    EXPECT_LT(mid_pt.at(2, 0), 1.0f);
+}
+
+TEST(UnitTestMatStandalone, MatPerspectiveRightHandedZeroToOne)
+{
+    const auto proj = mat_perspective_right_handed<float, MatStoreType::ROW_MAJOR, NDCDepthRange::ZERO_TO_ONE>(
+            90.f, 16.f / 9.f, 0.1f, 1000.f);
+
+    // Near plane point (negative z for right-handed) should map to z ~ 0
+    auto near_pt = proj * mat_column_from_vector<float>({0, 0, -0.1f});
+    near_pt /= near_pt.at(3, 0);
+    EXPECT_NEAR(near_pt.at(2, 0), 0.0f, 1e-4f);
+
+    // Far plane point should map to z ~ 1
+    auto far_pt = proj * mat_column_from_vector<float>({0, 0, -1000.f});
+    far_pt /= far_pt.at(3, 0);
+    EXPECT_NEAR(far_pt.at(2, 0), 1.0f, 1e-4f);
+
+    // Mid-range point should be in [0, 1]
+    auto mid_pt = proj * mat_column_from_vector<float>({0, 0, -500.f});
+    mid_pt /= mid_pt.at(3, 0);
+    EXPECT_GT(mid_pt.at(2, 0), 0.0f);
+    EXPECT_LT(mid_pt.at(2, 0), 1.0f);
+}
+
+TEST(UnitTestMatStandalone, MatPerspectiveNegativeOneToOneRange)
+{
+    // Verify existing [-1, 1] behavior with explicit template arg matches default
+    const auto proj_default = mat_perspective_left_handed(90.f, 16.f / 9.f, 0.1f, 1000.f);
+    const auto proj_explicit = mat_perspective_left_handed<float, MatStoreType::ROW_MAJOR,
+            NDCDepthRange::NEGATIVE_ONE_TO_ONE>(90.f, 16.f / 9.f, 0.1f, 1000.f);
+
+    EXPECT_EQ(proj_default, proj_explicit);
+
+    // Near plane should map to z ~ -1
+    auto near_pt = proj_default * mat_column_from_vector<float>({0, 0, 0.1f});
+    near_pt /= near_pt.at(3, 0);
+    EXPECT_NEAR(near_pt.at(2, 0), -1.0f, 1e-3f);
+
+    // Far plane should map to z ~ 1
+    auto far_pt = proj_default * mat_column_from_vector<float>({0, 0, 1000.f});
+    far_pt /= far_pt.at(3, 0);
+    EXPECT_NEAR(far_pt.at(2, 0), 1.0f, 1e-3f);
+}
+
+TEST(UnitTestMatStandalone, MatPerspectiveZeroToOneEquanity)
+{
+    // LH and RH should produce same NDC for mirrored z
+    constexpr omath::Vector3<float> left_handed = {0, 2, 10};
+    constexpr omath::Vector3<float> right_handed = {0, 2, -10};
+
+    const auto proj_lh = mat_perspective_left_handed<float, MatStoreType::ROW_MAJOR, NDCDepthRange::ZERO_TO_ONE>(
+            90.f, 16.f / 9.f, 0.1f, 1000.f);
+    const auto proj_rh = mat_perspective_right_handed<float, MatStoreType::ROW_MAJOR, NDCDepthRange::ZERO_TO_ONE>(
+            90.f, 16.f / 9.f, 0.1f, 1000.f);
+
+    auto ndc_lh = proj_lh * mat_column_from_vector(left_handed);
+    auto ndc_rh = proj_rh * mat_column_from_vector(right_handed);
+
+    ndc_lh /= ndc_lh.at(3, 0);
+    ndc_rh /= ndc_rh.at(3, 0);
+
+    EXPECT_EQ(ndc_lh, ndc_rh);
+}
+
+TEST(UnitTestMatStandalone, MatOrthoLeftHandedZeroToOne)
+{
+    const auto ortho = mat_ortho_left_handed<float, MatStoreType::ROW_MAJOR, NDCDepthRange::ZERO_TO_ONE>(
+            -1.f, 1.f, -1.f, 1.f, 0.1f, 100.f);
+
+    // Near plane should map to z ~ 0
+    auto near_pt = ortho * mat_column_from_vector<float>({0, 0, 0.1f});
+    EXPECT_NEAR(near_pt.at(2, 0), 0.0f, 1e-4f);
+
+    // Far plane should map to z ~ 1
+    auto far_pt = ortho * mat_column_from_vector<float>({0, 0, 100.f});
+    EXPECT_NEAR(far_pt.at(2, 0), 1.0f, 1e-4f);
+}
+
+TEST(UnitTestMatStandalone, MatOrthoRightHandedZeroToOne)
+{
+    const auto ortho = mat_ortho_right_handed<float, MatStoreType::ROW_MAJOR, NDCDepthRange::ZERO_TO_ONE>(
+            -1.f, 1.f, -1.f, 1.f, 0.1f, 100.f);
+
+    // Near plane (negative z for RH) should map to z ~ 0
+    auto near_pt = ortho * mat_column_from_vector<float>({0, 0, -0.1f});
+    EXPECT_NEAR(near_pt.at(2, 0), 0.0f, 1e-4f);
+
+    // Far plane should map to z ~ 1
+    auto far_pt = ortho * mat_column_from_vector<float>({0, 0, -100.f});
+    EXPECT_NEAR(far_pt.at(2, 0), 1.0f, 1e-4f);
+}
+
+TEST(UnitTestMatStandalone, MatOrthoNegativeOneToOneDefault)
+{
+    // Verify explicit [-1, 1] matches default
+    const auto ortho_default = mat_ortho_left_handed(-1.f, 1.f, -1.f, 1.f, 0.1f, 100.f);
+    const auto ortho_explicit = mat_ortho_left_handed<float, MatStoreType::ROW_MAJOR,
+            NDCDepthRange::NEGATIVE_ONE_TO_ONE>(-1.f, 1.f, -1.f, 1.f, 0.1f, 100.f);
+
+    EXPECT_EQ(ortho_default, ortho_explicit);
 }

tests/general/unit_test_projection.cpp

@@ -4,6 +4,8 @@
 #include "omath/engines/unity_engine/camera.hpp"
 #include <complex>
 #include <gtest/gtest.h>
+#include <omath/3d_primitives/aabb.hpp>
+#include <omath/engines/opengl_engine/camera.hpp>
 #include <omath/engines/source_engine/camera.hpp>
 #include <omath/projection/camera.hpp>
 #include <print>
@@ -216,4 +218,296 @@ TEST(UnitTestProjection, ScreenToWorldBottomLeftCorner)
         EXPECT_NEAR(screen_cords->x, initial_screen_cords.x, 0.001f);
         EXPECT_NEAR(screen_cords->y, initial_screen_cords.y, 0.001f);
     }
+}
+
+TEST(UnitTestProjection, AabbInsideFrustumNotCulled)
+{
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(90.f);
+    const auto cam = omath::source_engine::Camera({0, 0, 0}, omath::source_engine::ViewAngles{}, {1920.f, 1080.f}, fov,
+                                                  0.01f, 1000.f);
+
+    // Small box directly in front of camera (Source Engine: +X forward, +Y left, +Z up)
+    const omath::primitives::Aabb<float> aabb{{90.f, -1.f, -1.f}, {110.f, 1.f, 1.f}};
+    EXPECT_FALSE(cam.is_aabb_culled_by_frustum(aabb));
+}
+
+TEST(UnitTestProjection, AabbBehindCameraCulled)
+{
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(90.f);
+    const auto cam = omath::source_engine::Camera({0, 0, 0}, omath::source_engine::ViewAngles{}, {1920.f, 1080.f}, fov,
+                                                  0.01f, 1000.f);
+
+    // Box entirely behind the camera
+    const omath::primitives::Aabb<float> aabb{{-200.f, -1.f, -1.f}, {-100.f, 1.f, 1.f}};
+    EXPECT_TRUE(cam.is_aabb_culled_by_frustum(aabb));
+}
+
+TEST(UnitTestProjection, AabbBeyondFarPlaneCulled)
+{
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(90.f);
+    const auto cam = omath::source_engine::Camera({0, 0, 0}, omath::source_engine::ViewAngles{}, {1920.f, 1080.f}, fov,
+                                                  0.01f, 1000.f);
+
+    // Box beyond far plane (1000)
+    const omath::primitives::Aabb<float> aabb{{1500.f, -1.f, -1.f}, {2000.f, 1.f, 1.f}};
+    EXPECT_TRUE(cam.is_aabb_culled_by_frustum(aabb));
+}
+
+TEST(UnitTestProjection, AabbFarLeftCulled)
+{
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(90.f);
+    const auto cam = omath::source_engine::Camera({0, 0, 0}, omath::source_engine::ViewAngles{}, {1920.f, 1080.f}, fov,
+                                                  0.01f, 1000.f);
+
+    // Box far to the side, outside the frustum
+    const omath::primitives::Aabb<float> aabb{{90.f, 4000.f, -1.f}, {110.f, 5000.f, 1.f}};
+    EXPECT_TRUE(cam.is_aabb_culled_by_frustum(aabb));
+}
+
+TEST(UnitTestProjection, AabbFarRightCulled)
+{
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(90.f);
+    const auto cam = omath::source_engine::Camera({0, 0, 0}, omath::source_engine::ViewAngles{}, {1920.f, 1080.f}, fov,
+                                                  0.01f, 1000.f);
+
+    // Box far to the other side, outside the frustum
+    const omath::primitives::Aabb<float> aabb{{90.f, -5000.f, -1.f}, {110.f, -4000.f, 1.f}};
+    EXPECT_TRUE(cam.is_aabb_culled_by_frustum(aabb));
+}
+
+TEST(UnitTestProjection, AabbAboveCulled)
+{
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(90.f);
+    const auto cam = omath::source_engine::Camera({0, 0, 0}, omath::source_engine::ViewAngles{}, {1920.f, 1080.f}, fov,
+                                                  0.01f, 1000.f);
+
+    // Box far above the frustum
+    const omath::primitives::Aabb<float> aabb{{90.f, -1.f, 5000.f}, {110.f, 1.f, 6000.f}};
+    EXPECT_TRUE(cam.is_aabb_culled_by_frustum(aabb));
+}
+
+TEST(UnitTestProjection, AabbPartiallyInsideNotCulled)
+{
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(90.f);
+    const auto cam = omath::source_engine::Camera({0, 0, 0}, omath::source_engine::ViewAngles{}, {1920.f, 1080.f}, fov,
+                                                  0.01f, 1000.f);
+
+    // Large box that straddles the frustum boundary — partially inside
+    const omath::primitives::Aabb<float> aabb{{50.f, -5000.f, -5000.f}, {500.f, 5000.f, 5000.f}};
+    EXPECT_FALSE(cam.is_aabb_culled_by_frustum(aabb));
+}
+
+TEST(UnitTestProjection, AabbStraddlesNearPlaneNotCulled)
+{
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(90.f);
+    const auto cam = omath::source_engine::Camera({0, 0, 0}, omath::source_engine::ViewAngles{}, {1920.f, 1080.f}, fov,
+                                                  0.01f, 1000.f);
+
+    // Box that straddles the near plane — partially in front
+    const omath::primitives::Aabb<float> aabb{{-5.f, -1.f, -1.f}, {5.f, 1.f, 1.f}};
+    EXPECT_FALSE(cam.is_aabb_culled_by_frustum(aabb));
+}
+
+TEST(UnitTestProjection, AabbStraddlesFarPlaneNotCulled)
+{
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(90.f);
+    const auto cam = omath::source_engine::Camera({0, 0, 0}, omath::source_engine::ViewAngles{}, {1920.f, 1080.f}, fov,
+                                                  0.01f, 1000.f);
+
+    // Box that straddles the far plane
+    const omath::primitives::Aabb<float> aabb{{900.f, -1.f, -1.f}, {1100.f, 1.f, 1.f}};
+    EXPECT_FALSE(cam.is_aabb_culled_by_frustum(aabb));
+}
+
+TEST(UnitTestProjection, AabbCulledUnityEngine)
+{
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(60.f);
+    const auto cam = omath::unity_engine::Camera({0, 0, 0}, {}, {1280.f, 720.f}, fov, 0.03f, 1000.f);
+
+    // Box in front — not culled
+    const omath::primitives::Aabb<float> inside{{-1.f, -1.f, 50.f}, {1.f, 1.f, 100.f}};
+    EXPECT_FALSE(cam.is_aabb_culled_by_frustum(inside));
+
+    // Box behind — culled
+    const omath::primitives::Aabb<float> behind{{-1.f, -1.f, -200.f}, {1.f, 1.f, -100.f}};
+    EXPECT_TRUE(cam.is_aabb_culled_by_frustum(behind));
+}
+
+TEST(UnitTestProjection, AabbBelowCulled)
+{
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(90.f);
+    const auto cam = omath::source_engine::Camera({0, 0, 0}, omath::source_engine::ViewAngles{}, {1920.f, 1080.f}, fov,
+                                                  0.01f, 1000.f);
+
+    // Box far below the frustum (Source Engine: +Z up)
+    const omath::primitives::Aabb<float> aabb{{90.f, -1.f, -6000.f}, {110.f, 1.f, -5000.f}};
+    EXPECT_TRUE(cam.is_aabb_culled_by_frustum(aabb));
+}
+
+TEST(UnitTestProjection, AabbEnclosesCameraNotCulled)
+{
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(90.f);
+    const auto cam = omath::source_engine::Camera({0, 0, 0}, omath::source_engine::ViewAngles{}, {1920.f, 1080.f}, fov,
+                                                  0.01f, 1000.f);
+
+    // Huge box that fully encloses the camera
+    const omath::primitives::Aabb<float> aabb{{-500.f, -500.f, -500.f}, {500.f, 500.f, 500.f}};
+    EXPECT_FALSE(cam.is_aabb_culled_by_frustum(aabb));
+}
+
+TEST(UnitTestProjection, AabbExactlyAtNearPlaneNotCulled)
+{
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(90.f);
+    const auto cam = omath::source_engine::Camera({0, 0, 0}, omath::source_engine::ViewAngles{}, {1920.f, 1080.f}, fov,
+                                                  0.01f, 1000.f);
+
+    // Box starting exactly at the near plane distance
+    const omath::primitives::Aabb<float> aabb{{0.01f, -1.f, -1.f}, {10.f, 1.f, 1.f}};
+    EXPECT_FALSE(cam.is_aabb_culled_by_frustum(aabb));
+}
+
+TEST(UnitTestProjection, AabbExactlyAtFarPlaneNotCulled)
+{
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(90.f);
+    const auto cam = omath::source_engine::Camera({0, 0, 0}, omath::source_engine::ViewAngles{}, {1920.f, 1080.f}, fov,
+                                                  0.01f, 1000.f);
+
+    // Box ending exactly at the far plane distance
+    const omath::primitives::Aabb<float> aabb{{990.f, -1.f, -1.f}, {1000.f, 1.f, 1.f}};
+    EXPECT_FALSE(cam.is_aabb_culled_by_frustum(aabb));
+}
+
+TEST(UnitTestProjection, AabbNarrowFovCulledAtEdge)
+{
+    // Narrow FOV — box that would be visible at 90 is culled at 30
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(30.f);
+    const auto cam = omath::source_engine::Camera({0, 0, 0}, omath::source_engine::ViewAngles{}, {1920.f, 1080.f}, fov,
+                                                  0.01f, 1000.f);
+
+    const omath::primitives::Aabb<float> aabb{{100.f, 200.f, -1.f}, {110.f, 210.f, 1.f}};
+    EXPECT_TRUE(cam.is_aabb_culled_by_frustum(aabb));
+}
+
+TEST(UnitTestProjection, AabbWideFovNotCulledAtEdge)
+{
+    // Wide FOV — same box is visible
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(120.f);
+    const auto cam = omath::source_engine::Camera({0, 0, 0}, omath::source_engine::ViewAngles{}, {1920.f, 1080.f}, fov,
+                                                  0.01f, 1000.f);
+
+    const omath::primitives::Aabb<float> aabb{{100.f, 200.f, -1.f}, {110.f, 210.f, 1.f}};
+    EXPECT_FALSE(cam.is_aabb_culled_by_frustum(aabb));
+}
+
+TEST(UnitTestProjection, AabbCameraOffOrigin)
+{
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(90.f);
+    const auto cam = omath::source_engine::Camera({500.f, 0, 0}, omath::source_engine::ViewAngles{}, {1920.f, 1080.f},
+                                                  fov, 0.01f, 1000.f);
+
+    // Box in front of shifted camera
+    const omath::primitives::Aabb<float> in_front{{600.f, -1.f, -1.f}, {700.f, 1.f, 1.f}};
+    EXPECT_FALSE(cam.is_aabb_culled_by_frustum(in_front));
+
+    // Box behind shifted camera
+    const omath::primitives::Aabb<float> behind{{0.f, -1.f, -1.f}, {100.f, 1.f, 1.f}};
+    EXPECT_TRUE(cam.is_aabb_culled_by_frustum(behind));
+}
+
+TEST(UnitTestProjection, AabbShortFarPlaneCulled)
+{
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(90.f);
+    // Very short far plane
+    const auto cam = omath::source_engine::Camera({0, 0, 0}, omath::source_engine::ViewAngles{}, {1920.f, 1080.f}, fov,
+                                                  0.01f, 50.f);
+
+    // Box at distance 100 — beyond the 50-unit far plane
+    const omath::primitives::Aabb<float> aabb{{90.f, -1.f, -1.f}, {110.f, 1.f, 1.f}};
+    EXPECT_TRUE(cam.is_aabb_culled_by_frustum(aabb));
+
+    // Box at distance 30 — within range
+    const omath::primitives::Aabb<float> near_box{{25.f, -1.f, -1.f}, {35.f, 1.f, 1.f}};
+    EXPECT_FALSE(cam.is_aabb_culled_by_frustum(near_box));
+}
+
+TEST(UnitTestProjection, AabbPointSizedInsideNotCulled)
+{
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(90.f);
+    const auto cam = omath::source_engine::Camera({0, 0, 0}, omath::source_engine::ViewAngles{}, {1920.f, 1080.f}, fov,
+                                                  0.01f, 1000.f);
+
+    // Degenerate zero-volume AABB (a point) inside the frustum
+    const omath::primitives::Aabb<float> aabb{{100.f, 0.f, 0.f}, {100.f, 0.f, 0.f}};
+    EXPECT_FALSE(cam.is_aabb_culled_by_frustum(aabb));
+}
+
+TEST(UnitTestProjection, AabbOpenGlEngineInsideNotCulled)
+{
+    // OpenGL: COLUMN_MAJOR, NEGATIVE_ONE_TO_ONE, inverted_z, forward = -Z
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(90.f);
+    const auto cam = omath::opengl_engine::Camera({0, 0, 0}, {}, {1920.f, 1080.f}, fov, 0.01f, 1000.f);
+
+    // Box in front of camera (OpenGL: -Z forward)
+    const omath::primitives::Aabb<float> aabb{{-1.f, -1.f, -110.f}, {1.f, 1.f, -90.f}};
+    EXPECT_FALSE(cam.is_aabb_culled_by_frustum(aabb));
+}
+
+TEST(UnitTestProjection, AabbOpenGlEngineBehindCulled)
+{
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(90.f);
+    const auto cam = omath::opengl_engine::Camera({0, 0, 0}, {}, {1920.f, 1080.f}, fov, 0.01f, 1000.f);
+
+    // Box behind (OpenGL: +Z is behind the camera)
+    const omath::primitives::Aabb<float> aabb{{-1.f, -1.f, 100.f}, {1.f, 1.f, 200.f}};
+    EXPECT_TRUE(cam.is_aabb_culled_by_frustum(aabb));
+}
+
+TEST(UnitTestProjection, AabbOpenGlEngineBeyondFarCulled)
+{
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(90.f);
+    const auto cam = omath::opengl_engine::Camera({0, 0, 0}, {}, {1920.f, 1080.f}, fov, 0.01f, 1000.f);
+
+    // Box beyond far plane along -Z
+    const omath::primitives::Aabb<float> aabb{{-1.f, -1.f, -2000.f}, {1.f, 1.f, -1500.f}};
+    EXPECT_TRUE(cam.is_aabb_culled_by_frustum(aabb));
+}
+
+TEST(UnitTestProjection, AabbOpenGlEngineSideCulled)
+{
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(90.f);
+    const auto cam = omath::opengl_engine::Camera({0, 0, 0}, {}, {1920.f, 1080.f}, fov, 0.01f, 1000.f);
+
+    // Box far to the right (OpenGL: +X right)
+    const omath::primitives::Aabb<float> aabb{{4000.f, -1.f, -110.f}, {5000.f, 1.f, -90.f}};
+    EXPECT_TRUE(cam.is_aabb_culled_by_frustum(aabb));
+}
+
+TEST(UnitTestProjection, AabbUnityEngineBeyondFarCulled)
+{
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(60.f);
+    const auto cam = omath::unity_engine::Camera({0, 0, 0}, {}, {1280.f, 720.f}, fov, 0.03f, 500.f);
+
+    // Box beyond 500-unit far plane (Unity: +Z forward)
+    const omath::primitives::Aabb<float> aabb{{-1.f, -1.f, 600.f}, {1.f, 1.f, 700.f}};
+    EXPECT_TRUE(cam.is_aabb_culled_by_frustum(aabb));
+}
+
+TEST(UnitTestProjection, AabbUnityEngineSideCulled)
+{
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(60.f);
+    const auto cam = omath::unity_engine::Camera({0, 0, 0}, {}, {1280.f, 720.f}, fov, 0.03f, 1000.f);
+
+    // Box far above (Unity: +Y up)
+    const omath::primitives::Aabb<float> aabb{{-1.f, 5000.f, 50.f}, {1.f, 6000.f, 100.f}};
+    EXPECT_TRUE(cam.is_aabb_culled_by_frustum(aabb));
+}
+
+TEST(UnitTestProjection, AabbUnityEngineStraddlesNearNotCulled)
+{
+    constexpr auto fov = omath::projection::FieldOfView::from_degrees(60.f);
+    const auto cam = omath::unity_engine::Camera({0, 0, 0}, {}, {1280.f, 720.f}, fov, 0.03f, 1000.f);
+
+    // Box straddles near plane (Unity: +Z forward)
+    const omath::primitives::Aabb<float> aabb{{-1.f, -1.f, -5.f}, {1.f, 1.f, 5.f}};
+    EXPECT_FALSE(cam.is_aabb_culled_by_frustum(aabb));
 }