Picking and seeing through a tree of boxes per object (made again only when it changes)

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

src/CMakeLists.txt

@@ -1,6 +1,6 @@
 # Fry3DX: the app (app.cpp) on the engine, and its modelling core (mesh, scene) as a library the
 # tests use too.
-add_library(fry_core STATIC mesh.cpp tools.cpp cut.cpp uv.cpp image.cpp sticker.cpp shapes.cpp scene.cpp package.cpp export.cpp print.cpp boolean.cpp modify.cpp sculpt.cpp)
+add_library(fry_core STATIC mesh.cpp bvh.cpp tools.cpp cut.cpp uv.cpp image.cpp sticker.cpp shapes.cpp scene.cpp package.cpp export.cpp print.cpp boolean.cpp modify.cpp sculpt.cpp)
 target_include_directories(fry_core PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
 target_link_libraries(fry_core PUBLIC glm::glm libzstd_static)
 target_compile_definitions(fry_core PUBLIC GLM_ENABLE_EXPERIMENTAL)

src/app.cpp

@@ -3,6 +3,7 @@
 // and clear while you model (Solid), or properly lit (Rendered). The panels are Dear ImGui.
 //
 //   fry3dx [file.fry]
+#include "bvh.h"
 #include "commands.h"
 #include "icons.h"
 #include "mesh.h"
@@ -3966,15 +3967,27 @@ private:
         for (const Object& o : scene_.objects) {
             if (!o.visible || (only && o.id != only)) continue;
             const Mesh& mesh = base || (mode_ == Mode::Edit && o.id == active_) ? o.mesh : shownMesh(o);
-            const glm::mat4 m = objectMatrix(o);
+            if (auto h = treeOf(o, mesh).nearest(ray.origin, ray.dir, best.t)) best = {o.id, h->face, h->t};
-            std::vector<glm::vec3> w(mesh.verts.size());
-            for (size_t i = 0; i < w.size(); ++i) w[i] = glm::vec3(m * glm::vec4(mesh.verts[i], 1.0f));
-            for (size_t f = 0; f < mesh.faces.size(); ++f)
-                for (const auto& tri : triangulate(mesh, int(f)))
-                    if (auto t = hitTriangle(ray, w[size_t(tri[0])], w[size_t(tri[1])], w[size_t(tri[2])]); t && *t < best.t) best = {o.id, int(f), *t};
         }
         return best;
     }
+    // Each object's triangles in the world, in a tree of boxes: made again only when it changes
+    // (its shape, where it is), not each time something's picked.
+    struct Tree {
+        uint64_t version = 0;
+        glm::mat4 matrix{0.0f};
+        size_t faces = 0;
+        fry::TriangleTree tree;
+    };
+    std::map<std::pair<uint32_t, bool>, Tree> trees_; // (an object's own mesh, or the one shown)
+    const fry::TriangleTree& treeOf(const Object& o, const Mesh& mesh) {
+        Tree& t = trees_[{o.id, &mesh != &o.mesh}];
+        const glm::mat4 m = objectMatrix(o);
+        const uint64_t version = versions_[o.id];
+        if (t.version != version || !version || t.matrix != m || t.faces != mesh.faces.size() || t.tree.points().size() != mesh.verts.size())
+            t.version = version, t.matrix = m, t.faces = mesh.faces.size(), t.tree = fry::TriangleTree(mesh, m);
+        return t.tree;
+    }
     // Can the eye see this point of the object (nothing of it in the way)?
     bool visible(const Object& o, const glm::vec3& world) {
         if (xray_) return true;
@@ -4000,48 +4013,16 @@ private:
         if (v.version == version && glm::distance(v.eye, eye_) < 1e-6f && v.vp == vp_ && v.point.size() == mesh.verts.size()) return v.point;
         v.version = version, v.eye = eye_, v.vp = vp_;
         v.point.assign(mesh.verts.size(), 1);
-        const glm::mat4 m = objectMatrix(o);
+        // (each point: is anything of it between the point and the eye)
-        std::vector<glm::vec3> w(mesh.verts.size());
+        const fry::TriangleTree& tree = treeOf(o, mesh);
-        std::vector<glm::vec4> clip(w.size());
+        const std::vector<glm::vec3>& w = tree.points();
-        for (size_t i = 0; i < w.size(); ++i) w[i] = glm::vec3(m * glm::vec4(mesh.verts[i], 1.0f)), clip[i] = vp_ * glm::vec4(w[i], 1.0f);
-        std::vector<std::array<int, 3>> tris;
-        for (size_t f = 0; f < mesh.faces.size(); ++f)
-            for (const auto& t : triangulate(mesh, int(f))) tris.push_back(t);
-        constexpr int G = 64;
-        auto cell = [](float ndc) { return std::clamp(int((ndc * 0.5f + 0.5f) * G), 0, G - 1); };
-        std::vector<std::vector<int>> grid(G * G);
-        std::vector<int> everywhere; // (reaching behind the eye: in every square)
-        for (size_t t = 0; t < tris.size(); ++t) {
-            glm::vec2 lo(1e30f), hi(-1e30f);
-            bool behind = false;
-            for (int k : tris[t]) {
-                const glm::vec4 c = clip[size_t(k)];
-                if (c.w < 1e-5f) { behind = true; break; }
-                lo = glm::min(lo, glm::vec2(c) / c.w), hi = glm::max(hi, glm::vec2(c) / c.w);
-            }
-            if (behind) { everywhere.push_back(int(t)); continue; }
-            if (hi.x < -1 || hi.y < -1 || lo.x > 1 || lo.y > 1) continue; // (off the screen: covers nothing on it)
-            for (int y = cell(lo.y); y <= cell(hi.y); ++y)
-                for (int x = cell(lo.x); x <= cell(hi.x); ++x) grid[size_t(y * G + x)].push_back(int(t));
-        }
         for (size_t i = 0; i < w.size(); ++i) {
-            const glm::vec4 c = clip[i];
+            const glm::vec4 c = vp_ * glm::vec4(w[i], 1.0f);
-            if (c.w < 1e-5f) continue; // (behind the eye: not drawn anyway)
+            if (c.w < 1e-5f || std::abs(c.x) > c.w || std::abs(c.y) > c.w) continue; // (off the screen: not drawn anyway)
-            const glm::vec2 n = glm::vec2(c) / c.w;
-            if (std::abs(n.x) > 1 || std::abs(n.y) > 1) continue;
             const glm::vec3 d = w[i] - eye_;
             const float len = glm::length(d);
-            const Ray ray{eye_, d / len};
+            if (len < 1e-6f) continue;
-            const float limit = len - std::max(len * 0.002f, 0.0005f);
+            if (tree.blocked(eye_, d / len, len - std::max(len * 0.002f, 0.0005f), int(i))) v.point[i] = 0;
-            auto hides = [&](int t) {
-                const auto& tri = tris[size_t(t)];
-                if (tri[0] == int(i) || tri[1] == int(i) || tri[2] == int(i)) return false;
-                auto hit = hitTriangle(ray, w[size_t(tri[0])], w[size_t(tri[1])], w[size_t(tri[2])]);
-                return hit && *hit < limit;
-            };
-            if (std::any_of(everywhere.begin(), everywhere.end(), hides)) { v.point[i] = 0; continue; }
-            const auto& here = grid[size_t(cell(n.y) * G + cell(n.x))];
-            if (std::any_of(here.begin(), here.end(), hides)) v.point[i] = 0;
         }
         return v.point;
     }

src/bvh.cpp

@@ -0,0 +1,119 @@
+#include "bvh.h"
+
+#include <algorithm>
+#include <cmath>
+
+namespace fry {
+
+std::optional<float> rayTriangle(const glm::vec3& origin, const glm::vec3& dir, const glm::vec3& a, const glm::vec3& b, const glm::vec3& c) {
+    const glm::vec3 e1 = b - a, e2 = c - a, p = glm::cross(dir, e2);
+    const float det = glm::dot(e1, p);
+    if (std::abs(det) < 1e-12f) return std::nullopt;
+    const float inv = 1.0f / det;
+    const glm::vec3 s = origin - a;
+    const float u = glm::dot(s, p) * inv;
+    if (u < 0 || u > 1) return std::nullopt;
+    const glm::vec3 q = glm::cross(s, e1);
+    const float v = glm::dot(dir, q) * inv;
+    if (v < 0 || u + v > 1) return std::nullopt;
+    const float t = glm::dot(e2, q) * inv;
+    return t > 1e-6f ? std::optional<float>(t) : std::nullopt;
+}
+
+namespace {
+
+// Where the ray goes into the box and out, or nothing (in the ray's lengths).
+bool rayBox(const glm::vec3& o, const glm::vec3& inv, const glm::vec3& lo, const glm::vec3& hi, float maxT) {
+    const glm::vec3 a = (lo - o) * inv, b = (hi - o) * inv;
+    const glm::vec3 n = glm::min(a, b), f = glm::max(a, b);
+    const float in = std::max(std::max(n.x, n.y), std::max(n.z, 0.0f)), out = std::min(std::min(f.x, f.y), std::min(f.z, maxT));
+    return in <= out;
+}
+
+} // namespace
+
+TriangleTree::TriangleTree(const Mesh& m, const glm::mat4& matrix) {
+    points_.resize(m.verts.size());
+    for (size_t i = 0; i < points_.size(); ++i) points_[i] = glm::vec3(matrix * glm::vec4(m.verts[i], 1.0f));
+    for (size_t f = 0; f < m.faces.size(); ++f) {
+        const auto& fv = m.faces[f];
+        if (fv.size() == 3) tris_.push_back({{fv[0], fv[1], fv[2]}, int(f)}); // (the usual ones, without triangulating)
+        else if (fv.size() == 4 && glm::dot(glm::cross(m.verts[size_t(fv[1])] - m.verts[size_t(fv[0])], m.verts[size_t(fv[2])] - m.verts[size_t(fv[0])]),
+                                            glm::cross(m.verts[size_t(fv[2])] - m.verts[size_t(fv[0])], m.verts[size_t(fv[3])] - m.verts[size_t(fv[0])])) > 0)
+            tris_.push_back({{fv[0], fv[1], fv[2]}, int(f)}), tris_.push_back({{fv[0], fv[2], fv[3]}, int(f)}); // (a quad with no dent)
+        else
+            for (const auto& t : triangulate(m, int(f))) tris_.push_back({t, int(f)});
+    }
+    if (!tris_.empty()) nodes_.reserve(tris_.size() * 2 / 3 + 1), build(0, int(tris_.size()));
+}
+
+int TriangleTree::build(int first, int count) {
+    const int at = int(nodes_.size());
+    nodes_.push_back({});
+    glm::vec3 lo(1e30f), hi(-1e30f), clo(1e30f), chi(-1e30f);
+    for (int i = first; i < first + count; ++i) {
+        glm::vec3 c(0.0f);
+        for (int k : tris_[size_t(i)].v) lo = glm::min(lo, points_[size_t(k)]), hi = glm::max(hi, points_[size_t(k)]), c += points_[size_t(k)];
+        clo = glm::min(clo, c / 3.0f), chi = glm::max(chi, c / 3.0f);
+    }
+    nodes_[size_t(at)].lo = lo, nodes_[size_t(at)].hi = hi;
+    if (count <= 4) {
+        nodes_[size_t(at)].first = first, nodes_[size_t(at)].count = count;
+        return at;
+    }
+    // Split down the middle of the widest way the triangles' middles spread.
+    const glm::vec3 size = chi - clo;
+    const int axis = size.x > size.y && size.x > size.z ? 0 : size.y > size.z ? 1 : 2;
+    auto mid = [&](const Tri& t) { return points_[size_t(t.v[0])][axis] + points_[size_t(t.v[1])][axis] + points_[size_t(t.v[2])][axis]; };
+    const int half = count / 2;
+    std::nth_element(tris_.begin() + first, tris_.begin() + first + half, tris_.begin() + first + count, [&](const Tri& a, const Tri& b) { return mid(a) < mid(b); });
+    build(first, half);
+    const int right = build(first + half, count - half);
+    nodes_[size_t(at)].first = right, nodes_[size_t(at)].count = 0;
+    return at;
+}
+
+std::optional<TriangleTree::Hit> TriangleTree::nearest(const glm::vec3& o, const glm::vec3& d, float maxT) const {
+    if (nodes_.empty()) return std::nullopt;
+    const glm::vec3 inv(1.0f / d.x, 1.0f / d.y, 1.0f / d.z);
+    std::optional<Hit> best;
+    int stack[128], n = 0;
+    stack[n++] = 0;
+    while (n) {
+        const Node& node = nodes_[size_t(stack[--n])];
+        if (!rayBox(o, inv, node.lo, node.hi, best ? best->t : maxT)) continue;
+        if (node.count) {
+            for (int i = node.first; i < node.first + node.count; ++i) {
+                const Tri& t = tris_[size_t(i)];
+                if (auto h = rayTriangle(o, d, points_[size_t(t.v[0])], points_[size_t(t.v[1])], points_[size_t(t.v[2])]); h && *h < (best ? best->t : maxT))
+                    best = Hit{*h, t.face, t.v};
+            }
+        } else if (n + 2 <= 128) {
+            stack[n++] = node.first, stack[n++] = int(&node - nodes_.data()) + 1;
+        }
+    }
+    return best;
+}
+
+bool TriangleTree::blocked(const glm::vec3& o, const glm::vec3& d, float maxT, int skip) const {
+    if (nodes_.empty()) return false;
+    const glm::vec3 inv(1.0f / d.x, 1.0f / d.y, 1.0f / d.z);
+    int stack[128], n = 0;
+    stack[n++] = 0;
+    while (n) {
+        const Node& node = nodes_[size_t(stack[--n])];
+        if (!rayBox(o, inv, node.lo, node.hi, maxT)) continue;
+        if (node.count) {
+            for (int i = node.first; i < node.first + node.count; ++i) {
+                const Tri& t = tris_[size_t(i)];
+                if (t.v[0] == skip || t.v[1] == skip || t.v[2] == skip) continue;
+                if (auto h = rayTriangle(o, d, points_[size_t(t.v[0])], points_[size_t(t.v[1])], points_[size_t(t.v[2])]); h && *h < maxT) return true;
+            }
+        } else if (n + 2 <= 128) {
+            stack[n++] = node.first, stack[n++] = int(&node - nodes_.data()) + 1;
+        }
+    }
+    return false;
+}
+
+} // namespace fry

src/bvh.h

@@ -0,0 +1,51 @@
+#pragma once
+// A tree of boxes round a mesh's triangles (a bounding volume hierarchy): what a ray hits first,
+// or whether anything's in its way, without trying every triangle.
+#include "mesh.h"
+
+#include <glm/glm.hpp>
+
+#include <array>
+#include <cstdint>
+#include <optional>
+#include <vector>
+
+namespace fry {
+
+class TriangleTree {
+public:
+    TriangleTree() = default;
+    // The mesh's triangles, its points moved by `matrix` (into the world).
+    TriangleTree(const Mesh& m, const glm::mat4& matrix);
+
+    struct Hit {
+        float t = 0;
+        int face = -1;           // (the mesh's face it's part of)
+        std::array<int, 3> tri{}; // (its points)
+    };
+    // The nearest hit along the ray (dir needn't be unit length: t is in its lengths), before maxT.
+    std::optional<Hit> nearest(const glm::vec3& origin, const glm::vec3& dir, float maxT = 1e30f) const;
+    // Anything in the way before maxT, but triangles with point `skip` (the point being looked at).
+    bool blocked(const glm::vec3& origin, const glm::vec3& dir, float maxT, int skip = -1) const;
+    bool empty() const { return tris_.empty(); }
+    const std::vector<glm::vec3>& points() const { return points_; } // (in the world)
+
+private:
+    struct Tri {
+        std::array<int, 3> v;
+        int face;
+    };
+    struct Node {
+        glm::vec3 lo, hi;
+        int first, count; // (a leaf: its triangles; else count 0 and `first` is the right child, the left one next)
+    };
+    std::vector<glm::vec3> points_;
+    std::vector<Tri> tris_;
+    std::vector<Node> nodes_;
+    int build(int first, int count);
+};
+
+// Where a ray meets a triangle (Möller–Trumbore), in the ray's lengths.
+std::optional<float> rayTriangle(const glm::vec3& origin, const glm::vec3& dir, const glm::vec3& a, const glm::vec3& b, const glm::vec3& c);
+
+} // namespace fry

tests/mesh_test.cpp

@@ -1,5 +1,6 @@
 // The modelling core without a window: shapes are closed and face out, faces with dents
 // triangulate right, extrude / inset / merge keep the mesh whole, and files round-trip.
+#include "bvh.h"
 #include "mesh.h"
 #include "image.h"
 #include "scene.h"
@@ -362,6 +363,25 @@ int main() {
         CHECK(combine(a, drill, Combine::Subtract) && checkForPrint(a).closed() && checkForPrint(a).volume < 1.0f && checkForPrint(a).volume > 0.85f,
               "a cylinder drilled through a cube: closed, %g m3", checkForPrint(a).volume);
     }
+    // The tree of boxes finds the same nearest hit as trying every triangle (and the same "in the way").
+    {
+        Mesh m = makeSphere(32, 16, 1.0f);
+        const glm::mat4 M = glm::translate(glm::mat4(1), glm::vec3(0.3f, -0.2f, 0.5f)) * glm::rotate(glm::mat4(1), 0.7f, glm::vec3(0, 1, 0));
+        const TriangleTree tree(m, M);
+        uint32_t seed = 7;
+        auto rnd = [&] { return (seed = seed * 1664525u + 1013904223u) / 4294967296.0f * 2 - 1; };
+        int same = 0, tries = 300;
+        for (int i = 0; i < tries; ++i) {
+            const glm::vec3 o(rnd() * 3, rnd() * 3, rnd() * 3 + 3), d = glm::normalize(glm::vec3(rnd(), rnd(), rnd() - 1.5f));
+            float best = 1e30f;
+            for (size_t f = 0; f < m.faces.size(); ++f)
+                for (const auto& t : triangulate(m, int(f)))
+                    if (auto h = rayTriangle(o, d, glm::vec3(M * glm::vec4(m.verts[size_t(t[0])], 1)), glm::vec3(M * glm::vec4(m.verts[size_t(t[1])], 1)), glm::vec3(M * glm::vec4(m.verts[size_t(t[2])], 1)))) best = std::min(best, *h);
+            const auto h = tree.nearest(o, d);
+            if ((best > 1e29f && !h) || (h && std::abs(h->t - best) < 1e-4f && tree.blocked(o, d, best + 1e-3f) && !tree.blocked(o, d, best - 1e-3f))) ++same;
+        }
+        CHECK(same == tries, "the tree of boxes finds what trying every triangle does (%d of %d)", same, tries);
+    }
     // Shapes drawn and cut in, again and again inside the last cut, on a turned object (the shape
     // carried into its own space, as the Draw tool does). This once built its tree forever and
     // crashed; and the cutting's rounding left small gaps.