Parallel CPU LBVH Implementation

.clang-format

@@ -1,5 +1,4 @@
 ---
-Language: Cpp
 BasedOnStyle: WebKit
 AlignAfterOpenBracket: AlwaysBreak
 AlignTrailingComments:
@@ -38,6 +37,5 @@ IncludeCategories:
     SortPriority:    1
     CaseSensitive:   true
 PackConstructorInitializers: CurrentLine
-RemoveEmptyLinesInUnwrappedLines: true
 SortIncludes: CaseInsensitive
 ...

CMakeLists.txt

@@ -82,15 +82,16 @@ if(CMAKE_CUDA_COMPILER)
 else()
   option(IPC_TOOLKIT_WITH_CUDA "Enable CUDA CCD" OFF)
 endif()
-option(IPC_TOOLKIT_WITH_RATIONAL_INTERSECTION "Use rational edge-triangle intersection check" OFF)
-option(IPC_TOOLKIT_WITH_ROBIN_MAP             "Use Tessil's robin-map rather than std maps"    ON)
-option(IPC_TOOLKIT_WITH_ABSEIL                "Use Abseil's hash functions"                    ON)
-option(IPC_TOOLKIT_WITH_FILIB                 "Use filib for interval arithmetic"              ON)
-option(IPC_TOOLKIT_WITH_INEXACT_CCD           "Use the original inexact CCD method of IPC"    OFF)
+option(IPC_TOOLKIT_WITH_RATIONAL_INTERSECTION "Use rational edge-triangle intersection check"           OFF)
+option(IPC_TOOLKIT_WITH_ROBIN_MAP             "Use Tessil's robin-map rather than std maps"              ON)
+option(IPC_TOOLKIT_WITH_ABSEIL                "Use Abseil's hash functions"                              ON)
+option(IPC_TOOLKIT_WITH_FILIB                 "Use filib for interval arithmetic"                        ON)
+option(IPC_TOOLKIT_WITH_INEXACT_CCD           "Use the original inexact CCD method of IPC"              OFF)
+option(IPC_TOOLKIT_WITH_PROFILER              "Enable performance profiler" ${IPC_TOOLKIT_TOPLEVEL_PROJECT})
 
 # Advanced options
-option(IPC_TOOLKIT_WITH_SIMD                  "Enable SIMD"                                   OFF)
-option(IPC_TOOLKIT_WITH_CODE_COVERAGE         "Enable coverage reporting"                     OFF)
+option(IPC_TOOLKIT_WITH_SIMD                  "Enable SIMD"                                             OFF)
+option(IPC_TOOLKIT_WITH_CODE_COVERAGE         "Enable coverage reporting"                               OFF)
 
 mark_as_advanced(IPC_TOOLKIT_WITH_SIMD)          # This does not work reliably
 mark_as_advanced(IPC_TOOLKIT_WITH_CODE_COVERAGE) # This is used in GitHub Actions
@@ -207,6 +208,12 @@ if(IPC_TOOLKIT_WITH_FILIB)
   target_link_libraries(ipc_toolkit PUBLIC filib::filib)
 endif()
 
+if(IPC_TOOLKIT_WITH_PROFILER)
+  # Add nlohmann/json for the profiler
+  include(json)
+  target_link_libraries(ipc_toolkit PRIVATE nlohmann_json::nlohmann_json)
+endif()
+
 # Extra warnings (link last for highest priority)
 include(ipc_toolkit_warnings)
 target_link_libraries(ipc_toolkit PRIVATE ipc::toolkit::warnings)

python/examples/find_ipctk.py

@@ -0,0 +1,4 @@
+import sys
+import pathlib
+sys.path.append(str(pathlib.Path(__file__).parents[1]))
+from _find_ipctk import ipctk  # noqa

python/examples/lbvh.py

@@ -0,0 +1,77 @@
+from find_ipctk import ipctk
+import meshio
+import polyscope as ps
+from polyscope import imgui
+import numpy as np
+
+import pathlib
+
+mesh = meshio.read(pathlib.Path(
+    __file__).parents[2] / "tests/data/puffer-ball/20.ply")
+
+lbvh = ipctk.LBVH()
+lbvh.build(mesh.points, np.array([], dtype=int), mesh.cells_dict["triangle"])
+
+ps.init()
+
+ps.set_give_focus_on_show(True)
+
+ps_mesh = ps.register_surface_mesh(
+    "bunny",
+    mesh.points,
+    mesh.cells_dict["triangle"]
+)
+
+nodes = lbvh.vertex_nodes
+
+
+def traverse_lbvh(node, max_depth):
+    if node.is_inner and max_depth > 0:
+        V_left, E_left = traverse_lbvh(nodes[node.left], max_depth - 1)
+        V_right, E_right = traverse_lbvh(nodes[node.right], max_depth - 1)
+        return np.vstack([V_left, V_right]), np.vstack([E_left, E_right + V_left.shape[0]])
+
+    E = np.array([
+        [0, 1],
+        [0, 2],
+        [0, 3],
+        [1, 5],
+        [1, 4],
+        [2, 4],
+        [2, 6],
+        [3, 5],
+        [3, 6],
+        [7, 4],
+        [7, 5],
+        [7, 6],
+    ])
+    V = np.array([
+        node.aabb_min,
+        [node.aabb_min[0], node.aabb_min[1], node.aabb_max[2]],
+        [node.aabb_min[0], node.aabb_max[1], node.aabb_min[2]],
+        [node.aabb_max[0], node.aabb_min[1], node.aabb_min[2]],
+        [node.aabb_min[0], node.aabb_max[1], node.aabb_max[2]],
+        [node.aabb_max[0], node.aabb_min[1], node.aabb_max[2]],
+        [node.aabb_max[0], node.aabb_max[1], node.aabb_min[2]],
+        node.aabb_max
+    ])
+    return V, E
+
+
+max_depth = 0
+bvh_nodes, bvh_edges = traverse_lbvh(nodes[0], max_depth=max_depth)
+
+ps.register_curve_network("bvh", bvh_nodes, bvh_edges)
+
+
+def foo():
+    global max_depth
+    changed, max_depth = imgui.SliderInt("max depth", max_depth, 0, 20)
+    if changed:
+        bvh_nodes, bvh_edges = traverse_lbvh(nodes[0], max_depth=max_depth)
+        ps.register_curve_network("bvh", bvh_nodes, bvh_edges)
+
+
+ps.set_user_callback(foo)
+
+ps.show()

python/src/bindings.cpp

@@ -27,6 +27,7 @@ PYBIND11_MODULE(ipctk, m)
     define_brute_force(m);
     define_bvh(m);
     define_hash_grid(m);
+    define_lbvh(m);
     define_spatial_hash(m);
     define_sweep_and_prune(m);
     define_sweep_and_tiniest_queue(m);

python/src/broad_phase/CMakeLists.txt

@@ -4,6 +4,7 @@ set(SOURCES
   brute_force.cpp
   bvh.cpp
   hash_grid.cpp
+  lbvh.cpp
   spatial_hash.cpp
   sweep_and_prune.cpp
   sweep_and_tiniest_queue.cpp

python/src/broad_phase/bindings.hpp

@@ -7,6 +7,7 @@ void define_broad_phase(py::module_& m);
 void define_brute_force(py::module_& m);
 void define_bvh(py::module_& m);
 void define_hash_grid(py::module_& m);
+void define_lbvh(py::module_& m);
 void define_spatial_hash(py::module_& m);
 void define_sweep_and_prune(py::module_& m);
 void define_sweep_and_tiniest_queue(py::module_& m);

python/src/broad_phase/lbvh.cpp

@@ -0,0 +1,23 @@
+#include <common.hpp>
+
+#include <ipc/broad_phase/lbvh.hpp>
+
+using namespace ipc;
+
+void define_lbvh(py::module_& m)
+{
+    py::class_<LBVH::Node>(m, "LBVH_Node")
+        .def_readonly("aabb_min", &LBVH::Node::aabb_min)
+        .def_readonly("aabb_max", &LBVH::Node::aabb_max)
+        .def_readonly("left", &LBVH::Node::left)
+        .def_readonly("right", &LBVH::Node::right)
+        .def_property_readonly("is_leaf", &LBVH::Node::is_leaf)
+        .def_property_readonly("is_inner", &LBVH::Node::is_inner)
+        .def_property_readonly("is_valid", &LBVH::Node::is_valid);
+
+    py::class_<LBVH, BroadPhase, std::shared_ptr<LBVH>>(m, "LBVH")
+        .def(py::init())
+        .def_property_readonly("vertex_nodes", &LBVH::vertex_nodes)
+        .def_property_readonly("edge_nodes", &LBVH::edge_nodes)
+        .def_property_readonly("face_nodes", &LBVH::face_nodes);
+}

src/ipc/broad_phase/CMakeLists.txt

@@ -10,6 +10,8 @@ set(SOURCES
   default_broad_phase.hpp
   hash_grid.cpp
   hash_grid.hpp
+  lbvh.cpp
+  lbvh.hpp
   spatial_hash.cpp
   spatial_hash.hpp
   sweep_and_prune.cpp

src/ipc/broad_phase/bvh.cpp

@@ -1,6 +1,7 @@
 #include "bvh.hpp"
 
 #include <ipc/utils/merge_thread_local.hpp>
+#include <ipc/utils/profiler.hpp>
 
 #include <SimpleBVH/BVH.hpp>
 #include <tbb/blocked_range.h>
@@ -37,6 +38,8 @@ void BVH::init_bvh(const std::vector<AABB>& boxes, SimpleBVH::BVH& bvh)
         return;
     }
 
+    IPC_TOOLKIT_PROFILE_BLOCK("BVH::init_bvh");
+
     std::vector<std::array<Eigen::Vector3d, 2>> vector_boxes(boxes.size());
     for (int i = 0; i < boxes.size(); i++) {
         vector_boxes[i] = { { to_3D(boxes[i].min), to_3D(boxes[i].max) } };
@@ -47,6 +50,7 @@ void BVH::init_bvh(const std::vector<AABB>& boxes, SimpleBVH::BVH& bvh)
 
 void BVH::clear()
 {
+    BroadPhase::clear();
     vertex_bvh->clear();
     edge_bvh->clear();
     face_bvh->clear();
@@ -105,6 +109,8 @@ void BVH::detect_vertex_vertex_candidates(
         return;
     }
 
+    IPC_TOOLKIT_PROFILE_BLOCK("BVH::detect_vertex_vertex_candidates");
+
     detect_candidates<
         VertexVertexCandidate, /*swap_order=*/false, /*triangular=*/true>(
         vertex_boxes, *vertex_bvh, can_vertices_collide, candidates);
@@ -117,6 +123,8 @@ void BVH::detect_edge_vertex_candidates(
         return;
     }
 
+    IPC_TOOLKIT_PROFILE_BLOCK("BVH::detect_edge_vertex_candidates");
+
     // In 2D and for codimensional edge-vertex collisions, there are more
     // vertices than edges, so we want to iterate over the edges.
     detect_candidates(
@@ -131,6 +139,8 @@ void BVH::detect_edge_edge_candidates(
         return;
     }
 
+    IPC_TOOLKIT_PROFILE_BLOCK("BVH::detect_edge_edge_candidates");
+
     detect_candidates<
         EdgeEdgeCandidate, /*swap_order=*/false, /*triangular=*/true>(
         edge_boxes, *edge_bvh, std::bind(&BVH::can_edges_collide, this, _1, _2),
@@ -144,6 +154,8 @@ void BVH::detect_face_vertex_candidates(
         return;
     }
 
+    IPC_TOOLKIT_PROFILE_BLOCK("BVH::detect_face_vertex_candidates");
+
     // The ratio vertices:faces is 1:2, so we want to iterate over the vertices.
     detect_candidates<FaceVertexCandidate, /*swap_order=*/true>(
         vertex_boxes, *face_bvh,
@@ -157,6 +169,8 @@ void BVH::detect_edge_face_candidates(
         return;
     }
 
+    IPC_TOOLKIT_PROFILE_BLOCK("BVH::detect_edge_face_candidates");
+
     // The ratio edges:faces is 3:2, so we want to iterate over the faces.
     detect_candidates<EdgeFaceCandidate, /*swap_order=*/true>(
         face_boxes, *edge_bvh,
@@ -170,6 +184,8 @@ void BVH::detect_face_face_candidates(
         return;
     }
 
+    IPC_TOOLKIT_PROFILE_BLOCK("BVH::detect_face_face_candidates");
+
     detect_candidates<
         FaceFaceCandidate, /*swap_order=*/false, /*triangular=*/true>(
         face_boxes, *face_bvh, std::bind(&BVH::can_faces_collide, this, _1, _2),